11 FieldSpace space, boost::shared_ptr<MatrixDouble> inv_jac_ptr)
27 "This operator can be used only with element which faces");
30 return applyTransform<2, 2, 2, 2>(diff_n);
33 if (!(
type == MBVERTEX && sPace ==
L2)) {
43 CHKERR apply_transform(*(
m.second));
48 CHKERR apply_transform(*ptr);
62 "This operator can be used only with element which face");
65 return applyTransform<2, 3, 3, 3>(diff_n);
68 if (!(
type == MBVERTEX && sPace ==
L2)) {
77 CHKERR apply_transform(*(
m.second));
82 CHKERR apply_transform(*ptr);
96 "This operator can be used only with element which faces");
100 size_t nb_functions = diff2_n.size2() / 4;
102 size_t nb_gauss_pts = diff2_n.size1();
105 if (nb_gauss_pts != getGaussPts().size2())
107 "Wrong number of Gauss Pts");
110 diffNinvJac.resize(diff2_n.size1(), diff2_n.size2(),
false);
116 auto t_diff2_n = getFTensor2FromPtr<2, 2>(&*diffNinvJac.data().begin());
117 auto t_diff2_n_ref = getFTensor2FromPtr<2, 2>(&*diff2_n.data().begin());
118 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
119 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
120 for (
size_t dd = 0;
dd != nb_functions; ++
dd) {
122 t_inv_jac(
k,
i) * (t_inv_jac(
l,
j) * t_diff2_n_ref(
k,
l));
128 diff2_n.swap(diffNinvJac);
133 if (!(
type == MBVERTEX && sPace ==
L2)) {
138 data.
getN(base, BaseDerivatives::SecondDerivative));
141 auto error = [&](
auto &
m) {
145 "Operator do not work for Bernstein-Bezier. This "
155 CHKERR apply_transform(*(
m.second));
162 CHKERR apply_transform(*ptr);
180 "This operator can be used only with element which is triangle");
190 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
191 if (nb_base_functions) {
192 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
194 diffHcurlInvJac.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
197 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
199 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
205 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
206 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
207 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
208 t_inv_diff_n(
i,
j) = t_diff_n(
i,
k) * t_inv_jac(
k,
j);
214 data.
getDiffN(base).swap(diffHcurlInvJac);
231 "This operator can be used only with element which is triangle");
241 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 6;
242 if (nb_base_functions) {
243 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
245 diffHcurlInvJac.resize(nb_gauss_pts, nb_base_functions * 9,
false);
248 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
250 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
259 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
260 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
261 for (
unsigned int bb = 0; bb != nb_base_functions; bb++) {
262 t_inv_diff_n(
i,
j) = t_diff_n(
i,
K) * t_inv_jac(
K,
j);
268 data.
getDiffN(base).swap(diffHcurlInvJac);
284 "This operator can be used only with element which is face");
290 const size_t nb_base_functions = data.
getN(base).size2() / 3;
291 if (nb_base_functions) {
293 const size_t nb_gauss_pts = data.
getN(base).size1();
297 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
298 for (
size_t bb = 0; bb != nb_base_functions; ++bb) {
300 const double a = t_n(0);
304 for (
auto n : {0, 1}) {
305 const double b = t_diff_n(0,
n);
306 t_diff_n(0,
n) = -t_diff_n(1,
n);
324 invJacPtr(inv_jac_ptr) {}
330 const auto type_dim = moab::CN::Dimension(
type);
331 if (type_dim != 1 && type_dim != 2)
342 auto &baseN = data.
getN(base);
343 auto &diffBaseN = data.
getDiffN(base);
345 int nb_dofs = baseN.size2() / 3;
346 int nb_gauss_pts = baseN.size1();
348 piolaN.resize(baseN.size1(), baseN.size2());
350 if (nb_dofs > 0 && nb_gauss_pts > 0) {
353 auto t_transformed_h_curl =
354 getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
355 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
356 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
357 for (
int ll = 0; ll != nb_dofs; ll++) {
358 t_transformed_h_curl(
i) = t_inv_jac(
j,
i) * t_h_curl(
j);
360 ++t_transformed_h_curl;
366 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
367 if (diffBaseN.data().size() > 0) {
369 auto t_transformed_diff_h_curl =
371 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
372 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
373 for (
int ll = 0; ll != nb_dofs; ll++) {
374 t_transformed_diff_h_curl(
i,
k) =
375 t_inv_jac(
j,
i) * t_diff_h_curl(
j,
k);
377 ++t_transformed_diff_h_curl;
381 diffBaseN.swap(diffPiolaN);
405 const size_t nb_base_functions = data.
getN(base).size2() / 3;
406 if (nb_base_functions) {
408 const size_t nb_gauss_pts = data.
getN(base).size1();
409 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
410 if (data.
getN(base).size2() > 0) {
412 double *t_transformed_n_ptr = &*piolaN.data().begin();
415 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
416 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
417 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
420 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
421 t_transformed_n(
i) = t_jac(
i,
k) * t_n(
k) / det;
426 data.
getN(base).swap(piolaN);
429 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
430 if (data.
getDiffN(base).data().size() > 0) {
432 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
434 t_transformed_diff_n(t_transformed_diff_n_ptr,
435 &t_transformed_diff_n_ptr[
HVEC0_1],
437 &t_transformed_diff_n_ptr[
HVEC1_0],
438 &t_transformed_diff_n_ptr[
HVEC1_1],
440 &t_transformed_diff_n_ptr[
HVEC2_0],
441 &t_transformed_diff_n_ptr[
HVEC2_1]);
442 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
443 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
446 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
447 t_transformed_diff_n(
i,
k) = t_jac(
i,
j) * t_diff_n(
j,
k) / det;
449 ++t_transformed_diff_n;
452 data.
getDiffN(base).swap(piolaDiffN);
476 const size_t nb_base_functions = data.
getN(base).size2() / 3;
477 if (nb_base_functions) {
479 const size_t nb_gauss_pts = data.
getN(base).size1();
480 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
481 if (data.
getN(base).size2() > 0) {
483 double *t_transformed_n_ptr = &*piolaN.data().begin();
486 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
487 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
488 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
491 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
492 t_transformed_n(
i) = t_jac(
i,
j) * t_n(
j) / det;
497 data.
getN(base).swap(piolaN);
500 piolaDiffN.resize(nb_gauss_pts, data.
getDiffN(base).size2(),
false);
501 if (data.
getDiffN(base).size2() > 0) {
503 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
505 t_transformed_diff_n(t_transformed_diff_n_ptr,
506 &t_transformed_diff_n_ptr[
HVEC0_1],
508 &t_transformed_diff_n_ptr[
HVEC1_0],
509 &t_transformed_diff_n_ptr[
HVEC1_1],
511 &t_transformed_diff_n_ptr[
HVEC2_0],
512 &t_transformed_diff_n_ptr[
HVEC2_1]);
514 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
515 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
518 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
519 t_transformed_diff_n(
i,
K) = t_jac(
i,
j) * t_diff_n(
j,
K) / det;
521 ++t_transformed_diff_n;
524 data.
getDiffN(base).swap(piolaDiffN);
546 size_t nb_gauss_pts = data.
getN(base).size1();
547 size_t nb_dofs = data.
getN(base).size2() / 3;
554 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
556 auto t_dir = getFTensor1TangentAtGaussPts<3>();
558 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
563 const auto l2 = t_normal(
i) * t_normal(
i);
564 for (
int ll = 0; ll != nb_dofs; ++ll) {
565 const double val = t_h_div(0);
566 const double a = val / l2;
567 t_h_div(
i) = t_normal(
i) *
a;
575 if (cc != nb_gauss_pts * nb_dofs)
587 if (
type == MBVERTEX) {
590 double *coords_ptr = &*coords.data().begin();
592 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
600 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
605 for (
int bb = 0; bb != 6; bb++) {
606 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
620 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
637 if (
type == MBVERTEX) {
640 double *coords_ptr = &*coords.data().begin();
642 const int nb_gauss_pts = data.
getN(
NOBASE).size1();
644 invJac.resize(9, nb_gauss_pts,
false);
646 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
652 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
657 for (
int bb = 0; bb != 6; bb++) {
658 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
684 if (data.
getN(base).size2() == 0)
687 const int nb_gauss_pts = data.
getN(base).size1();
697 auto t_inv_jac = getFTensor2FromMat<3, 3>(
invJac);
699 const int nb_dofs = data.
getN(base).size2();
700 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
701 for (
int bb = 0; bb != nb_dofs; bb++) {
702 t_inv_diff_n(
i) = t_diff_n(
j) * t_inv_jac(
j,
i);
721 if (
type == MBVERTEX) {
724 double *coords_ptr = &*coords.data().begin();
727 double j00_f3, j01_f3, j10_f3, j11_f3;
728 for (
int gg = 0; gg < 1; gg++) {
731 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
732 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
733 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
734 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
736 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
759 unsigned int nb_dofs = data.
getN(base).size2();
762 unsigned int nb_gauss_pts = data.
getN(base).size1();
763 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs,
false);
766 if (
type != MBVERTEX) {
767 if (nb_dofs != data.
getDiffN(base).size2() / 2) {
769 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
771 nb_dofs, data.
getDiffN(base).size2());
780 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
781 for (
unsigned int dd = 0;
dd < nb_dofs;
dd++) {
782 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
799 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
800 const EntityType zero_type,
const int zero_side)
803 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
812 const auto nb_integration_points = getGaussPts().size2();
813 if (
type == zeroType && side == zeroSide) {
814 dataPtr->resize(3, nb_integration_points,
false);
823 const auto nb_base_functions = data.
getN().size2() / 3;
825 auto t_data = getFTensor1FromMat<3>(*dataPtr);
826 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
829 for (; bb != nb_dofs; ++bb) {
834 for (; bb < nb_base_functions; ++bb)
843 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
844 const EntityType zero_type,
const int zero_side)
847 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
856 const auto nb_integration_points = getGaussPts().size2();
857 if (
type == zeroType && side == zeroSide) {
858 dataPtr->resize(2, nb_integration_points,
false);
868 const auto nb_base_functions = data.
getN().size2();
870 auto t_data = getFTensor1FromMat<2>(*dataPtr);
871 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
874 for (; bb != nb_dofs; ++bb) {
879 for (; bb < nb_base_functions; ++bb)
888 const std::string
field_name, boost::shared_ptr<MatrixDouble> data_ptr,
889 const EntityType zero_type,
const int zero_side)
892 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
901 const auto nb_integration_points = getGaussPts().size2();
902 if (
type == zeroType && side == zeroSide) {
903 dataPtr->resize(3, nb_integration_points,
false);
914 const auto nb_base_functions = data.
getN().size2();
916 auto t_data = getFTensor1FromMat<3>(*dataPtr);
917 auto t_normal = getFTensor1NormalsAtGaussPts();
918 for (
auto gg = 0; gg != nb_integration_points; ++gg) {
920 auto nrm = t_normal.l2();
922 for (; bb != nb_dofs; ++bb) {
923 t_data(
k) += (t_dof / nrm) *
928 for (; bb < nb_base_functions; ++bb)