6#ifndef __HO_DATA_OPERATORS_HPP__
7#define __HO_DATA_OPERATORS_HPP__
29 boost::shared_ptr<MatrixDouble>
jacPtr;
40template <
int FIELD_DIM = 3>
74 using OpSetHOInvJacToScalarBasesImpl::OpSetHOInvJacToScalarBasesImpl;
94 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
99 "Pointer for invJacPtr not allocated");
140 "Pointer for detPtr not allocated");
158 boost::shared_ptr<VectorDouble> det_ptr,
159 boost::shared_ptr<MatrixDouble> jac_ptr)
184 boost::shared_ptr<MatrixDouble> jac_inv_ptr)
192 "Pointer for jacPtr not allocated");
254template <
int FIELD_DIM = 3>
277 boost::shared_ptr<MatrixDouble> normals_at_gauss_pts =
nullptr)
296 boost::shared_ptr<MatrixDouble> tangent_at_pts =
nullptr)
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)
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");
342template <
int FIELD_DIM = 3>
348 boost::shared_ptr<MatrixDouble> tangents_at_pts =
nullptr)
367 boost::shared_ptr<VectorDouble> det_jac_ptr,
const FieldSpace space =
L2)
381 auto &base_fun = data.
getN(base);
382 auto &diff_base_fun = data.
getDiffN(base);
386 const auto nb_int_pts = base_fun.size1();
388 if (nb_int_pts != det_vec.size())
390 "Number of integration pts in detJacPtr does not mush "
391 "number of integration pts in base function");
393 auto get_row = [](
auto &
m,
auto gg) {
394 return ublas::matrix_row<decltype(m)>(
m, gg);
397 for (
auto gg = 0; gg != nb_int_pts; ++gg)
398 get_row(base_fun) /= det_vec[gg];
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];
408 if (this->getFEDim() == 3) {
428 }
else if (this->getFEDim() == 2) {
445 }
else if (this->getFEDim() == 1) {
477template <
int FE_DIM,
int PROBLEM_DIM,
int SPACE_DIM>
struct AddHOOps;
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);
492 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
493 std::vector<FieldSpace> spaces, std::string geom_field_name =
"");
498 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
499 std::vector<FieldSpace> space, std::string geom_field_name =
"");
505 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
506 std::vector<FieldSpace> space, std::string geom_field_name =
"");
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);
519template <
int FIELD_DIM>
527 if (type == MBVERTEX)
528 getCoordsAtGaussPts().clear();
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) {
536 for (; bb != nb_dofs; ++bb) {
537 t_coords(
i) += t_base * t_dof(
i);
541 for (; bb < nb_base_functions; ++bb)
549template <
int FIELD_DIM>
560 &
m(0, 0), &
m(0, 1), &
m(0, 2));
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);
574 tangent1_at_gauss_pts.clear();
575 tangent2_at_gauss_pts.clear();
582 if (2 * data.
getN().size2() != data.
getDiffN().size2()) {
584 "expected two direcatives in local element coordinates");
588 "expected that number of dofs is multiplicative of field "
595 for (; nn != nb_dofs; nn++) {
601 "Data inconsistency for base %s",
609 const int nb_base_functions = data.
getN().size2();
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) {
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);
625 for (; bb != nb_base_functions; ++bb) {
638 if (moab::CN::Dimension(type) == 2) {
640 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
641 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
642 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
644 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
645 normal_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
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) {
664template <
int FIELD_DIM>
672 return *tangentsAtPts;
674 return getTangentAtGaussPts();
677 auto &tangent = get_tangent();
678 int nb_gauss_pts = getGaussPts().size2();
680 if (type == MBVERTEX) {
681 tangent.resize(nb_gauss_pts, 3,
false);
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());
696 "Inconsistent number of dofs and filed dimension");
698 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
701 for (; bb < size; ++bb) {
702 tangent_at_gauss_pts(
i) += field_data(
i) * (*diff_base_function);
704 ++diff_base_function;
706 for (; bb != nb_base_functions; ++bb) {
707 ++diff_base_function;
709 ++tangent_at_gauss_pts;
729 bool hcurl,
bool hdiv,
bool l2) {
730 std::vector<FieldSpace> spaces;
732 spaces.push_back(
H1);
734 spaces.push_back(
HCURL);
736 spaces.push_back(
HDIV);
738 spaces.push_back(
L2);
754 bool hcurl,
bool hdiv) {
755 std::vector<FieldSpace> spaces;
757 spaces.push_back(
HCURL);
759 spaces.push_back(
HDIV);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#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
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#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 > >::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.
implementation of Data Operators for Forces and Sources
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)
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
default operator for EDGE element
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
default operator for TRI element
@ 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)
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
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)
MatrixDouble diffHdivInvJac
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)