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_base_fun = base_fun.size2();
387 const auto nb_int_pts = base_fun.size1();
388
389 if (nb_int_pts != det_vec.size())
390 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
391 "Number of integration pts in detJacPtr does not mush "
392 "number of integration pts in base function");
393
394 auto get_row = [](auto &m, auto gg) {
395 return ublas::matrix_row<decltype(m)>(m, gg);
396 };
397
398 for (auto gg = 0; gg != nb_int_pts; ++gg)
399 get_row(base_fun) /= det_vec[gg];
400
401 if (diff_base_fun.size1() == nb_int_pts) {
402 for (auto gg = 0; gg != nb_int_pts; ++gg)
403 get_row(diff_base_fun) /= det_vec[gg];
404 }
405 }
406 }
407 };
408
409 if (this->getFEDim() == 3) {
410 switch (fieldSpace) {
411 case H1:
412 scale();
413 break;
414 case HCURL:
415 if (type >= MBEDGE)
416 scale();
417 break;
418 case HDIV:
419 if (type >= MBTRI)
420 scale();
421 break;
422 case L2:
423 if (type >= MBTET)
424 scale();
425 break;
426 default:
427 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
428 }
429 } else if (this->getFEDim() == 2) {
430 switch (fieldSpace) {
431 case H1:
432 scale();
433 break;
434 case HCURL:
435 if (type >= MBEDGE)
436 scale();
437 break;
438 case HDIV:
439 case L2:
440 if (type >= MBTRI)
441 scale();
442 break;
443 default:
444 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
445 }
446 } else if (this->getFEDim() == 1) {
447 switch (fieldSpace) {
448 case H1:
449 scale();
450 break;
451 case HCURL:
452 case L2:
453 if (type >= MBEDGE)
454 scale();
455 break;
456 default:
457 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
458 }
459 } else {
460 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "impossible case");
461 }
462
464 }
465
466private:
468 boost::shared_ptr<VectorDouble> detJacPtr;
469};
470
471/**
472 * @brief Add operators pushing bases from local to physical configuration
473 *
474 * @tparam FE_DIM dimension of element
475 * @tparam PROBLEM_DIM problem dimension
476 * @tparam SPACE_DIM space dimension
477 */
478template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
479
480template <> struct AddHOOps<2, 2, 2> {
481 AddHOOps() = delete;
482 static MoFEMErrorCode
483 add(boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
484 std::vector<FieldSpace> spaces, std::string geom_field_name = "",
485 boost::shared_ptr<MatrixDouble> jac = nullptr,
486 boost::shared_ptr<VectorDouble> det = nullptr,
487 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
488};
489
490template <> struct AddHOOps<1, 2, 2> {
491 AddHOOps() = delete;
492 static MoFEMErrorCode
493 add(boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
494 std::vector<FieldSpace> spaces, std::string geom_field_name = "");
495};
496template <> struct AddHOOps<1, 3, 3> {
497 AddHOOps() = delete;
498 static MoFEMErrorCode
499 add(boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
500 std::vector<FieldSpace> space, std::string geom_field_name = "");
501};
502
503template <> struct AddHOOps<2, 3, 3> {
504 AddHOOps() = delete;
505 static MoFEMErrorCode
506 add(boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
507 std::vector<FieldSpace> space, std::string geom_field_name = "");
508};
509
510template <> struct AddHOOps<3, 3, 3> {
511 AddHOOps() = delete;
512 static MoFEMErrorCode
513 add(boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
514 std::vector<FieldSpace> space, std::string geom_field_name = "",
515 boost::shared_ptr<MatrixDouble> jac = nullptr,
516 boost::shared_ptr<VectorDouble> det = nullptr,
517 boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
518};
519
520template <int FIELD_DIM>
526 const auto nb_dofs = data.getFieldData().size() / FIELD_DIM;
527 if (nb_dofs) {
528 if (type == MBVERTEX)
529 getCoordsAtGaussPts().clear();
530 auto t_base = data.getFTensor0N();
531 auto t_coords = getFTensor1CoordsAtGaussPts();
532 const auto nb_integration_pts = data.getN().size1();
533 const auto nb_base_functions = data.getN().size2();
534 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
535 auto t_dof = data.getFTensor1FieldData<FIELD_DIM>();
536 size_t bb = 0;
537 for (; bb != nb_dofs; ++bb) {
538 t_coords(i) += t_base * t_dof(i);
539 ++t_dof;
540 ++t_base;
541 }
542 for (; bb < nb_base_functions; ++bb)
543 ++t_base;
544 ++t_coords;
545 }
546 }
548}
549
550template <int FIELD_DIM>
555
558
559 auto get_ftensor1 = [](MatrixDouble &m) {
561 &m(0, 0), &m(0, 1), &m(0, 2));
562 };
563
564 unsigned int nb_dofs = data.getFieldData().size();
565 if (nb_dofs != 0) {
566
567 int nb_gauss_pts = data.getN().size1();
568 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
569 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
570 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3, false);
571 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3, false);
572
573 switch (type) {
574 case MBVERTEX: {
575 tangent1_at_gauss_pts.clear();
576 tangent2_at_gauss_pts.clear();
577 }
578 case MBEDGE:
579 case MBTRI:
580 case MBQUAD: {
581
582#ifndef NDEBUG
583 if (2 * data.getN().size2() != data.getDiffN().size2()) {
584 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
585 "expected two direcatives in local element coordinates");
586 }
587 if (nb_dofs % FIELD_DIM != 0) {
588 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
589 "expected that number of dofs is multiplicative of field "
590 "dimension");
591 }
592#endif
593
594 if (nb_dofs > FIELD_DIM * data.getN().size2()) {
595 unsigned int nn = 0;
596 for (; nn != nb_dofs; nn++) {
597 if (!data.getFieldDofs()[nn]->getActive())
598 break;
599 }
600 if (nn > FIELD_DIM * data.getN().size2()) {
601 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
602 "Data inconsistency for base %s",
604 } else {
605 nb_dofs = nn;
606 if (!nb_dofs)
608 }
609 }
610 const int nb_base_functions = data.getN().size2();
611 auto t_base = data.getFTensor0N();
612 auto t_diff_base = data.getFTensor1DiffN<2>();
613 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
614 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
615 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
616 auto t_data = data.getFTensor1FieldData<FIELD_DIM>();
617 int bb = 0;
618 for (; bb != nb_dofs / FIELD_DIM; ++bb) {
620 t_t1(i) += t_data(i) * t_diff_base(N0);
621 t_t2(i) += t_data(i) * t_diff_base(N1);
622 ++t_data;
623 ++t_base;
624 ++t_diff_base;
625 }
626 for (; bb != nb_base_functions; ++bb) {
627 ++t_base;
628 ++t_diff_base;
629 }
630 ++t_t1;
631 ++t_t2;
632 }
633 } break;
634 default:
635 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
636 }
637 }
638
639 if (moab::CN::Dimension(type) == 2) {
640
641 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
642 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
643 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
644
645 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
646 normal_at_gauss_pts.resize(nb_gauss_pts, 3, false);
647
648 auto t_normal = get_ftensor1(normal_at_gauss_pts);
649 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
650 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
651 for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
655 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
656 ++t_normal;
657 ++t_t1;
658 ++t_t2;
659 }
660 }
661
663}
664
665template <int FIELD_DIM>
670
671 auto get_tangent = [&]() -> MatrixDouble & {
672 if (tangentsAtPts)
673 return *tangentsAtPts;
674 else
675 return getTangentAtGaussPts();
676 };
677
678 auto &tangent = get_tangent();
679 int nb_gauss_pts = getGaussPts().size2();
680
681 if (type == MBVERTEX) {
682 tangent.resize(nb_gauss_pts, 3, false);
683 tangent.clear();
684 }
685
686 int nb_dofs = data.getFieldData().size();
687 if (nb_dofs != 0) {
688 const int nb_base_functions = data.getN().size2();
689 double *diff_base_function = &*data.getDiffN().data().begin();
690 auto tangent_at_gauss_pts =
691 getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
692
694 int size = nb_dofs / FIELD_DIM;
695 if (nb_dofs % FIELD_DIM) {
696 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
697 "Inconsistent number of dofs and filed dimension");
698 }
699 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
700 auto field_data = data.getFTensor1FieldData<FIELD_DIM>();
701 int bb = 0;
702 tangent_at_gauss_pts(i) = 0;
703 for (; bb < size; ++bb) {
704 tangent_at_gauss_pts(i) += field_data(i) * (*diff_base_function);
705 ++field_data;
706 ++diff_base_function;
707 }
708 for (; bb != nb_base_functions; ++bb) {
709 ++diff_base_function;
710 }
711 ++tangent_at_gauss_pts;
712 }
713 }
715}
716
717/**
718 * @deprecated do not use this function, instead use AddHOOps.
719 *
720 * @tparam E
721 * @param field
722 * @param e
723 * @param h1
724 * @param hcurl
725 * @param hdiv
726 * @param l2
727 * @return MoFEMErrorCode
728 */
729template <typename E>
730MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1,
731 bool hcurl, bool hdiv, bool l2) {
732 std::vector<FieldSpace> spaces;
733 if (h1)
734 spaces.push_back(H1);
735 if (hcurl)
736 spaces.push_back(HCURL);
737 if (hdiv)
738 spaces.push_back(HDIV);
739 if (l2)
740 spaces.push_back(L2);
741 return AddHOOps<3, 3, 3>::add(e.getOpPtrVector(), spaces, field);
742}
743
744/**
745 * @deprecated do not use this function, instead use AddHOOps.
746 *
747 * @tparam E
748 * @param field
749 * @param e
750 * @param hcurl
751 * @param hdiv
752 * @return MoFEMErrorCode
753 */
754template <typename E>
755MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e,
756 bool hcurl, bool hdiv) {
757 std::vector<FieldSpace> spaces;
758 if (hcurl)
759 spaces.push_back(HCURL);
760 if (hdiv)
761 spaces.push_back(HDIV);
762 return AddHOOps<2, 3, 3>::add(e.getOpPtrVector(), spaces, field);
763}
764
765}; // namespace MoFEM
766
767#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:94
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)