v0.13.2
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 MoFEMErrorCode doWork(int side, EntityType type,
108
109private:
110 boost::shared_ptr<MatrixDouble> invJacPtr;
113};
114
115/**
116 * @brief Modify integration weights on face to take in account higher-order
117 * geometry
118 * @ingroup mofem_forces_and_sources_tri_element
119 *
120 */
125 MoFEMErrorCode doWork(int side, EntityType type,
127};
128
129/**
130 * \brief Set inverse jacobian to base functions
131 *
132 */
134
135 OpSetHOWeights(boost::shared_ptr<VectorDouble> det_ptr)
137 detPtr(det_ptr) {
138 if (!detPtr)
140 "Pointer for detPtr not allocated");
141 }
142
143 MoFEMErrorCode doWork(int side, EntityType type,
145
146private:
147 boost::shared_ptr<VectorDouble> detPtr;
148};
149
150/** \brief Apply contravariant (Piola) transfer to Hdiv space for HO geometry
151
152* \ingroup mofem_forces_and_sources
153*/
156
158 boost::shared_ptr<VectorDouble> det_ptr,
159 boost::shared_ptr<MatrixDouble> jac_ptr)
161 jacPtr(jac_ptr) {
162 doVertices = false;
163 if (space == HDIV)
164 doEdges = false;
165 }
166
167 MoFEMErrorCode doWork(int side, EntityType type,
169
170private:
171 boost::shared_ptr<VectorDouble> detPtr;
172 boost::shared_ptr<MatrixDouble> jacPtr;
173
176};
177
178/** \brief Apply covariant (Piola) transfer to Hcurl space for HO geometry
179 */
182
184 boost::shared_ptr<MatrixDouble> jac_inv_ptr)
186 jacInvPtr(jac_inv_ptr) {
187 doVertices = false;
188 if (space == HDIV)
189 doEdges = false;
190 if (!jacInvPtr)
192 "Pointer for jacPtr not allocated");
193 }
194
195 MoFEMErrorCode doWork(int side, EntityType type,
197
198private:
199 boost::shared_ptr<MatrixDouble> jacInvPtr;
200
203};
204
205/** \brief Calculate jacobian for face element
206
207 It is assumed that face element is XY plane. Applied
208 only for 2d problems and 2d problems embedded in 3d space
209
210 \note If you push operators for HO normal befor this operator, HO geometry is
211 taken into account when you calculate jacobian.
212
213*/
214template <int DIM> struct OpCalculateHOJacForFaceImpl;
215
216template <>
219
220 OpCalculateHOJacForFaceImpl(boost::shared_ptr<MatrixDouble> jac_ptr);
221
222 MoFEMErrorCode doWork(int side, EntityType type,
224
225protected:
226 boost::shared_ptr<MatrixDouble> jacPtr;
227};
228
229template <>
231
232 using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
233
234 MoFEMErrorCode doWork(int side, EntityType type,
236};
237
240
241template <int DIM> struct OpCalculateHOJac;
242
243template <> struct OpCalculateHOJac<3> : public OpCalculateHOJacForVolume {
245};
246
247template <> struct OpCalculateHOJac<2> : public OpCalculateHOJacForFaceImpl<2> {
248 using OpCalculateHOJacForFaceImpl<2>::OpCalculateHOJacForFaceImpl;
249};
250
251/** \brief Calculate normals at Gauss points of triangle element
252 * \ingroup mofem_forces_and_source
253 */
254template <int FIELD_DIM = 3>
257
260
261 MoFEMErrorCode doWork(int side, EntityType type,
263};
264
265/** \brief transform Hdiv base fluxes from reference element to physical
266 * triangle \ingroup mofem_forces_and_sources
267 *
268 * \note Matrix which keeps normal is assumed to have three columns, and number
269 * of rows should be equal to number of integration points.
270 *
271 */
274
276 const FieldSpace space,
277 boost::shared_ptr<MatrixDouble> normals_at_gauss_pts = nullptr)
279 normalsAtGaussPts(normals_at_gauss_pts) {}
280
281 MoFEMErrorCode doWork(int side, EntityType type,
283
284private:
285 boost::shared_ptr<MatrixDouble> normalsAtGaussPts;
286};
287
288/** \brief transform Hcurl base fluxes from reference element to physical edge
289 * \ingroup mofem_forces_and_sources
290 */
293
295 const FieldSpace space = HCURL,
296 boost::shared_ptr<MatrixDouble> tangent_at_pts = nullptr)
298 tangentAtGaussPts(tangent_at_pts) {}
299
300 MoFEMErrorCode doWork(int side, EntityType type,
302
303private:
304 boost::shared_ptr<MatrixDouble> tangentAtGaussPts;
305};
306
307/** \brief transform Hcurl base fluxes from reference element to physical
308 * triangle \ingroup mofem_forces_and_sources
309 */
312
314 const FieldSpace space,
315 boost::shared_ptr<MatrixDouble> normals_at_pts = nullptr,
316 boost::shared_ptr<MatrixDouble> tangent1_at_pts = nullptr,
317 boost::shared_ptr<MatrixDouble> tangent2_at_pts = nullptr)
319 normalsAtPts(normals_at_pts), tangent1AtPts(tangent1_at_pts),
320 tangent2AtPts(tangent2_at_pts) {
321 if (normals_at_pts || tangent1_at_pts || tangent2_at_pts)
322 if (normals_at_pts && tangent1_at_pts && tangent2_at_pts)
324 "All elements in constructor have to set pointer");
325 }
326
327 MoFEMErrorCode doWork(int side, EntityType type,
329
330private:
331 boost::shared_ptr<MatrixDouble> normalsAtPts;
332 boost::shared_ptr<MatrixDouble> tangent1AtPts;
333 boost::shared_ptr<MatrixDouble> tangent2AtPts;
334
337};
338
339/** \brief Calculate tangent vector on edge form HO geometry approximation
340 * \ingroup mofem_forces_and_sources
341 */
342template <int FIELD_DIM = 3>
345
347 std::string field_name,
348 boost::shared_ptr<MatrixDouble> tangents_at_pts = nullptr)
350 tangentsAtPts(tangents_at_pts) {}
351
352 MoFEMErrorCode doWork(int side, EntityType type,
354
355private:
356 boost::shared_ptr<MatrixDouble> tangentsAtPts;
357};
358
359/**
360 * @brief Scale base functions by inverses of measure of element
361 *
362 * @tparam OP
363 */
364template <typename OP> struct OpScaleBaseBySpaceInverseOfMeasure : public OP {
365
367 boost::shared_ptr<VectorDouble> det_jac_ptr, const FieldSpace space = L2)
368 : OP(space), fieldSpace(space), detJacPtr(det_jac_ptr) {}
369
373
374 if (!detJacPtr) {
375 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
376 }
377
378 auto scale = [&]() {
379 for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
380 FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
381 auto &base_fun = data.getN(base);
382 auto &diff_base_fun = data.getDiffN(base);
383 if (detJacPtr) {
384
385 auto &det_vec = *detJacPtr;
386 const auto nb_int_pts = base_fun.size1();
387
388 if (nb_int_pts != det_vec.size())
389 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
390 "Number of integration pts in detJacPtr does not mush "
391 "number of integration pts in base function");
392
393 auto get_row = [](auto &m, auto gg) {
394 return ublas::matrix_row<decltype(m)>(m, gg);
395 };
396
397 for (auto gg = 0; gg != nb_int_pts; ++gg)
398 get_row(base_fun) /= det_vec[gg];
399
400 if (diff_base_fun.size1() == nb_int_pts) {
401 for (auto gg = 0; gg != nb_int_pts; ++gg)
402 get_row(diff_base_fun) /= det_vec[gg];
403 }
404 }
405 }
406 };
407
408 if (this->getFEDim() == 3) {
409 switch (fieldSpace) {
410 case H1:
411 scale();
412 break;
413 case HCURL:
414 if (type >= MBEDGE)
415 scale();
416 break;
417 case HDIV:
418 if (type >= MBTRI)
419 scale();
420 break;
421 case L2:
422 if (type >= MBTET)
423 scale();
424 break;
425 default:
426 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
427 }
428 } else if (this->getFEDim() == 2) {
429 switch (fieldSpace) {
430 case H1:
431 scale();
432 break;
433 case HCURL:
434 if (type >= MBEDGE)
435 scale();
436 break;
437 case HDIV:
438 case L2:
439 if (type >= MBTRI)
440 scale();
441 break;
442 default:
443 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
444 }
445 } else if (this->getFEDim() == 1) {
446 switch (fieldSpace) {
447 case H1:
448 scale();
449 break;
450 case HCURL:
451 case L2:
452 if (type >= MBEDGE)
453 scale();
454 break;
455 default:
456 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
457 }
458 } else {
459 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
460 }
461
463 }
464
465private:
467 boost::shared_ptr<VectorDouble> detJacPtr;
468};
469
470/**
471 * @brief Add operators pushing bases from local to physical configuration
472 *
473 * @tparam FE_DIM dimension of element
474 * @tparam PROBLEM_DIM problem dimension
475 * @tparam SPACE_DIM space dimension
476 */
477template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
478
479template <> struct AddHOOps<2, 2, 2> {
480 AddHOOps() = delete;
481 static MoFEMErrorCode
482 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
483 std::vector<FieldSpace> spaces, std::string geom_field_name = "",
484 boost::shared_ptr<MatrixDouble> jac = nullptr,
485 boost::shared_ptr<VectorDouble> det = nullptr,
486 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
487};
488
489template <> struct AddHOOps<1, 2, 2> {
490 AddHOOps() = delete;
491 static MoFEMErrorCode
492 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
493 std::vector<FieldSpace> spaces, std::string geom_field_name = "");
494};
495template <> struct AddHOOps<1, 3, 3> {
496 AddHOOps() = delete;
497 static MoFEMErrorCode
498 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
499 std::vector<FieldSpace> space, std::string geom_field_name = "");
500};
501
502template <> struct AddHOOps<2, 3, 3> {
503 AddHOOps() = delete;
504 static MoFEMErrorCode
505 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
506 std::vector<FieldSpace> space, std::string geom_field_name = "");
507};
508
509template <> struct AddHOOps<3, 3, 3> {
510 AddHOOps() = delete;
511 static MoFEMErrorCode
512 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
513 std::vector<FieldSpace> space, std::string geom_field_name = "",
514 boost::shared_ptr<MatrixDouble> jac = nullptr,
515 boost::shared_ptr<VectorDouble> det = nullptr,
516 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
517};
518
519template <int FIELD_DIM>
525 const auto nb_dofs = data.getFieldData().size() / FIELD_DIM;
526 if (nb_dofs) {
527 if (type == MBVERTEX)
528 getCoordsAtGaussPts().clear();
529 auto t_base = data.getFTensor0N();
530 auto t_coords = getFTensor1CoordsAtGaussPts();
531 const auto nb_integration_pts = data.getN().size1();
532 const auto nb_base_functions = data.getN().size2();
533 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
534 auto t_dof = data.getFTensor1FieldData<FIELD_DIM>();
535 size_t bb = 0;
536 for (; bb != nb_dofs; ++bb) {
537 t_coords(i) += t_base * t_dof(i);
538 ++t_dof;
539 ++t_base;
540 }
541 for (; bb < nb_base_functions; ++bb)
542 ++t_base;
543 ++t_coords;
544 }
545 }
547}
548
549template <int FIELD_DIM>
554
557
558 auto get_ftensor1 = [](MatrixDouble &m) {
560 &m(0, 0), &m(0, 1), &m(0, 2));
561 };
562
563 unsigned int nb_dofs = data.getFieldData().size();
564 if (nb_dofs != 0) {
565
566 int nb_gauss_pts = data.getN().size1();
567 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
568 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
569 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
570 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
571
572 switch (type) {
573 case MBVERTEX: {
574 tangent1_at_gauss_pts.clear();
575 tangent2_at_gauss_pts.clear();
576 }
577 case MBEDGE:
578 case MBTRI:
579 case MBQUAD: {
580
581#ifndef NDEBUG
582 if (2 * data.getN().size2() != data.getDiffN().size2()) {
583 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
584 "expected two direcatives in local element coordinates");
585 }
586 if (nb_dofs % FIELD_DIM != 0) {
587 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
588 "expected that number of dofs is multiplicative of field "
589 "dimension");
590 }
591#endif
592
593 if (nb_dofs > FIELD_DIM * data.getN().size2()) {
594 unsigned int nn = 0;
595 for (; nn != nb_dofs; nn++) {
596 if (!data.getFieldDofs()[nn]->getActive())
597 break;
598 }
599 if (nn > FIELD_DIM * data.getN().size2()) {
600 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
601 "Data inconsistency for base %s",
603 } else {
604 nb_dofs = nn;
605 if (!nb_dofs)
607 }
608 }
609 const int nb_base_functions = data.getN().size2();
610 auto t_base = data.getFTensor0N();
611 auto t_diff_base = data.getFTensor1DiffN<2>();
612 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
613 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
614 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
615 auto t_data = data.getFTensor1FieldData<FIELD_DIM>();
616 int bb = 0;
617 for (; bb != nb_dofs / FIELD_DIM; ++bb) {
619 t_t1(i) += t_data(i) * t_diff_base(N0);
620 t_t2(i) += t_data(i) * t_diff_base(N1);
621 ++t_data;
622 ++t_base;
623 ++t_diff_base;
624 }
625 for (; bb != nb_base_functions; ++bb) {
626 ++t_base;
627 ++t_diff_base;
628 }
629 ++t_t1;
630 ++t_t2;
631 }
632 } break;
633 default:
634 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
635 }
636 }
637
638 if (moab::CN::Dimension(type) == 2) {
639
640 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
641 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
642 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
643
644 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
645 normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
646
647 auto t_normal = get_ftensor1(normal_at_gauss_pts);
648 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
649 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
650 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
654 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
655 ++t_normal;
656 ++t_t1;
657 ++t_t2;
658 }
659 }
660
662}
663
664template <int FIELD_DIM>
669
670 auto get_tangent = [&]() -> MatrixDouble & {
671 if (tangentsAtPts)
672 return *tangentsAtPts;
673 else
674 return getTangentAtGaussPts();
675 };
676
677 auto &tangent = get_tangent();
678 int nb_gauss_pts = getGaussPts().size2();
679
680 if (type == MBVERTEX) {
681 tangent.resize(nb_gauss_pts, 3, false);
682 tangent.clear();
683 }
684
685 int nb_dofs = data.getFieldData().size();
686 if (nb_dofs != 0) {
687 const int nb_base_functions = data.getN().size2();
688 double *diff_base_function = &*data.getDiffN().data().begin();
689 auto tangent_at_gauss_pts =
690 getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
691
693 int size = nb_dofs / FIELD_DIM;
694 if (nb_dofs % FIELD_DIM) {
695 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
696 "Inconsistent number of dofs and filed dimension");
697 }
698 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
699 auto field_data = data.getFTensor1FieldData<FIELD_DIM>();
700 int bb = 0;
701 for (; bb < size; ++bb) {
702 tangent_at_gauss_pts(i) += field_data(i) * (*diff_base_function);
703 ++field_data;
704 ++diff_base_function;
705 }
706 for (; bb != nb_base_functions; ++bb) {
707 ++diff_base_function;
708 }
709 ++tangent_at_gauss_pts;
710 }
711 }
713}
714
715/**
716 * @deprecated do not use this function, instead use AddHOOps.
717 *
718 * @tparam E
719 * @param field
720 * @param e
721 * @param h1
722 * @param hcurl
723 * @param hdiv
724 * @param l2
725 * @return MoFEMErrorCode
726 */
727template <typename E>
728MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1,
729 bool hcurl, bool hdiv, bool l2) {
730 std::vector<FieldSpace> spaces;
731 if (h1)
732 spaces.push_back(H1);
733 if (hcurl)
734 spaces.push_back(HCURL);
735 if (hdiv)
736 spaces.push_back(HDIV);
737 if (l2)
738 spaces.push_back(L2);
739 return AddHOOps<3, 3, 3>::add(e.getOpPtrVector(), spaces, field);
740}
741
742/**
743 * @deprecated do not use this function, instead use AddHOOps.
744 *
745 * @tparam E
746 * @param field
747 * @param e
748 * @param hcurl
749 * @param hdiv
750 * @return MoFEMErrorCode
751 */
752template <typename E>
753MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e,
754 bool hcurl, bool hdiv) {
755 std::vector<FieldSpace> spaces;
756 if (hcurl)
757 spaces.push_back(HCURL);
758 if (hdiv)
759 spaces.push_back(HDIV);
760 return AddHOOps<2, 3, 3>::add(e.getOpPtrVector(), spaces, field);
761}
762
763}; // namespace MoFEM
764
765#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:97
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.
Transform local reference derivatives of shape functions to global derivatives.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)