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
8namespace 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++) {
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++) {
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
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
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition: Templates.hpp:1524
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:921
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:1559
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
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)