v0.14.0
UserDataOperators.cpp
Go to the documentation of this file.
1 /** \file UserDataOperators.cpp
2 
3 \brief Generic user data operators for evaluate fields, and other common
4 purposes.
5 
6 */
7 
8 namespace MoFEM {
9 
11  FieldSpace space, boost::shared_ptr<MatrixDouble> inv_jac_ptr)
12  : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {
13  if (!inv_jac_ptr)
14  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "invJacPtr not allocated");
15  if (space == L2) {
16  doVertices = false;
17  }
18 }
19 
24 
25  if (getFEDim() != 2)
26  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
27  "This operator can be used only with element which faces");
28 
29  auto apply_transform = [&](MatrixDouble &diff_n) {
30  return applyTransform<2, 2, 2, 2>(diff_n);
31  };
32 
33  if (!(type == MBVERTEX && sPace == L2)) {
34 
35  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
36  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
37  CHKERR apply_transform(data.getDiffN(base));
38  }
39 
40  switch (type) {
41  case MBVERTEX:
42  for (auto &m : data.getBBDiffNMap())
43  CHKERR apply_transform(*(m.second));
44  break;
45  default:
46  for (auto &ptr : data.getBBDiffNByOrderArray())
47  if (ptr)
48  CHKERR apply_transform(*ptr);
49  }
50  }
51 
53 }
54 
59 
60  if (getFEDim() != 2)
61  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
62  "This operator can be used only with element which face");
63 
64  auto apply_transform = [&](MatrixDouble &diff_n) {
65  return applyTransform<2, 3, 3, 3>(diff_n);
66  };
67 
68  if (!(type == MBVERTEX && sPace == L2)) {
69  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
70  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
71  CHKERR apply_transform(data.getDiffN(base));
72  }
73 
74  switch (type) {
75  case MBVERTEX:
76  for (auto &m : data.getBBDiffNMap())
77  CHKERR apply_transform(*(m.second));
78  break;
79  default:
80  for (auto &ptr : data.getBBDiffNByOrderArray())
81  if (ptr)
82  CHKERR apply_transform(*ptr);
83  }
84  }
85 
87 }
88 
93 
94  if (getFEDim() != 2)
95  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
96  "This operator can be used only with element which faces");
97 
98  auto apply_transform = [&](MatrixDouble &diff2_n) {
100  size_t nb_functions = diff2_n.size2() / 4;
101  if (nb_functions) {
102  size_t nb_gauss_pts = diff2_n.size1();
103 
104 #ifndef NDEBUG
105  if (nb_gauss_pts != getGaussPts().size2())
106  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
107  "Wrong number of Gauss Pts");
108 #endif
109 
110  diffNinvJac.resize(diff2_n.size1(), diff2_n.size2(), false);
111 
116  auto t_diff2_n = getFTensor2FromPtr<2, 2>(&*diffNinvJac.data().begin());
117  auto t_diff2_n_ref = getFTensor2FromPtr<2, 2>(&*diff2_n.data().begin());
118  auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
119  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
120  for (size_t dd = 0; dd != nb_functions; ++dd) {
121  t_diff2_n(i, j) =
122  t_inv_jac(k, i) * (t_inv_jac(l, j) * t_diff2_n_ref(k, l));
123  ++t_diff2_n;
124  ++t_diff2_n_ref;
125  }
126  }
127 
128  diff2_n.swap(diffNinvJac);
129  }
131  };
132 
133  if (!(type == MBVERTEX && sPace == L2)) {
134 
135  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
136  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
137  CHKERR apply_transform(
138  data.getN(base, BaseDerivatives::SecondDerivative));
139  }
140 
141  auto error = [&](auto &m) {
143  if (m.size2())
144  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
145  "Operator do not work for Bernstein-Bezier. This "
146  "functionality is "
147  "not implemented.");
149  };
150 
151  switch (type) {
152  case MBVERTEX:
153  for (auto &m : data.getBBDiffNMap()) {
154  CHKERR error(*(m.second));
155  CHKERR apply_transform(*(m.second));
156  }
157  break;
158  default:
159  for (auto &ptr : data.getBBDiffNByOrderArray())
160  if (ptr) {
161  CHKERR error(*ptr);
162  CHKERR apply_transform(*ptr);
163  }
164  }
165  }
166 
168 }
169 
174 
175  if (type != MBEDGE && type != MBTRI && type != MBQUAD)
177 
178  if (getFEDim() != 2)
179  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
180  "This operator can be used only with element which is triangle");
181 
185 
186  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
187 
188  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
189 
190  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
191  if (nb_base_functions) {
192  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
193 
194  diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
195 
196  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
197  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
199  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
200 
201  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
202 
203  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
204 
205  auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
206  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
207  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
208  t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
209  ++t_diff_n;
210  ++t_inv_diff_n;
211  }
212  }
213 
214  data.getDiffN(base).swap(diffHcurlInvJac);
215  }
216  }
217 
219 }
220 
225 
226  if (type != MBEDGE && type != MBTRI && type != MBQUAD)
228 
229  if (getFEDim() != 2)
230  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
231  "This operator can be used only with element which is triangle");
232 
236 
237  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
238 
239  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
240 
241  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
242  if (nb_base_functions) {
243  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
244 
245  diffHcurlInvJac.resize(nb_gauss_pts, nb_base_functions * 9, false);
246 
247  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
248  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
250  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
251  &t_inv_diff_n_ptr[HVEC0_2],
252 
253  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
254  &t_inv_diff_n_ptr[HVEC1_2],
255 
256  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
257  &t_inv_diff_n_ptr[HVEC2_2]);
258 
259  auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
260  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
261  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
262  t_inv_diff_n(i, j) = t_diff_n(i, K) * t_inv_jac(K, j);
263  ++t_diff_n;
264  ++t_inv_diff_n;
265  }
266  }
267 
268  data.getDiffN(base).swap(diffHcurlInvJac);
269  }
270  }
271 
273 }
274 
278 
279  if (type != MBEDGE && type != MBTRI && type != MBQUAD)
281 
282  if (getFEDim() != 2)
283  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
284  "This operator can be used only with element which is face");
285 
286  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
287 
288  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
289 
290  const size_t nb_base_functions = data.getN(base).size2() / 3;
291  if (nb_base_functions) {
292 
293  const size_t nb_gauss_pts = data.getN(base).size1();
294 
295  auto t_n = data.getFTensor1N<3>(base);
296  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
297  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
298  for (size_t bb = 0; bb != nb_base_functions; ++bb) {
299 
300  const double a = t_n(0);
301  t_n(0) = -t_n(1);
302  t_n(1) = a;
303 
304  for (auto n : {0, 1}) {
305  const double b = t_diff_n(0, n);
306  t_diff_n(0, n) = -t_diff_n(1, n);
307  t_diff_n(1, n) = b;
308  }
309 
310  ++t_n;
311  ++t_diff_n;
312  }
313  }
314  }
315  }
316 
318 }
319 
321  2>::OpSetCovariantPiolaTransformOnFace2DImpl(boost::shared_ptr<MatrixDouble>
322  inv_jac_ptr)
324  invJacPtr(inv_jac_ptr) {}
325 
327  int side, EntityType type, EntitiesFieldData::EntData &data) {
329 
330  const auto type_dim = moab::CN::Dimension(type);
331  if (type_dim != 1 && type_dim != 2)
333 
337 
338  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
339 
340  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
341 
342  auto &baseN = data.getN(base);
343  auto &diffBaseN = data.getDiffN(base);
344 
345  int nb_dofs = baseN.size2() / 3;
346  int nb_gauss_pts = baseN.size1();
347 
348  piolaN.resize(baseN.size1(), baseN.size2());
349 
350  if (nb_dofs > 0 && nb_gauss_pts > 0) {
351 
352  auto t_h_curl = data.getFTensor1N<3>(base);
353  auto t_transformed_h_curl =
354  getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
355  auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
356  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
357  for (int ll = 0; ll != nb_dofs; ll++) {
358  t_transformed_h_curl(i) = t_inv_jac(j, i) * t_h_curl(j);
359  ++t_h_curl;
360  ++t_transformed_h_curl;
361  }
362  ++t_inv_jac;
363  }
364  baseN.swap(piolaN);
365 
366  diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
367  if (diffBaseN.data().size() > 0) {
368  auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
369  auto t_transformed_diff_h_curl =
370  getFTensor2HVecFromPtr<3, 2>(&*diffPiolaN.data().begin());
371  auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
372  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
373  for (int ll = 0; ll != nb_dofs; ll++) {
374  t_transformed_diff_h_curl(i, k) =
375  t_inv_jac(j, i) * t_diff_h_curl(j, k);
376  ++t_diff_h_curl;
377  ++t_transformed_diff_h_curl;
378  }
379  ++t_inv_jac;
380  }
381  diffBaseN.swap(diffPiolaN);
382  }
383  }
384  }
385 
387 }
388 
390  int side, EntityType type, EntitiesFieldData::EntData &data) {
391 
393 
394  if (type != MBEDGE && type != MBTRI && type != MBQUAD)
396 
400 
401  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
402 
403  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
404 
405  const size_t nb_base_functions = data.getN(base).size2() / 3;
406  if (nb_base_functions) {
407 
408  const size_t nb_gauss_pts = data.getN(base).size1();
409  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
410  if (data.getN(base).size2() > 0) {
411  auto t_n = data.getFTensor1N<3>(base);
412  double *t_transformed_n_ptr = &*piolaN.data().begin();
414  t_transformed_n_ptr, // HVEC0
415  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
416  auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
417  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
418  double det;
419  CHKERR determinantTensor2by2(t_jac, det);
420  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
421  t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
422  ++t_n;
423  ++t_transformed_n;
424  }
425  }
426  data.getN(base).swap(piolaN);
427  }
428 
429  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
430  if (data.getDiffN(base).data().size() > 0) {
431  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
432  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
434  t_transformed_diff_n(t_transformed_diff_n_ptr,
435  &t_transformed_diff_n_ptr[HVEC0_1],
436 
437  &t_transformed_diff_n_ptr[HVEC1_0],
438  &t_transformed_diff_n_ptr[HVEC1_1],
439 
440  &t_transformed_diff_n_ptr[HVEC2_0],
441  &t_transformed_diff_n_ptr[HVEC2_1]);
442  auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
443  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
444  double det;
445  CHKERR determinantTensor2by2(t_jac, det);
446  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
447  t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
448  ++t_diff_n;
449  ++t_transformed_diff_n;
450  }
451  }
452  data.getDiffN(base).swap(piolaDiffN);
453  }
454  }
455  }
456 
458 }
459 
461  int side, EntityType type, EntitiesFieldData::EntData &data) {
462 
464 
465  if (type != MBEDGE && type != MBTRI && type != MBQUAD)
467 
471 
472  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
473 
474  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
475 
476  const size_t nb_base_functions = data.getN(base).size2() / 3;
477  if (nb_base_functions) {
478 
479  const size_t nb_gauss_pts = data.getN(base).size1();
480  piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
481  if (data.getN(base).size2() > 0) {
482  auto t_n = data.getFTensor1N<3>(base);
483  double *t_transformed_n_ptr = &*piolaN.data().begin();
485  t_transformed_n_ptr, // HVEC0
486  &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
487  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
488  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
489  double det;
490  CHKERR determinantTensor3by3(t_jac, det);
491  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
492  t_transformed_n(i) = t_jac(i, j) * t_n(j) / det;
493  ++t_n;
494  ++t_transformed_n;
495  }
496  }
497  data.getN(base).swap(piolaN);
498  }
499 
500  piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
501  if (data.getDiffN(base).size2() > 0) {
502  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
503  double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
505  t_transformed_diff_n(t_transformed_diff_n_ptr,
506  &t_transformed_diff_n_ptr[HVEC0_1],
507 
508  &t_transformed_diff_n_ptr[HVEC1_0],
509  &t_transformed_diff_n_ptr[HVEC1_1],
510 
511  &t_transformed_diff_n_ptr[HVEC2_0],
512  &t_transformed_diff_n_ptr[HVEC2_1]);
513 
514  auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
515  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
516  double det;
517  CHKERR determinantTensor3by3(t_jac, det);
518  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
519  t_transformed_diff_n(i, K) = t_jac(i, j) * t_diff_n(j, K) / det;
520  ++t_diff_n;
521  ++t_transformed_diff_n;
522  }
523  }
524  data.getDiffN(base).swap(piolaDiffN);
525  }
526  }
527  }
528 
530 }
531 
533  int side, EntityType type, EntitiesFieldData::EntData &data) {
535 
536  if (type != MBEDGE)
538 
542 
543  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
544 
545  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
546  size_t nb_gauss_pts = data.getN(base).size1();
547  size_t nb_dofs = data.getN(base).size2() / 3;
548  if (nb_dofs > 0) {
549 
550  auto t_h_div = data.getFTensor1N<3>(base);
551 
552  size_t cc = 0;
553  const auto &edge_direction_at_gauss_pts = getTangentAtGaussPts();
554  if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
555 
556  auto t_dir = getFTensor1TangentAtGaussPts<3>();
557 
558  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
560  t_normal(i) = (FTensor::levi_civita(i, j, k) *
562  t_dir(k);
563  const auto l2 = t_normal(i) * t_normal(i);
564  for (int ll = 0; ll != nb_dofs; ++ll) {
565  const double val = t_h_div(0);
566  const double a = val / l2;
567  t_h_div(i) = t_normal(i) * a;
568  ++t_h_div;
569  ++cc;
570  }
571  ++t_dir;
572  }
573  }
574 
575  if (cc != nb_gauss_pts * nb_dofs)
576  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
577  }
578  }
579 
581 }
582 
584  int side, EntityType type, EntitiesFieldData::EntData &data) {
586 
587  if (type == MBVERTEX) {
588 
589  VectorDouble &coords = getCoords();
590  double *coords_ptr = &*coords.data().begin();
591 
592  const int nb_gauss_pts = data.getN(NOBASE).size1();
593  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
594 
598 
599  auto t_w = getFTensor0IntegrationWeight();
600  for (int gg = 0; gg != nb_gauss_pts; gg++) {
601 
602  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
603  &coords_ptr[2], 3);
604  t_jac(i, j) = 0;
605  for (int bb = 0; bb != 6; bb++) {
606  t_jac(i, j) += t_coords(i) * t_diff_n(j);
607  ++t_diff_n;
608  ++t_coords;
609  }
610 
611  double det;
612  CHKERR determinantTensor3by3(t_jac, det);
613  t_w *= det / 2.;
614 
615  ++t_w;
616  }
617 
618  double &vol = getVolume();
619  auto t_w_scaled = getFTensor0IntegrationWeight();
620  for (int gg = 0; gg != nb_gauss_pts; gg++) {
621  t_w_scaled /= vol;
622  ++t_w_scaled;
623  }
624  }
625 
626  doEntities[MBVERTEX] = true;
627  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
628 
630 }
631 
636 
637  if (type == MBVERTEX) {
638 
639  VectorDouble &coords = getCoords();
640  double *coords_ptr = &*coords.data().begin();
641 
642  const int nb_gauss_pts = data.getN(NOBASE).size1();
643  auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
644  invJac.resize(9, nb_gauss_pts, false);
645  invJac.clear();
646  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
647 
651 
652  for (int gg = 0; gg != nb_gauss_pts; gg++) {
653 
654  FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
655  &coords_ptr[2], 3);
656  t_jac(i, j) = 0;
657  for (int bb = 0; bb != 6; bb++) {
658  t_jac(i, j) += t_coords(i) * t_diff_n(j);
659  ++t_diff_n;
660  ++t_coords;
661  }
662 
663  double det;
664  CHKERR determinantTensor3by3(t_jac, det);
665  CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
666  ++t_inv_jac;
667  }
668  }
669 
670  doEntities[MBVERTEX] = true;
671  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
672 
674 }
675 
680 
681  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
682 
683  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
684  if (data.getN(base).size2() == 0)
685  continue;
686 
687  const int nb_gauss_pts = data.getN(base).size1();
688  auto t_diff_n = data.getFTensor1DiffN<3>(base);
689  diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
690  false);
691 
694 
696  &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
697  auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
698 
699  const int nb_dofs = data.getN(base).size2();
700  for (int gg = 0; gg != nb_gauss_pts; gg++) {
701  for (int bb = 0; bb != nb_dofs; bb++) {
702  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
703  ++t_inv_diff_n;
704  ++t_diff_n;
705  }
706  ++t_inv_jac;
707  }
708 
709  data.getDiffN(base).swap(diffNinvJac);
710  }
711 
713 }
714 
718 
720 
721  if (type == MBVERTEX) {
722 
723  VectorDouble &coords = getCoords();
724  double *coords_ptr = &*coords.data().begin();
725  double diff_n[6];
726  CHKERR ShapeDiffMBTRI(diff_n);
727  double j00_f3, j01_f3, j10_f3, j11_f3;
728  for (int gg = 0; gg < 1; gg++) {
729  // this is triangle, derivative of nodal shape functions is constant.
730  // So only need to do one node.
731  j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
732  j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
733  j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
734  j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
735  }
736  double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
737  invJacF3.resize(2, 2, false);
738  invJacF3(0, 0) = j11_f3 / det_f3;
739  invJacF3(0, 1) = -j01_f3 / det_f3;
740  invJacF3(1, 0) = -j10_f3 / det_f3;
741  invJacF3(1, 1) = j00_f3 / det_f3;
742  }
743 
744  doEntities[MBVERTEX] = true;
745  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
746 
748 }
749 
754 
755  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
756 
757  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
758 
759  unsigned int nb_dofs = data.getN(base).size2();
760  if (nb_dofs == 0)
762  unsigned int nb_gauss_pts = data.getN(base).size1();
763  diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
764 
765 #ifndef NDEBUG
766  if (type != MBVERTEX) {
767  if (nb_dofs != data.getDiffN(base).size2() / 2) {
768  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
769  "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
770  "%u/2 )",
771  nb_dofs, data.getDiffN(base).size2());
772  }
773  }
774 #endif
775 
776  switch (type) {
777  case MBVERTEX:
778  case MBEDGE:
779  case MBTRI: {
780  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
781  for (unsigned int dd = 0; dd < nb_dofs; dd++) {
782  cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
783  &*invJacF3.data().begin(), 2,
784  &data.getDiffN(base)(gg, 2 * dd), 1, 0,
785  &diffNinvJac(gg, 2 * dd), 1);
786  }
787  }
788  data.getDiffN(base).swap(diffNinvJac);
789  } break;
790  default:
791  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
792  }
793  }
794 
796 }
797 
799  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
800  const EntityType zero_type, const int zero_side)
803  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
804  if (!dataPtr)
805  THROW_MESSAGE("Pointer is not set");
806 }
807 
812  const auto nb_integration_points = getGaussPts().size2();
813  if (type == zeroType && side == zeroSide) {
814  dataPtr->resize(3, nb_integration_points, false);
815  dataPtr->clear();
816  }
817  const auto nb_dofs = data.getFieldData().size();
818  if (!nb_dofs)
823  const auto nb_base_functions = data.getN().size2() / 3;
824  auto t_n_diff_hcurl = data.getFTensor2DiffN<3, 3>();
825  auto t_data = getFTensor1FromMat<3>(*dataPtr);
826  for (auto gg = 0; gg != nb_integration_points; ++gg) {
827  auto t_dof = data.getFTensor0FieldData();
828  int bb = 0;
829  for (; bb != nb_dofs; ++bb) {
830  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
831  ++t_n_diff_hcurl;
832  ++t_dof;
833  }
834  for (; bb < nb_base_functions; ++bb)
835  ++t_n_diff_hcurl;
836  ++t_data;
837  }
838 
840 }
841 
843  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
844  const EntityType zero_type, const int zero_side)
847  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
848  if (!dataPtr)
849  THROW_MESSAGE("Pointer is not set");
850 }
851 
856  const auto nb_integration_points = getGaussPts().size2();
857  if (type == zeroType && side == zeroSide) {
858  dataPtr->resize(2, nb_integration_points, false);
859  dataPtr->clear();
860  }
861  const auto nb_dofs = data.getFieldData().size();
862  if (!nb_dofs)
864 
867 
868  const auto nb_base_functions = data.getN().size2();
869  auto t_n_diff_hcurl = data.getFTensor1DiffN<2>();
870  auto t_data = getFTensor1FromMat<2>(*dataPtr);
871  for (auto gg = 0; gg != nb_integration_points; ++gg) {
872  auto t_dof = data.getFTensor0FieldData();
873  int bb = 0;
874  for (; bb != nb_dofs; ++bb) {
875  t_data(i) += t_dof * (levi_civita(i, j) * t_n_diff_hcurl(j));
876  ++t_n_diff_hcurl;
877  ++t_dof;
878  }
879  for (; bb < nb_base_functions; ++bb)
880  ++t_n_diff_hcurl;
881  ++t_data;
882  }
883 
885 }
886 
888  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
889  const EntityType zero_type, const int zero_side)
892  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
893  if (!dataPtr)
894  THROW_MESSAGE("Pointer is not set");
895 }
896 
901  const auto nb_integration_points = getGaussPts().size2();
902  if (type == zeroType && side == zeroSide) {
903  dataPtr->resize(3, nb_integration_points, false);
904  dataPtr->clear();
905  }
906  const auto nb_dofs = data.getFieldData().size();
907  if (!nb_dofs)
909 
913 
914  const auto nb_base_functions = data.getN().size2();
915  auto t_n_diff_hcurl = data.getFTensor1DiffN<3>();
916  auto t_data = getFTensor1FromMat<3>(*dataPtr);
917  auto t_normal = getFTensor1NormalsAtGaussPts();
918  for (auto gg = 0; gg != nb_integration_points; ++gg) {
919  auto t_dof = data.getFTensor0FieldData();
920  auto nrm = t_normal.l2();
921  int bb = 0;
922  for (; bb != nb_dofs; ++bb) {
923  t_data(k) += (t_dof / nrm) *
924  ((levi_civita(i, j, k) * t_normal(i)) * t_n_diff_hcurl(j));
925  ++t_n_diff_hcurl;
926  ++t_dof;
927  }
928  for (; bb < nb_base_functions; ++bb)
929  ++t_n_diff_hcurl;
930  ++t_data;
931  ++t_normal;
932  }
933 
935 }
936 
937 } // namespace MoFEM
UBlasMatrix< double >
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl
Transform Hcurl base fluxes from reference element to physical triangle.
Definition: UserDataOperators.hpp:2998
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::OpCalculateInvJacForFatPrism::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:633
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
HVEC0_2
@ HVEC0_2
Definition: definitions.h:211
MoFEM::OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:583
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::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
HVEC0_1
@ HVEC0_1
Definition: definitions.h:208
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getCoords
VectorDouble & getCoords()
nodal coordinates
Definition: VolumeElementForcesAndSourcesCore.hpp:180
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
MoFEM::OpSetInvJacH1ForFatPrism::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:677
MoFEM::EntitiesFieldData::EntData::getBBDiffNByOrderArray
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
Definition: EntitiesFieldData.cpp:902
MoFEM::OpSetInvJacH1ForFatPrism::diffNinvJac
MatrixDouble diffNinvJac
Definition: UserDataOperators.hpp:3189
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
HVEC1_1
@ HVEC1_1
Definition: definitions.h:209
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MoFEM::DataOperator::doVertices
bool & doVertices
\deprectaed If false skip vertices
Definition: DataOperators.hpp:90
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
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::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1240
MoFEM::OpSetInvJacToScalarBasesBasic::OpSetInvJacToScalarBasesBasic
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.cpp:10
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
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
MoFEM::determinantTensor2by2
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition: Templates.hpp:1553
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
HVEC2_2
@ HVEC2_2
Definition: definitions.h:213
HVEC1_2
@ HVEC1_2
Definition: definitions.h:212
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::OpCalculateInvJacForFlatPrism::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:716
convert.type
type
Definition: convert.py:64
MoFEM::OpSetInvJacH1ForFlatPrism::invJacF3
MatrixDouble & invJacF3
Definition: UserDataOperators.hpp:3238
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpSetInvJacHcurlFaceImpl
Transform local reference derivatives of shape function to global derivatives for face.
Definition: UserDataOperators.hpp:2953
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 2 >
convert.n
n
Definition: convert.py:82
MoFEM::EdgeElementForcesAndSourcesCore::tFaceOrientation
static FTensor::Tensor1< double, 3 > tFaceOrientation
Definition: EdgeElementForcesAndSourcesCore.hpp:44
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
MoFEM::EntitiesFieldData::EntData::getFTensor0FieldData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Definition: EntitiesFieldData.cpp:506
MoFEM::DataOperator::doEntities
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Definition: DataOperators.hpp:85
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
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::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::OpCalculateInvJacForFatPrism::invJac
MatrixDouble & invJac
Definition: UserDataOperators.hpp:3160
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEDim
int getFEDim() const
Get dimension of finite element.
Definition: ForcesAndSourcesCore.hpp:1008
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:957
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
UBlasVector< double >
MoFEM::OpSetInvJacH1ForFlatPrism::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:751
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpMakeHdivFromHcurl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:275
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::OpCalculateInvJacForFlatPrism::invJacF3
MatrixDouble & invJacF3
Definition: UserDataOperators.hpp:3215
MoFEM::FlatPrismElementForcesAndSourcesCore::UserDataOperator::getCoords
VectorDouble & getCoords()
get triangle coordinates
Definition: FlatPrismElementForcesAndSourcesCore.hpp:208
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
HVEC2_0
@ HVEC2_0
Definition: definitions.h:207
MoFEM::OpSetContravariantPiolaTransformOnEdge2D::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: UserDataOperators.cpp:532
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::OpSetInvJacH1ForFatPrism::invJac
MatrixDouble & invJac
Definition: UserDataOperators.hpp:3188
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
HVEC2
@ HVEC2
Definition: definitions.h:199
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::OpSetInvJacH1ForFlatPrism::diffNinvJac
MatrixDouble diffNinvJac
Definition: UserDataOperators.hpp:3239
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl
Apply contravariant (Piola) transfer to Hdiv space on face.
Definition: UserDataOperators.hpp:3045
MoFEM::OpCalculateHcurlVectorCurl
Calculate curl of vector field.
Definition: UserDataOperators.hpp:2488
MoFEM::OpSetInvJacSpaceForFaceImpl
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:2833
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
HVEC1_0
@ HVEC1_0
Definition: definitions.h:206