v0.13.1
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
7/* This file is part of MoFEM.
8 * MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20
21namespace MoFEM {
22
24 boost::shared_ptr<MatrixDouble> jac_ptr)
26 jacPtr(jac_ptr) {
27
28 for (auto t = MBEDGE; t != MBMAXTYPE; ++t)
29 doEntities[t] = false;
30
31 if (!jacPtr)
32 THROW_MESSAGE("Jac pointer not allocated");
33}
34
39
40 const auto nb_base_functions = data.getN(NOBASE).size2();
41 if (nb_base_functions) {
42
43 const auto nb_gauss_pts = getGaussPts().size2();
44 auto t_diff_base = data.getFTensor1DiffN<3>(NOBASE);
45 auto &coords = getCoords();
46
47#ifndef NDEBUG
48 if (nb_gauss_pts != data.getDiffN().size1())
49 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
50 "Inconsistent number base functions and gauss points");
51 if (nb_base_functions != data.getDiffN().size2() / 3)
52 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
53 "Inconsistent number of base functions");
54 if (coords.size() != 3 * nb_base_functions)
55 SETERRQ2(
56 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
57 "Number of vertex coordinates and base functions is inconsistent "
58 "%d != %d",
59 coords.size() / 3, nb_base_functions);
60#endif
61
62 jacPtr->resize(9, nb_gauss_pts, false);
63 jacPtr->clear();
64 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
65 auto t_vol_inv_jac = getInvJac();
66
70
71 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
73 &coords[0], &coords[1], &coords[2]);
74 for (size_t bb = 0; bb != nb_base_functions; ++bb) {
75 t_jac(i, j) += t_coords(i) * (t_vol_inv_jac(k, j) * t_diff_base(k));
76 ++t_diff_base;
77 ++t_coords;
78 }
79 ++t_jac;
80 }
81 }
82
84}
85
91 const auto nb_dofs = data.getFieldData().size() / 3;
92 if (nb_dofs) {
93 if (type == MBVERTEX)
94 getCoordsAtGaussPts().clear();
95 auto t_base = data.getFTensor0N();
96 auto t_coords = getFTensor1CoordsAtGaussPts();
97 const auto nb_integration_pts = data.getN().size1();
98 const auto nb_base_functions = data.getN().size2();
99 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
100 auto t_dof = data.getFTensor1FieldData<3>();
101 size_t bb = 0;
102 for (; bb != nb_dofs; ++bb) {
103 t_coords(i) += t_base * t_dof(i);
104 ++t_dof;
105 ++t_base;
106 }
107 for (; bb < nb_base_functions; ++bb)
108 ++t_base;
109 ++t_coords;
110 }
111 }
113};
114
119
120 if (getFEDim() == 3) {
121
122 auto transform_base = [&](MatrixDouble &diff_n) {
124 return applyTransform<3, 3, 3, 3>(diff_n);
126 };
127
128 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
130 CHKERR transform_base(data.getDiffN(base));
131 }
132
133 switch (type) {
134 case MBVERTEX:
135 for (auto &m : data.getBBDiffNMap())
136 CHKERR transform_base(*(m.second));
137 break;
138 default:
139 for (auto &ptr : data.getBBDiffNByOrderArray())
140 if (ptr)
141 CHKERR transform_base(*ptr);
142 }
143
144 } else {
145 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Use different operator");
146 }
147
149}
150
155
159
160 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
161
163
164 diffHdivInvJac.resize(data.getDiffN(base).size1(),
165 data.getDiffN(base).size2(), false);
166
167 unsigned int nb_gauss_pts = data.getDiffN(base).size1();
168 unsigned int nb_base_functions = data.getDiffN(base).size2() / 9;
169 if (nb_base_functions == 0)
170 continue;
171
172 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
173 double *t_inv_diff_n_ptr = &*diffHdivInvJac.data().begin();
175 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
176 &t_inv_diff_n_ptr[HVEC0_2], &t_inv_diff_n_ptr[HVEC1_0],
177 &t_inv_diff_n_ptr[HVEC1_1], &t_inv_diff_n_ptr[HVEC1_2],
178 &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1],
179 &t_inv_diff_n_ptr[HVEC2_2]);
180 auto t_inv_jac = getFTensor2FromMat<3, 3>(*invJacPtr);
181 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
182 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
183 t_inv_diff_n(i, j) = t_inv_jac(k, j) * t_diff_n(i, k);
184 ++t_diff_n;
185 ++t_inv_diff_n;
186 }
187 ++t_inv_jac;
188 }
189
190 data.getDiffN(base).swap(diffHdivInvJac);
191 }
192
194}
195
200 const size_t nb_int_pts = getGaussPts().size2();
201 if (getNormalsAtGaussPts().size1()) {
202 if (getNormalsAtGaussPts().size1() == nb_int_pts) {
203 double a = getMeasure();
204 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI)
205 a *= 2;
206 auto t_w = getFTensor0IntegrationWeight();
207 auto t_normal = getFTensor1NormalsAtGaussPts();
209 for (size_t gg = 0; gg != nb_int_pts; ++gg) {
210 t_w *= sqrt(t_normal(i) * t_normal(i)) / a;
211 ++t_w;
212 ++t_normal;
213 }
214 } else {
215 SETERRQ2(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
216 "Number of rows in getNormalsAtGaussPts should be equal to "
217 "number of integration points, but is not, i.e. %d != %d",
218 getNormalsAtGaussPts().size1(), nb_int_pts);
219 }
220 }
222}
223
227
228 const auto nb_integration_pts = detPtr->size();
229
230#ifndef NDEBUG
231 if (nb_integration_pts != getGaussPts().size2())
232 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
233 "Inconsistent number of data points");
234#endif
235
236 auto t_w = getFTensor0IntegrationWeight();
237 auto t_det = getFTensor0FromVec(*detPtr);
238 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
239 t_w *= t_det;
240 ++t_w;
241 ++t_det;
242 }
243
245}
246
248 int side, EntityType type, EntitiesFieldData::EntData &data) {
250
254
255#ifndef NDEBUG
256 if (!detPtr)
257 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
258 "Pointer for detPtr not allocated");
259 if (!jacPtr)
260 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
261 "Pointer for jacPtr not allocated");
262#endif
263
264 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
265
267
268 auto nb_gauss_pts = data.getN(base).size1();
269 auto nb_base_functions = data.getN(base).size2() / 3;
270
271#ifndef NDEBUG
272 if (data.getDiffN(base).size1() != nb_gauss_pts)
273 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
274 "Wrong number integration points");
275
276 if (data.getDiffN(base).size2() / 9 != nb_base_functions)
277 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
278 "Wrong number base functions %d != %d",
279 data.getDiffN(base).size2(), nb_base_functions);
280#endif
281
282 if (nb_gauss_pts && nb_base_functions) {
283
284 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
285 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
286
287 auto t_n = data.getFTensor1N<3>(base);
288 double *t_transformed_n_ptr = &*piolaN.data().begin();
290 t_transformed_n_ptr, // HVEC0
291 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
292 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
293 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
295 t_transformed_diff_n(t_transformed_diff_n_ptr,
296 &t_transformed_diff_n_ptr[HVEC0_1],
297 &t_transformed_diff_n_ptr[HVEC0_2],
298 &t_transformed_diff_n_ptr[HVEC1_0],
299 &t_transformed_diff_n_ptr[HVEC1_1],
300 &t_transformed_diff_n_ptr[HVEC1_2],
301 &t_transformed_diff_n_ptr[HVEC2_0],
302 &t_transformed_diff_n_ptr[HVEC2_1],
303 &t_transformed_diff_n_ptr[HVEC2_2]);
304
305 auto t_det = getFTensor0FromVec(*detPtr);
306 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
307
308 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
309 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
310 const double a = 1. / t_det;
311 t_transformed_n(i) = a * t_jac(i, k) * t_n(k);
312 t_transformed_diff_n(i, k) = a * t_jac(i, j) * t_diff_n(j, k);
313 ++t_n;
314 ++t_transformed_n;
315 ++t_diff_n;
316 ++t_transformed_diff_n;
317 }
318 ++t_det;
319 ++t_jac;
320 }
321
322 data.getN(base).swap(piolaN);
323 data.getDiffN(base).swap(piolaDiffN);
324 }
325 }
326
328}
329
331 int side, EntityType type, EntitiesFieldData::EntData &data) {
333
337
338 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
339
341
342 unsigned int nb_gauss_pts = data.getN(base).size1();
343 unsigned int nb_base_functions = data.getN(base).size2() / 3;
344 piolaN.resize(nb_gauss_pts, data.getN(base).size2(), false);
345 piolaDiffN.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
346
347 auto t_n = data.getFTensor1N<3>(base);
348 double *t_transformed_n_ptr = &*piolaN.data().begin();
350 t_transformed_n_ptr, // HVEC0
351 &t_transformed_n_ptr[HVEC1], &t_transformed_n_ptr[HVEC2]);
352 auto t_diff_n = data.getFTensor2DiffN<3, 3>(base);
353 double *t_transformed_diff_n_ptr = &*piolaDiffN.data().begin();
354 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_transformed_diff_n(
355 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[HVEC0_1],
356 &t_transformed_diff_n_ptr[HVEC0_2], &t_transformed_diff_n_ptr[HVEC1_0],
357 &t_transformed_diff_n_ptr[HVEC1_1], &t_transformed_diff_n_ptr[HVEC1_2],
358 &t_transformed_diff_n_ptr[HVEC2_0], &t_transformed_diff_n_ptr[HVEC2_1],
359 &t_transformed_diff_n_ptr[HVEC2_2]);
360
361 auto t_inv_jac = getFTensor2FromMat<3, 3>(*jacInvPtr);
362
363 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
364 for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
365 t_transformed_n(i) = t_inv_jac(k, i) * t_n(k);
366 t_transformed_diff_n(i, k) = t_inv_jac(j, i) * t_diff_n(j, k);
367 ++t_n;
368 ++t_transformed_n;
369 ++t_diff_n;
370 ++t_transformed_diff_n;
371 }
372 ++t_inv_jac;
373 }
374
375 data.getN(base).swap(piolaN);
376 data.getDiffN(base).swap(piolaDiffN);
377 }
378
380}
381
383 boost::shared_ptr<MatrixDouble> jac_ptr)
385 jacPtr(jac_ptr) {}
386
388 int side, EntityType type, EntitiesFieldData::EntData &data) {
389
391
392 const auto nb_gauss_pts = getGaussPts().size2();
393
394 jacPtr->resize(4, nb_gauss_pts, false);
395 auto t_jac = getFTensor2FromMat<2, 2>(*jacPtr);
396
399
400 auto t_t1 = getFTensor1Tangent1AtGaussPts();
401 auto t_t2 = getFTensor1Tangent2AtGaussPts();
404
405 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
406 t_jac(i, N0) = t_t1(i);
407 t_jac(i, N1) = t_t2(i);
408 ++t_t1;
409 ++t_t2;
410 ++t_jac;
411 }
412
414};
415
417 int side, EntityType type, EntitiesFieldData::EntData &data) {
419
423
424 size_t nb_gauss_pts = getGaussPts().size2();
425 jacPtr->resize(9, nb_gauss_pts, false);
426
427 auto t_t1 = getFTensor1Tangent1AtGaussPts();
428 auto t_t2 = getFTensor1Tangent2AtGaussPts();
429 auto t_normal = getFTensor1NormalsAtGaussPts();
430
434
435 auto t_jac = getFTensor2FromMat<3, 3>(*jacPtr);
436 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_jac) {
437
438 t_jac(i, N0) = t_t1(i);
439 t_jac(i, N1) = t_t2(i);
440
441 const double l = sqrt(t_normal(j) * t_normal(j));
442
443 t_jac(i, N2) = t_normal(i) / l;
444
445 ++t_t1;
446 ++t_t2;
447 ++t_normal;
448 }
449
451}
452
455}
456
461
465
468
469 auto get_ftensor1 = [](MatrixDouble &m) {
471 &m(0, 0), &m(0, 1), &m(0, 2));
472 };
473
474 unsigned int nb_dofs = data.getFieldData().size();
475 if (nb_dofs != 0) {
476
477 int nb_gauss_pts = data.getN().size1();
478 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
479 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
480 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
481 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
482
483 switch (type) {
484 case MBVERTEX: {
485 tangent1_at_gauss_pts.clear();
486 tangent2_at_gauss_pts.clear();
487 }
488 case MBEDGE:
489 case MBTRI:
490 case MBQUAD: {
491
492#ifndef NDEBUG
493 if (2 * data.getN().size2() != data.getDiffN().size2()) {
494 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
495 "data inconsistency");
496 }
497 if (nb_dofs % 3 != 0) {
498 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
499 "data inconsistency");
500 }
501#endif
502
503 if (nb_dofs > 3 * data.getN().size2()) {
504 unsigned int nn = 0;
505 for (; nn != nb_dofs; nn++) {
506 if (!data.getFieldDofs()[nn]->getActive())
507 break;
508 }
509 if (nn > 3 * data.getN().size2()) {
510 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511 "data inconsistency for base %s",
513 } else {
514 nb_dofs = nn;
515 if (!nb_dofs)
517 }
518 }
519 const int nb_base_functions = data.getN().size2();
520 auto t_base = data.getFTensor0N();
521 auto t_diff_base = data.getFTensor1DiffN<2>();
522 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
523 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
524 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
525 auto t_data = data.getFTensor1FieldData<3>();
526 int bb = 0;
527 for (; bb != nb_dofs / 3; ++bb) {
528 t_t1(i) += t_data(i) * t_diff_base(N0);
529 t_t2(i) += t_data(i) * t_diff_base(N1);
530 ++t_data;
531 ++t_base;
532 ++t_diff_base;
533 }
534 for (; bb != nb_base_functions; ++bb) {
535 ++t_base;
536 ++t_diff_base;
537 }
538 ++t_t1;
539 ++t_t2;
540 }
541 } break;
542 default:
543 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
544 }
545 }
546
547 if (moab::CN::Dimension(type) == 2) {
548
549 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
550 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
551 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
552
553 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
554 normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
555
556 auto t_normal = get_ftensor1(normal_at_gauss_pts);
557 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
558 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
559 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
560 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
561 ++t_normal;
562 ++t_t1;
563 ++t_t2;
564 }
565 }
566
568}
569
571 int side, EntityType type, EntitiesFieldData::EntData &data) {
574
575 if (moab::CN::Dimension(type) != 2)
577
578 auto get_normals_ptr = [&]() {
580 return &*normalsAtGaussPts->data().begin();
581 else
582 return &*getNormalsAtGaussPts().data().begin();
583 };
584
585 auto apply_transform_nonlinear_geometry = [&](auto base, auto nb_gauss_pts,
586 auto nb_base_functions) {
588
589 auto ptr = get_normals_ptr();
591 &ptr[0], &ptr[1], &ptr[2]);
592
593 auto t_base = data.getFTensor1N<3>(base);
594 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
595 const auto l2 = t_normal(i) * t_normal(i);
596 for (int bb = 0; bb != nb_base_functions; ++bb) {
597 const auto v = t_base(0);
598 t_base(i) = (v / l2) * t_normal(i);
599 ++t_base;
600 }
601 ++t_normal;
602 }
604 };
605
606 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
607
609 const auto &base_functions = data.getN(base);
610 const auto nb_gauss_pts = base_functions.size1();
611
612 if (nb_gauss_pts) {
613
614 const auto nb_base_functions = base_functions.size2() / 3;
615 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
616 nb_base_functions);
617 }
618 }
619
621}
622
624 int side, EntityType type, EntitiesFieldData::EntData &data) {
626
627 const auto type_dim = moab::CN::Dimension(type);
628 if (type_dim != 1 && type_dim != 2)
630
634
635 auto get_jac = [&]() {
637 double *ptr_n = &*normalsAtPts->data().begin();
638 double *ptr_t1 = &*tangent1AtPts->data().begin();
639 double *ptr_t2 = &*tangent2AtPts->data().begin();
641 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
642
643 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
644
645 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
646 } else {
647 double *ptr_n = &*getNormalsAtGaussPts().data().begin();
648 double *ptr_t1 = &*getTangent1AtGaussPts().data().begin();
649 double *ptr_t2 = &*getTangent2AtGaussPts().data().begin();
651 &ptr_t1[0], &ptr_t2[0], &ptr_n[0],
652
653 &ptr_t1[1], &ptr_t2[1], &ptr_n[1],
654
655 &ptr_t1[2], &ptr_t2[2], &ptr_n[2]);
656 }
657 };
658
659 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; ++b) {
660
662
663 auto &baseN = data.getN(base);
664 auto &diffBaseN = data.getDiffN(base);
665
666 int nb_dofs = baseN.size2() / 3;
667 int nb_gauss_pts = baseN.size1();
668
669 piolaN.resize(baseN.size1(), baseN.size2());
670 diffPiolaN.resize(diffBaseN.size1(), diffBaseN.size2());
671
672 if (nb_dofs > 0 && nb_gauss_pts > 0) {
673
674 auto t_h_curl = data.getFTensor1N<3>(base);
675 auto t_diff_h_curl = data.getFTensor2DiffN<3, 2>(base);
676
677 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_transformed_h_curl(
678 &piolaN(0, HVEC0), &piolaN(0, HVEC1), &piolaN(0, HVEC2));
679
681 t_transformed_diff_h_curl(
685
686 int cc = 0;
687 double det;
689
690 // HO geometry is set, so jacobian is different at each gauss point
691 auto t_m_at_pts = get_jac();
692 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
693 CHKERR determinantTensor3by3(t_m_at_pts, det);
694 CHKERR invertTensor3by3(t_m_at_pts, det, t_inv_m);
695 for (int ll = 0; ll != nb_dofs; ll++) {
696 t_transformed_h_curl(i) = t_inv_m(j, i) * t_h_curl(j);
697 t_transformed_diff_h_curl(i, k) = t_inv_m(j, i) * t_diff_h_curl(j, k);
698 ++t_h_curl;
699 ++t_transformed_h_curl;
700 ++t_diff_h_curl;
701 ++t_transformed_diff_h_curl;
702 ++cc;
703 }
704 ++t_m_at_pts;
705 }
706
707#ifndef NDEBUG
708 if (cc != nb_gauss_pts * nb_dofs)
709 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
710#endif
711
712 baseN.swap(piolaN);
713 diffBaseN.swap(diffPiolaN);
714 }
715 }
716
718}
719
721 int side, EntityType type, EntitiesFieldData::EntData &data) {
724
725 if (type != MBEDGE)
727
728 const auto nb_gauss_pts = getGaussPts().size2();
729
731 if (tangentAtGaussPts->size1() != nb_gauss_pts)
732 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
733 "Wrong number of integration points %d!=%d",
734 tangentAtGaussPts->size1(), nb_gauss_pts);
735
736 auto get_base_at_pts = [&](auto base) {
738 &data.getN(base)(0, HVEC0), &data.getN(base)(0, HVEC1),
739 &data.getN(base)(0, HVEC2));
740 return t_h_curl;
741 };
742
743 auto get_tangent_ptr = [&]() {
744 if (tangentAtGaussPts) {
745 return &*tangentAtGaussPts->data().begin();
746 } else {
747 return &*getTangetAtGaussPts().data().begin();
748 }
749 };
750
751 auto get_tangent_at_pts = [&]() {
752 auto ptr = get_tangent_ptr();
754 &ptr[0], &ptr[1], &ptr[2]);
755 return t_m_at_pts;
756 };
757
758 auto calculate_squared_edge_length = [&]() {
759 std::vector<double> l1;
760 if (nb_gauss_pts) {
761 l1.resize(nb_gauss_pts);
762 auto t_m_at_pts = get_tangent_at_pts();
763 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
764 l1[gg] = t_m_at_pts(i) * t_m_at_pts(i);
765 ++t_m_at_pts;
766 }
767 }
768 return l1;
769 };
770
771 auto l1 = calculate_squared_edge_length();
772
773 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
774
776 const auto nb_dofs = data.getN(base).size2() / 3;
777
778 if (nb_gauss_pts && nb_dofs) {
779 auto t_h_curl = get_base_at_pts(base);
780 int cc = 0;
781 auto t_m_at_pts = get_tangent_at_pts();
782 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
783 const double l0 = l1[gg];
784 for (int ll = 0; ll != nb_dofs; ++ll) {
785 const double val = t_h_curl(0);
786 const double a = val / l0;
787 t_h_curl(i) = t_m_at_pts(i) * a;
788 ++t_h_curl;
789 ++cc;
790 }
791 ++t_m_at_pts;
792 }
793
794#ifndef NDEBUG
795 if (cc != nb_gauss_pts * nb_dofs)
796 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "Data inconsistency");
797#endif
798 }
799 }
800
802}
803
808
809 int nb_dofs = data.getFieldData().size();
810 if (nb_dofs == 0)
812
813 auto get_tangent = [&]() -> MatrixDouble & {
814 if (tangentsAtPts)
815 return *tangentsAtPts;
816 else
817 return getTangetAtGaussPts();
818 };
819
820 auto &tangent = get_tangent();
821
822 int nb_gauss_pts = data.getN().size1();
823 tangent.resize(nb_gauss_pts, 3, false);
824
825 int nb_approx_fun = data.getN().size2();
826 double *diff = &*data.getDiffN().data().begin();
827 double *dofs[] = {&data.getFieldData()[0], &data.getFieldData()[1],
828 &data.getFieldData()[2]};
829
830 tangent.resize(nb_gauss_pts, 3, false);
831
832 switch (type) {
833 case MBVERTEX:
834 for (int dd = 0; dd != 3; dd++) {
835 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
836 tangent(gg, dd) = cblas_ddot(2, diff, 1, dofs[dd], 3);
837 }
838 }
839 break;
840 case MBEDGE:
841#ifndef NDEBUG
842 if (nb_dofs % 3) {
843 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
844 "Approximated field should be rank 3, i.e. vector in 3d space");
845 }
846#endif
847 for (int dd = 0; dd != 3; dd++) {
848 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
849 tangent(gg, dd) +=
850 cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[dd], 3);
851 }
852 }
853 break;
854 default:
855 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE,
856 "This operator can calculate tangent vector only on edge");
857 }
858
860}
861
862} // 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:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ NOBASE
Definition: definitions.h:72
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
@ NOSPACE
Definition: definitions.h:96
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ HVEC0
Definition: definitions.h:199
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
@ HVEC1_1
Definition: definitions.h:209
@ HVEC0_1
Definition: definitions.h:208
@ HVEC1_0
Definition: definitions.h:206
@ HVEC2_1
Definition: definitions.h:210
@ HVEC1_2
Definition: definitions.h:212
@ HVEC2_2
Definition: definitions.h:213
@ HVEC2_0
Definition: definitions.h:207
@ HVEC0_2
Definition: definitions.h:211
@ HVEC0_0
Definition: definitions.h:205
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
FTensor::Index< 'i', SPACE_DIM > i
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
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:1218
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr auto field_name
FTensor::Index< 'm', 3 > m
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
MatrixDouble & getTangetAtGaussPts()
get tangent vector to edge curve at integration points
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FieldApproximationBase & getBase()
Get approximation base.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
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.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
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 getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
int getFEDim() const
Get dimension of finite element.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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)
OpGetHONormalsOnFace(std::string field_name)
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 > tangentsAtPts
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 > tangentAtGaussPts
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 > normalsAtGaussPts
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 > 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
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
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.
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.
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.
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian