v0.15.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
22OpCalculateHOJacForVolume::doWork(int side, EntityType type,
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 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
42 "Number of vertex coordinates and base functions is inconsistent "
43 "%ld != %ld",
44 coords.size() / 3, nb_base_functions);
45#endif
46
47 jacPtr->resize(9, nb_gauss_pts, false);
48 jacPtr->clear();
49 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
50 auto t_vol_inv_jac = getInvJac();
51
52 FTensor::Index<'i', 3> i;
53 FTensor::Index<'j', 3> j;
54 FTensor::Index<'k', 3> k;
55
56 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
58 &coords[0], &coords[1], &coords[2]);
59 for (size_t bb = 0; bb != nb_base_functions; ++bb) {
60 t_jac(i, j) += t_coords(i) * (t_vol_inv_jac(k, j) * t_diff_base(k));
61 ++t_diff_base;
62 ++t_coords;
63 }
64 ++t_jac;
65 }
66 }
67
69}
70
75
76 if (getFEDim() == 3) {
77
78 auto transform_base = [&](MatrixDouble &diff_n) {
79 return applyTransform<3, 3, 3, 3>(diff_n);
80 };
81
82 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
84 CHKERR transform_base(data.getDiffN(base));
85 }
86
87 switch (type) {
88 case MBVERTEX:
89 for (auto &m : data.getBBDiffNMap())
90 CHKERR transform_base(*(m.second));
91 break;
92 default:
93 for (auto &ptr : data.getBBDiffNByOrderArray())
94 if (ptr)
95 CHKERR transform_base(*ptr);
96 }
97
98 } else {
99 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Use different operator");
100 }
101
103}
104
106OpSetHOInvJacVectorBase::doWork(int side, EntityType type,
109
110 FTensor::Index<'i', 3> i;
111 FTensor::Index<'j', 3> j;
112 FTensor::Index<'k', 3> k;
113
114 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
115
116 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
117
118 diffHdivInvJac.resize(data.getDiffN(base).size1(),
119 data.getDiffN(base).size2(), false);
120
121 unsigned int nb_gauss_pts = data.getDiffN(base).size1();
122 unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
123 if (nb_base_functions == 0)
124 continue;
125
126 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
127 double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
129 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
130 &t_inv_diff_n_ptr[HVEC0_2], &t_inv_diff_n_ptr[HVEC1_0],
131 &t_inv_diff_n_ptr[HVEC1_1], &t_inv_diff_n_ptr[HVEC1_2],
132 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
133 &t_inv_diff_n_ptr[HVEC2_2]);
134 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
135 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
136 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
137 t_inv_diff_n(i, j) = t_inv_jac(k, j) * t_diff_n(i, k);
138 ++t_diff_n;
139 ++t_inv_diff_n;
140 }
141 ++t_inv_jac;
142 }
143
144 data.getDiffN(base).swap(diffHdivInvJac);
145 }
146
148}
149
153
154 /* Note that is different treatment to OpSetHOWeight, that is because volume
155 * element by default calculate HO jacobian on physical element, whereas face
156 * element uses reference element. */
157
158 const size_t nb_int_pts = getGaussPts().size2();
159 if (getNormalsAtGaussPts().size1()) {
160 if (getNormalsAtGaussPts().size1() == nb_int_pts) {
161 double a = getMeasure();
162 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI)
163 a *= 2;
164 auto t_w = getFTensor0IntegrationWeight();
165 auto t_normal = getFTensor1NormalsAtGaussPts();
166 FTensor::Index<'i', 3> i;
167 for (size_t gg = 0; gg != nb_int_pts; ++gg) {
168 t_w *= sqrt(t_normal(i) * t_normal(i)) / a;
169 ++t_w;
170 ++t_normal;
171 }
172 } else {
173 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
174 "Number of rows in getNormalsAtGaussPts should be equal to "
175 "number of integration points, but is not, i.e. %ld != %ld",
176 getNormalsAtGaussPts().size1(), nb_int_pts);
177 }
178 }
180}
181
185 const size_t nb_int_pts = getGaussPts().size2();
186 if (getTangentAtGaussPts().size1()) {
187 if (getTangentAtGaussPts().size1() == nb_int_pts) {
188 double a = getMeasure();
189 auto t_w = getFTensor0IntegrationWeight();
190 auto t_tangent = getFTensor1TangentAtGaussPts();
191 FTensor::Index<'i', 3> i;
192 for (size_t gg = 0; gg != nb_int_pts; ++gg) {
193 t_w *= sqrt(t_tangent(i) * t_tangent(i)) / a;
194 ++t_w;
195 ++t_tangent;
196 }
197 } else {
198 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
199 "Number of rows in getTangentAtGaussPts should be equal to "
200 "number of integration points, but is not, i.e. %ld != %ld",
201 getTangentAtGaussPts().size1(), nb_int_pts);
202 }
203 }
205}
206
207MoFEMErrorCode OpSetHOWeights::doWork(int side, EntityType type,
210
211 const auto nb_integration_pts = detPtr->size();
212
213#ifndef NDEBUG
214 if (nb_integration_pts != getGaussPts().size2())
215 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
216 "Inconsistent number of data points");
217#endif
218
219 auto t_w = getFTensor0IntegrationWeight();
220 auto t_det = getFTensor0FromVec(*detPtr);
221 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
222 t_w *= t_det;
223 ++t_w;
224 ++t_det;
225 }
226
228}
229
234
235 FTensor::Index<'i', 3> i;
236 FTensor::Index<'j', 3> j;
237 FTensor::Index<'k', 3> k;
238
239#ifndef NDEBUG
240 if (!detPtr)
241 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
242 "Pointer for detPtr not allocated");
243 if (!jacPtr)
244 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
245 "Pointer for jacPtr not allocated");
246#endif
247
248 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
249
250 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
251
252 auto nb_gauss_pts = data.getNSharedPtr(base) ? data.getN(base).size1() : 0;
253 auto nb_base_functions =
254 data.getNSharedPtr(base) ? data.getN(base).size2() / 3 : 0;
255 auto nb_diff_base_functions =
256 data.getDiffNSharedPtr(base) ? data.getDiffN(base).size2() / 9 : 0;
257
258#ifndef NDEBUG
259 if (nb_diff_base_functions) {
260 if (nb_diff_base_functions != nb_base_functions)
261 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
262 "Wrong number base functions %ld != %ld",
263 nb_diff_base_functions, nb_base_functions);
264 if (data.getDiffN(base).size1() != nb_gauss_pts)
265 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
266 "Wrong number integration points");
267 }
268#endif
269
270 if (nb_gauss_pts && nb_base_functions) {
271
272 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
273
274 auto t_n = data.getFTensor1N<3>(base);
275 double *t_transformed_n_ptr = &*piolaN.data().begin();
277 t_transformed_n_ptr, // HVEC0
278 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
279
280 auto t_det = getFTensor0FromVec(*detPtr);
281 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
282
283 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
284 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
285 const double a = 1. / t_det;
286 t_transformed_n(i) = a * (t_jac(i, k) * t_n(k));
287 ++t_n;
288 ++t_transformed_n;
289 }
290 ++t_det;
291 ++t_jac;
292 }
293
294 data.getN(base).swap(piolaN);
295 }
296
297 if (nb_gauss_pts && nb_diff_base_functions) {
298
299 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
300
301 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
302 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
304 t_transformed_diff_n(t_transformed_diff_n_ptr,
305 &t_transformed_diff_n_ptr[HVEC0_1],
306 &t_transformed_diff_n_ptr[HVEC0_2],
307 &t_transformed_diff_n_ptr[HVEC1_0],
308 &t_transformed_diff_n_ptr[HVEC1_1],
309 &t_transformed_diff_n_ptr[HVEC1_2],
310 &t_transformed_diff_n_ptr[HVEC2_0],
311 &t_transformed_diff_n_ptr[HVEC2_1],
312 &t_transformed_diff_n_ptr[HVEC2_2]);
313
314 auto t_det = getFTensor0FromVec(*detPtr);
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_diff_n(i, k);
319 ++t_diff_n;
320 ++t_transformed_diff_n;
321 }
322 ++t_det;
323 }
324
325 data.getDiffN(base).swap(piolaDiffN);
326 }
327 }
328
330}
331
336
337 FTensor::Index<'i', 3> i;
338 FTensor::Index<'j', 3> j;
339 FTensor::Index<'k', 3> k;
340
341 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
342
343 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
344
345 unsigned int nb_gauss_pts = data.getN(base).size1();
346 unsigned int nb_base_functions = data.getN(base).size2() / 3;
347 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
348 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
349
350 auto t_n = data.getFTensor1N<3>(base);
351 double *t_transformed_n_ptr = &*piolaN.data().begin();
353 t_transformed_n_ptr, // HVEC0
354 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
355 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
356 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
357 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
358 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
359 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
360 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
361 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
362 &t_transformed_diff_n_ptr[HVEC2_2]);
363
364 auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
365
366 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
367 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
368 t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
369 t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
370 ++t_n;
371 ++t_transformed_n;
372 ++t_diff_n;
373 ++t_transformed_diff_n;
374 }
375 ++t_inv_jac;
376 }
377
378 data.getN(base).swap(piolaN);
379 data.getDiffN(base).swap(piolaDiffN);
380 }
381
383}
384
386 boost::shared_ptr<MatrixDouble> jac_ptr)
388 jacPtr(jac_ptr) {}
389
393
395
396 const auto nb_gauss_pts = getGaussPts().size2();
397
398 jacPtr->resize(4, nb_gauss_pts, false);
399 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
400
401 FTensor::Index<'i', 2> i;
402 FTensor::Index<'j', 2> j;
403
404 auto t_t1 = getFTensor1Tangent1AtGaussPts();
405 auto t_t2 = getFTensor1Tangent2AtGaussPts();
408
409 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
410 t_jac(i, N0) = t_t1(i);
411 t_jac(i, N1) = t_t2(i);
412 ++t_t1;
413 ++t_t2;
414 ++t_jac;
415 }
416
418};
419
424
425 FTensor::Index<'i', 3> i;
426 FTensor::Index<'j', 3> j;
427 FTensor::Index<'k', 3> k;
428
429 size_t nb_gauss_pts = getGaussPts().size2();
430 jacPtr->resize(9, nb_gauss_pts, false);
431
432 auto t_t1 = getFTensor1Tangent1AtGaussPts();
433 auto t_t2 = getFTensor1Tangent2AtGaussPts();
434 auto t_normal = getFTensor1NormalsAtGaussPts();
435
439
440 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
441 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
442
443 t_jac(i, N0) = t_t1(i);
444 t_jac(i, N1) = t_t2(i);
445
446 const double l = sqrt(t_normal(j) * t_normal(j));
447
448 t_jac(i, N2) = t_normal(i) / l;
449
450 ++t_t1;
451 ++t_t2;
452 ++t_normal;
453 }
454
456}
457
459 int side, EntityType type, EntitiesFieldData::EntData &data) {
460 FTensor::Index<'i', 3> i;
462
463 if (moab::CN::Dimension(type) != 2)
465
466 auto get_normals_ptr = [&]() {
468 return &*normalsAtGaussPts->data().begin();
469 else
470 return &*getNormalsAtGaussPts().data().begin();
471 };
472
473 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
474 auto nb_base_functions) {
476
477 auto ptr = get_normals_ptr();
479 &ptr[0], &ptr[1], &ptr[2]);
480
481 auto t_base = data.getFTensor1N<3>(base);
482 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
483 const auto l2 = t_normal(i) * t_normal(i);
484 for (int bb = 0; bb != nb_base_functions; ++bb) {
485 const auto v = t_base(0);
486 t_base(i) = (v / l2) * t_normal(i);
487 ++t_base;
488 }
489 ++t_normal;
490 }
492 };
493
494 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
495
496 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
497 const auto &base_functions = data.getN(base);
498 const auto nb_gauss_pts = base_functions.size1();
499
500 if (nb_gauss_pts) {
501
502 const auto nb_base_functions = base_functions.size2() / 3;
503 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
504 nb_base_functions);
505 }
506 }
507
509}
510
512 int side, EntityType type, EntitiesFieldData::EntData &data) {
514
515 const auto type_dim = moab::CN::Dimension(type);
516 if (type_dim != 1 && type_dim != 2)
518
519 FTensor::Index<'i', 3> i;
520 FTensor::Index<'j', 3> j;
521 FTensor::Index<'k', 2> k;
522
523 auto get_jac = [&]() {
525 double *ptr_n = &*normalsAtPts->data().begin();
526 double *ptr_t1 = &*tangent1AtPts->data().begin();
527 double *ptr_t2 = &*tangent2AtPts->data().begin();
529 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
530
531 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
532
533 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
534 } else {
535 double *ptr_n = &*getNormalsAtGaussPts().data().begin();
536 double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
537 double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
539 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
540
541 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
542
543 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
544 }
545 };
546
547 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
548
549 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
550
551 auto &baseN = data.getN(base);
552 auto &diffBaseN = data.getDiffN(base);
553
554 int nb_dofs = baseN.size2() / 3;
555 int nb_gauss_pts = baseN.size1();
556
557 piolaN.resize(baseN.size1(), baseN.size2());
558 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
559
560 if (nb_dofs > 0 && nb_gauss_pts > 0) {
561
562 auto t_h_curl = data.getFTensor1N<3>(base);
563 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
564
565 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
566 &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
567
569 t_transformed_diff_h_curl(
573
574 int cc = 0;
575 double det;
577
578 // HO geometry is set, so jacobian is different at each gauss point
579 auto t_m_at_pts = get_jac();
580 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
581 CHKERR determinantTensor3by3(t_m_at_pts, det);
582 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
583 for (int ll = 0; ll != nb_dofs; ll++) {
584 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
585 t_transformed_diff_h_curl(i, k) = t_inv_m(j, i) * t_diff_h_curl(j, k);
586 ++t_h_curl;
587 ++t_transformed_h_curl;
588 ++t_diff_h_curl;
589 ++t_transformed_diff_h_curl;
590 ++cc;
591 }
592 ++t_m_at_pts;
593 }
594
595#ifndef NDEBUG
596 if (cc != nb_gauss_pts * nb_dofs)
597 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
598#endif
599
600 baseN.swap(piolaN);
601 diffBaseN.swap(diffPiolaN);
602 }
603 }
604
606}
607
609 int side, EntityType type, EntitiesFieldData::EntData &data) {
610 FTensor::Index<'i', 3> i;
612
613 if (type != MBEDGE)
615
616 const auto nb_gauss_pts = getGaussPts().size2();
617
619 if (tangentAtGaussPts->size1() != nb_gauss_pts)
620 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
621 "Wrong number of integration points %ld!=%ld",
622 tangentAtGaussPts->size1(), nb_gauss_pts);
623
624 auto get_base_at_pts = [&](auto base) {
626 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
627 &data.getN(base)(0, HVEC2));
628 return t_h_curl;
629 };
630
631 auto get_tangent_ptr = [&]() {
632 if (tangentAtGaussPts) {
633 return &*tangentAtGaussPts->data().begin();
634 } else {
635 return &*getTangentAtGaussPts().data().begin();
636 }
637 };
638
639 auto get_tangent_at_pts = [&]() {
640 auto ptr = get_tangent_ptr();
642 &ptr[0], &ptr[1], &ptr[2]);
643 return t_m_at_pts;
644 };
645
646 auto calculate_squared_edge_length = [&]() {
647 std::vector<double> l1;
648 if (nb_gauss_pts) {
649 l1.resize(nb_gauss_pts);
650 auto t_m_at_pts = get_tangent_at_pts();
651 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
652 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
653 ++t_m_at_pts;
654 }
655 }
656 return l1;
657 };
658
659 auto l1 = calculate_squared_edge_length();
660
661 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
662
663 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
664 const auto nb_dofs = data.getN(base).size2() / 3;
665
666 if (nb_gauss_pts && nb_dofs) {
667 auto t_h_curl = get_base_at_pts(base);
668 int cc = 0;
669 auto t_m_at_pts = get_tangent_at_pts();
670 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
671 const double l0 = l1[gg];
672 for (int ll = 0; ll != nb_dofs; ++ll) {
673 const double val = t_h_curl(0);
674 const double a = val / l0;
675 t_h_curl(i) = t_m_at_pts(i) * a;
676 ++t_h_curl;
677 ++cc;
678 }
679 ++t_m_at_pts;
680 }
681
682#ifndef NDEBUG
683 if (cc != nb_gauss_pts * nb_dofs)
684 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Data inconsistency");
685#endif
686 }
687 }
688
690}
691
693 const FieldSpace space, const FieldApproximationBase base,
694 boost::shared_ptr<VectorDouble> det_jac_ptr,
695 boost::shared_ptr<double> scale_det_ptr)
696 : OP(space), fieldSpace(space), fieldBase(base), detJacPtr(det_jac_ptr),
697 scaleDetPtr(scale_det_ptr) {
698 if (!detJacPtr) {
699 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
700 }
701}
702
707
708 double scale_det = 1.0;
709 if (scaleDetPtr)
710 scale_det = *scaleDetPtr;
711
712 auto scale = [&](auto base) {
714 auto &base_fun = data.getN(base);
715 if (base_fun.size2()) {
716
717 auto &det_vec = *detJacPtr;
718 const auto nb_int_pts = base_fun.size1();
719
720 if (nb_int_pts != det_vec.size())
721 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
722 "Number of integration pts in detJacPtr does not mush "
723 "number of integration pts in base function");
724
725 for (auto gg = 0; gg != nb_int_pts; ++gg) {
726 auto row = ublas::matrix_row<MatrixDouble>(base_fun, gg);
727 row *= scale_det / det_vec[gg];
728 }
729
730 auto &diff_base_fun = data.getDiffN(base);
731 if (diff_base_fun.size2()) {
732 for (auto gg = 0; gg != nb_int_pts; ++gg) {
733 auto row = ublas::matrix_row<MatrixDouble>(diff_base_fun, gg);
734 row *= scale_det / det_vec[gg];
735 }
736 }
737 }
739 };
740
741 if (this->getFEDim() == 3) {
742 switch (fieldSpace) {
743 case L2:
744 if (type >= MBTET)
746 break;
747 default:
748 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
749 }
750 } else if (this->getFEDim() == 2) {
751 switch (fieldSpace) {
752 case L2:
753 if (type >= MBTRI)
755 break;
756 default:
757 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
758 }
759 } else if (this->getFEDim() == 1) {
760 switch (fieldSpace) {
761 case L2:
762 if (type >= MBEDGE)
764 break;
765 default:
766 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not handled");
767 }
768 } else {
769 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "dimension not handled");
770 }
771
773}
774
776 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
777 std::vector<FieldSpace> spaces, std::string geom_field_name,
778 boost::shared_ptr<MatrixDouble> jac_ptr,
779 boost::shared_ptr<VectorDouble> det_ptr,
780 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
782
783 if (!jac_ptr)
784 jac_ptr = boost::make_shared<MatrixDouble>();
785 if (!det_ptr)
786 det_ptr = boost::make_shared<VectorDouble>();
787 if (!inv_jac_ptr)
788 inv_jac_ptr = boost::make_shared<MatrixDouble>();
789
790 if (geom_field_name.empty()) {
791
792 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
793
794 } else {
795
796 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
797 pipeline.push_back(
798 new OpCalculateVectorFieldGradient<2, 2>(geom_field_name, jac_ptr));
799 pipeline.push_back(new OpGetHONormalsOnFace<2>(geom_field_name));
800 }
801
802 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
803 pipeline.push_back(new OpSetHOWeightsOnFace());
804
805 for (auto s : spaces) {
806 switch (s) {
807 case NOSPACE:
808 break;
809 case H1:
810 case L2:
811 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(s, inv_jac_ptr));
812 break;
813 case HCURL:
814 pipeline.push_back(new OpSetCovariantPiolaTransformOnFace2D(inv_jac_ptr));
815 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
816 break;
817 case HDIV:
818 pipeline.push_back(new OpMakeHdivFromHcurl());
819 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
820 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
821 break;
822 default:
823 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
824 "Space %s not yet implemented", FieldSpaceNames[s]);
825 }
826 }
827
829}
830
832 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
833 std::vector<FieldSpace> spaces, std::string geom_field_name,
834 boost::shared_ptr<MatrixDouble> jac_ptr,
835 boost::shared_ptr<VectorDouble> det_ptr,
836 boost::shared_ptr<MatrixDouble> inv_jac_ptr) {
838
839 if (!jac_ptr)
840 jac_ptr = boost::make_shared<MatrixDouble>();
841 if (!det_ptr)
842 det_ptr = boost::make_shared<VectorDouble>();
843 if (!inv_jac_ptr)
844 inv_jac_ptr = boost::make_shared<MatrixDouble>();
845
846 if (geom_field_name.empty()) {
847
848 } else {
849
850 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
851 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
852 }
853
854 pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
855 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
856 pipeline.push_back(new OpSetHOWeightsOnFace());
857
858 for (auto s : spaces) {
859 switch (s) {
860 case NOSPACE:
861 break;
862 case H1:
863 pipeline.push_back(
865 break;
866 case HDIV:
867 pipeline.push_back(new OpMakeHdivFromHcurl());
868 pipeline.push_back(
870 jac_ptr));
871 pipeline.push_back(
873 break;
874 case L2:
875 pipeline.push_back(
877 break;
878 default:
879 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
880 "Space %s not yet implemented", FieldSpaceNames[s]);
881 }
882 }
883
885}
886
888 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
889 std::vector<FieldSpace> spaces, std::string geom_field_name) {
891
892 if (geom_field_name.empty()) {
893
894 } else {
895
896 pipeline.push_back(new OpCalculateHOCoords<2>(geom_field_name));
897 pipeline.push_back(new OpGetHOTangentsOnEdge<2>(geom_field_name));
898 }
899
900 for (auto s : spaces) {
901 switch (s) {
902 case NOSPACE:
903 break;
904 case HCURL:
905 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
906 break;
907 case HDIV:
908 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
909 break;
910 default:
911 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
912 "Space %s not yet implemented", FieldSpaceNames[s]);
913 }
914 }
915
917}
918
920 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
921 std::vector<FieldSpace> spaces, std::string geom_field_name) {
923
924 if (geom_field_name.empty()) {
925
926 } else {
927
928 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
929 pipeline.push_back(new OpGetHOTangentsOnEdge<3>(geom_field_name));
930 }
931
932 for (auto s : spaces) {
933 switch (s) {
934 case HCURL:
935 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
936 break;
937 default:
938 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
939 "Space %s not yet implemented", FieldSpaceNames[s]);
940 }
941 }
942
944}
945
947 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
948 std::vector<FieldSpace> spaces, std::string geom_field_name) {
950
951 if (geom_field_name.empty()) {
952 } else {
953
954 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
955 pipeline.push_back(new OpGetHONormalsOnFace<3>(geom_field_name));
956 }
957
958 for (auto s : spaces) {
959 switch (s) {
960 case NOSPACE:
961 break;
962 case HCURL:
963 pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
964 break;
965 case HDIV:
966 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
967 break;
968 case L2:
969 break;
970 default:
971 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
972 "Space %s not yet implemented", FieldSpaceNames[s]);
973 break;
974 }
975 }
976
978}
979
981 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
982 std::vector<FieldSpace> spaces, std::string geom_field_name,
983 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
984 boost::shared_ptr<MatrixDouble> inv_jac) {
986
987 if (!jac)
988 jac = boost::make_shared<MatrixDouble>();
989 if (!det)
990 det = boost::make_shared<VectorDouble>();
991 if (!inv_jac)
992 inv_jac = boost::make_shared<MatrixDouble>();
993
994 if (geom_field_name.empty()) {
995
996 pipeline.push_back(new OpCalculateHOJac<3>(jac));
997
998 } else {
999
1000 pipeline.push_back(new OpCalculateHOCoords<3>(geom_field_name));
1001 pipeline.push_back(
1002 new OpCalculateVectorFieldGradient<3, 3>(geom_field_name, jac));
1003 }
1004
1005 pipeline.push_back(new OpInvertMatrix<3>(jac, det, inv_jac));
1006 pipeline.push_back(new OpSetHOWeights(det));
1007
1008 for (auto s : spaces) {
1009 switch (s) {
1010 case NOSPACE:
1011 break;
1012 case H1:
1013 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac));
1014 break;
1015 case HCURL:
1016 pipeline.push_back(new OpSetHOCovariantPiolaTransform(HCURL, inv_jac));
1017 pipeline.push_back(new OpSetHOInvJacVectorBase(HCURL, inv_jac));
1018 break;
1019 case HDIV:
1020 pipeline.push_back(
1022 break;
1023 case L2:
1024 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(L2, inv_jac));
1025 break;
1026 default:
1027 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1028 "Space %s not yet implemented", FieldSpaceNames[s]);
1029 break;
1030 }
1031 }
1032
1034}
1035
1036} // namespace MoFEM
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 nme:nme847.
Definition definitions.h:60
@ NOBASE
Definition definitions.h:59
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
Definition definitions.h:82
@ 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< '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.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
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
double scale
Definition plastic.cpp:123
constexpr double t
plate stiffness
Definition plate.cpp:58
FTensor::Index< 'm', 3 > m
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
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
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 derivatives 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)
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 field 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.
boost::shared_ptr< VectorDouble > detJacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpScaleBaseBySpaceInverseOfMeasure(const FieldSpace space, const FieldApproximationBase base, boost::shared_ptr< VectorDouble > det_jac_ptr, boost::shared_ptr< double > scale_det_ptr=nullptr)
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
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.
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.
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.
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 applyTransform(MatrixDouble &diff_n)
Apply transformation to the input matrix.
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian