v0.14.0
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 const size_t nb_int_pts = getGaussPts().size2();
184 if (getTangentAtGaussPts().size1()) {
185 if (getTangentAtGaussPts().size1() == nb_int_pts) {
186 double a = getMeasure();
187 auto t_w = getFTensor0IntegrationWeight();
188 auto t_tangent = getFTensor1TangentAtGaussPts();
190 for (size_t gg = 0; gg != nb_int_pts; ++gg) {
191 t_w *= sqrt(t_tangent(i) * t_tangent(i)) / a;
192 ++t_w;
193 ++t_tangent;
194 }
195 } else {
196 SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
197 "Number of rows in getTangentAtGaussPts should be equal to "
198 "number of integration points, but is not, i.e. %d != %d",
199 getTangentAtGaussPts().size1(), nb_int_pts);
200 }
201 }
203}
204
208
209 const auto nb_integration_pts = detPtr->size();
210
211#ifndef NDEBUG
212 if (nb_integration_pts != getGaussPts().size2())
213 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
214 "Inconsistent number of data points");
215#endif
216
217 auto t_w = getFTensor0IntegrationWeight();
218 auto t_det = getFTensor0FromVec(*detPtr);
219 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
220 t_w *= t_det;
221 ++t_w;
222 ++t_det;
223 }
224
226}
227
232
236
237#ifndef NDEBUG
238 if (!detPtr)
239 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
240 "Pointer for detPtr not allocated");
241 if (!jacPtr)
242 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
243 "Pointer for jacPtr not allocated");
244#endif
245
246 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
247
248 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
249
250 auto nb_gauss_pts = data.getNSharedPtr(base) ? data.getN(base).size1() : 0;
251 auto nb_base_functions =
252 data.getNSharedPtr(base) ? data.getN(base).size2() / 3 : 0;
253 auto nb_diff_base_functions =
254 data.getDiffNSharedPtr(base) ? data.getDiffN(base).size2() / 9 : 0;
255
256#ifndef NDEBUG
257 if (nb_diff_base_functions) {
258 if (nb_diff_base_functions != nb_base_functions)
259 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
260 "Wrong number base functions %d != %d", nb_diff_base_functions,
261 nb_base_functions);
262 if (data.getDiffN(base).size1() != nb_gauss_pts)
263 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
264 "Wrong number integration points");
265 }
266#endif
267
268 if (nb_gauss_pts && nb_base_functions) {
269
270 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
271
272 auto t_n = data.getFTensor1N<3>(base);
273 double *t_transformed_n_ptr = &*piolaN.data().begin();
275 t_transformed_n_ptr, // HVEC0
276 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
277
278 auto t_det = getFTensor0FromVec(*detPtr);
279 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
280
281 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
282 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
283 const double a = 1. / t_det;
284 t_transformed_n(i) = a * t_jac(i, k) * t_n(k);
285 ++t_n;
286 ++t_transformed_n;
287 }
288 ++t_det;
289 ++t_jac;
290 }
291
292 data.getN(base).swap(piolaN);
293 }
294
295 if (nb_gauss_pts && nb_diff_base_functions) {
296
297 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
298
299 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
300 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
302 t_transformed_diff_n(t_transformed_diff_n_ptr,
303 &t_transformed_diff_n_ptr[HVEC0_1],
304 &t_transformed_diff_n_ptr[HVEC0_2],
305 &t_transformed_diff_n_ptr[HVEC1_0],
306 &t_transformed_diff_n_ptr[HVEC1_1],
307 &t_transformed_diff_n_ptr[HVEC1_2],
308 &t_transformed_diff_n_ptr[HVEC2_0],
309 &t_transformed_diff_n_ptr[HVEC2_1],
310 &t_transformed_diff_n_ptr[HVEC2_2]);
311
312 auto t_det = getFTensor0FromVec(*detPtr);
313 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
314
315 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
316 for (unsigned int bb = 0; bb != nb_diff_base_functions; ++bb) {
317 const double a = 1. / t_det;
318 t_transformed_diff_n(i, k) = a * t_jac(i, j) * t_diff_n(j, k);
319 ++t_diff_n;
320 ++t_transformed_diff_n;
321 }
322 ++t_det;
323 ++t_jac;
324 }
325
326 data.getDiffN(base).swap(piolaDiffN);
327 }
328 }
329
331}
332
337
341
342 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
343
344 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
345
346 unsigned int nb_gauss_pts = data.getN(base).size1();
347 unsigned int nb_base_functions = data.getN(base).size2() / 3;
348 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
349 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
350
351 auto t_n = data.getFTensor1N<3>(base);
352 double *t_transformed_n_ptr = &*piolaN.data().begin();
354 t_transformed_n_ptr, // HVEC0
355 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
356 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
357 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
358 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
359 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
360 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
361 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
362 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
363 &t_transformed_diff_n_ptr[HVEC2_2]);
364
365 auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
366
367 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
368 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
369 t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
370 t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
371 ++t_n;
372 ++t_transformed_n;
373 ++t_diff_n;
374 ++t_transformed_diff_n;
375 }
376 ++t_inv_jac;
377 }
378
379 data.getN(base).swap(piolaN);
380 data.getDiffN(base).swap(piolaDiffN);
381 }
382
384}
385
387 boost::shared_ptr<MatrixDouble> jac_ptr)
389 jacPtr(jac_ptr) {}
390
394
396
397 const auto nb_gauss_pts = getGaussPts().size2();
398
399 jacPtr->resize(4, nb_gauss_pts, false);
400 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
401
404
405 auto t_t1 = getFTensor1Tangent1AtGaussPts();
406 auto t_t2 = getFTensor1Tangent2AtGaussPts();
409
410 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
411 t_jac(i, N0) = t_t1(i);
412 t_jac(i, N1) = t_t2(i);
413 ++t_t1;
414 ++t_t2;
415 ++t_jac;
416 }
417
419};
420
425
429
430 size_t nb_gauss_pts = getGaussPts().size2();
431 jacPtr->resize(9, nb_gauss_pts, false);
432
433 auto t_t1 = getFTensor1Tangent1AtGaussPts();
434 auto t_t2 = getFTensor1Tangent2AtGaussPts();
435 auto t_normal = getFTensor1NormalsAtGaussPts();
436
440
441 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
442 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
443
444 t_jac(i, N0) = t_t1(i);
445 t_jac(i, N1) = t_t2(i);
446
447 const double l = sqrt(t_normal(j) * t_normal(j));
448
449 t_jac(i, N2) = t_normal(i) / l;
450
451 ++t_t1;
452 ++t_t2;
453 ++t_normal;
454 }
455
457}
458
460 int side, EntityType type, EntitiesFieldData::EntData &data) {
463
464 if (moab::CN::Dimension(type) != 2)
466
467 auto get_normals_ptr = [&]() {
469 return &*normalsAtGaussPts->data().begin();
470 else
471 return &*getNormalsAtGaussPts().data().begin();
472 };
473
474 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
475 auto nb_base_functions) {
477
478 auto ptr = get_normals_ptr();
480 &ptr[0], &ptr[1], &ptr[2]);
481
482 auto t_base = data.getFTensor1N<3>(base);
483 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
484 const auto l2 = t_normal(i) * t_normal(i);
485 for (int bb = 0; bb != nb_base_functions; ++bb) {
486 const auto v = t_base(0);
487 t_base(i) = (v / l2) * t_normal(i);
488 ++t_base;
489 }
490 ++t_normal;
491 }
493 };
494
495 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
496
497 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
498 const auto &base_functions = data.getN(base);
499 const auto nb_gauss_pts = base_functions.size1();
500
501 if (nb_gauss_pts) {
502
503 const auto nb_base_functions = base_functions.size2() / 3;
504 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
505 nb_base_functions);
506 }
507 }
508
510}
511
513 int side, EntityType type, EntitiesFieldData::EntData &data) {
515
516 const auto type_dim = moab::CN::Dimension(type);
517 if (type_dim != 1 && type_dim != 2)
519
523
524 auto get_jac = [&]() {
526 double *ptr_n = &*normalsAtPts->data().begin();
527 double *ptr_t1 = &*tangent1AtPts->data().begin();
528 double *ptr_t2 = &*tangent2AtPts->data().begin();
530 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
531
532 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
533
534 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
535 } else {
536 double *ptr_n = &*getNormalsAtGaussPts().data().begin();
537 double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
538 double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
540 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
541
542 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
543
544 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
545 }
546 };
547
548 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
549
550 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
551
552 auto &baseN = data.getN(base);
553 auto &diffBaseN = data.getDiffN(base);
554
555 int nb_dofs = baseN.size2() / 3;
556 int nb_gauss_pts = baseN.size1();
557
558 piolaN.resize(baseN.size1(), baseN.size2());
559 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
560
561 if (nb_dofs > 0 && nb_gauss_pts > 0) {
562
563 auto t_h_curl = data.getFTensor1N<3>(base);
564 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
565
566 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
567 &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
568
570 t_transformed_diff_h_curl(
574
575 int cc = 0;
576 double det;
578
579 // HO geometry is set, so jacobian is different at each gauss point
580 auto t_m_at_pts = get_jac();
581 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
582 CHKERR determinantTensor3by3(t_m_at_pts, det);
583 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
584 for (int ll = 0; ll != nb_dofs; ll++) {
585 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
586 t_transformed_diff_h_curl(i, k) = t_inv_m(j, i) * t_diff_h_curl(j, k);
587 ++t_h_curl;
588 ++t_transformed_h_curl;
589 ++t_diff_h_curl;
590 ++t_transformed_diff_h_curl;
591 ++cc;
592 }
593 ++t_m_at_pts;
594 }
595
596#ifndef NDEBUG
597 if (cc != nb_gauss_pts * nb_dofs)
598 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
599#endif
600
601 baseN.swap(piolaN);
602 diffBaseN.swap(diffPiolaN);
603 }
604 }
605
607}
608
610 int side, EntityType type, EntitiesFieldData::EntData &data) {
613
614 if (type != MBEDGE)
616
617 const auto nb_gauss_pts = getGaussPts().size2();
618
620 if (tangentAtGaussPts->size1() != nb_gauss_pts)
621 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
622 "Wrong number of integration points %d!=%d",
623 tangentAtGaussPts->size1(), nb_gauss_pts);
624
625 auto get_base_at_pts = [&](auto base) {
627 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
628 &data.getN(base)(0, HVEC2));
629 return t_h_curl;
630 };
631
632 auto get_tangent_ptr = [&]() {
633 if (tangentAtGaussPts) {
634 return &*tangentAtGaussPts->data().begin();
635 } else {
636 return &*getTangentAtGaussPts().data().begin();
637 }
638 };
639
640 auto get_tangent_at_pts = [&]() {
641 auto ptr = get_tangent_ptr();
643 &ptr[0], &ptr[1], &ptr[2]);
644 return t_m_at_pts;
645 };
646
647 auto calculate_squared_edge_length = [&]() {
648 std::vector<double> l1;
649 if (nb_gauss_pts) {
650 l1.resize(nb_gauss_pts);
651 auto t_m_at_pts = get_tangent_at_pts();
652 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
653 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
654 ++t_m_at_pts;
655 }
656 }
657 return l1;
658 };
659
660 auto l1 = calculate_squared_edge_length();
661
662 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
663
664 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
665 const auto nb_dofs = data.getN(base).size2() / 3;
666
667 if (nb_gauss_pts && nb_dofs) {
668 auto t_h_curl = get_base_at_pts(base);
669 int cc = 0;
670 auto t_m_at_pts = get_tangent_at_pts();
671 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
672 const double l0 = l1[gg];
673 for (int ll = 0; ll != nb_dofs; ++ll) {
674 const double val = t_h_curl(0);
675 const double a = val / l0;
676 t_h_curl(i) = t_m_at_pts(i) * a;
677 ++t_h_curl;
678 ++cc;
679 }
680 ++t_m_at_pts;
681 }
682
683#ifndef NDEBUG
684 if (cc != nb_gauss_pts * nb_dofs)
685 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
686#endif
687 }
688 }
689
691}
692
694 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
695 std::vector<FieldSpace> spaces, std::string geom_field_name,
696 boost::shared_ptr<MatrixDouble> jac_ptr,
697 boost::shared_ptr<VectorDouble> det_ptr,
698 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
700
701 if (!jac_ptr)
702 jac_ptr = boost::make_shared<MatrixDouble>();
703 if (!det_ptr)
704 det_ptr = boost::make_shared<VectorDouble>();
705 if (!inv_jac_ptr)
706 inv_jac_ptr = boost::make_shared<MatrixDouble>();
707
708 if (geom_field_name.empty()) {
709
710 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
711
712 } else {
713
714 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
715 pipeline.push_back(
716 new OpCalculateVectorFieldGradient<2, 2>(geom_field_name, jac_ptr));
717 pipeline.push_back(new OpGetHONormalsOnFace<2>(geom_field_name));
718 }
719
720 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
721 pipeline.push_back(new OpSetHOWeightsOnFace());
722
723 for (auto s : spaces) {
724 switch (s) {
725 case NOSPACE:
726 break;
727 case H1:
728 case L2:
729 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(s, inv_jac_ptr));
730 break;
731 case HCURL:
732 pipeline.push_back(new OpSetCovariantPiolaTransformOnFace2D(inv_jac_ptr));
733 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
734 break;
735 case HDIV:
736 pipeline.push_back(new OpMakeHdivFromHcurl());
737 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
738 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
739 break;
740 default:
741 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
742 "Space %s not yet implemented", FieldSpaceNames[s]);
743 }
744 }
745
747}
748
750 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
751 std::vector<FieldSpace> spaces, std::string geom_field_name) {
753
754 if (geom_field_name.empty()) {
755
756 } else {
757
758 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
759 pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
760 }
761
762 for (auto s : spaces) {
763 switch (s) {
764 case NOSPACE:
765 break;
766 case HCURL:
767 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
768 break;
769 case HDIV:
770 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
771 break;
772 default:
773 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
774 "Space %s not yet implemented", FieldSpaceNames[s]);
775 }
776 }
777
779}
780
782 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
783 std::vector<FieldSpace> spaces, std::string geom_field_name) {
785
786 if (geom_field_name.empty()) {
787
788 } else {
789
790 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
791 pipeline.push_back(new OpGetHOTangentsOnEdge<3>(geom_field_name));
792 }
793
794 for (auto s : spaces) {
795 switch (s) {
796 case HCURL:
797 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
798 break;
799 default:
800 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
801 "Space %s not yet implemented", FieldSpaceNames[s]);
802 }
803 }
804
806}
807
809 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
810 std::vector<FieldSpace> spaces, std::string geom_field_name) {
812
813 if (geom_field_name.empty()) {
814 } else {
815
816 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
817 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
818 }
819
820 for (auto s : spaces) {
821 switch (s) {
822 case NOSPACE:
823 break;
824 case HCURL:
825 pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
826 break;
827 case HDIV:
828 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
829 break;
830 default:
831 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
832 "Space %s not yet implemented", FieldSpaceNames[s]);
833 break;
834 }
835 }
836
838}
839
841 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
842 std::vector<FieldSpace> spaces, std::string geom_field_name,
843 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
844 boost::shared_ptr<MatrixDouble> inv_jac) {
846
847 if (!jac)
848 jac = boost::make_shared<MatrixDouble>();
849 if (!det)
850 det = boost::make_shared<VectorDouble>();
851 if (!inv_jac)
852 inv_jac = boost::make_shared<MatrixDouble>();
853
854 if (geom_field_name.empty()) {
855
856 pipeline.push_back(new OpCalculateHOJac<3>(jac));
857
858 } else {
859
860 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
861 pipeline.push_back(
862 new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
863 }
864
865 pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
866 pipeline.push_back(new OpSetHOWeights(det));
867
868 for (auto s : spaces) {
869 switch (s) {
870 case NOSPACE:
871 break;
872 case H1:
873 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
874 break;
875 case HCURL:
876 pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
877 pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
878 break;
879 case HDIV:
880 pipeline.push_back(
882 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac));
883 break;
884 case L2:
885 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
886 break;
887 default:
888 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
889 "Space %s not yet implemented", FieldSpaceNames[s]);
890 break;
891 }
892 }
893
895}
896
897} // 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:1608
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:1557
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
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, DIM > getFTensor1TangentAtGaussPts()
Get tangent vector 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.
virtual boost::shared_ptr< MatrixDouble > & getDiffNSharedPtr(const FieldApproximationBase base)
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
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.
double getMeasure() const
get measure of element
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
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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