v0.14.0
Loading...
Searching...
No Matches
HODataOperators.hpp
Go to the documentation of this file.
1/** \file HODataOperators.hpp
2 * \brief Operators managing HO geometry
3
4*/
5
6#ifndef __HO_DATA_OPERATORS_HPP__
7#define __HO_DATA_OPERATORS_HPP__
8
9namespace MoFEM {
10
11/**
12 * @brief Calculate jacobian on Hex or other volume which is not simplex.
13 *
14 * Using not a field data, but geometrical coordinates of element.
15 *
16 * \note Use OpCalculateVectorFieldGradient to calculate Jacobian from field
17 * data.
18 *
19 */
22
23 OpCalculateHOJacForVolume(boost::shared_ptr<MatrixDouble> jac_ptr);
24
25 MoFEMErrorCode doWork(int side, EntityType type,
27
28private:
29 boost::shared_ptr<MatrixDouble> jacPtr;
30};
31
32/** \deprecated use OpCalculateHOJacForVolume
33 */
35
36/**
37 * @brief Calculate HO coordinates at gauss points
38 *
39 */
40template <int FIELD_DIM = 3>
49
50/**
51 * @deprecated This class should be DIM = 3 specialization when default
52 * parameter is removed
53 *
54 */
62
63/**
64 * \brief Set inverse jacobian to base functions
65 *
66 * \deprecated Version with default variant DIM = 3 is deprecated, keep for
67 * back compatibility. Should be removed in future. Use of default variant make
68 * code implicit, what can be source of some hidden error. Explict interface is
69 * better, when user see and control outcome, and is aware of existing variants.
70 *
71 */
72template <int DIM = 3>
74 using OpSetHOInvJacToScalarBasesImpl::OpSetHOInvJacToScalarBasesImpl;
75};
76
77template <>
79 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
80
81 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
82};
83
84/**
85 * \brief transform local reference derivatives of shape function to global
86 derivatives if higher order geometry is given
87 *
88
89 * \ingroup mofem_forces_and_sources
90*/
92
94 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
95 : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {
96
97 if (!invJacPtr)
99 "Pointer for invJacPtr not allocated");
100
101 doVertices = false;
102 if (space == HDIV)
103 doEdges = false;
104 }
105
106
107 MoFEMErrorCode doWork(int side, EntityType type,
109
110private:
111 boost::shared_ptr<MatrixDouble> invJacPtr;
114};
115
116/**
117 * @brief Modify integration weights on face to take in account higher-order
118 * geometry
119 * @ingroup mofem_forces_and_sources_tri_element
120 *
121 */
129
130/**
131 * @brief Modify integration weights on face to take in account higher-order
132 * geometry
133 * @ingroup mofem_forces_and_sources_tri_element
134 *
135 */
143
144template <int SPACE_DIM>
146
150
154
155/**
156 * \brief Set inverse jacobian to base functions
157 *
158 */
160
161 OpSetHOWeights(boost::shared_ptr<VectorDouble> det_ptr)
163 detPtr(det_ptr) {
164 if (!detPtr)
166 "Pointer for detPtr not allocated");
167 }
168
169 MoFEMErrorCode doWork(int side, EntityType type,
171
172private:
173 boost::shared_ptr<VectorDouble> detPtr;
174};
175
176/** \brief Apply contravariant (Piola) transfer to Hdiv space for HO geometry
177
178* \ingroup mofem_forces_and_sources
179*/
182
184 boost::shared_ptr<VectorDouble> det_ptr,
185 boost::shared_ptr<MatrixDouble> jac_ptr)
187 jacPtr(jac_ptr) {
188 doVertices = false;
189 if (space == HDIV)
190 doEdges = false;
191 }
192
193 MoFEMErrorCode doWork(int side, EntityType type,
195
196private:
197 boost::shared_ptr<VectorDouble> detPtr;
198 boost::shared_ptr<MatrixDouble> jacPtr;
199
202};
203
204/** \brief Apply covariant (Piola) transfer to Hcurl space for HO geometry
205 */
208
210 boost::shared_ptr<MatrixDouble> jac_inv_ptr)
212 jacInvPtr(jac_inv_ptr) {
213 doVertices = false;
214 if (space == HDIV)
215 doEdges = false;
216 if (!jacInvPtr)
218 "Pointer for jacPtr not allocated");
219 }
220
221 MoFEMErrorCode doWork(int side, EntityType type,
223
224private:
225 boost::shared_ptr<MatrixDouble> jacInvPtr;
226
229};
230
231/** \brief Calculate jacobian for face element
232
233 It is assumed that face element is XY plane. Applied
234 only for 2d problems and 2d problems embedded in 3d space
235
236 \note If you push operators for HO normal befor this operator, HO geometry is
237 taken into account when you calculate jacobian.
238
239*/
240template <int DIM> struct OpCalculateHOJacForFaceImpl;
241
242template <>
245
246 OpCalculateHOJacForFaceImpl(boost::shared_ptr<MatrixDouble> jac_ptr);
247
248 MoFEMErrorCode doWork(int side, EntityType type,
250
251protected:
252 boost::shared_ptr<MatrixDouble> jacPtr;
253};
254
255template <>
257
258 using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
259
260 MoFEMErrorCode doWork(int side, EntityType type,
262};
263
266
267template <int DIM> struct OpCalculateHOJac;
268
272
273template <> struct OpCalculateHOJac<2> : public OpCalculateHOJacForFaceImpl<2> {
274 using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
275};
276
277/** \brief Calculate normals at Gauss points of triangle element
278 * \ingroup mofem_forces_and_source
279 */
280template <int FIELD_DIM = 3>
290
291/** \brief transform Hdiv base fluxes from reference element to physical
292 * triangle \ingroup mofem_forces_and_sources
293 *
294 * \note Matrix which keeps normal is assumed to have three columns, and number
295 * of rows should be equal to number of integration points.
296 *
297 */
300
302 const FieldSpace space,
303 boost::shared_ptr<MatrixDouble> normals_at_gauss_pts = nullptr)
305 normalsAtGaussPts(normals_at_gauss_pts) {}
306
307 MoFEMErrorCode doWork(int side, EntityType type,
309
310private:
311 boost::shared_ptr<MatrixDouble> normalsAtGaussPts;
312};
313
314/** \brief transform Hcurl base fluxes from reference element to physical edge
315 * \ingroup mofem_forces_and_sources
316 */
319
321 const FieldSpace space = HCURL,
322 boost::shared_ptr<MatrixDouble> tangent_at_pts = nullptr)
324 tangentAtGaussPts(tangent_at_pts) {}
325
326 MoFEMErrorCode doWork(int side, EntityType type,
328
329private:
330 boost::shared_ptr<MatrixDouble> tangentAtGaussPts;
331};
332
333/** \brief transform Hcurl base fluxes from reference element to physical
334 * triangle \ingroup mofem_forces_and_sources
335 */
338
340 const FieldSpace space,
341 boost::shared_ptr<MatrixDouble> normals_at_pts = nullptr,
342 boost::shared_ptr<MatrixDouble> tangent1_at_pts = nullptr,
343 boost::shared_ptr<MatrixDouble> tangent2_at_pts = nullptr)
345 normalsAtPts(normals_at_pts), tangent1AtPts(tangent1_at_pts),
346 tangent2AtPts(tangent2_at_pts) {
347 if (normals_at_pts || tangent1_at_pts || tangent2_at_pts)
348 if (normals_at_pts && tangent1_at_pts && tangent2_at_pts)
350 "All elements in constructor have to set pointer");
351 }
352
353 MoFEMErrorCode doWork(int side, EntityType type,
355
356private:
357 boost::shared_ptr<MatrixDouble> normalsAtPts;
358 boost::shared_ptr<MatrixDouble> tangent1AtPts;
359 boost::shared_ptr<MatrixDouble> tangent2AtPts;
360
363};
364
365/** \brief Calculate tangent vector on edge form HO geometry approximation
366 * \ingroup mofem_forces_and_sources
367 */
368template <int FIELD_DIM = 3>
371
373 std::string field_name,
374 boost::shared_ptr<MatrixDouble> tangents_at_pts = nullptr)
376 tangentsAtPts(tangents_at_pts) {}
377
378 MoFEMErrorCode doWork(int side, EntityType type,
380
381private:
382 boost::shared_ptr<MatrixDouble> tangentsAtPts;
383};
384
385/**
386 * @brief Scale base functions by inverses of measure of element
387 *
388 * @tparam OP
389 */
390template <typename OP> struct OpScaleBaseBySpaceInverseOfMeasure : public OP {
391
393 boost::shared_ptr<VectorDouble> det_jac_ptr, const FieldSpace space = L2)
394 : OP(space), fieldSpace(space), detJacPtr(det_jac_ptr) {}
395
399
400 if (!detJacPtr) {
401 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
402 }
403
404 auto scale = [&]() {
405 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
406 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
407 auto &base_fun = data.getN(base);
408 auto &diff_base_fun = data.getDiffN(base);
409 if (detJacPtr) {
410
411 auto &det_vec = *detJacPtr;
412 const auto nb_int_pts = base_fun.size1();
413
414 if (nb_int_pts != det_vec.size())
415 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
416 "Number of integration pts in detJacPtr does not mush "
417 "number of integration pts in base function");
418
419 auto get_row = [](auto &m, auto gg) {
420 return ublas::matrix_row<decltype(m)>(m, gg);
421 };
422
423 for (auto gg = 0; gg != nb_int_pts; ++gg)
424 get_row(base_fun) /= det_vec[gg];
425
426 if (diff_base_fun.size1() == nb_int_pts) {
427 for (auto gg = 0; gg != nb_int_pts; ++gg)
428 get_row(diff_base_fun) /= det_vec[gg];
429 }
430 }
431 }
432 };
433
434 if (this->getFEDim() == 3) {
435 switch (fieldSpace) {
436 case H1:
437 scale();
438 break;
439 case HCURL:
440 if (type >= MBEDGE)
441 scale();
442 break;
443 case HDIV:
444 if (type >= MBTRI)
445 scale();
446 break;
447 case L2:
448 if (type >= MBTET)
449 scale();
450 break;
451 default:
452 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
453 }
454 } else if (this->getFEDim() == 2) {
455 switch (fieldSpace) {
456 case H1:
457 scale();
458 break;
459 case HCURL:
460 if (type >= MBEDGE)
461 scale();
462 break;
463 case HDIV:
464 case L2:
465 if (type >= MBTRI)
466 scale();
467 break;
468 default:
469 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
470 }
471 } else if (this->getFEDim() == 1) {
472 switch (fieldSpace) {
473 case H1:
474 scale();
475 break;
476 case HCURL:
477 case L2:
478 if (type >= MBEDGE)
479 scale();
480 break;
481 default:
482 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
483 }
484 } else {
485 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
486 }
487
489 }
490
491private:
493 boost::shared_ptr<VectorDouble> detJacPtr;
494};
495
496/**
497 * @brief Add operators pushing bases from local to physical configuration
498 *
499 * @tparam FE_DIM dimension of element
500 * @tparam PROBLEM_DIM problem dimension
501 * @tparam SPACE_DIM space dimension
502 */
503template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
504
505template <> struct AddHOOps<2, 2, 2> {
506 AddHOOps() = delete;
507 static MoFEMErrorCode
508 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
509 std::vector<FieldSpace> spaces, std::string geom_field_name = "",
510 boost::shared_ptr<MatrixDouble> jac = nullptr,
511 boost::shared_ptr<VectorDouble> det = nullptr,
512 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
513};
514
515template <> struct AddHOOps<2, 2, 3> {
516 AddHOOps() = delete;
517 static MoFEMErrorCode
518 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
519 std::vector<FieldSpace> spaces, std::string geom_field_name = "",
520 boost::shared_ptr<MatrixDouble> jac = nullptr,
521 boost::shared_ptr<VectorDouble> det = nullptr,
522 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
523};
524
525template <> struct AddHOOps<1, 2, 2> {
526 AddHOOps() = delete;
527 static MoFEMErrorCode
528 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
529 std::vector<FieldSpace> spaces, std::string geom_field_name = "");
530};
531template <> struct AddHOOps<1, 3, 3> {
532 AddHOOps() = delete;
533 static MoFEMErrorCode
534 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
535 std::vector<FieldSpace> space, std::string geom_field_name = "");
536};
537
538template <> struct AddHOOps<2, 3, 3> {
539 AddHOOps() = delete;
540 static MoFEMErrorCode
541 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
542 std::vector<FieldSpace> space, std::string geom_field_name = "");
543};
544
545template <> struct AddHOOps<3, 3, 3> {
546 AddHOOps() = delete;
547 static MoFEMErrorCode
548 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
549 std::vector<FieldSpace> space, std::string geom_field_name = "",
550 boost::shared_ptr<MatrixDouble> jac = nullptr,
551 boost::shared_ptr<VectorDouble> det = nullptr,
552 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
553};
554
555template <int FIELD_DIM>
561 const auto nb_dofs = data.getFieldData().size() / FIELD_DIM;
562 if (nb_dofs) {
563 if (type == MBVERTEX)
564 getCoordsAtGaussPts().clear();
565 auto t_base = data.getFTensor0N();
566 auto t_coords = getFTensor1CoordsAtGaussPts();
567 const auto nb_integration_pts = data.getN().size1();
568 const auto nb_base_functions = data.getN().size2();
569 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
570 auto t_dof = data.getFTensor1FieldData<FIELD_DIM>();
571 size_t bb = 0;
572 for (; bb != nb_dofs; ++bb) {
573 t_coords(i) += t_base * t_dof(i);
574 ++t_dof;
575 ++t_base;
576 }
577 for (; bb < nb_base_functions; ++bb)
578 ++t_base;
579 ++t_coords;
580 }
581 }
583}
584
585template <int FIELD_DIM>
590
593
594 auto get_ftensor1 = [](MatrixDouble &m) {
596 &m(0, 0), &m(0, 1), &m(0, 2));
597 };
598
599 unsigned int nb_dofs = data.getFieldData().size();
600 if (nb_dofs != 0) {
601
602 int nb_gauss_pts = data.getN().size1();
603 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
604 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
605 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
606 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
607
608 switch (type) {
609 case MBVERTEX: {
610 tangent1_at_gauss_pts.clear();
611 tangent2_at_gauss_pts.clear();
612 }
613 case MBEDGE:
614 case MBTRI:
615 case MBQUAD: {
616
617#ifndef NDEBUG
618 if (2 * data.getN().size2() != data.getDiffN().size2()) {
619 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
620 "expected two direcatives in local element coordinates");
621 }
622 if (nb_dofs % FIELD_DIM != 0) {
623 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
624 "expected that number of dofs is multiplicative of field "
625 "dimension");
626 }
627#endif
628
629 if (nb_dofs > FIELD_DIM * data.getN().size2()) {
630 unsigned int nn = 0;
631 for (; nn != nb_dofs; nn++) {
632 if (!data.getFieldDofs()[nn]->getActive())
633 break;
634 }
635 if (nn > FIELD_DIM * data.getN().size2()) {
636 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
637 "Data inconsistency for base %s",
639 } else {
640 nb_dofs = nn;
641 if (!nb_dofs)
643 }
644 }
645 const int nb_base_functions = data.getN().size2();
646 auto t_base = data.getFTensor0N();
647 auto t_diff_base = data.getFTensor1DiffN<2>();
648 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
649 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
650 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
651 auto t_data = data.getFTensor1FieldData<FIELD_DIM>();
652 int bb = 0;
653 for (; bb != nb_dofs / FIELD_DIM; ++bb) {
655 t_t1(i) += t_data(i) * t_diff_base(N0);
656 t_t2(i) += t_data(i) * t_diff_base(N1);
657 ++t_data;
658 ++t_base;
659 ++t_diff_base;
660 }
661 for (; bb != nb_base_functions; ++bb) {
662 ++t_base;
663 ++t_diff_base;
664 }
665 ++t_t1;
666 ++t_t2;
667 }
668 } break;
669 default:
670 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
671 }
672 }
673
674 if (moab::CN::Dimension(type) == 2) {
675
676 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
677 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
678 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
679
680 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
681 normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
682
683 auto t_normal = get_ftensor1(normal_at_gauss_pts);
684 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
685 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
686 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
687 FTensor::Index<'i', 3> i;
688 FTensor::Index<'j', 3> j;
689 FTensor::Index<'k', 3> k;
690 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
691 ++t_normal;
692 ++t_t1;
693 ++t_t2;
694 }
695 }
696
698}
699
700template <int FIELD_DIM>
705
706 auto get_tangent = [&]() -> MatrixDouble & {
707 if (tangentsAtPts)
708 return *tangentsAtPts;
709 else
710 return getTangentAtGaussPts();
711 };
712
713 auto &tangent = get_tangent();
714 int nb_gauss_pts = getGaussPts().size2();
715
716 if (type == MBVERTEX) {
717 tangent.resize(nb_gauss_pts, 3, false);
718 tangent.clear();
719 }
720
721 int nb_dofs = data.getFieldData().size();
722 if (nb_dofs != 0) {
723 const int nb_base_functions = data.getN().size2();
724 double *diff_base_function = &*data.getDiffN().data().begin();
725 auto tangent_at_gauss_pts =
726 getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
727
729 int size = nb_dofs / FIELD_DIM;
730 if (nb_dofs % FIELD_DIM) {
731 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
732 "Inconsistent number of dofs and filed dimension");
733 }
734 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
735 auto field_data = data.getFTensor1FieldData<FIELD_DIM>();
736 int bb = 0;
737 for (; bb < size; ++bb) {
738 tangent_at_gauss_pts(i) += field_data(i) * (*diff_base_function);
739 ++field_data;
740 ++diff_base_function;
741 }
742 for (; bb != nb_base_functions; ++bb) {
743 ++diff_base_function;
744 }
745 ++tangent_at_gauss_pts;
746 }
747 }
749}
750
751/**
752 * @deprecated do not use this function, instead use AddHOOps.
753 *
754 * @tparam E
755 * @param field
756 * @param e
757 * @param h1
758 * @param hcurl
759 * @param hdiv
760 * @param l2
761 * @return MoFEMErrorCode
762 */
763template <typename E>
764MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1,
765 bool hcurl, bool hdiv, bool l2) {
766 std::vector<FieldSpace> spaces;
767 if (h1)
768 spaces.push_back(H1);
769 if (hcurl)
770 spaces.push_back(HCURL);
771 if (hdiv)
772 spaces.push_back(HDIV);
773 if (l2)
774 spaces.push_back(L2);
775 return AddHOOps<3, 3, 3>::add(e.getOpPtrVector(), spaces, field);
776}
777
778/**
779 * @deprecated do not use this function, instead use AddHOOps.
780 *
781 * @tparam E
782 * @param field
783 * @param e
784 * @param hcurl
785 * @param hdiv
786 * @return MoFEMErrorCode
787 */
788template <typename E>
789MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e,
790 bool hcurl, bool hdiv) {
791 std::vector<FieldSpace> spaces;
792 if (hcurl)
793 spaces.push_back(HCURL);
794 if (hdiv)
795 spaces.push_back(HDIV);
796 return AddHOOps<2, 3, 3>::add(e.getOpPtrVector(), spaces, field);
797}
798
799}; // namespace MoFEM
800
801#endif
static Number< 1 > N1
static Number< 0 > N0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int FIELD_DIM
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ USER_BASE
user implemented approximation base
Definition definitions.h:68
#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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
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 > >::typ levi_civita)(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
double scale
Definition plastic.cpp:170
constexpr auto field_name
Add operators pushing bases from local to physical configuration.
bool & doVertices
\deprectaed If false skip vertices
bool & doEdges
\deprectaed If false skip edges
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
structure to get information form mofem into EntitiesFieldData
Calculate HO coordinates at gauss points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHOCoords(const std::string field_name)
boost::shared_ptr< MatrixDouble > jacPtr
Calculate jacobian for face element.
Calculate jacobian on Hex or other volume which is not simplex.
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)
Calculate normals at Gauss points of triangle element.
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.
Calculate tangent vector on edge form HO geometry approximation.
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.
OpGetHOTangentsOnEdge(std::string field_name, boost::shared_ptr< MatrixDouble > tangents_at_pts=nullptr)
transform Hcurl base fluxes from reference element to physical edge
boost::shared_ptr< MatrixDouble > tangentAtGaussPts
OpHOSetContravariantPiolaTransformOnEdge3D(const FieldSpace space=HCURL, boost::shared_ptr< MatrixDouble > tangent_at_pts=nullptr)
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.
OpHOSetContravariantPiolaTransformOnFace3D(const FieldSpace space, boost::shared_ptr< MatrixDouble > normals_at_gauss_pts=nullptr)
transform Hcurl base fluxes from reference element to physical triangle
OpHOSetCovariantPiolaTransformOnFace3D(const FieldSpace space, boost::shared_ptr< MatrixDouble > normals_at_pts=nullptr, boost::shared_ptr< MatrixDouble > tangent1_at_pts=nullptr, boost::shared_ptr< MatrixDouble > tangent2_at_pts=nullptr)
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
Scale base functions by inverses of measure of element.
OpScaleBaseBySpaceInverseOfMeasure(boost::shared_ptr< VectorDouble > det_jac_ptr, const FieldSpace space=L2)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > detJacPtr
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
OpSetHOContravariantPiolaTransform(const FieldSpace space, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > jac_ptr)
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
OpSetHOCovariantPiolaTransform(const FieldSpace space, boost::shared_ptr< MatrixDouble > jac_inv_ptr)
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.
OpSetHOInvJacVectorBase(const FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Set inverse jacobian to base functions.
OpSetHOWeights(boost::shared_ptr< VectorDouble > det_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > detPtr
Modify integration weights on face to take in account higher-order geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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.
Transform local reference derivatives of shape functions to global derivatives.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)