v0.13.2
Loading...
Searching...
No Matches
HODataOperators.cpp
Go to the documentation of this file.
1/** \file HODataOperators.cpp
2
3\brief Set of operators for high-order geometry approximation.
4
5*/
6
7namespace MoFEM {
8
10 boost::shared_ptr<MatrixDouble> jac_ptr)
12 jacPtr(jac_ptr) {
13
14 for (auto t = MBEDGE; t != MBMAXTYPE; ++t)
15 doEntities[t] = false;
16
17 if (!jacPtr)
18 THROW_MESSAGE("Jac pointer not allocated");
19}
20
25
26 const auto nb_base_functions = data.getN(NOBASE).size2();
27 if (nb_base_functions) {
28
29 const auto nb_gauss_pts = getGaussPts().size2();
30 auto t_diff_base = data.getFTensor1DiffN<3>(NOBASE);
31 auto &coords = getCoords();
32
33#ifndef NDEBUG
34 if (nb_gauss_pts != data.getDiffN().size1())
35 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
36 "Inconsistent number base functions and gauss points");
37 if (nb_base_functions != data.getDiffN().size2() / 3)
38 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
39 "Inconsistent number of base functions");
40 if (coords.size() != 3 * nb_base_functions)
41 SETERRQ2(
42 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
43 "Number of vertex coordinates and base functions is inconsistent "
44 "%d != %d",
45 coords.size() / 3, nb_base_functions);
46#endif
47
48 jacPtr->resize(9, nb_gauss_pts, false);
49 jacPtr->clear();
50 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
51 auto t_vol_inv_jac = getInvJac();
52
56
57 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
59 &coords[0], &coords[1], &coords[2]);
60 for (size_t bb = 0; bb != nb_base_functions; ++bb) {
61 t_jac(i, j) += t_coords(i) * (t_vol_inv_jac(k, j) * t_diff_base(k));
62 ++t_diff_base;
63 ++t_coords;
64 }
65 ++t_jac;
66 }
67 }
68
70}
71
76
77 if (getFEDim() == 3) {
78
79 auto transform_base = [&](MatrixDouble &diff_n) {
81 return applyTransform<3, 3, 3, 3>(diff_n);
83 };
84
85 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
87 CHKERR transform_base(data.getDiffN(base));
88 }
89
90 switch (type) {
91 case MBVERTEX:
92 for (auto &m : data.getBBDiffNMap())
93 CHKERR transform_base(*(m.second));
94 break;
95 default:
96 for (auto &ptr : data.getBBDiffNByOrderArray())
97 if (ptr)
98 CHKERR transform_base(*ptr);
99 }
100
101 } else {
102 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Use different operator");
103 }
104
106}
107
112
116
117 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
118
119 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
120
121 diffHdivInvJac.resize(data.getDiffN(base).size1(),
122 data.getDiffN(base).size2(), false);
123
124 unsigned int nb_gauss_pts = data.getDiffN(base).size1();
125 unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
126 if (nb_base_functions == 0)
127 continue;
128
129 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
130 double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
132 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
133 &t_inv_diff_n_ptr[HVEC0_2], &t_inv_diff_n_ptr[HVEC1_0],
134 &t_inv_diff_n_ptr[HVEC1_1], &t_inv_diff_n_ptr[HVEC1_2],
135 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
136 &t_inv_diff_n_ptr[HVEC2_2]);
137 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
138 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
139 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
140 t_inv_diff_n(i, j) = t_inv_jac(k, j) * t_diff_n(i, k);
141 ++t_diff_n;
142 ++t_inv_diff_n;
143 }
144 ++t_inv_jac;
145 }
146
147 data.getDiffN(base).swap(diffHdivInvJac);
148 }
149
151}
152
156 const size_t nb_int_pts = getGaussPts().size2();
157 if (getNormalsAtGaussPts().size1()) {
158 if (getNormalsAtGaussPts().size1() == nb_int_pts) {
159 double a = getMeasure();
160 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI)
161 a *= 2;
162 auto t_w = getFTensor0IntegrationWeight();
163 auto t_normal = getFTensor1NormalsAtGaussPts();
165 for (size_t gg = 0; gg != nb_int_pts; ++gg) {
166 t_w *= sqrt(t_normal(i) * t_normal(i)) / a;
167 ++t_w;
168 ++t_normal;
169 }
170 } else {
171 SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
172 "Number of rows in getNormalsAtGaussPts should be equal to "
173 "number of integration points, but is not, i.e. %d != %d",
174 getNormalsAtGaussPts().size1(), nb_int_pts);
175 }
176 }
178}
179
183
184 const auto nb_integration_pts = detPtr->size();
185
186#ifndef NDEBUG
187 if (nb_integration_pts != getGaussPts().size2())
188 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
189 "Inconsistent number of data points");
190#endif
191
192 auto t_w = getFTensor0IntegrationWeight();
193 auto t_det = getFTensor0FromVec(*detPtr);
194 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
195 t_w *= t_det;
196 ++t_w;
197 ++t_det;
198 }
199
201}
202
207
211
212#ifndef NDEBUG
213 if (!detPtr)
214 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
215 "Pointer for detPtr not allocated");
216 if (!jacPtr)
217 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
218 "Pointer for jacPtr not allocated");
219#endif
220
221 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
222
223 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
224
225 auto nb_gauss_pts = data.getN(base).size1();
226 auto nb_base_functions = data.getN(base).size2() / 3;
227
228#ifndef NDEBUG
229 if (data.getDiffN(base).size1() != nb_gauss_pts)
230 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
231 "Wrong number integration points");
232
233 if (data.getDiffN(base).size2() / 9 != nb_base_functions)
234 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
235 "Wrong number base functions %d != %d",
236 data.getDiffN(base).size2(), nb_base_functions);
237#endif
238
239 if (nb_gauss_pts && nb_base_functions) {
240
241 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
242 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
243
244 auto t_n = data.getFTensor1N<3>(base);
245 double *t_transformed_n_ptr = &*piolaN.data().begin();
247 t_transformed_n_ptr, // HVEC0
248 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
249 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
250 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
252 t_transformed_diff_n(t_transformed_diff_n_ptr,
253 &t_transformed_diff_n_ptr[HVEC0_1],
254 &t_transformed_diff_n_ptr[HVEC0_2],
255 &t_transformed_diff_n_ptr[HVEC1_0],
256 &t_transformed_diff_n_ptr[HVEC1_1],
257 &t_transformed_diff_n_ptr[HVEC1_2],
258 &t_transformed_diff_n_ptr[HVEC2_0],
259 &t_transformed_diff_n_ptr[HVEC2_1],
260 &t_transformed_diff_n_ptr[HVEC2_2]);
261
262 auto t_det = getFTensor0FromVec(*detPtr);
263 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
264
265 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
266 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
267 const double a = 1. / t_det;
268 t_transformed_n(i) = a * t_jac(i, k) * t_n(k);
269 t_transformed_diff_n(i, k) = a * t_jac(i, j) * t_diff_n(j, k);
270 ++t_n;
271 ++t_transformed_n;
272 ++t_diff_n;
273 ++t_transformed_diff_n;
274 }
275 ++t_det;
276 ++t_jac;
277 }
278
279 data.getN(base).swap(piolaN);
280 data.getDiffN(base).swap(piolaDiffN);
281 }
282 }
283
285}
286
291
295
296 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
297
298 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
299
300 unsigned int nb_gauss_pts = data.getN(base).size1();
301 unsigned int nb_base_functions = data.getN(base).size2() / 3;
302 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
303 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
304
305 auto t_n = data.getFTensor1N<3>(base);
306 double *t_transformed_n_ptr = &*piolaN.data().begin();
308 t_transformed_n_ptr, // HVEC0
309 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
310 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
311 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
312 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
313 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
314 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
315 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
316 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
317 &t_transformed_diff_n_ptr[HVEC2_2]);
318
319 auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
320
321 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
322 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
323 t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
324 t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
325 ++t_n;
326 ++t_transformed_n;
327 ++t_diff_n;
328 ++t_transformed_diff_n;
329 }
330 ++t_inv_jac;
331 }
332
333 data.getN(base).swap(piolaN);
334 data.getDiffN(base).swap(piolaDiffN);
335 }
336
338}
339
341 boost::shared_ptr<MatrixDouble> jac_ptr)
343 jacPtr(jac_ptr) {}
344
348
350
351 const auto nb_gauss_pts = getGaussPts().size2();
352
353 jacPtr->resize(4, nb_gauss_pts, false);
354 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
355
358
359 auto t_t1 = getFTensor1Tangent1AtGaussPts();
360 auto t_t2 = getFTensor1Tangent2AtGaussPts();
363
364 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
365 t_jac(i, N0) = t_t1(i);
366 t_jac(i, N1) = t_t2(i);
367 ++t_t1;
368 ++t_t2;
369 ++t_jac;
370 }
371
373};
374
379
383
384 size_t nb_gauss_pts = getGaussPts().size2();
385 jacPtr->resize(9, nb_gauss_pts, false);
386
387 auto t_t1 = getFTensor1Tangent1AtGaussPts();
388 auto t_t2 = getFTensor1Tangent2AtGaussPts();
389 auto t_normal = getFTensor1NormalsAtGaussPts();
390
394
395 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
396 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
397
398 t_jac(i, N0) = t_t1(i);
399 t_jac(i, N1) = t_t2(i);
400
401 const double l = sqrt(t_normal(j) * t_normal(j));
402
403 t_jac(i, N2) = t_normal(i) / l;
404
405 ++t_t1;
406 ++t_t2;
407 ++t_normal;
408 }
409
411}
412
414 int side, EntityType type, EntitiesFieldData::EntData &data) {
417
418 if (moab::CN::Dimension(type) != 2)
420
421 auto get_normals_ptr = [&]() {
423 return &*normalsAtGaussPts->data().begin();
424 else
425 return &*getNormalsAtGaussPts().data().begin();
426 };
427
428 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
429 auto nb_base_functions) {
431
432 auto ptr = get_normals_ptr();
434 &ptr[0], &ptr[1], &ptr[2]);
435
436 auto t_base = data.getFTensor1N<3>(base);
437 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
438 const auto l2 = t_normal(i) * t_normal(i);
439 for (int bb = 0; bb != nb_base_functions; ++bb) {
440 const auto v = t_base(0);
441 t_base(i) = (v / l2) * t_normal(i);
442 ++t_base;
443 }
444 ++t_normal;
445 }
447 };
448
449 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
450
451 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
452 const auto &base_functions = data.getN(base);
453 const auto nb_gauss_pts = base_functions.size1();
454
455 if (nb_gauss_pts) {
456
457 const auto nb_base_functions = base_functions.size2() / 3;
458 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
459 nb_base_functions);
460 }
461 }
462
464}
465
467 int side, EntityType type, EntitiesFieldData::EntData &data) {
469
470 const auto type_dim = moab::CN::Dimension(type);
471 if (type_dim != 1 && type_dim != 2)
473
477
478 auto get_jac = [&]() {
480 double *ptr_n = &*normalsAtPts->data().begin();
481 double *ptr_t1 = &*tangent1AtPts->data().begin();
482 double *ptr_t2 = &*tangent2AtPts->data().begin();
484 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
485
486 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
487
488 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
489 } else {
490 double *ptr_n = &*getNormalsAtGaussPts().data().begin();
491 double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
492 double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
494 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
495
496 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
497
498 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
499 }
500 };
501
502 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
503
504 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
505
506 auto &baseN = data.getN(base);
507 auto &diffBaseN = data.getDiffN(base);
508
509 int nb_dofs = baseN.size2() / 3;
510 int nb_gauss_pts = baseN.size1();
511
512 piolaN.resize(baseN.size1(), baseN.size2());
513 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
514
515 if (nb_dofs > 0 && nb_gauss_pts > 0) {
516
517 auto t_h_curl = data.getFTensor1N<3>(base);
518 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
519
520 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
521 &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
522
524 t_transformed_diff_h_curl(
528
529 int cc = 0;
530 double det;
532
533 // HO geometry is set, so jacobian is different at each gauss point
534 auto t_m_at_pts = get_jac();
535 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
536 CHKERR determinantTensor3by3(t_m_at_pts, det);
537 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
538 for (int ll = 0; ll != nb_dofs; ll++) {
539 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
540 t_transformed_diff_h_curl(i, k) = t_inv_m(j, i) * t_diff_h_curl(j, k);
541 ++t_h_curl;
542 ++t_transformed_h_curl;
543 ++t_diff_h_curl;
544 ++t_transformed_diff_h_curl;
545 ++cc;
546 }
547 ++t_m_at_pts;
548 }
549
550#ifndef NDEBUG
551 if (cc != nb_gauss_pts * nb_dofs)
552 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
553#endif
554
555 baseN.swap(piolaN);
556 diffBaseN.swap(diffPiolaN);
557 }
558 }
559
561}
562
564 int side, EntityType type, EntitiesFieldData::EntData &data) {
567
568 if (type != MBEDGE)
570
571 const auto nb_gauss_pts = getGaussPts().size2();
572
574 if (tangentAtGaussPts->size1() != nb_gauss_pts)
575 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
576 "Wrong number of integration points %d!=%d",
577 tangentAtGaussPts->size1(), nb_gauss_pts);
578
579 auto get_base_at_pts = [&](auto base) {
581 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
582 &data.getN(base)(0, HVEC2));
583 return t_h_curl;
584 };
585
586 auto get_tangent_ptr = [&]() {
587 if (tangentAtGaussPts) {
588 return &*tangentAtGaussPts->data().begin();
589 } else {
590 return &*getTangentAtGaussPts().data().begin();
591 }
592 };
593
594 auto get_tangent_at_pts = [&]() {
595 auto ptr = get_tangent_ptr();
597 &ptr[0], &ptr[1], &ptr[2]);
598 return t_m_at_pts;
599 };
600
601 auto calculate_squared_edge_length = [&]() {
602 std::vector<double> l1;
603 if (nb_gauss_pts) {
604 l1.resize(nb_gauss_pts);
605 auto t_m_at_pts = get_tangent_at_pts();
606 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
607 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
608 ++t_m_at_pts;
609 }
610 }
611 return l1;
612 };
613
614 auto l1 = calculate_squared_edge_length();
615
616 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
617
618 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
619 const auto nb_dofs = data.getN(base).size2() / 3;
620
621 if (nb_gauss_pts && nb_dofs) {
622 auto t_h_curl = get_base_at_pts(base);
623 int cc = 0;
624 auto t_m_at_pts = get_tangent_at_pts();
625 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
626 const double l0 = l1[gg];
627 for (int ll = 0; ll != nb_dofs; ++ll) {
628 const double val = t_h_curl(0);
629 const double a = val / l0;
630 t_h_curl(i) = t_m_at_pts(i) * a;
631 ++t_h_curl;
632 ++cc;
633 }
634 ++t_m_at_pts;
635 }
636
637#ifndef NDEBUG
638 if (cc != nb_gauss_pts * nb_dofs)
639 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
640#endif
641 }
642 }
643
645}
646
648 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
649 std::vector<FieldSpace> spaces, std::string geom_field_name,
650 boost::shared_ptr<MatrixDouble> jac_ptr,
651 boost::shared_ptr<VectorDouble> det_ptr,
652 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
654
655 if (!jac_ptr)
656 jac_ptr = boost::make_shared<MatrixDouble>();
657 if (!det_ptr)
658 det_ptr = boost::make_shared<VectorDouble>();
659 if (!inv_jac_ptr)
660 inv_jac_ptr = boost::make_shared<MatrixDouble>();
661
662 if (geom_field_name.empty()) {
663
664 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
665
666 } else {
667
668 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
669 pipeline.push_back(
670 new OpCalculateVectorFieldGradient<2, 2>(geom_field_name, jac_ptr));
671 pipeline.push_back(new OpGetHONormalsOnFace<2>(geom_field_name));
672 }
673
674 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
675 pipeline.push_back(new OpSetHOWeightsOnFace());
676
677 for (auto s : spaces) {
678 switch (s) {
679 case NOSPACE:
680 break;
681 case H1:
682 case L2:
683 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(s, inv_jac_ptr));
684 break;
685 case HCURL:
686 pipeline.push_back(new OpSetCovariantPiolaTransformOnFace2D(inv_jac_ptr));
687 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
688 break;
689 case HDIV:
690 pipeline.push_back(new OpMakeHdivFromHcurl());
691 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
692 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
693 break;
694 default:
695 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
696 "Space %s not yet implemented", FieldSpaceNames[s]);
697 }
698 }
699
701}
702
704 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
705 std::vector<FieldSpace> spaces, std::string geom_field_name) {
707
708 if (geom_field_name.empty()) {
709
710 } else {
711
712 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
713 pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
714 }
715
716 for (auto s : spaces) {
717 switch (s) {
718 case NOSPACE:
719 break;
720 case HCURL:
721 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
722 break;
723 case HDIV:
724 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
725 break;
726 default:
727 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
728 "Space %s not yet implemented", FieldSpaceNames[s]);
729 }
730 }
731
733}
734
736 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
737 std::vector<FieldSpace> spaces, std::string geom_field_name) {
739
740 if (geom_field_name.empty()) {
741
742 } else {
743
744 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
745 pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
746 }
747
748 for (auto s : spaces) {
749 switch (s) {
750 case HCURL:
751 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
752 break;
753 default:
754 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
755 "Space %s not yet implemented", FieldSpaceNames[s]);
756 }
757 }
758
760}
761
763 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
764 std::vector<FieldSpace> spaces, std::string geom_field_name) {
766
767 if (geom_field_name.empty()) {
768 } else {
769
770 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
771 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
772 }
773
774 for (auto s : spaces) {
775 switch (s) {
776 case NOSPACE:
777 break;
778 case HCURL:
779 pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
780 break;
781 case HDIV:
782 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
783 break;
784 default:
785 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
786 "Space %s not yet implemented", FieldSpaceNames[s]);
787 break;
788 }
789 }
790
792}
793
795 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
796 std::vector<FieldSpace> spaces, std::string geom_field_name,
797 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
798 boost::shared_ptr<MatrixDouble> inv_jac) {
800
801 if (!jac)
802 jac = boost::make_shared<MatrixDouble>();
803 if (!det)
804 det = boost::make_shared<VectorDouble>();
805 if (!inv_jac)
806 inv_jac = boost::make_shared<MatrixDouble>();
807
808 if (geom_field_name.empty()) {
809
810 pipeline.push_back(new OpCalculateHOJac<3>(jac));
811
812 } else {
813
814 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
815 pipeline.push_back(
816 new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
817 }
818
819 pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
820 pipeline.push_back(new OpSetHOWeights(det));
821
822 for (auto s : spaces) {
823 switch (s) {
824 case NOSPACE:
825 break;
826 case H1:
827 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
828 break;
829 case HCURL:
830 pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
831 pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
832 break;
833 case HDIV:
834 pipeline.push_back(
836 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac));
837 break;
838 case L2:
839 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
840 break;
841 default:
842 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
843 "Space %s not yet implemented", FieldSpaceNames[s]);
844 break;
845 }
846 }
847
849}
850
851} // namespace MoFEM
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
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
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
static const char *const FieldSpaceNames[]
Definition: definitions.h:92
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ HVEC0
Definition: definitions.h:186
@ 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
@ HVEC0_0
Definition: definitions.h:192
#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
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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
OpSetCovariantPiolaTransformOnFace2DImpl< 2 > OpSetCovariantPiolaTransformOnFace2D
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
constexpr double t
plate stiffness
Definition: plate.cpp:59
Add operators pushing bases from local to physical configuration.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
MatrixDouble & getTangentAtGaussPts()
get tangent vector to edge curve at integration points
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....
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::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MatrixDouble & getTangent2AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
MatrixDouble & getTangent1AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
int getFEDim() const
Get dimension of finite element.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Calculate HO coordinates at gauss points.
Calculate jacobian for face element.
boost::shared_ptr< MatrixDouble > jacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHOJacForVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Calculate tangent vector on edge form HO geometry approximation.
transform Hcurl base fluxes from reference element to physical edge
boost::shared_ptr< MatrixDouble > tangentAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
transform Hdiv base fluxes from reference element to physical triangle
boost::shared_ptr< MatrixDouble > normalsAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
transform Hcurl base fluxes from reference element to physical triangle
boost::shared_ptr< MatrixDouble > normalsAtPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > tangent2AtPts
boost::shared_ptr< MatrixDouble > tangent1AtPts
Make Hdiv space from Hcurl space in 2d.
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
boost::shared_ptr< VectorDouble > detPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > jacPtr
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > jacInvPtr
Set inverse jacobian to base functions.
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 if higher order geometr...
boost::shared_ptr< MatrixDouble > invJacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Set inverse jacobian to base functions.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > detPtr
Modify integration weights on face to take in account higher-order geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian