12 meshPositionsFieldName(
"MESH_NODE_POSITIONS"), aRea(elementMeasure) {}
29 auto get_ftensor_n_diff = [&]() {
37 &
m(0, 0), &
m(0, 1), &
m(0, 2));
42 const size_t nb_gauss_pts =
gaussPts.size2();
58 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
67 }
else if (
type == MBQUAD) {
75 const size_t nb_gauss_pts =
gaussPts.size2();
90 auto t_diff = get_ftensor_n_diff();
91 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
92 auto t_coords = get_ftensor_from_vec_3d(
coords);
94 t_t1(
i) += t_coords(
i) * t_diff(N0);
95 t_t2(
i) += t_coords(
i) * t_diff(N1);
107 "Element type not implemented");
124 auto calc_normal = [&](
const double *diff_ptr) {
135 diff_ptr, &diff_ptr[1]);
146 for (
int nn = 0; nn !=
num_nodes; ++nn) {
147 t_t1(
i) += t_coords(
i) * t_diff(N0);
148 t_t2(
i) += t_coords(
i) * t_diff(N1);
153 aRea = sqrt(t_normal(
i) * t_normal(
i));
157 const double *diff_ptr;
161 CHKERR calc_normal(diff_ptr);
167 CHKERR calc_normal(diff_ptr);
171 "Element type not implemented");
186 auto get_rule_by_type = [&]() {
189 return getRule(order_row + 1, order_col + 1, order_data + 1);
191 return getRule(order_row, order_col, order_data);
195 const int rule = get_rule_by_type();
197 auto set_integration_pts_for_tri = [&]() {
208 gaussPts.resize(3, nb_gauss_pts,
false);
209 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
211 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
219 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
232 auto calc_base_for_tri = [&]() {
234 const size_t nb_gauss_pts =
gaussPts.size2();
237 base.resize(nb_gauss_pts, 3,
false);
238 diff_base.resize(3, 2,
false);
247 auto calc_base_for_quad = [&]() {
249 const size_t nb_gauss_pts =
gaussPts.size2();
252 base.resize(nb_gauss_pts, 4,
false);
253 diff_base.resize(nb_gauss_pts, 8,
false);
254 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
276 CHKERR set_integration_pts_for_tri();
281 CHKERR calc_base_for_quad();
285 "Element type not implemented: %d",
type);
291 const size_t nb_gauss_pts =
gaussPts.size2();
295 CHKERR calc_base_for_tri();
298 CHKERR calc_base_for_quad();
302 "Element type not implemented: %d",
type);
378 const size_t nb_nodes =
380 double *shape_functions =
382 const size_t nb_gauss_pts =
gaussPts.size2();
384 for (
int gg = 0; gg != nb_gauss_pts; ++gg)
385 for (
int dd = 0;
dd != 3; ++
dd)
387 nb_nodes, &shape_functions[nb_nodes * gg], 1, &
coords[
dd], 3);
397 "User operator and finite element do not work together");
444 return loopSide(fe_name, &fe_method, 3);
452 if (toElePtr->gaussPts.size1() != getGaussPts().size1()) {
454 "Inconsistent numer of weights %d != %d",
455 toElePtr->gaussPts.size1(), getGaussPts().size1());
457 if (toElePtr->gaussPts.size2() != getGaussPts().size2()) {
459 "Inconsistent numer of integtaion pts %d != %d",
460 toElePtr->gaussPts.size2(), getGaussPts().size2());
465 switch (getFEType()) {
470 "Element type not implemented");
475 &
m(0, 0), &
m(0, 1), &
m(0, 2));
480 auto get_local_coords_triangle = [&]() {
481 double ref_gauss_pts[2][3] = {{0, 1, 0}, {0, 0, 1}};
482 std::array<double, 3> ksi0 = {0, 1, 0};
483 std::array<double, 3> ksi1 = {0, 0, 1};
484 std::array<double, 9> ref_shapes;
485 CHKERR Tools::shapeFunMBTRI<1>(ref_shapes.data(), ksi0.data(), ksi1.data(),
487 auto &node_coords = getCoords();
488 auto &glob_coords = toElePtr->coords;
489 std::array<double, 6> local_coords;
491 &*node_coords.begin(), &*glob_coords.begin(), 3, local_coords.data());
496 auto get_diff_triangle = [&]() {
499 diff_ptr, &diff_ptr[1]);
503 auto get_jac = [&](
auto &&local_coords,
auto &&t_diff) {
507 auto t_local_coords = getFTensor1FromPtr<2>(local_coords.data());
509 for (
int nn = 0; nn != 3; ++nn) {
510 t_jac(
I,
J) += t_local_coords(
I) * t_diff(
J);
518 auto t_mat_tangent = [&](
auto &t1,
auto &t2) {
520 &t1(0), &t1(1), &t1(2), &t2(0), &t2(1), &t2(2)};
524 auto transfrom = [&](
auto &&t_mat_t,
auto &&t_mat_out_t,
auto &&t_inv_jac) {
530 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
532 t_mat_out_t(
J,
i) = t_mat_t(
I,
i) * t_inv_jac(
I,
J);
539 auto calc_normal = [&](
auto &
n,
auto &t1,
auto &t2) {
543 auto t_t1 = get_ftensor_from_mat_3d(t1);
544 auto t_t2 = get_ftensor_from_mat_3d(t2);
545 auto t_n = get_ftensor_from_mat_3d(
n);
546 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
556 t_mat_tangent(getTangent1AtGaussPts(), getTangent2AtGaussPts()),
557 t_mat_tangent(toElePtr->tangentOneAtGaussPts,
558 toElePtr->tangentTwoAtGaussPts),
559 get_jac(get_local_coords_triangle(), get_diff_triangle())
562 calc_normal(toElePtr->normalsAtGaussPts, toElePtr->tangentOneAtGaussPts,
563 toElePtr->tangentTwoAtGaussPts);