v0.14.0
HODataOperators.cpp
Go to the documentation of this file.
1 /** \file HODataOperators.cpp
2 
3 \brief Set of operators for high-order geometry approximation.
4 
5 */
6 
7 namespace MoFEM {
8 
10  boost::shared_ptr<MatrixDouble> jac_ptr)
12  jacPtr(jac_ptr) {
13 
14  for (auto t = MBEDGE; t != MBMAXTYPE; ++t)
15  doEntities[t] = false;
16 
17  if (!jacPtr)
18  THROW_MESSAGE("Jac pointer not allocated");
19 }
20 
25 
26  const auto nb_base_functions = data.getN(NOBASE).size2();
27  if (nb_base_functions) {
28 
29  const auto nb_gauss_pts = getGaussPts().size2();
30  auto t_diff_base = data.getFTensor1DiffN<3>(NOBASE);
31  auto &coords = getCoords();
32 
33 #ifndef NDEBUG
34  if (nb_gauss_pts != data.getDiffN().size1())
35  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
36  "Inconsistent number base functions and gauss points");
37  if (nb_base_functions != data.getDiffN().size2() / 3)
38  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
39  "Inconsistent number of base functions");
40  if (coords.size() != 3 * nb_base_functions)
41  SETERRQ2(
42  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
43  "Number of vertex coordinates and base functions is inconsistent "
44  "%d != %d",
45  coords.size() / 3, nb_base_functions);
46 #endif
47 
48  jacPtr->resize(9, nb_gauss_pts, false);
49  jacPtr->clear();
50  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
51  auto t_vol_inv_jac = getInvJac();
52 
56 
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));
62  ++t_diff_base;
63  ++t_coords;
64  }
65  ++t_jac;
66  }
67  }
68 
70 }
71 
76 
77  if (getFEDim() == 3) {
78 
79  auto transform_base = [&](MatrixDouble &diff_n) {
81  return applyTransform<3, 3, 3, 3>(diff_n);
83  };
84 
85  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
86  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
87  CHKERR transform_base(data.getDiffN(base));
88  }
89 
90  switch (type) {
91  case MBVERTEX:
92  for (auto &m : data.getBBDiffNMap())
93  CHKERR transform_base(*(m.second));
94  break;
95  default:
96  for (auto &ptr : data.getBBDiffNByOrderArray())
97  if (ptr)
98  CHKERR transform_base(*ptr);
99  }
100 
101  } else {
102  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Use different operator");
103  }
104 
106 }
107 
112 
116 
117  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
118 
119  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
120 
121  diffHdivInvJac.resize(data.getDiffN(base).size1(),
122  data.getDiffN(base).size2(), false);
123 
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)
127  continue;
128 
129  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
130  double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
132  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
133  &t_inv_diff_n_ptr[HVEC0_2], &t_inv_diff_n_ptr[HVEC1_0],
134  &t_inv_diff_n_ptr[HVEC1_1], &t_inv_diff_n_ptr[HVEC1_2],
135  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
136  &t_inv_diff_n_ptr[HVEC2_2]);
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);
141  ++t_diff_n;
142  ++t_inv_diff_n;
143  }
144  ++t_inv_jac;
145  }
146 
147  data.getDiffN(base).swap(diffHdivInvJac);
148  }
149 
151 }
152 
156  const size_t nb_int_pts = getGaussPts().size2();
157  if (getNormalsAtGaussPts().size1()) {
158  if (getNormalsAtGaussPts().size1() == nb_int_pts) {
159  double a = getMeasure();
160  if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI)
161  a *= 2;
162  auto t_w = getFTensor0IntegrationWeight();
163  auto t_normal = getFTensor1NormalsAtGaussPts();
165  for (size_t gg = 0; gg != nb_int_pts; ++gg) {
166  t_w *= sqrt(t_normal(i) * t_normal(i)) / a;
167  ++t_w;
168  ++t_normal;
169  }
170  } else {
171  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
172  "Number of rows in getNormalsAtGaussPts should be equal to "
173  "number of integration points, but is not, i.e. %d != %d",
174  getNormalsAtGaussPts().size1(), nb_int_pts);
175  }
176  }
178 }
179 
183  const size_t nb_int_pts = getGaussPts().size2();
184  if (getTangentAtGaussPts().size1()) {
185  if (getTangentAtGaussPts().size1() == nb_int_pts) {
186  double a = getMeasure();
187  auto t_w = getFTensor0IntegrationWeight();
188  auto t_tangent = getFTensor1TangentAtGaussPts();
190  for (size_t gg = 0; gg != nb_int_pts; ++gg) {
191  t_w *= sqrt(t_tangent(i) * t_tangent(i)) / a;
192  ++t_w;
193  ++t_tangent;
194  }
195  } else {
196  SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
197  "Number of rows in getTangentAtGaussPts should be equal to "
198  "number of integration points, but is not, i.e. %d != %d",
199  getTangentAtGaussPts().size1(), nb_int_pts);
200  }
201  }
203 }
204 
208 
209  const auto nb_integration_pts = detPtr->size();
210 
211 #ifndef NDEBUG
212  if (nb_integration_pts != getGaussPts().size2())
213  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
214  "Inconsistent number of data points");
215 #endif
216 
217  auto t_w = getFTensor0IntegrationWeight();
218  auto t_det = getFTensor0FromVec(*detPtr);
219  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
220  t_w *= t_det;
221  ++t_w;
222  ++t_det;
223  }
224 
226 }
227 
232 
236 
237 #ifndef NDEBUG
238  if (!detPtr)
239  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
240  "Pointer for detPtr not allocated");
241  if (!jacPtr)
242  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
243  "Pointer for jacPtr not allocated");
244 #endif
245 
246  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
247 
248  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
249 
250  auto nb_gauss_pts = data.getNSharedPtr(base) ? data.getN(base).size1() : 0;
251  auto nb_base_functions =
252  data.getNSharedPtr(base) ? data.getN(base).size2() / 3 : 0;
253  auto nb_diff_base_functions =
254  data.getDiffNSharedPtr(base) ? data.getDiffN(base).size2() / 9 : 0;
255 
256 #ifndef NDEBUG
257  if (nb_diff_base_functions) {
258  if (nb_diff_base_functions != nb_base_functions)
259  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
260  "Wrong number base functions %d != %d", nb_diff_base_functions,
261  nb_base_functions);
262  if (data.getDiffN(base).size1() != nb_gauss_pts)
263  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
264  "Wrong number integration points");
265  }
266 #endif
267 
268  if (nb_gauss_pts && nb_base_functions) {
269 
270  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
271 
272  auto t_n = data.getFTensor1N<3>(base);
273  double *t_transformed_n_ptr = &*piolaN.data().begin();
275  t_transformed_n_ptr, // HVEC0
276  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
277 
278  auto t_det = getFTensor0FromVec(*detPtr);
279  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
280 
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));
285  ++t_n;
286  ++t_transformed_n;
287  }
288  ++t_det;
289  ++t_jac;
290  }
291 
292  data.getN(base).swap(piolaN);
293  }
294 
295  if (nb_gauss_pts && nb_diff_base_functions) {
296 
297  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
298 
299  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
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]);
311 
312  auto t_det = getFTensor0FromVec(*detPtr);
313  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
314  for (unsigned int bb = 0; bb != nb_diff_base_functions; ++bb) {
315  const double a = 1. / t_det;
316  t_transformed_diff_n(i, k) = a * t_diff_n(i, k);
317  ++t_diff_n;
318  ++t_transformed_diff_n;
319  }
320  ++t_det;
321  }
322 
323  data.getDiffN(base).swap(piolaDiffN);
324  }
325  }
326 
328 }
329 
334 
338 
339  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
340 
341  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
342 
343  unsigned int nb_gauss_pts = data.getN(base).size1();
344  unsigned int nb_base_functions = data.getN(base).size2() / 3;
345  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
346  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
347 
348  auto t_n = data.getFTensor1N<3>(base);
349  double *t_transformed_n_ptr = &*piolaN.data().begin();
351  t_transformed_n_ptr, // HVEC0
352  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
353  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
354  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
355  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
356  t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
357  &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
358  &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
359  &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
360  &t_transformed_diff_n_ptr[HVEC2_2]);
361 
362  auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
363 
364  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
365  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
366  t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
367  t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
368  ++t_n;
369  ++t_transformed_n;
370  ++t_diff_n;
371  ++t_transformed_diff_n;
372  }
373  ++t_inv_jac;
374  }
375 
376  data.getN(base).swap(piolaN);
377  data.getDiffN(base).swap(piolaDiffN);
378  }
379 
381 }
382 
384  boost::shared_ptr<MatrixDouble> jac_ptr)
386  jacPtr(jac_ptr) {}
387 
391 
393 
394  const auto nb_gauss_pts = getGaussPts().size2();
395 
396  jacPtr->resize(4, nb_gauss_pts, false);
397  auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
398 
401 
402  auto t_t1 = getFTensor1Tangent1AtGaussPts();
403  auto t_t2 = getFTensor1Tangent2AtGaussPts();
406 
407  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
408  t_jac(i, N0) = t_t1(i);
409  t_jac(i, N1) = t_t2(i);
410  ++t_t1;
411  ++t_t2;
412  ++t_jac;
413  }
414 
416 };
417 
422 
426 
427  size_t nb_gauss_pts = getGaussPts().size2();
428  jacPtr->resize(9, nb_gauss_pts, false);
429 
430  auto t_t1 = getFTensor1Tangent1AtGaussPts();
431  auto t_t2 = getFTensor1Tangent2AtGaussPts();
432  auto t_normal = getFTensor1NormalsAtGaussPts();
433 
437 
438  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
439  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
440 
441  t_jac(i, N0) = t_t1(i);
442  t_jac(i, N1) = t_t2(i);
443 
444  const double l = sqrt(t_normal(j) * t_normal(j));
445 
446  t_jac(i, N2) = t_normal(i) / l;
447 
448  ++t_t1;
449  ++t_t2;
450  ++t_normal;
451  }
452 
454 }
455 
457  int side, EntityType type, EntitiesFieldData::EntData &data) {
460 
461  if (moab::CN::Dimension(type) != 2)
463 
464  auto get_normals_ptr = [&]() {
465  if (normalsAtGaussPts)
466  return &*normalsAtGaussPts->data().begin();
467  else
468  return &*getNormalsAtGaussPts().data().begin();
469  };
470 
471  auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
472  auto nb_base_functions) {
474 
475  auto ptr = get_normals_ptr();
477  &ptr[0], &ptr[1], &ptr[2]);
478 
479  auto t_base = data.getFTensor1N<3>(base);
480  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
481  const auto l2 = t_normal(i) * t_normal(i);
482  for (int bb = 0; bb != nb_base_functions; ++bb) {
483  const auto v = t_base(0);
484  t_base(i) = (v / l2) * t_normal(i);
485  ++t_base;
486  }
487  ++t_normal;
488  }
490  };
491 
492  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
493 
494  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
495  const auto &base_functions = data.getN(base);
496  const auto nb_gauss_pts = base_functions.size1();
497 
498  if (nb_gauss_pts) {
499 
500  const auto nb_base_functions = base_functions.size2() / 3;
501  CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
502  nb_base_functions);
503  }
504  }
505 
507 }
508 
510  int side, EntityType type, EntitiesFieldData::EntData &data) {
512 
513  const auto type_dim = moab::CN::Dimension(type);
514  if (type_dim != 1 && type_dim != 2)
516 
520 
521  auto get_jac = [&]() {
523  double *ptr_n = &*normalsAtPts->data().begin();
524  double *ptr_t1 = &*tangent1AtPts->data().begin();
525  double *ptr_t2 = &*tangent2AtPts->data().begin();
527  &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
528 
529  &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
530 
531  &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
532  } else {
533  double *ptr_n = &*getNormalsAtGaussPts().data().begin();
534  double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
535  double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
537  &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
538 
539  &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
540 
541  &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
542  }
543  };
544 
545  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
546 
547  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
548 
549  auto &baseN = data.getN(base);
550  auto &diffBaseN = data.getDiffN(base);
551 
552  int nb_dofs = baseN.size2() / 3;
553  int nb_gauss_pts = baseN.size1();
554 
555  piolaN.resize(baseN.size1(), baseN.size2());
556  diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
557 
558  if (nb_dofs > 0 && nb_gauss_pts > 0) {
559 
560  auto t_h_curl = data.getFTensor1N<3>(base);
561  auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
562 
563  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
564  &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
565 
567  t_transformed_diff_h_curl(
570  &diffPiolaN(0, HVEC2_0), &diffPiolaN(0, HVEC2_1));
571 
572  int cc = 0;
573  double det;
575 
576  // HO geometry is set, so jacobian is different at each gauss point
577  auto t_m_at_pts = get_jac();
578  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
579  CHKERR determinantTensor3by3(t_m_at_pts, det);
580  CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
581  for (int ll = 0; ll != nb_dofs; ll++) {
582  t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
583  t_transformed_diff_h_curl(i, k) = t_inv_m(j, i) * t_diff_h_curl(j, k);
584  ++t_h_curl;
585  ++t_transformed_h_curl;
586  ++t_diff_h_curl;
587  ++t_transformed_diff_h_curl;
588  ++cc;
589  }
590  ++t_m_at_pts;
591  }
592 
593 #ifndef NDEBUG
594  if (cc != nb_gauss_pts * nb_dofs)
595  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
596 #endif
597 
598  baseN.swap(piolaN);
599  diffBaseN.swap(diffPiolaN);
600  }
601  }
602 
604 }
605 
607  int side, EntityType type, EntitiesFieldData::EntData &data) {
610 
611  if (type != MBEDGE)
613 
614  const auto nb_gauss_pts = getGaussPts().size2();
615 
616  if (tangentAtGaussPts)
617  if (tangentAtGaussPts->size1() != nb_gauss_pts)
618  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
619  "Wrong number of integration points %d!=%d",
620  tangentAtGaussPts->size1(), nb_gauss_pts);
621 
622  auto get_base_at_pts = [&](auto base) {
624  &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
625  &data.getN(base)(0, HVEC2));
626  return t_h_curl;
627  };
628 
629  auto get_tangent_ptr = [&]() {
630  if (tangentAtGaussPts) {
631  return &*tangentAtGaussPts->data().begin();
632  } else {
633  return &*getTangentAtGaussPts().data().begin();
634  }
635  };
636 
637  auto get_tangent_at_pts = [&]() {
638  auto ptr = get_tangent_ptr();
640  &ptr[0], &ptr[1], &ptr[2]);
641  return t_m_at_pts;
642  };
643 
644  auto calculate_squared_edge_length = [&]() {
645  std::vector<double> l1;
646  if (nb_gauss_pts) {
647  l1.resize(nb_gauss_pts);
648  auto t_m_at_pts = get_tangent_at_pts();
649  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
650  l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
651  ++t_m_at_pts;
652  }
653  }
654  return l1;
655  };
656 
657  auto l1 = calculate_squared_edge_length();
658 
659  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
660 
661  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
662  const auto nb_dofs = data.getN(base).size2() / 3;
663 
664  if (nb_gauss_pts && nb_dofs) {
665  auto t_h_curl = get_base_at_pts(base);
666  int cc = 0;
667  auto t_m_at_pts = get_tangent_at_pts();
668  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
669  const double l0 = l1[gg];
670  for (int ll = 0; ll != nb_dofs; ++ll) {
671  const double val = t_h_curl(0);
672  const double a = val / l0;
673  t_h_curl(i) = t_m_at_pts(i) * a;
674  ++t_h_curl;
675  ++cc;
676  }
677  ++t_m_at_pts;
678  }
679 
680 #ifndef NDEBUG
681  if (cc != nb_gauss_pts * nb_dofs)
682  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
683 #endif
684  }
685  }
686 
688 }
689 
691  const FieldSpace space, const FieldApproximationBase base,
692  boost::shared_ptr<VectorDouble> det_jac_ptr,
693  boost::shared_ptr<double> scale_det_ptr)
694  : OP(space), fieldSpace(space), fieldBase(base), detJacPtr(det_jac_ptr),
695  scaleDetPtr(scale_det_ptr) {
696  if (!detJacPtr) {
697  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
698  }
699 }
700 
705 
706  double scale_det = 1.0;
707  if (scaleDetPtr)
708  scale_det = *scaleDetPtr;
709 
710  auto scale = [&](auto base) {
712  auto &base_fun = data.getN(base);
713  if (base_fun.size2()) {
714 
715  auto &det_vec = *detJacPtr;
716  const auto nb_int_pts = base_fun.size1();
717 
718  if (nb_int_pts != det_vec.size())
719  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
720  "Number of integration pts in detJacPtr does not mush "
721  "number of integration pts in base function");
722 
723  for (auto gg = 0; gg != nb_int_pts; ++gg) {
724  auto row = ublas::matrix_row<MatrixDouble>(base_fun, gg);
725  row *= scale_det / det_vec[gg];
726  }
727 
728  auto &diff_base_fun = data.getDiffN(base);
729  if (diff_base_fun.size2()) {
730  for (auto gg = 0; gg != nb_int_pts; ++gg) {
731  auto row = ublas::matrix_row<MatrixDouble>(diff_base_fun, gg);
732  row *= scale_det / det_vec[gg];
733  }
734  }
735  }
737  };
738 
739  if (this->getFEDim() == 3) {
740  switch (fieldSpace) {
741  case L2:
742  if (type >= MBTET)
744  break;
745  default:
746  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
747  }
748  } else if (this->getFEDim() == 2) {
749  switch (fieldSpace) {
750  case L2:
751  if (type >= MBTRI)
753  break;
754  default:
755  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
756  }
757  } else if (this->getFEDim() == 1) {
758  switch (fieldSpace) {
759  case L2:
760  if (type >= MBEDGE)
762  break;
763  default:
764  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
765  }
766  } else {
767  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "dimension not handled");
768  }
769 
771 }
772 
774  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
775  std::vector<FieldSpace> spaces, std::string geom_field_name,
776  boost::shared_ptr<MatrixDouble> jac_ptr,
777  boost::shared_ptr<VectorDouble> det_ptr,
778  boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
780 
781  if (!jac_ptr)
782  jac_ptr = boost::make_shared<MatrixDouble>();
783  if (!det_ptr)
784  det_ptr = boost::make_shared<VectorDouble>();
785  if (!inv_jac_ptr)
786  inv_jac_ptr = boost::make_shared<MatrixDouble>();
787 
788  if (geom_field_name.empty()) {
789 
790  pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
791 
792  } else {
793 
794  pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
795  pipeline.push_back(
796  new OpCalculateVectorFieldGradient<2, 2>(geom_field_name, jac_ptr));
797  pipeline.push_back(new OpGetHONormalsOnFace<2>(geom_field_name));
798  }
799 
800  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
801  pipeline.push_back(new OpSetHOWeightsOnFace());
802 
803  for (auto s : spaces) {
804  switch (s) {
805  case NOSPACE:
806  break;
807  case H1:
808  case L2:
809  pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(s, inv_jac_ptr));
810  break;
811  case HCURL:
812  pipeline.push_back(new OpSetCovariantPiolaTransformOnFace2D(inv_jac_ptr));
813  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
814  break;
815  case HDIV:
816  pipeline.push_back(new OpMakeHdivFromHcurl());
817  pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
818  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
819  break;
820  default:
821  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
822  "Space %s not yet implemented", FieldSpaceNames[s]);
823  }
824  }
825 
827 }
828 
830  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
831  std::vector<FieldSpace> spaces, std::string geom_field_name,
832  boost::shared_ptr<MatrixDouble> jac_ptr,
833  boost::shared_ptr<VectorDouble> det_ptr,
834  boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
836 
837  if (!jac_ptr)
838  jac_ptr = boost::make_shared<MatrixDouble>();
839  if (!det_ptr)
840  det_ptr = boost::make_shared<VectorDouble>();
841  if (!inv_jac_ptr)
842  inv_jac_ptr = boost::make_shared<MatrixDouble>();
843 
844  if (geom_field_name.empty()) {
845 
846  } else {
847 
848  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
849  pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
850  }
851 
852  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
853  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
854  pipeline.push_back(new OpSetHOWeightsOnFace());
855 
856  for (auto s : spaces) {
857  switch (s) {
858  case NOSPACE:
859  break;
860  case H1:
861  pipeline.push_back(
862  new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
863  break;
864  case HDIV:
865  pipeline.push_back(new OpMakeHdivFromHcurl());
866  pipeline.push_back(
868  jac_ptr));
869  pipeline.push_back(
870  new OpSetInvJacHcurlFaceEmbeddedIn3DSpace(inv_jac_ptr));
871  break;
872  case L2:
873  pipeline.push_back(
874  new OpSetInvJacL2ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
875  break;
876  default:
877  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
878  "Space %s not yet implemented", FieldSpaceNames[s]);
879  }
880  }
881 
883 }
884 
886  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
887  std::vector<FieldSpace> spaces, std::string geom_field_name) {
889 
890  if (geom_field_name.empty()) {
891 
892  } else {
893 
894  pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
895  pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
896  }
897 
898  for (auto s : spaces) {
899  switch (s) {
900  case NOSPACE:
901  break;
902  case HCURL:
903  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
904  break;
905  case HDIV:
906  pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
907  break;
908  default:
909  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
910  "Space %s not yet implemented", FieldSpaceNames[s]);
911  }
912  }
913 
915 }
916 
918  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
919  std::vector<FieldSpace> spaces, std::string geom_field_name) {
921 
922  if (geom_field_name.empty()) {
923 
924  } else {
925 
926  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
927  pipeline.push_back(new OpGetHOTangentsOnEdge<3>(geom_field_name));
928  }
929 
930  for (auto s : spaces) {
931  switch (s) {
932  case HCURL:
933  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
934  break;
935  default:
936  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
937  "Space %s not yet implemented", FieldSpaceNames[s]);
938  }
939  }
940 
942 }
943 
945  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
946  std::vector<FieldSpace> spaces, std::string geom_field_name) {
948 
949  if (geom_field_name.empty()) {
950  } else {
951 
952  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
953  pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
954  }
955 
956  for (auto s : spaces) {
957  switch (s) {
958  case NOSPACE:
959  break;
960  case HCURL:
961  pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
962  break;
963  case HDIV:
964  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
965  break;
966  case L2:
967  break;
968  default:
969  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
970  "Space %s not yet implemented", FieldSpaceNames[s]);
971  break;
972  }
973  }
974 
976 }
977 
979  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
980  std::vector<FieldSpace> spaces, std::string geom_field_name,
981  boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
982  boost::shared_ptr<MatrixDouble> inv_jac) {
984 
985  if (!jac)
986  jac = boost::make_shared<MatrixDouble>();
987  if (!det)
988  det = boost::make_shared<VectorDouble>();
989  if (!inv_jac)
990  inv_jac = boost::make_shared<MatrixDouble>();
991 
992  if (geom_field_name.empty()) {
993 
994  pipeline.push_back(new OpCalculateHOJac<3>(jac));
995 
996  } else {
997 
998  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
999  pipeline.push_back(
1000  new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
1001  }
1002 
1003  pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
1004  pipeline.push_back(new OpSetHOWeights(det));
1005 
1006  for (auto s : spaces) {
1007  switch (s) {
1008  case NOSPACE:
1009  break;
1010  case H1:
1011  pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
1012  break;
1013  case HCURL:
1014  pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
1015  pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
1016  break;
1017  case HDIV:
1018  pipeline.push_back(
1019  new OpSetHOContravariantPiolaTransform(HDIV, det, jac));
1020  break;
1021  case L2:
1022  pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
1023  break;
1024  default:
1025  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1026  "Space %s not yet implemented", FieldSpaceNames[s]);
1027  break;
1028  }
1029  }
1030 
1032 }
1033 
1034 } // namespace MoFEM
UBlasMatrix< double >
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::OpSetHOCovariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: HODataOperators.hpp:227
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::OpSetHOInvJacToScalarBasesImpl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:73
MoFEM::OpSetHOContravariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: HODataOperators.hpp:201
MoFEM::OpSetHOWeights::detPtr
boost::shared_ptr< VectorDouble > detPtr
Definition: HODataOperators.hpp:173
MoFEM::OpGetHOTangentsOnEdge
Calculate tangent vector on edge form HO geometry approximation.
Definition: HODataOperators.hpp:369
MoFEM::EntitiesFieldData::EntData::getBBDiffNMap
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of derivatives base function for BB base, key is a field name
Definition: EntitiesFieldData.cpp:869
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:3413
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::OpSetHOContravariantPiolaTransform::jacPtr
boost::shared_ptr< MatrixDouble > jacPtr
Definition: HODataOperators.hpp:198
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::OpCalculateHOJacForFaceEmbeddedIn3DSpace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
Definition: HODataOperators.hpp:265
HVEC0_2
@ HVEC0_2
Definition: definitions.h:211
MoFEM::OpSetContravariantPiolaTransformOnFace2D
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
Definition: UserDataOperators.hpp:3511
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::normalsAtPts
boost::shared_ptr< MatrixDouble > normalsAtPts
Definition: HODataOperators.hpp:357
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::OpHOSetContravariantPiolaTransformOnEdge3D::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:606
MoFEM::OpSetHOContravariantPiolaTransform
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
Definition: HODataOperators.hpp:180
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:702
MoFEM::OpSetHOCovariantPiolaTransform::piolaDiffN
MatrixDouble piolaDiffN
Definition: HODataOperators.hpp:228
HVEC0_1
@ HVEC0_1
Definition: definitions.h:208
MoFEM::OpHOSetContravariantPiolaTransformOnFace3D
transform Hdiv base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:298
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::OpSetHOContravariantPiolaTransform::piolaN
MatrixDouble piolaN
Definition: HODataOperators.hpp:200
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getCoords
VectorDouble & getCoords()
nodal coordinates
Definition: VolumeElementForcesAndSourcesCore.hpp:181
MoFEM::OpCalculateHOJacForVolume::OpCalculateHOJacForVolume
OpCalculateHOJacForVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
Definition: HODataOperators.cpp:9
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
MoFEM::OpSetCovariantPiolaTransformOnFace2D
OpSetCovariantPiolaTransformOnFace2DImpl< 2 > OpSetCovariantPiolaTransformOnFace2D
Definition: UserDataOperators.hpp:3463
MoFEM::EntitiesFieldData::EntData::getBBDiffNByOrderArray
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
Definition: EntitiesFieldData.cpp:902
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::detJacPtr
boost::shared_ptr< VectorDouble > detJacPtr
Definition: HODataOperators.hpp:406
HVEC0_0
@ HVEC0_0
Definition: definitions.h:205
scale
double scale
Definition: plastic.cpp:123
HVEC1_1
@ HVEC1_1
Definition: definitions.h:209
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MoFEM::OpCalculateHOCoords
Calculate HO coordinates at gauss points.
Definition: HODataOperators.hpp:41
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
MoFEM::OpCalculateHOJacForFaceImpl
Calculate jacobian for face element.
Definition: HODataOperators.hpp:240
HVEC1
@ HVEC1
Definition: definitions.h:199
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1241
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1276
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1588
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator::getTangentAtGaussPts
MatrixDouble & getTangentAtGaussPts()
get tangent vector to edge curve at integration points
Definition: EdgeElementForcesAndSourcesCore.hpp:204
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1237
MoFEM::OpCalculateHOJacForVolume::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:22
HVEC2_1
@ HVEC2_1
Definition: definitions.h:210
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
MoFEM::OpSetHOContravariantPiolaTransform::detPtr
boost::shared_ptr< VectorDouble > detPtr
Definition: HODataOperators.hpp:197
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPts
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:290
MoFEM::OpSetInvJacL2ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3375
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getTangent2AtGaussPts
MatrixDouble & getTangent2AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:310
HVEC2_2
@ HVEC2_2
Definition: definitions.h:213
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getTangent1AtGaussPts
MatrixDouble & getTangent1AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:304
HVEC1_2
@ HVEC1_2
Definition: definitions.h:212
MoFEM::OpSetContravariantPiolaTransformOnFace2DEmbeddedIn3DSpace
OpSetContravariantPiolaTransformOnFace2DImpl< 3 > OpSetContravariantPiolaTransformOnFace2DEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3513
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::OpSetHOCovariantPiolaTransform
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
Definition: HODataOperators.hpp:206
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::fieldSpace
FieldSpace fieldSpace
Definition: HODataOperators.hpp:404
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getInvJac
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian
Definition: VolumeElementForcesAndSourcesCore.hpp:176
convert.type
type
Definition: convert.py:64
MoFEM::OpSetHOWeights::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:205
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::piolaN
MatrixDouble piolaN
Definition: HODataOperators.hpp:361
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpSetHOContravariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:229
MoFEM::OpSetContravariantPiolaTransformOnEdge2D
Definition: UserDataOperators.hpp:3521
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1001
MoFEM::OpSetHOCovariantPiolaTransform::jacInvPtr
boost::shared_ptr< MatrixDouble > jacInvPtr
Definition: HODataOperators.hpp:225
MoFEM::OpSetHOInvJacVectorBase::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:109
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::EntitiesFieldData::EntData::getDiffNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
Definition: EntitiesFieldData.cpp:37
MoFEM::OpCalculateHOJac< 3 >
Definition: HODataOperators.hpp:269
get_jac
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
Definition: NonLinearElasticElement.cpp:731
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::diffPiolaN
MatrixDouble diffPiolaN
Definition: HODataOperators.hpp:362
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::tangent1AtPts
boost::shared_ptr< MatrixDouble > tangent1AtPts
Definition: HODataOperators.hpp:358
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpSetHOInvJacVectorBase::invJacPtr
boost::shared_ptr< MatrixDouble > invJacPtr
Definition: HODataOperators.hpp:111
FTensor::Index< 'i', 3 >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::OpCalculateHOJacForVolume::jacPtr
boost::shared_ptr< MatrixDouble > jacPtr
Definition: HODataOperators.hpp:29
MoFEM::OpSetHOWeightsOnEdge::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:180
MoFEM::DataOperator::doEntities
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Definition: DataOperators.hpp:85
FieldSpaceNames
const static char *const FieldSpaceNames[]
Definition: definitions.h:92
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D
transform Hcurl base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:336
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::OpScaleBaseBySpaceInverseOfMeasure
OpScaleBaseBySpaceInverseOfMeasure(const FieldSpace space, const FieldApproximationBase base, boost::shared_ptr< VectorDouble > det_jac_ptr, boost::shared_ptr< double > scale_det_ptr=nullptr)
Definition: HODataOperators.cpp:690
MoFEM::OpHOSetContravariantPiolaTransformOnEdge3D::tangentAtGaussPts
boost::shared_ptr< MatrixDouble > tangentAtGaussPts
Definition: HODataOperators.hpp:330
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::tangent2AtPts
boost::shared_ptr< MatrixDouble > tangent2AtPts
Definition: HODataOperators.hpp:359
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
MoFEM::OpSetHOCovariantPiolaTransform::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:331
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEDim
int getFEDim() const
Get dimension of finite element.
Definition: ForcesAndSourcesCore.hpp:1009
MoFEM::OpSetHOInvJacVectorBase::diffHdivInvJac
MatrixDouble diffHdivInvJac
Definition: HODataOperators.hpp:112
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:159
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3684
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
HVEC2_0
@ HVEC2_0
Definition: definitions.h:207
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getFTensor1NormalsAtGaussPts
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
Definition: FaceElementForcesAndSourcesCore.hpp:316
MoFEM::OpHOSetContravariantPiolaTransformOnFace3D::normalsAtGaussPts
boost::shared_ptr< MatrixDouble > normalsAtGaussPts
Definition: HODataOperators.hpp:311
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::OpSetHOInvJacVectorBase
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Definition: HODataOperators.hpp:91
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::OpMakeHdivFromHcurl
Make Hdiv space from Hcurl space in 2d.
Definition: UserDataOperators.hpp:3420
HVEC2
@ HVEC2
Definition: definitions.h:199
MoFEM::OpSetInvJacHcurlFaceEmbeddedIn3DSpace
OpSetInvJacHcurlFaceImpl< 3 > OpSetInvJacHcurlFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3414
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::OpHOSetContravariantPiolaTransformOnFace3D::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:456
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3367
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::OpSetHOWeightsOnFace::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:153
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator::getFTensor1TangentAtGaussPts
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, DIM > getFTensor1TangentAtGaussPts()
Get tangent vector at integration points.
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::fieldBase
FieldApproximationBase fieldBase
Definition: HODataOperators.hpp:405
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HODataOperators.cpp:509
HVEC1_0
@ HVEC1_0
Definition: definitions.h:206
HVEC0
@ HVEC0
Definition: definitions.h:199
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::scaleDetPtr
boost::shared_ptr< double > scaleDetPtr
Definition: HODataOperators.hpp:407
MoFEM::OpHOSetContravariantPiolaTransformOnEdge3D
transform Hcurl base fluxes from reference element to physical edge
Definition: HODataOperators.hpp:317
MoFEM::EntitiesFieldData::EntData::getNSharedPtr
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
Definition: EntitiesFieldData.cpp:27