v0.14.0
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
340
341 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
342
343 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
344
345 auto &baseN = data.getN(base);
346 auto &diffBaseN = data.getDiffN(base);
347
348 int nb_dofs = baseN.size2() / 3;
349 int nb_gauss_pts = baseN.size1();
350
351 piolaN.resize(baseN.size1(), baseN.size2());
352
353 if (nb_dofs > 0 && nb_gauss_pts > 0) {
354
355 auto t_h_curl = data.getFTensor1N<3>(base);
356 auto t_transformed_h_curl =
357 getFTensor1FromPtr<3, 3>(&*piolaN.data().begin());
358 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
359 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
360 for (int ll = 0; ll != nb_dofs; ll++) {
361 t_transformed_h_curl(i) = t_inv_jac(j, i) * t_h_curl(j);
362 ++t_h_curl;
363 ++t_transformed_h_curl;
364 }
365 ++t_inv_jac;
366 }
367 baseN.swap(piolaN);
368
369 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
370 if (diffBaseN.data().size() > 0) {
371 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
372 auto t_transformed_diff_h_curl =
373 getFTensor2FromPtr<3, 2>(&*diffPiolaN.data().begin());
374 auto t_inv_jac = getFTensor2FromMat<2, 2>(*invJacPtr);
375 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
376 for (int ll = 0; ll != nb_dofs; ll++) {
377 t_transformed_diff_h_curl(i, k) =
378 t_inv_jac(j, i) * t_diff_h_curl(j, k);
379 ++t_diff_h_curl;
380 ++t_transformed_diff_h_curl;
381 }
382 ++t_inv_jac;
383 }
384 diffBaseN.swap(diffPiolaN);
385 }
386 }
387 }
388
390}
391
393 int side, EntityType type, EntitiesFieldData::EntData &data) {
394
396
397 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
399
403
404 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
405
406 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
407
408 const size_t nb_base_functions = data.getN(base).size2() / 3;
409 if (nb_base_functions) {
410
411 const size_t nb_gauss_pts = data.getN(base).size1();
412 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
413 if (data.getN(base).size2() > 0) {
414 auto t_n = data.getFTensor1N<3>(base);
415 double *t_transformed_n_ptr = &*piolaN.data().begin();
417 t_transformed_n_ptr, // HVEC0
418 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
419 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
420 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
421 double det;
422 CHKERR determinantTensor2by2(t_jac, det);
423 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
424 t_transformed_n(i) = t_jac(i, k) * t_n(k) / det;
425 ++t_n;
426 ++t_transformed_n;
427 }
428 }
429 data.getN(base).swap(piolaN);
430 }
431
432 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
433 if (data.getDiffN(base).data().size() > 0) {
434 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
435 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
437 t_transformed_diff_n(t_transformed_diff_n_ptr,
438 &t_transformed_diff_n_ptr[HVEC0_1],
439
440 &t_transformed_diff_n_ptr[HVEC1_0],
441 &t_transformed_diff_n_ptr[HVEC1_1],
442
443 &t_transformed_diff_n_ptr[HVEC2_0],
444 &t_transformed_diff_n_ptr[HVEC2_1]);
445 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
446 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
447 double det;
448 CHKERR determinantTensor2by2(t_jac, det);
449 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
450 t_transformed_diff_n(i, k) = t_jac(i, j) * t_diff_n(j, k) / det;
451 ++t_diff_n;
452 ++t_transformed_diff_n;
453 }
454 }
455 data.getDiffN(base).swap(piolaDiffN);
456 }
457 }
458 }
459
461}
462
464 int side, EntityType type, EntitiesFieldData::EntData &data) {
465
467
468 if (type != MBEDGE && type != MBTRI && type != MBQUAD)
470
474
475 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
476
477 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
478
479 const size_t nb_base_functions = data.getN(base).size2() / 3;
480 if (nb_base_functions) {
481
482 const size_t nb_gauss_pts = data.getN(base).size1();
483 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
484 if (data.getN(base).size2() > 0) {
485 auto t_n = data.getFTensor1N<3>(base);
486 double *t_transformed_n_ptr = &*piolaN.data().begin();
488 t_transformed_n_ptr, // HVEC0
489 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
490 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
491 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
492 double det;
493 CHKERR determinantTensor3by3(t_jac, det);
494 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
495 t_transformed_n(i) = t_jac(i, j) * t_n(j) / det;
496 ++t_n;
497 ++t_transformed_n;
498 }
499 }
500 data.getN(base).swap(piolaN);
501 }
502
503 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
504 if (data.getDiffN(base).size2() > 0) {
505 auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
506 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
508 t_transformed_diff_n(t_transformed_diff_n_ptr,
509 &t_transformed_diff_n_ptr[HVEC0_1],
510
511 &t_transformed_diff_n_ptr[HVEC1_0],
512 &t_transformed_diff_n_ptr[HVEC1_1],
513
514 &t_transformed_diff_n_ptr[HVEC2_0],
515 &t_transformed_diff_n_ptr[HVEC2_1]);
516
517 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
518 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
519 double det;
520 CHKERR determinantTensor3by3(t_jac, det);
521 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
522 t_transformed_diff_n(i, K) = t_jac(i, j) * t_diff_n(j, K) / det;
523 ++t_diff_n;
524 ++t_transformed_diff_n;
525 }
526 }
527 data.getDiffN(base).swap(piolaDiffN);
528 }
529 }
530 }
531
533}
534
536 int side, EntityType type, EntitiesFieldData::EntData &data) {
538
539 if (type != MBEDGE)
541
545
546 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
547
548 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
549 size_t nb_gauss_pts = data.getN(base).size1();
550 size_t nb_dofs = data.getN(base).size2() / 3;
551 if (nb_dofs > 0) {
552
553 auto t_h_div = data.getFTensor1N<3>(base);
554
555 size_t cc = 0;
556 const auto &edge_direction_at_gauss_pts = getTangentAtGaussPts();
557 if (edge_direction_at_gauss_pts.size1() == nb_gauss_pts) {
558
559 auto t_dir = getFTensor1TangentAtGaussPts<3>();
560
561 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
563 t_normal(i) = (FTensor::levi_civita(i, j, k) *
565 t_dir(k);
566 const auto l2 = t_normal(i) * t_normal(i);
567 for (int ll = 0; ll != nb_dofs; ++ll) {
568 const double val = t_h_div(0);
569 const double a = val / l2;
570 t_h_div(i) = t_normal(i) * a;
571 ++t_h_div;
572 ++cc;
573 }
574 ++t_dir;
575 }
576 }
577
578 if (cc != nb_gauss_pts * nb_dofs)
579 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
580 }
581 }
582
584}
585
587 int side, EntityType type, EntitiesFieldData::EntData &data) {
589
590 if (type == MBVERTEX) {
591
592 VectorDouble &coords = getCoords();
593 double *coords_ptr = &*coords.data().begin();
594
595 const int nb_gauss_pts = data.getN(NOBASE).size1();
596 auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
597
601
602 auto t_w = getFTensor0IntegrationWeight();
603 for (int gg = 0; gg != nb_gauss_pts; gg++) {
604
605 FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
606 &coords_ptr[2], 3);
607 t_jac(i, j) = 0;
608 for (int bb = 0; bb != 6; bb++) {
609 t_jac(i, j) += t_coords(i) * t_diff_n(j);
610 ++t_diff_n;
611 ++t_coords;
612 }
613
614 double det;
615 CHKERR determinantTensor3by3(t_jac, det);
616 t_w *= det / 2.;
617
618 ++t_w;
619 }
620
621 double &vol = getVolume();
622 auto t_w_scaled = getFTensor0IntegrationWeight();
623 for (int gg = 0; gg != nb_gauss_pts; gg++) {
624 t_w_scaled /= vol;
625 ++t_w_scaled;
626 }
627 }
628
629 doEntities[MBVERTEX] = true;
630 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
631
633}
634
639
640 if (type == MBVERTEX) {
641
642 VectorDouble &coords = getCoords();
643 double *coords_ptr = &*coords.data().begin();
644
645 const int nb_gauss_pts = data.getN(NOBASE).size1();
646 auto t_diff_n = data.getFTensor1DiffN<3>(NOBASE);
647 invJac.resize(9, nb_gauss_pts, false);
648 invJac.clear();
649 auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
650
654
655 for (int gg = 0; gg != nb_gauss_pts; gg++) {
656
657 FTensor::Tensor1<double *, 3> t_coords(coords_ptr, &coords_ptr[1],
658 &coords_ptr[2], 3);
659 t_jac(i, j) = 0;
660 for (int bb = 0; bb != 6; bb++) {
661 t_jac(i, j) += t_coords(i) * t_diff_n(j);
662 ++t_diff_n;
663 ++t_coords;
664 }
665
666 double det;
667 CHKERR determinantTensor3by3(t_jac, det);
668 CHKERR invertTensor3by3(t_jac, det, t_inv_jac);
669 ++t_inv_jac;
670 }
671 }
672
673 doEntities[MBVERTEX] = true;
674 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
675
677}
678
683
684 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
685
686 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
687 if (data.getN(base).size2() == 0)
688 continue;
689
690 const int nb_gauss_pts = data.getN(base).size1();
691 auto t_diff_n = data.getFTensor1DiffN<3>(base);
692 diffNinvJac.resize(data.getDiffN(base).size1(), data.getDiffN(base).size2(),
693 false);
694
697
699 &diffNinvJac(0, 0), &diffNinvJac(0, 1), &diffNinvJac(0, 2));
700 auto t_inv_jac = getFTensor2FromMat<3, 3>(invJac);
701
702 const int nb_dofs = data.getN(base).size2();
703 for (int gg = 0; gg != nb_gauss_pts; gg++) {
704 for (int bb = 0; bb != nb_dofs; bb++) {
705 t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
706 ++t_inv_diff_n;
707 ++t_diff_n;
708 }
709 ++t_inv_jac;
710 }
711
712 data.getDiffN(base).swap(diffNinvJac);
713 }
714
716}
717
721
723
724 if (type == MBVERTEX) {
725
726 VectorDouble &coords = getCoords();
727 double *coords_ptr = &*coords.data().begin();
728 double diff_n[6];
729 CHKERR ShapeDiffMBTRI(diff_n);
730 double j00_f3, j01_f3, j10_f3, j11_f3;
731 for (int gg = 0; gg < 1; gg++) {
732 // this is triangle, derivative of nodal shape functions is constant.
733 // So only need to do one node.
734 j00_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[0], 2);
735 j01_f3 = cblas_ddot(3, &coords_ptr[0], 3, &diff_n[1], 2);
736 j10_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[0], 2);
737 j11_f3 = cblas_ddot(3, &coords_ptr[1], 3, &diff_n[1], 2);
738 }
739 double det_f3 = j00_f3 * j11_f3 - j01_f3 * j10_f3;
740 invJacF3.resize(2, 2, false);
741 invJacF3(0, 0) = j11_f3 / det_f3;
742 invJacF3(0, 1) = -j01_f3 / det_f3;
743 invJacF3(1, 0) = -j10_f3 / det_f3;
744 invJacF3(1, 1) = j00_f3 / det_f3;
745 }
746
747 doEntities[MBVERTEX] = true;
748 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
749
751}
752
757
758 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
759
760 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
761
762 unsigned int nb_dofs = data.getN(base).size2();
763 if (nb_dofs == 0)
765 unsigned int nb_gauss_pts = data.getN(base).size1();
766 diffNinvJac.resize(nb_gauss_pts, 2 * nb_dofs, false);
767
768#ifndef NDEBUG
769 if (type != MBVERTEX) {
770 if (nb_dofs != data.getDiffN(base).size2() / 2) {
771 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
772 "data inconsistency nb_dofs != data.diffN.size2()/2 ( %u != "
773 "%u/2 )",
774 nb_dofs, data.getDiffN(base).size2());
775 }
776 }
777#endif
778
779 switch (type) {
780 case MBVERTEX:
781 case MBEDGE:
782 case MBTRI: {
783 for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
784 for (unsigned int dd = 0; dd < nb_dofs; dd++) {
785 cblas_dgemv(CblasRowMajor, CblasTrans, 2, 2, 1,
786 &*invJacF3.data().begin(), 2,
787 &data.getDiffN(base)(gg, 2 * dd), 1, 0,
788 &diffNinvJac(gg, 2 * dd), 1);
789 }
790 }
791 data.getDiffN(base).swap(diffNinvJac);
792 } break;
793 default:
794 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
795 }
796 }
797
799}
800
802 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
803 const EntityType zero_type, const int zero_side)
806 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
807 if (!dataPtr)
808 THROW_MESSAGE("Pointer is not set");
809}
810
815 const auto nb_integration_points = getGaussPts().size2();
816 if (type == zeroType && side == zeroSide) {
817 dataPtr->resize(3, nb_integration_points, false);
818 dataPtr->clear();
819 }
820 const auto nb_dofs = data.getFieldData().size();
821 if (!nb_dofs)
826 const auto nb_base_functions = data.getN().size2() / 3;
827 auto t_n_diff_hcurl = data.getFTensor2DiffN<3, 3>();
828 auto t_data = getFTensor1FromMat<3>(*dataPtr);
829 for (auto gg = 0; gg != nb_integration_points; ++gg) {
830 auto t_dof = data.getFTensor0FieldData();
831 int bb = 0;
832 for (; bb != nb_dofs; ++bb) {
833 t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
834 ++t_n_diff_hcurl;
835 ++t_dof;
836 }
837 for (; bb < nb_base_functions; ++bb)
838 ++t_n_diff_hcurl;
839 ++t_data;
840 }
841
843}
844
846 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
847 const EntityType zero_type, const int zero_side)
850 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
851 if (!dataPtr)
852 THROW_MESSAGE("Pointer is not set");
853}
854
859 const auto nb_integration_points = getGaussPts().size2();
860 if (type == zeroType && side == zeroSide) {
861 dataPtr->resize(2, nb_integration_points, false);
862 dataPtr->clear();
863 }
864 const auto nb_dofs = data.getFieldData().size();
865 if (!nb_dofs)
867
870
871 const auto nb_base_functions = data.getN().size2();
872 auto t_n_diff_hcurl = data.getFTensor1DiffN<2>();
873 auto t_data = getFTensor1FromMat<2>(*dataPtr);
874 for (auto gg = 0; gg != nb_integration_points; ++gg) {
875 auto t_dof = data.getFTensor0FieldData();
876 int bb = 0;
877 for (; bb != nb_dofs; ++bb) {
878 t_data(i) += t_dof * (levi_civita(i, j) * t_n_diff_hcurl(j));
879 ++t_n_diff_hcurl;
880 ++t_dof;
881 }
882 for (; bb < nb_base_functions; ++bb)
883 ++t_n_diff_hcurl;
884 ++t_data;
885 }
886
888}
889
890} // 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:1534
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
Definition: Templates.hpp:970
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:1504
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1493
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:948
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)