10 boost::shared_ptr<MatrixDouble> jac_ptr)
14 for (
auto t = MBEDGE;
t != MBMAXTYPE; ++
t)
26 const auto nb_base_functions = data.
getN(
NOBASE).size2();
27 if (nb_base_functions) {
34 if (nb_gauss_pts != data.
getDiffN().size1())
36 "Inconsistent number base functions and gauss points");
37 if (nb_base_functions != data.
getDiffN().size2() / 3)
39 "Inconsistent number of base functions");
40 if (coords.size() != 3 * nb_base_functions)
43 "Number of vertex coordinates and base functions is inconsistent "
45 coords.size() / 3, nb_base_functions);
48 jacPtr->resize(9, nb_gauss_pts,
false);
50 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
57 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
59 &coords[0], &coords[1], &coords[2]);
60 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
61 t_jac(
i,
j) += t_coords(
i) * (t_vol_inv_jac(
k,
j) * t_diff_base(
k));
81 return applyTransform<3, 3, 3, 3>(diff_n);
93 CHKERR transform_base(*(
m.second));
98 CHKERR transform_base(*ptr);
122 data.
getDiffN(base).size2(),
false);
124 unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
125 unsigned int nb_base_functions = data.
getDiffN(base).size2() / 9;
126 if (nb_base_functions == 0)
132 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
137 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
invJacPtr);
138 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
139 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
140 t_inv_diff_n(
i,
j) = t_inv_jac(
k,
j) * t_diff_n(
i,
k);
165 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
166 t_w *= sqrt(t_normal(
i) * t_normal(
i)) /
a;
172 "Number of rows in getNormalsAtGaussPts should be equal to "
173 "number of integration points, but is not, i.e. %d != %d",
190 for (
size_t gg = 0; gg != nb_int_pts; ++gg) {
191 t_w *= sqrt(t_tangent(
i) * t_tangent(
i)) /
a;
197 "Number of rows in getTangentAtGaussPts should be equal to "
198 "number of integration points, but is not, i.e. %d != %d",
209 const auto nb_integration_pts =
detPtr->size();
214 "Inconsistent number of data points");
219 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
240 "Pointer for detPtr not allocated");
243 "Pointer for jacPtr not allocated");
251 auto nb_base_functions =
253 auto nb_diff_base_functions =
257 if (nb_diff_base_functions) {
258 if (nb_diff_base_functions != nb_base_functions)
260 "Wrong number base functions %d != %d", nb_diff_base_functions,
262 if (data.
getDiffN(base).size1() != nb_gauss_pts)
264 "Wrong number integration points");
268 if (nb_gauss_pts && nb_base_functions) {
270 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
273 double *t_transformed_n_ptr = &*
piolaN.data().begin();
276 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
279 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
281 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
282 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
283 const double a = 1. / t_det;
284 t_transformed_n(
i) =
a * (t_jac(
i,
k) * t_n(
k));
295 if (nb_gauss_pts && nb_diff_base_functions) {
300 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
302 t_transformed_diff_n(t_transformed_diff_n_ptr,
303 &t_transformed_diff_n_ptr[
HVEC0_1],
304 &t_transformed_diff_n_ptr[
HVEC0_2],
305 &t_transformed_diff_n_ptr[
HVEC1_0],
306 &t_transformed_diff_n_ptr[
HVEC1_1],
307 &t_transformed_diff_n_ptr[
HVEC1_2],
308 &t_transformed_diff_n_ptr[
HVEC2_0],
309 &t_transformed_diff_n_ptr[
HVEC2_1],
310 &t_transformed_diff_n_ptr[
HVEC2_2]);
313 auto t_jac = getFTensor2FromMat<3, 3>(*
jacPtr);
315 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
316 for (
unsigned int bb = 0; bb != nb_diff_base_functions; ++bb) {
317 const double a = 1. / t_det;
318 t_transformed_diff_n(
i,
k) =
a * (t_jac(
i,
j) * t_diff_n(
j,
k));
320 ++t_transformed_diff_n;
346 unsigned int nb_gauss_pts = data.
getN(base).size1();
347 unsigned int nb_base_functions = data.
getN(base).size2() / 3;
348 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
352 double *t_transformed_n_ptr = &*
piolaN.data().begin();
355 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
357 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
359 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
360 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
361 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
362 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
363 &t_transformed_diff_n_ptr[
HVEC2_2]);
365 auto t_inv_jac = getFTensor2FromMat<3, 3>(*
jacInvPtr);
367 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
368 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
369 t_transformed_n(
i) = t_inv_jac(
k,
i) * t_n(
k);
370 t_transformed_diff_n(
i,
k) = t_inv_jac(
j,
i) * t_diff_n(
j,
k);
374 ++t_transformed_diff_n;
387 boost::shared_ptr<MatrixDouble> jac_ptr)
397 const auto nb_gauss_pts = getGaussPts().size2();
399 jacPtr->resize(4, nb_gauss_pts,
false);
400 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
405 auto t_t1 = getFTensor1Tangent1AtGaussPts();
406 auto t_t2 = getFTensor1Tangent2AtGaussPts();
410 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
411 t_jac(
i, N0) = t_t1(
i);
412 t_jac(
i, N1) = t_t2(
i);
430 size_t nb_gauss_pts = getGaussPts().size2();
431 jacPtr->resize(9, nb_gauss_pts,
false);
433 auto t_t1 = getFTensor1Tangent1AtGaussPts();
434 auto t_t2 = getFTensor1Tangent2AtGaussPts();
435 auto t_normal = getFTensor1NormalsAtGaussPts();
441 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
442 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
444 t_jac(
i, N0) = t_t1(
i);
445 t_jac(
i, N1) = t_t2(
i);
447 const double l = sqrt(t_normal(
j) * t_normal(
j));
449 t_jac(
i, N2) = t_normal(
i) /
l;
464 if (moab::CN::Dimension(
type) != 2)
467 auto get_normals_ptr = [&]() {
474 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
475 auto nb_base_functions) {
478 auto ptr = get_normals_ptr();
480 &ptr[0], &ptr[1], &ptr[2]);
483 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
484 const auto l2 = t_normal(
i) * t_normal(
i);
485 for (
int bb = 0; bb != nb_base_functions; ++bb) {
486 const auto v = t_base(0);
487 t_base(
i) = (
v / l2) * t_normal(
i);
498 const auto &base_functions = data.
getN(base);
499 const auto nb_gauss_pts = base_functions.size1();
503 const auto nb_base_functions = base_functions.size2() / 3;
504 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
516 const auto type_dim = moab::CN::Dimension(
type);
517 if (type_dim != 1 && type_dim != 2)
530 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
532 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
534 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
540 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
542 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
544 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
552 auto &baseN = data.
getN(base);
553 auto &diffBaseN = data.
getDiffN(base);
555 int nb_dofs = baseN.size2() / 3;
556 int nb_gauss_pts = baseN.size1();
558 piolaN.resize(baseN.size1(), baseN.size2());
559 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
561 if (nb_dofs > 0 && nb_gauss_pts > 0) {
570 t_transformed_diff_h_curl(
581 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
584 for (
int ll = 0; ll != nb_dofs; ll++) {
585 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
586 t_transformed_diff_h_curl(
i,
k) = t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
588 ++t_transformed_h_curl;
590 ++t_transformed_diff_h_curl;
597 if (cc != nb_gauss_pts * nb_dofs)
622 "Wrong number of integration points %d!=%d",
625 auto get_base_at_pts = [&](
auto base) {
632 auto get_tangent_ptr = [&]() {
640 auto get_tangent_at_pts = [&]() {
641 auto ptr = get_tangent_ptr();
643 &ptr[0], &ptr[1], &ptr[2]);
647 auto calculate_squared_edge_length = [&]() {
648 std::vector<double> l1;
650 l1.resize(nb_gauss_pts);
651 auto t_m_at_pts = get_tangent_at_pts();
652 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
653 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
660 auto l1 = calculate_squared_edge_length();
665 const auto nb_dofs = data.
getN(base).size2() / 3;
667 if (nb_gauss_pts && nb_dofs) {
668 auto t_h_curl = get_base_at_pts(base);
670 auto t_m_at_pts = get_tangent_at_pts();
671 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
672 const double l0 = l1[gg];
673 for (
int ll = 0; ll != nb_dofs; ++ll) {
674 const double val = t_h_curl(0);
675 const double a = val / l0;
676 t_h_curl(
i) = t_m_at_pts(
i) *
a;
684 if (cc != nb_gauss_pts * nb_dofs)
694 const FieldSpace space, boost::shared_ptr<VectorDouble> det_jac_ptr)
695 :
OP(space), fieldSpace(space), detJacPtr(det_jac_ptr) {
710 auto &base_fun = data.
getN(base);
711 if (base_fun.size2()) {
714 const auto nb_int_pts = base_fun.size1();
716 if (nb_int_pts != det_vec.size())
718 "Number of integration pts in detJacPtr does not mush "
719 "number of integration pts in base function");
721 for (
auto gg = 0; gg != nb_int_pts; ++gg) {
722 auto row = ublas::matrix_row<MatrixDouble>(base_fun, gg);
726 auto &diff_base_fun = data.
getDiffN(base);
727 if (diff_base_fun.size2()) {
728 for (
auto gg = 0; gg != nb_int_pts; ++gg) {
729 auto row = ublas::matrix_row<MatrixDouble>(diff_base_fun, gg);
796 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
797 std::vector<FieldSpace> spaces, std::string geom_field_name,
798 boost::shared_ptr<MatrixDouble> jac_ptr,
799 boost::shared_ptr<VectorDouble> det_ptr,
800 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
804 jac_ptr = boost::make_shared<MatrixDouble>();
806 det_ptr = boost::make_shared<VectorDouble>();
808 inv_jac_ptr = boost::make_shared<MatrixDouble>();
810 if (geom_field_name.empty()) {
825 for (
auto s : spaces) {
852 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
853 std::vector<FieldSpace> spaces, std::string geom_field_name,
854 boost::shared_ptr<MatrixDouble> jac_ptr,
855 boost::shared_ptr<VectorDouble> det_ptr,
856 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
860 jac_ptr = boost::make_shared<MatrixDouble>();
862 det_ptr = boost::make_shared<VectorDouble>();
864 inv_jac_ptr = boost::make_shared<MatrixDouble>();
866 if (geom_field_name.empty()) {
878 for (
auto s : spaces) {
908 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
909 std::vector<FieldSpace> spaces, std::string geom_field_name) {
912 if (geom_field_name.empty()) {
920 for (
auto s : spaces) {
940 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
941 std::vector<FieldSpace> spaces, std::string geom_field_name) {
944 if (geom_field_name.empty()) {
952 for (
auto s : spaces) {
967 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
968 std::vector<FieldSpace> spaces, std::string geom_field_name) {
971 if (geom_field_name.empty()) {
978 for (
auto s : spaces) {
1001 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1002 std::vector<FieldSpace> spaces, std::string geom_field_name,
1003 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
1004 boost::shared_ptr<MatrixDouble> inv_jac) {
1008 jac = boost::make_shared<MatrixDouble>();
1010 det = boost::make_shared<VectorDouble>();
1012 inv_jac = boost::make_shared<MatrixDouble>();
1014 if (geom_field_name.empty()) {
1028 for (
auto s : spaces) {