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
53 FTensor::Index<'i', 3> i;
54 FTensor::Index<'j', 3> j;
55 FTensor::Index<'k', 3> k;
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
113 FTensor::Index<'i', 3> i;
114 FTensor::Index<'j', 3> j;
115 FTensor::Index<'k', 3> k;
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();
164 FTensor::Index<'i', 3> i;
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();
189 FTensor::Index<'i', 3> i;
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
233 FTensor::Index<'i', 3> i;
234 FTensor::Index<'j', 3> j;
235 FTensor::Index<'k', 3> k;
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
338 FTensor::Index<'i', 3> i;
339 FTensor::Index<'j', 3> j;
340 FTensor::Index<'k', 3> k;
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
402 FTensor::Index<'i', 2> i;
403 FTensor::Index<'j', 2> j;
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
426 FTensor::Index<'i', 3> i;
427 FTensor::Index<'j', 3> j;
428 FTensor::Index<'k', 3> k;
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) {
461 FTensor::Index<'i', 3> i;
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
520 FTensor::Index<'i', 3> i;
521 FTensor::Index<'j', 3> j;
522 FTensor::Index<'k', 2> k;
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) {
611 FTensor::Index<'i', 3> i;
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,
752 boost::shared_ptr<MatrixDouble> jac_ptr,
753 boost::shared_ptr<VectorDouble> det_ptr,
754 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
756
757 if (!jac_ptr)
758 jac_ptr = boost::make_shared<MatrixDouble>();
759 if (!det_ptr)
760 det_ptr = boost::make_shared<VectorDouble>();
761 if (!inv_jac_ptr)
762 inv_jac_ptr = boost::make_shared<MatrixDouble>();
763
764 if (geom_field_name.empty()) {
765
766 } else {
767
768 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
769 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
770 }
771
772 pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
773 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
774 pipeline.push_back(new OpSetHOWeightsOnFace());
775
776 for (auto s : spaces) {
777 switch (s) {
778 case NOSPACE:
779 break;
780 case H1:
781 pipeline.push_back(
783 break;
784 case HDIV:
785 pipeline.push_back(new OpMakeHdivFromHcurl());
786 pipeline.push_back(
788 jac_ptr));
789 pipeline.push_back(
791 break;
792 case L2:
793 pipeline.push_back(
795 break;
796 default:
797 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
798 "Space %s not yet implemented", FieldSpaceNames[s]);
799 }
800 }
801
803}
804
806 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
807 std::vector<FieldSpace> spaces, std::string geom_field_name) {
809
810 if (geom_field_name.empty()) {
811
812 } else {
813
814 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
815 pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
816 }
817
818 for (auto s : spaces) {
819 switch (s) {
820 case NOSPACE:
821 break;
822 case HCURL:
823 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
824 break;
825 case HDIV:
826 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
827 break;
828 default:
829 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
830 "Space %s not yet implemented", FieldSpaceNames[s]);
831 }
832 }
833
835}
836
838 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
839 std::vector<FieldSpace> spaces, std::string geom_field_name) {
841
842 if (geom_field_name.empty()) {
843
844 } else {
845
846 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
847 pipeline.push_back(new OpGetHOTangentsOnEdge<3>(geom_field_name));
848 }
849
850 for (auto s : spaces) {
851 switch (s) {
852 case HCURL:
853 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
854 break;
855 default:
856 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
857 "Space %s not yet implemented", FieldSpaceNames[s]);
858 }
859 }
860
862}
863
865 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
866 std::vector<FieldSpace> spaces, std::string geom_field_name) {
868
869 if (geom_field_name.empty()) {
870 } else {
871
872 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
873 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
874 }
875
876 for (auto s : spaces) {
877 switch (s) {
878 case NOSPACE:
879 break;
880 case HCURL:
881 pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
882 break;
883 case HDIV:
884 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
885 break;
886 default:
887 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
888 "Space %s not yet implemented", FieldSpaceNames[s]);
889 break;
890 }
891 }
892
894}
895
897 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
898 std::vector<FieldSpace> spaces, std::string geom_field_name,
899 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
900 boost::shared_ptr<MatrixDouble> inv_jac) {
902
903 if (!jac)
904 jac = boost::make_shared<MatrixDouble>();
905 if (!det)
906 det = boost::make_shared<VectorDouble>();
907 if (!inv_jac)
908 inv_jac = boost::make_shared<MatrixDouble>();
909
910 if (geom_field_name.empty()) {
911
912 pipeline.push_back(new OpCalculateHOJac<3>(jac));
913
914 } else {
915
916 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
917 pipeline.push_back(
918 new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
919 }
920
921 pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
922 pipeline.push_back(new OpSetHOWeights(det));
923
924 for (auto s : spaces) {
925 switch (s) {
926 case NOSPACE:
927 break;
928 case H1:
929 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
930 break;
931 case HCURL:
932 pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
933 pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
934 break;
935 case HDIV:
936 pipeline.push_back(
938 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac));
939 break;
940 case L2:
941 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
942 break;
943 default:
944 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
945 "Space %s not yet implemented", FieldSpaceNames[s]);
946 break;
947 }
948 }
949
951}
952
953} // 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
T data[Tensor_Dim]
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()
@ 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 ...
@ HVEC0
@ HVEC1
@ HVEC2
@ 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()
@ HVEC1_1
@ HVEC0_1
@ HVEC1_0
@ HVEC2_1
@ HVEC1_2
@ HVEC2_2
@ HVEC2_0
@ HVEC0_2
@ HVEC0_0
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
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.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
OpSetInvJacHcurlFaceImpl< 3 > OpSetInvJacHcurlFaceEmbeddedIn3DSpace
OpSetCovariantPiolaTransformOnFace2DImpl< 2 > OpSetCovariantPiolaTransformOnFace2D
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.
OpSetContravariantPiolaTransformOnFace2DImpl< 3 > OpSetContravariantPiolaTransformOnFace2DEmbeddedIn3DSpace
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
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