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  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
314 
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));
319  ++t_diff_n;
320  ++t_transformed_diff_n;
321  }
322  ++t_det;
323  ++t_jac;
324  }
325 
326  data.getDiffN(base).swap(piolaDiffN);
327  }
328  }
329 
331 }
332 
337 
341 
342  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
343 
344  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
345 
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);
349  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
350 
351  auto t_n = data.getFTensor1N<3>(base);
352  double *t_transformed_n_ptr = &*piolaN.data().begin();
354  t_transformed_n_ptr, // HVEC0
355  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
356  auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
357  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
358  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
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]);
364 
365  auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
366 
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);
371  ++t_n;
372  ++t_transformed_n;
373  ++t_diff_n;
374  ++t_transformed_diff_n;
375  }
376  ++t_inv_jac;
377  }
378 
379  data.getN(base).swap(piolaN);
380  data.getDiffN(base).swap(piolaDiffN);
381  }
382 
384 }
385 
387  boost::shared_ptr<MatrixDouble> jac_ptr)
389  jacPtr(jac_ptr) {}
390 
394 
396 
397  const auto nb_gauss_pts = getGaussPts().size2();
398 
399  jacPtr->resize(4, nb_gauss_pts, false);
400  auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
401 
404 
405  auto t_t1 = getFTensor1Tangent1AtGaussPts();
406  auto t_t2 = getFTensor1Tangent2AtGaussPts();
409 
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);
413  ++t_t1;
414  ++t_t2;
415  ++t_jac;
416  }
417 
419 };
420 
425 
429 
430  size_t nb_gauss_pts = getGaussPts().size2();
431  jacPtr->resize(9, nb_gauss_pts, false);
432 
433  auto t_t1 = getFTensor1Tangent1AtGaussPts();
434  auto t_t2 = getFTensor1Tangent2AtGaussPts();
435  auto t_normal = getFTensor1NormalsAtGaussPts();
436 
440 
441  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
442  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
443 
444  t_jac(i, N0) = t_t1(i);
445  t_jac(i, N1) = t_t2(i);
446 
447  const double l = sqrt(t_normal(j) * t_normal(j));
448 
449  t_jac(i, N2) = t_normal(i) / l;
450 
451  ++t_t1;
452  ++t_t2;
453  ++t_normal;
454  }
455 
457 }
458 
460  int side, EntityType type, EntitiesFieldData::EntData &data) {
463 
464  if (moab::CN::Dimension(type) != 2)
466 
467  auto get_normals_ptr = [&]() {
468  if (normalsAtGaussPts)
469  return &*normalsAtGaussPts->data().begin();
470  else
471  return &*getNormalsAtGaussPts().data().begin();
472  };
473 
474  auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
475  auto nb_base_functions) {
477 
478  auto ptr = get_normals_ptr();
480  &ptr[0], &ptr[1], &ptr[2]);
481 
482  auto t_base = data.getFTensor1N<3>(base);
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);
488  ++t_base;
489  }
490  ++t_normal;
491  }
493  };
494 
495  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
496 
497  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
498  const auto &base_functions = data.getN(base);
499  const auto nb_gauss_pts = base_functions.size1();
500 
501  if (nb_gauss_pts) {
502 
503  const auto nb_base_functions = base_functions.size2() / 3;
504  CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
505  nb_base_functions);
506  }
507  }
508 
510 }
511 
513  int side, EntityType type, EntitiesFieldData::EntData &data) {
515 
516  const auto type_dim = moab::CN::Dimension(type);
517  if (type_dim != 1 && type_dim != 2)
519 
523 
524  auto get_jac = [&]() {
526  double *ptr_n = &*normalsAtPts->data().begin();
527  double *ptr_t1 = &*tangent1AtPts->data().begin();
528  double *ptr_t2 = &*tangent2AtPts->data().begin();
530  &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
531 
532  &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
533 
534  &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
535  } else {
536  double *ptr_n = &*getNormalsAtGaussPts().data().begin();
537  double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
538  double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
540  &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
541 
542  &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
543 
544  &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
545  }
546  };
547 
548  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
549 
550  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
551 
552  auto &baseN = data.getN(base);
553  auto &diffBaseN = data.getDiffN(base);
554 
555  int nb_dofs = baseN.size2() / 3;
556  int nb_gauss_pts = baseN.size1();
557 
558  piolaN.resize(baseN.size1(), baseN.size2());
559  diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
560 
561  if (nb_dofs > 0 && nb_gauss_pts > 0) {
562 
563  auto t_h_curl = data.getFTensor1N<3>(base);
564  auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
565 
566  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
567  &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
568 
570  t_transformed_diff_h_curl(
573  &diffPiolaN(0, HVEC2_0), &diffPiolaN(0, HVEC2_1));
574 
575  int cc = 0;
576  double det;
578 
579  // HO geometry is set, so jacobian is different at each gauss point
580  auto t_m_at_pts = get_jac();
581  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
582  CHKERR determinantTensor3by3(t_m_at_pts, det);
583  CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
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);
587  ++t_h_curl;
588  ++t_transformed_h_curl;
589  ++t_diff_h_curl;
590  ++t_transformed_diff_h_curl;
591  ++cc;
592  }
593  ++t_m_at_pts;
594  }
595 
596 #ifndef NDEBUG
597  if (cc != nb_gauss_pts * nb_dofs)
598  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
599 #endif
600 
601  baseN.swap(piolaN);
602  diffBaseN.swap(diffPiolaN);
603  }
604  }
605 
607 }
608 
610  int side, EntityType type, EntitiesFieldData::EntData &data) {
613 
614  if (type != MBEDGE)
616 
617  const auto nb_gauss_pts = getGaussPts().size2();
618 
619  if (tangentAtGaussPts)
620  if (tangentAtGaussPts->size1() != nb_gauss_pts)
621  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
622  "Wrong number of integration points %d!=%d",
623  tangentAtGaussPts->size1(), nb_gauss_pts);
624 
625  auto get_base_at_pts = [&](auto base) {
627  &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
628  &data.getN(base)(0, HVEC2));
629  return t_h_curl;
630  };
631 
632  auto get_tangent_ptr = [&]() {
633  if (tangentAtGaussPts) {
634  return &*tangentAtGaussPts->data().begin();
635  } else {
636  return &*getTangentAtGaussPts().data().begin();
637  }
638  };
639 
640  auto get_tangent_at_pts = [&]() {
641  auto ptr = get_tangent_ptr();
643  &ptr[0], &ptr[1], &ptr[2]);
644  return t_m_at_pts;
645  };
646 
647  auto calculate_squared_edge_length = [&]() {
648  std::vector<double> l1;
649  if (nb_gauss_pts) {
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);
654  ++t_m_at_pts;
655  }
656  }
657  return l1;
658  };
659 
660  auto l1 = calculate_squared_edge_length();
661 
662  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
663 
664  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
665  const auto nb_dofs = data.getN(base).size2() / 3;
666 
667  if (nb_gauss_pts && nb_dofs) {
668  auto t_h_curl = get_base_at_pts(base);
669  int cc = 0;
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;
677  ++t_h_curl;
678  ++cc;
679  }
680  ++t_m_at_pts;
681  }
682 
683 #ifndef NDEBUG
684  if (cc != nb_gauss_pts * nb_dofs)
685  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
686 #endif
687  }
688  }
689 
691 }
692 
694  const FieldSpace space, boost::shared_ptr<VectorDouble> det_jac_ptr)
695  : OP(space), fieldSpace(space), detJacPtr(det_jac_ptr) {
696  if (!detJacPtr) {
697  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
698  }
699 }
700 
705 
706  auto scale = [&]() {
708  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
709  auto base = static_cast<FieldApproximationBase>(b);
710  auto &base_fun = data.getN(base);
711  if (base_fun.size2()) {
712 
713  auto &det_vec = *detJacPtr;
714  const auto nb_int_pts = base_fun.size1();
715 
716  if (nb_int_pts != det_vec.size())
717  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
718  "Number of integration pts in detJacPtr does not mush "
719  "number of integration pts in base function");
720 
721  for (auto gg = 0; gg != nb_int_pts; ++gg) {
722  auto row = ublas::matrix_row<MatrixDouble>(base_fun, gg);
723  row /= det_vec[gg];
724  }
725 
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);
730  row /= det_vec[gg];
731  }
732  }
733  }
734  }
736  };
737 
738  if (this->getFEDim() == 3) {
739  switch (fieldSpace) {
740  case H1:
741  CHKERR scale();
742  break;
743  case HCURL:
744  if (type >= MBEDGE)
745  CHKERR scale();
746  break;
747  case HDIV:
748  if (type >= MBTRI)
749  CHKERR scale();
750  break;
751  case L2:
752  if (type >= MBTET)
753  CHKERR scale();
754  break;
755  default:
756  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
757  }
758  } else if (this->getFEDim() == 2) {
759  switch (fieldSpace) {
760  case H1:
761  CHKERR scale();
762  break;
763  case HCURL:
764  if (type >= MBEDGE)
765  CHKERR scale();
766  break;
767  case HDIV:
768  case L2:
769  if (type >= MBTRI)
770  CHKERR scale();
771  break;
772  default:
773  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
774  }
775  } else if (this->getFEDim() == 1) {
776  switch (fieldSpace) {
777  case H1:
778  CHKERR scale();
779  break;
780  case HCURL:
781  case L2:
782  if (type >= MBEDGE)
783  CHKERR scale();
784  break;
785  default:
786  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
787  }
788  } else {
789  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
790  }
791 
793 }
794 
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) {
802 
803  if (!jac_ptr)
804  jac_ptr = boost::make_shared<MatrixDouble>();
805  if (!det_ptr)
806  det_ptr = boost::make_shared<VectorDouble>();
807  if (!inv_jac_ptr)
808  inv_jac_ptr = boost::make_shared<MatrixDouble>();
809 
810  if (geom_field_name.empty()) {
811 
812  pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
813 
814  } else {
815 
816  pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
817  pipeline.push_back(
818  new OpCalculateVectorFieldGradient<2, 2>(geom_field_name, jac_ptr));
819  pipeline.push_back(new OpGetHONormalsOnFace<2>(geom_field_name));
820  }
821 
822  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
823  pipeline.push_back(new OpSetHOWeightsOnFace());
824 
825  for (auto s : spaces) {
826  switch (s) {
827  case NOSPACE:
828  break;
829  case H1:
830  case L2:
831  pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(s, inv_jac_ptr));
832  break;
833  case HCURL:
834  pipeline.push_back(new OpSetCovariantPiolaTransformOnFace2D(inv_jac_ptr));
835  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
836  break;
837  case HDIV:
838  pipeline.push_back(new OpMakeHdivFromHcurl());
839  pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
840  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
841  break;
842  default:
843  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
844  "Space %s not yet implemented", FieldSpaceNames[s]);
845  }
846  }
847 
849 }
850 
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) {
858 
859  if (!jac_ptr)
860  jac_ptr = boost::make_shared<MatrixDouble>();
861  if (!det_ptr)
862  det_ptr = boost::make_shared<VectorDouble>();
863  if (!inv_jac_ptr)
864  inv_jac_ptr = boost::make_shared<MatrixDouble>();
865 
866  if (geom_field_name.empty()) {
867 
868  } else {
869 
870  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
871  pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
872  }
873 
874  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
875  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
876  pipeline.push_back(new OpSetHOWeightsOnFace());
877 
878  for (auto s : spaces) {
879  switch (s) {
880  case NOSPACE:
881  break;
882  case H1:
883  pipeline.push_back(
884  new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
885  break;
886  case HDIV:
887  pipeline.push_back(new OpMakeHdivFromHcurl());
888  pipeline.push_back(
890  jac_ptr));
891  pipeline.push_back(
892  new OpSetInvJacHcurlFaceEmbeddedIn3DSpace(inv_jac_ptr));
893  break;
894  case L2:
895  pipeline.push_back(
896  new OpSetInvJacL2ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
897  break;
898  default:
899  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
900  "Space %s not yet implemented", FieldSpaceNames[s]);
901  }
902  }
903 
905 }
906 
908  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
909  std::vector<FieldSpace> spaces, std::string geom_field_name) {
911 
912  if (geom_field_name.empty()) {
913 
914  } else {
915 
916  pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
917  pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
918  }
919 
920  for (auto s : spaces) {
921  switch (s) {
922  case NOSPACE:
923  break;
924  case HCURL:
925  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
926  break;
927  case HDIV:
928  pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
929  break;
930  default:
931  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
932  "Space %s not yet implemented", FieldSpaceNames[s]);
933  }
934  }
935 
937 }
938 
940  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
941  std::vector<FieldSpace> spaces, std::string geom_field_name) {
943 
944  if (geom_field_name.empty()) {
945 
946  } else {
947 
948  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
949  pipeline.push_back(new OpGetHOTangentsOnEdge<3>(geom_field_name));
950  }
951 
952  for (auto s : spaces) {
953  switch (s) {
954  case HCURL:
955  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
956  break;
957  default:
958  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
959  "Space %s not yet implemented", FieldSpaceNames[s]);
960  }
961  }
962 
964 }
965 
967  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
968  std::vector<FieldSpace> spaces, std::string geom_field_name) {
970 
971  if (geom_field_name.empty()) {
972  } else {
973 
974  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
975  pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
976  }
977 
978  for (auto s : spaces) {
979  switch (s) {
980  case NOSPACE:
981  break;
982  case HCURL:
983  pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
984  break;
985  case HDIV:
986  pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
987  break;
988  case L2:
989  break;
990  default:
991  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
992  "Space %s not yet implemented", FieldSpaceNames[s]);
993  break;
994  }
995  }
996 
998 }
999 
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) {
1006 
1007  if (!jac)
1008  jac = boost::make_shared<MatrixDouble>();
1009  if (!det)
1010  det = boost::make_shared<VectorDouble>();
1011  if (!inv_jac)
1012  inv_jac = boost::make_shared<MatrixDouble>();
1013 
1014  if (geom_field_name.empty()) {
1015 
1016  pipeline.push_back(new OpCalculateHOJac<3>(jac));
1017 
1018  } else {
1019 
1020  pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
1021  pipeline.push_back(
1022  new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
1023  }
1024 
1025  pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
1026  pipeline.push_back(new OpSetHOWeights(det));
1027 
1028  for (auto s : spaces) {
1029  switch (s) {
1030  case NOSPACE:
1031  break;
1032  case H1:
1033  pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
1034  break;
1035  case HCURL:
1036  pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
1037  pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
1038  break;
1039  case HDIV:
1040  pipeline.push_back(
1041  new OpSetHOContravariantPiolaTransform(HDIV, det, jac));
1042  pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac));
1043  break;
1044  case L2:
1045  pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
1046  break;
1047  default:
1048  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1049  "Space %s not yet implemented", FieldSpaceNames[s]);
1050  break;
1051  }
1052  }
1053 
1055 }
1056 
1057 } // 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:2978
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:3076
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:609
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:180
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:3028
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:403
HVEC0_0
@ HVEC0_0
Definition: definitions.h:205
scale
double scale
Definition: plastic.cpp:119
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
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1240
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1275
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:1236
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:2940
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:3078
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:402
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:175
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:3086
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
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:413
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
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::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:334
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEDim
int getFEDim() const
Get dimension of finite element.
Definition: ForcesAndSourcesCore.hpp:1008
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:3249
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:2985
HVEC2
@ HVEC2
Definition: definitions.h:199
MoFEM::OpSetInvJacHcurlFaceEmbeddedIn3DSpace
OpSetInvJacHcurlFaceImpl< 3 > OpSetInvJacHcurlFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:2979
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:459
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:2932
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.
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:512
HVEC1_0
@ HVEC1_0
Definition: definitions.h:206
MoFEM::OpScaleBaseBySpaceInverseOfMeasure::OpScaleBaseBySpaceInverseOfMeasure
OpScaleBaseBySpaceInverseOfMeasure(const FieldSpace space, boost::shared_ptr< VectorDouble > det_jac_ptr)
Definition: HODataOperators.cpp:693
HVEC0
@ HVEC0
Definition: definitions.h:199
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