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>
42
45
46 MoFEMErrorCode doWork(int side, EntityType type,
48};
49
50/**
51 * @deprecated This class should be DIM = 3 specialization when default
52 * parameter is removed
53 *
54 */
56
58
59 MoFEMErrorCode doWork(int side, EntityType type,
61};
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 */
126 MoFEMErrorCode doWork(int side, EntityType type,
128};
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 */
140 MoFEMErrorCode doWork(int side, EntityType type,
142};
143
144template <int SPACE_DIM>
146
147template <> struct OpSetHOWeightsOnSubDim<2> : public OpSetHOWeightsOnEdge {
149};
150
151template <> struct OpSetHOWeightsOnSubDim<3> : public OpSetHOWeightsOnFace {
153};
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
269template <> struct OpCalculateHOJac<3> : public OpCalculateHOJacForVolume {
271};
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>
283
286
287 MoFEMErrorCode doWork(int side, EntityType type,
289};
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<1, 2, 2> {
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};
521template <> struct AddHOOps<1, 3, 3> {
522 AddHOOps() = delete;
523 static MoFEMErrorCode
524 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
525 std::vector<FieldSpace> space, std::string geom_field_name = "");
526};
527
528template <> struct AddHOOps<2, 3, 3> {
529 AddHOOps() = delete;
530 static MoFEMErrorCode
531 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
532 std::vector<FieldSpace> space, std::string geom_field_name = "");
533};
534
535template <> struct AddHOOps<3, 3, 3> {
536 AddHOOps() = delete;
537 static MoFEMErrorCode
538 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
539 std::vector<FieldSpace> space, std::string geom_field_name = "",
540 boost::shared_ptr<MatrixDouble> jac = nullptr,
541 boost::shared_ptr<VectorDouble> det = nullptr,
542 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
543};
544
545template <int FIELD_DIM>
551 const auto nb_dofs = data.getFieldData().size() / FIELD_DIM;
552 if (nb_dofs) {
553 if (type == MBVERTEX)
554 getCoordsAtGaussPts().clear();
555 auto t_base = data.getFTensor0N();
556 auto t_coords = getFTensor1CoordsAtGaussPts();
557 const auto nb_integration_pts = data.getN().size1();
558 const auto nb_base_functions = data.getN().size2();
559 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
560 auto t_dof = data.getFTensor1FieldData<FIELD_DIM>();
561 size_t bb = 0;
562 for (; bb != nb_dofs; ++bb) {
563 t_coords(i) += t_base * t_dof(i);
564 ++t_dof;
565 ++t_base;
566 }
567 for (; bb < nb_base_functions; ++bb)
568 ++t_base;
569 ++t_coords;
570 }
571 }
573}
574
575template <int FIELD_DIM>
580
583
584 auto get_ftensor1 = [](MatrixDouble &m) {
586 &m(0, 0), &m(0, 1), &m(0, 2));
587 };
588
589 unsigned int nb_dofs = data.getFieldData().size();
590 if (nb_dofs != 0) {
591
592 int nb_gauss_pts = data.getN().size1();
593 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
594 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
595 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
596 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
597
598 switch (type) {
599 case MBVERTEX: {
600 tangent1_at_gauss_pts.clear();
601 tangent2_at_gauss_pts.clear();
602 }
603 case MBEDGE:
604 case MBTRI:
605 case MBQUAD: {
606
607#ifndef NDEBUG
608 if (2 * data.getN().size2() != data.getDiffN().size2()) {
609 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
610 "expected two direcatives in local element coordinates");
611 }
612 if (nb_dofs % FIELD_DIM != 0) {
613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614 "expected that number of dofs is multiplicative of field "
615 "dimension");
616 }
617#endif
618
619 if (nb_dofs > FIELD_DIM * data.getN().size2()) {
620 unsigned int nn = 0;
621 for (; nn != nb_dofs; nn++) {
622 if (!data.getFieldDofs()[nn]->getActive())
623 break;
624 }
625 if (nn > FIELD_DIM * data.getN().size2()) {
626 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
627 "Data inconsistency for base %s",
629 } else {
630 nb_dofs = nn;
631 if (!nb_dofs)
633 }
634 }
635 const int nb_base_functions = data.getN().size2();
636 auto t_base = data.getFTensor0N();
637 auto t_diff_base = data.getFTensor1DiffN<2>();
638 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
639 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
640 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
641 auto t_data = data.getFTensor1FieldData<FIELD_DIM>();
642 int bb = 0;
643 for (; bb != nb_dofs / FIELD_DIM; ++bb) {
645 t_t1(i) += t_data(i) * t_diff_base(N0);
646 t_t2(i) += t_data(i) * t_diff_base(N1);
647 ++t_data;
648 ++t_base;
649 ++t_diff_base;
650 }
651 for (; bb != nb_base_functions; ++bb) {
652 ++t_base;
653 ++t_diff_base;
654 }
655 ++t_t1;
656 ++t_t2;
657 }
658 } break;
659 default:
660 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
661 }
662 }
663
664 if (moab::CN::Dimension(type) == 2) {
665
666 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
667 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
668 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
669
670 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
671 normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
672
673 auto t_normal = get_ftensor1(normal_at_gauss_pts);
674 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
675 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
676 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
680 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
681 ++t_normal;
682 ++t_t1;
683 ++t_t2;
684 }
685 }
686
688}
689
690template <int FIELD_DIM>
695
696 auto get_tangent = [&]() -> MatrixDouble & {
697 if (tangentsAtPts)
698 return *tangentsAtPts;
699 else
700 return getTangentAtGaussPts();
701 };
702
703 auto &tangent = get_tangent();
704 int nb_gauss_pts = getGaussPts().size2();
705
706 if (type == MBVERTEX) {
707 tangent.resize(nb_gauss_pts, 3, false);
708 tangent.clear();
709 }
710
711 int nb_dofs = data.getFieldData().size();
712 if (nb_dofs != 0) {
713 const int nb_base_functions = data.getN().size2();
714 double *diff_base_function = &*data.getDiffN().data().begin();
715 auto tangent_at_gauss_pts =
716 getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
717
719 int size = nb_dofs / FIELD_DIM;
720 if (nb_dofs % FIELD_DIM) {
721 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
722 "Inconsistent number of dofs and filed dimension");
723 }
724 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
725 auto field_data = data.getFTensor1FieldData<FIELD_DIM>();
726 int bb = 0;
727 for (; bb < size; ++bb) {
728 tangent_at_gauss_pts(i) += field_data(i) * (*diff_base_function);
729 ++field_data;
730 ++diff_base_function;
731 }
732 for (; bb != nb_base_functions; ++bb) {
733 ++diff_base_function;
734 }
735 ++tangent_at_gauss_pts;
736 }
737 }
739}
740
741/**
742 * @deprecated do not use this function, instead use AddHOOps.
743 *
744 * @tparam E
745 * @param field
746 * @param e
747 * @param h1
748 * @param hcurl
749 * @param hdiv
750 * @param l2
751 * @return MoFEMErrorCode
752 */
753template <typename E>
754MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1,
755 bool hcurl, bool hdiv, bool l2) {
756 std::vector<FieldSpace> spaces;
757 if (h1)
758 spaces.push_back(H1);
759 if (hcurl)
760 spaces.push_back(HCURL);
761 if (hdiv)
762 spaces.push_back(HDIV);
763 if (l2)
764 spaces.push_back(L2);
765 return AddHOOps<3, 3, 3>::add(e.getOpPtrVector(), spaces, field);
766}
767
768/**
769 * @deprecated do not use this function, instead use AddHOOps.
770 *
771 * @tparam E
772 * @param field
773 * @param e
774 * @param hcurl
775 * @param hdiv
776 * @return MoFEMErrorCode
777 */
778template <typename E>
779MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e,
780 bool hcurl, bool hdiv) {
781 std::vector<FieldSpace> spaces;
782 if (hcurl)
783 spaces.push_back(HCURL);
784 if (hdiv)
785 spaces.push_back(HDIV);
786 return AddHOOps<2, 3, 3>::add(e.getOpPtrVector(), spaces, field);
787}
788
789}; // namespace MoFEM
790
791#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.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
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 ...
Definition: definitions.h:346
@ 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()
Definition: definitions.h:416
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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 > >::type 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.
Definition: Exceptions.hpp:56
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)