v0.13.2
Loading...
Searching...
No Matches
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
4purposes.
5
6*/
7
8
9
10namespace MoFEM {
11
13 FieldSpace space, boost::shared_ptr<MatrixDouble> inv_jac_ptr)
14 : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {
15 if (!inv_jac_ptr)
16 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "invJacPtr not allocated");
17 if (space == L2) {
18 doVertices = false;
19 }
20}
21
26
27 if (getFEDim() != 2)
28 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
29 "This operator can be used only with element which faces");
30
31 auto apply_transform = [&](MatrixDouble &diff_n) {
32 return applyTransform<2, 2, 2, 2>(diff_n);
33 };
34
35 if (!(type == MBVERTEX && sPace == L2)) {
36
37 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
39 CHKERR apply_transform(data.getDiffN(base));
40 }
41
42 switch (type) {
43 case MBVERTEX:
44 for (auto &m : data.getBBDiffNMap())
45 CHKERR apply_transform(*(m.second));
46 break;
47 default:
48 for (auto &ptr : data.getBBDiffNByOrderArray())
49 if (ptr)
50 CHKERR apply_transform(*ptr);
51 }
52 }
53
55}
56
61
62 if (getFEDim() != 2)
63 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
64 "This operator can be used only with element which face");
65
66 auto apply_transform = [&](MatrixDouble &diff_n) {
67 return applyTransform<2, 3, 3, 3>(diff_n);
68 };
69
70 if (!(type == MBVERTEX && sPace == L2)) {
71 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
73 CHKERR apply_transform(data.getDiffN(base));
74 }
75
76 switch (type) {
77 case MBVERTEX:
78 for (auto &m : data.getBBDiffNMap())
79 CHKERR apply_transform(*(m.second));
80 break;
81 default:
82 for (auto &ptr : data.getBBDiffNByOrderArray())
83 if (ptr)
84 CHKERR apply_transform(*ptr);
85 }
86 }
87
89}
90
95
96 if (getFEDim() != 2)
97 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
98 "This operator can be used only with element which faces");
99
100 auto apply_transform = [&](MatrixDouble &diff2_n) {
102 size_t nb_functions = diff2_n.size2() / 4;
103 if (nb_functions) {
104 size_t nb_gauss_pts = diff2_n.size1();
105
106#ifndef NDEBUG
107 if (nb_gauss_pts != getGaussPts().size2())
108 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
109 "Wrong number of Gauss Pts");
110#endif
111
112 diffNinvJac.resize(diff2_n.size1(), diff2_n.size2(), false);
113
118 auto t_diff2_n = getFTensor2FromPtr<2, 2>(&*diffNinvJac.data().begin());
119 auto t_diff2_n_ref = getFTensor2FromPtr<2, 2>(&*diff2_n.data().begin());
120 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
121 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
122 for (size_t dd = 0; dd != nb_functions; ++dd) {
123 t_diff2_n(i, j) =
124 t_inv_jac(k, i) * (t_inv_jac(l, j) * t_diff2_n_ref(k, l));
125 ++t_diff2_n;
126 ++t_diff2_n_ref;
127 }
128 }
129
130 diff2_n.swap(diffNinvJac);
131 }
133 };
134
135 if (!(type == MBVERTEX && sPace == L2)) {
136
137 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
138 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
139 CHKERR apply_transform(
140 data.getN(base, BaseDerivatives::SecondDerivative));
141 }
142
143 auto error = [&](auto &m) {
145 if (m.size2())
146 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
147 "Operator do not work for Bernstein-Bezier. This "
148 "functionality is "
149 "not implemented.");
151 };
152
153 switch (type) {
154 case MBVERTEX:
155 for (auto &m : data.getBBDiffNMap()) {
156 CHKERR error(*(m.second));
157 CHKERR apply_transform(*(m.second));
158 }
159 break;
160 default:
161 for (auto &ptr : data.getBBDiffNByOrderArray())
162 if (ptr) {
163 CHKERR error(*ptr);
164 CHKERR apply_transform(*ptr);
165 }
166 }
167 }
168
170}
171
176
177 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
179
180 if (getFEDim() != 2)
181 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
182 "This operator can be used only with element which is triangle");
183
187
188 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
189
190 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
191
192 const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
193 if (nb_base_functions) {
194 const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
195
196 diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
197
198 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
199 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
201 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
202
203 &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
204
205 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
206
207 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
208 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
209 for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
210 t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
211 ++t_diff_n;
212 ++t_inv_diff_n;
213 }
214 }
215
216 data.getDiffN(base).swap(diffHcurlInvJac);
217 }
218 }
219
221}
222
227
228 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
230
231 if (getFEDim() != 2)
232 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
233 "This operator can be used only with element which is triangle");
234
238
239 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
240
241 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
242
243 const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
244 if (nb_base_functions) {
245 const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
246
247 diffHcurlInvJac.resize(nb_gauss_pts, nb_base_functions * 9, false);
248
249 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
250 double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
252 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
253 &t_inv_diff_n_ptr[HVEC0_2],
254
255 &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
256 &t_inv_diff_n_ptr[HVEC1_2],
257
258 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
259 &t_inv_diff_n_ptr[HVEC2_2]);
260
261 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
262 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
263 for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
264 t_inv_diff_n(i, j) = t_diff_n(i, K) * t_inv_jac(K, j);
265 ++t_diff_n;
266 ++t_inv_diff_n;
267 }
268 }
269
270 data.getDiffN(base).swap(diffHcurlInvJac);
271 }
272 }
273
275}
276
281
282 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
284
285 if (getFEDim() != 2)
286 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
287 "This operator can be used only with element which is face");
288
289 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
290
291 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
292
293 const size_t nb_base_functions = data.getN(base).size2() / 3;
294 if (nb_base_functions) {
295
296 const size_t nb_gauss_pts = data.getN(base).size1();
297
298 auto t_n = data.getFTensor1N<3>(base);
299 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
300 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
301 for (size_t bb = 0; bb != nb_base_functions; ++bb) {
302
303 const double a = t_n(0);
304 t_n(0) = -t_n(1);
305 t_n(1) = a;
306
307 for (auto n : {0, 1}) {
308 const double b = t_diff_n(0, n);
309 t_diff_n(0, n) = -t_diff_n(1, n);
310 t_diff_n(1, n) = b;
311 }
312
313 ++t_n;
314 ++t_diff_n;
315 }
316 }
317 }
318 }
319
321}
322
324 2>::OpSetCovariantPiolaTransformOnFace2DImpl(boost::shared_ptr<MatrixDouble>
325 inv_jac_ptr)
327 invJacPtr(inv_jac_ptr) {}
328
330 int side, EntityType type, EntitiesFieldData::EntData &data) {
332
333 const auto type_dim = moab::CN::Dimension(type);
334 if (type_dim != 1 && type_dim != 2)
336
337 const auto nb_gauss_pts = getGaussPts().size2();
338
342
343 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
344
345 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
346
347 auto &baseN = data.getN(base);
348 auto &diffBaseN = data.getDiffN(base);
349
350 int nb_dofs = baseN.size2() / 3;
351 int nb_gauss_pts = baseN.size1();
352
353 piolaN.resize(baseN.size1(), baseN.size2());
354
355 if (nb_dofs > 0 && nb_gauss_pts > 0) {
356
357 auto t_h_curl = data.getFTensor1N<3>(base);
358 auto t_transformed_h_curl =
359 getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
360 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
361 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
362 for (int ll = 0; ll != nb_dofs; ll++) {
363 t_transformed_h_curl(i) = t_inv_jac(j, i) * t_h_curl(j);
364 ++t_h_curl;
365 ++t_transformed_h_curl;
366 }
367 ++t_inv_jac;
368 }
369 baseN.swap(piolaN);
370
371 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
372 if (diffBaseN.data().size() > 0) {
373 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
374 auto t_transformed_diff_h_curl =
375 getFTensor2FromPtr<3, 2>(&*diffPiolaN.data().begin());
376 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
377 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
378 for (int ll = 0; ll != nb_dofs; ll++) {
379 t_transformed_diff_h_curl(i, k) =
380 t_inv_jac(j, i) * t_diff_h_curl(j, k);
381 ++t_diff_h_curl;
382 ++t_transformed_diff_h_curl;
383 }
384 ++t_inv_jac;
385 }
386 diffBaseN.swap(diffPiolaN);
387 }
388 }
389 }
390
392}
393
395 int side, EntityType type, EntitiesFieldData::EntData &data) {
396
398
399 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
401
405
406 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
407
408 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
409
410 const size_t nb_base_functions = data.getN(base).size2() / 3;
411 if (nb_base_functions) {
412
413 const size_t nb_gauss_pts = data.getN(base).size1();
414 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
415 if (data.getN(base).size2() > 0) {
416 auto t_n = data.getFTensor1N<3>(base);
417 double *t_transformed_n_ptr = &*piolaN.data().begin();
419 t_transformed_n_ptr, // HVEC0
420 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
421 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
422 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
423 double det;
424 CHKERR determinantTensor2by2(t_jac, det);
425 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
426 t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
427 ++t_n;
428 ++t_transformed_n;
429 }
430 }
431 data.getN(base).swap(piolaN);
432 }
433
434 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
435 if (data.getDiffN(base).data().size() > 0) {
436 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
437 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
439 t_transformed_diff_n(t_transformed_diff_n_ptr,
440 &t_transformed_diff_n_ptr[HVEC0_1],
441
442 &t_transformed_diff_n_ptr[HVEC1_0],
443 &t_transformed_diff_n_ptr[HVEC1_1],
444
445 &t_transformed_diff_n_ptr[HVEC2_0],
446 &t_transformed_diff_n_ptr[HVEC2_1]);
447 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
448 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
449 double det;
450 CHKERR determinantTensor2by2(t_jac, det);
451 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
452 t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
453 ++t_diff_n;
454 ++t_transformed_diff_n;
455 }
456 }
457 data.getDiffN(base).swap(piolaDiffN);
458 }
459 }
460 }
461
463}
464
466 int side, EntityType type, EntitiesFieldData::EntData &data) {
467
469
470 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
472
476
477 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
478
479 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
480
481 const size_t nb_base_functions = data.getN(base).size2() / 3;
482 if (nb_base_functions) {
483
484 const size_t nb_gauss_pts = data.getN(base).size1();
485 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
486 if (data.getN(base).size2() > 0) {
487 auto t_n = data.getFTensor1N<3>(base);
488 double *t_transformed_n_ptr = &*piolaN.data().begin();
490 t_transformed_n_ptr, // HVEC0
491 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
492 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
493 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
494 double det;
495 CHKERR determinantTensor3by3(t_jac, det);
496 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
497 t_transformed_n(i) = t_jac(i, j) * t_n(j) / det;
498 ++t_n;
499 ++t_transformed_n;
500 }
501 }
502 data.getN(base).swap(piolaN);
503 }
504
505 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
506 if (data.getDiffN(base).size2() > 0) {
507 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
508 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
510 t_transformed_diff_n(t_transformed_diff_n_ptr,
511 &t_transformed_diff_n_ptr[HVEC0_1],
512
513 &t_transformed_diff_n_ptr[HVEC1_0],
514 &t_transformed_diff_n_ptr[HVEC1_1],
515
516 &t_transformed_diff_n_ptr[HVEC2_0],
517 &t_transformed_diff_n_ptr[HVEC2_1]);
518
519 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
520 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
521 double det;
522 CHKERR determinantTensor3by3(t_jac, det);
523 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
524 t_transformed_diff_n(i, K) = t_jac(i, j) * t_diff_n(j, K) / det;
525 ++t_diff_n;
526 ++t_transformed_diff_n;
527 }
528 }
529 data.getDiffN(base).swap(piolaDiffN);
530 }
531 }
532 }
533
535}
536
538 int side, EntityType type, EntitiesFieldData::EntData &data) {
540
541 if (type != MBEDGE)
543
547
548 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
549
550 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
551 size_t nb_gauss_pts = data.getN(base).size1();
552 size_t nb_dofs = data.getN(base).size2() / 3;
553 if (nb_dofs > 0) {
554
555 auto t_h_div = data.getFTensor1N<3>(base);
556
557 size_t cc = 0;
558 const auto &edge_direction_at_gauss_pts = getTangentAtGaussPts();
559 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
560
561 auto t_dir = getFTensor1TangentAtGaussPts<3>();
562
563 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
565 t_normal(i) = (FTensor::levi_civita(i, j, k) *
567 t_dir(k);
568 const auto l2 = t_normal(i) * t_normal(i);
569 for (int ll = 0; ll != nb_dofs; ++ll) {
570 const double val = t_h_div(0);
571 const double a = val / l2;
572 t_h_div(i) = t_normal(i) * a;
573 ++t_h_div;
574 ++cc;
575 }
576 ++t_dir;
577 }
578 }
579
580 if (cc != nb_gauss_pts * nb_dofs)
581 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
582 }
583 }
584
586}
587
589 int side, EntityType type, EntitiesFieldData::EntData &data) {
591
592 if (type == MBVERTEX) {
593
594 VectorDouble &coords = getCoords();
595 double *coords_ptr = &*coords.data().begin();
596
597 const int nb_gauss_pts = data.getN(NOBASE).size1();
598 auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
599
603
604 auto t_w = getFTensor0IntegrationWeight();
605 for (int gg = 0; gg != nb_gauss_pts; gg++) {
606
607 FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
608 &coords_ptr[2], 3);
609 t_jac(i, j) = 0;
610 for (int bb = 0; bb != 6; bb++) {
611 t_jac(i, j) += t_coords(i) * t_diff_n(j);
612 ++t_diff_n;
613 ++t_coords;
614 }
615
616 double det;
617 CHKERR determinantTensor3by3(t_jac, det);
618 t_w *= det / 2.;
619
620 ++t_w;
621 }
622
623 double &vol = getVolume();
624 auto t_w_scaled = getFTensor0IntegrationWeight();
625 for (int gg = 0; gg != nb_gauss_pts; gg++) {
626 t_w_scaled /= vol;
627 ++t_w_scaled;
628 }
629 }
630
631 doEntities[MBVERTEX] = true;
632 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
633
635}
636
641
642 if (type == MBVERTEX) {
643
644 VectorDouble &coords = getCoords();
645 double *coords_ptr = &*coords.data().begin();
646
647 const int nb_gauss_pts = data.getN(NOBASE).size1();
648 auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
649 invJac.resize(9, nb_gauss_pts, false);
650 invJac.clear();
651 auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
652
656
657 auto t_w = getFTensor0IntegrationWeight();
658 for (int gg = 0; gg != nb_gauss_pts; gg++) {
659
660 FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
661 &coords_ptr[2], 3);
662 t_jac(i, j) = 0;
663 for (int bb = 0; bb != 6; bb++) {
664 t_jac(i, j) += t_coords(i) * t_diff_n(j);
665 ++t_diff_n;
666 ++t_coords;
667 }
668
669 double det;
670 CHKERR determinantTensor3by3(t_jac, det);
671 CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
672 ++t_inv_jac;
673 }
674 }
675
676 doEntities[MBVERTEX] = true;
677 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
678
680}
681
686
687 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
688
689 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
690 if (data.getN(base).size2() == 0)
691 continue;
692
693 const int nb_gauss_pts = data.getN(base).size1();
694 auto t_diff_n = data.getFTensor1DiffN<3>(base);
695 diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
696 false);
697
700
702 &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
703 auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
704
705 const int nb_dofs = data.getN(base).size2();
706 for (int gg = 0; gg != nb_gauss_pts; gg++) {
707 for (int bb = 0; bb != nb_dofs; bb++) {
708 t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
709 ++t_inv_diff_n;
710 ++t_diff_n;
711 }
712 ++t_inv_jac;
713 }
714
715 data.getDiffN(base).swap(diffNinvJac);
716 }
717
719}
720
724
726
727 if (type == MBVERTEX) {
728
729 VectorDouble &coords = getCoords();
730 double *coords_ptr = &*coords.data().begin();
731 double diff_n[6];
732 CHKERR ShapeDiffMBTRI(diff_n);
733 double j00_f3, j01_f3, j10_f3, j11_f3;
734 for (int gg = 0; gg < 1; gg++) {
735 // this is triangle, derivative of nodal shape functions is constant.
736 // So only need to do one node.
737 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
738 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
739 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
740 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
741 }
742 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
743 invJacF3.resize(2, 2, false);
744 invJacF3(0, 0) = j11_f3 / det_f3;
745 invJacF3(0, 1) = -j01_f3 / det_f3;
746 invJacF3(1, 0) = -j10_f3 / det_f3;
747 invJacF3(1, 1) = j00_f3 / det_f3;
748 }
749
750 doEntities[MBVERTEX] = true;
751 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
752
754}
755
760
761 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
762
763 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
764
765 unsigned int nb_dofs = data.getN(base).size2();
766 if (nb_dofs == 0)
768 unsigned int nb_gauss_pts = data.getN(base).size1();
769 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
770
771#ifndef NDEBUG
772 if (type != MBVERTEX) {
773 if (nb_dofs != data.getDiffN(base).size2() / 2) {
774 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
775 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
776 "%u/2 )",
777 nb_dofs, data.getDiffN(base).size2());
778 }
779 }
780#endif
781
782 switch (type) {
783 case MBVERTEX:
784 case MBEDGE:
785 case MBTRI: {
786 for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
787 for (unsigned int dd = 0; dd < nb_dofs; dd++) {
788 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
789 &*invJacF3.data().begin(), 2,
790 &data.getDiffN(base)(gg, 2 * dd), 1, 0,
791 &diffNinvJac(gg, 2 * dd), 1);
792 }
793 }
794 data.getDiffN(base).swap(diffNinvJac);
795 } break;
796 default:
797 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
798 }
799 }
800
802}
803
805 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
806 const EntityType zero_type, const int zero_side)
809 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
810 if (!dataPtr)
811 THROW_MESSAGE("Pointer is not set");
812}
813
818 const auto nb_integration_points = getGaussPts().size2();
819 if (type == zeroType && side == zeroSide) {
820 dataPtr->resize(3, nb_integration_points, false);
821 dataPtr->clear();
822 }
823 const auto nb_dofs = data.getFieldData().size();
824 if (!nb_dofs)
829 const auto nb_base_functions = data.getN().size2() / 3;
830 auto t_n_diff_hcurl = data.getFTensor2DiffN<3, 3>();
831 auto t_data = getFTensor1FromMat<3>(*dataPtr);
832 for (auto gg = 0; gg != nb_integration_points; ++gg) {
833 auto t_dof = data.getFTensor0FieldData();
834 int bb = 0;
835 for (; bb != nb_dofs; ++bb) {
836 t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
837 ++t_n_diff_hcurl;
838 ++t_dof;
839 }
840 for (; bb < nb_base_functions; ++bb)
841 ++t_n_diff_hcurl;
842 ++t_data;
843 }
844
846}
847
849 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
850 const EntityType zero_type, const int zero_side)
853 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
854 if (!dataPtr)
855 THROW_MESSAGE("Pointer is not set");
856}
857
862 const auto nb_integration_points = getGaussPts().size2();
863 if (type == zeroType && side == zeroSide) {
864 dataPtr->resize(2, nb_integration_points, false);
865 dataPtr->clear();
866 }
867 const auto nb_dofs = data.getFieldData().size();
868 if (!nb_dofs)
870
873
874 const auto nb_base_functions = data.getN().size2();
875 auto t_n_diff_hcurl = data.getFTensor1DiffN<2>();
876 auto t_data = getFTensor1FromMat<2>(*dataPtr);
877 for (auto gg = 0; gg != nb_integration_points; ++gg) {
878 auto t_dof = data.getFTensor0FieldData();
879 int bb = 0;
880 for (; bb != nb_dofs; ++bb) {
881 t_data(i) += t_dof * (levi_civita(i, j) * t_n_diff_hcurl(j));
882 ++t_n_diff_hcurl;
883 ++t_dof;
884 }
885 for (; bb < nb_base_functions; ++bb)
886 ++t_n_diff_hcurl;
887 ++t_data;
888 }
889
891}
892
893} // namespace MoFEM
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ NOBASE
Definition: definitions.h:59
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
@ HVEC1_1
Definition: definitions.h:196
@ HVEC0_1
Definition: definitions.h:195
@ HVEC1_0
Definition: definitions.h:193
@ HVEC2_1
Definition: definitions.h:197
@ HVEC1_2
Definition: definitions.h:199
@ HVEC2_2
Definition: definitions.h:200
@ HVEC2_0
Definition: definitions.h:194
@ HVEC0_2
Definition: definitions.h:198
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
EntityType zero_type
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1393
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
Definition: Templates.hpp:884
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:1363
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:862
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
bool & doVertices
\deprectaed If false skip vertices
MatrixDouble & getTangentAtGaussPts()
get tangent vector to edge curve at integration points
static FTensor::Tensor1< double, 3 > tFaceOrientation
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of direvarives base function for BB base, key is a field name
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
auto getFTensor0IntegrationWeight()
Get integration weights.
int getFEDim() const
Get dimension of finite element.
structure to get information form mofem into EntitiesFieldData
Calculate curl of vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Apply contravariant (Piola) transfer to Hdiv space on face.
Transform Hcurl base fluxes from reference element to physical triangle.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Transform local reference derivatives of shape function to global derivatives for face.
Transform local reference derivatives of shape functions to global derivatives.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)