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");
144template <
int SPACE_DIM>
166 "Pointer for detPtr not allocated");
184 boost::shared_ptr<VectorDouble> det_ptr,
185 boost::shared_ptr<MatrixDouble> jac_ptr)
210 boost::shared_ptr<MatrixDouble> jac_inv_ptr)
218 "Pointer for jacPtr not allocated");
280template <
int FIELD_DIM = 3>
303 boost::shared_ptr<MatrixDouble> normals_at_gauss_pts =
nullptr)
322 boost::shared_ptr<MatrixDouble> tangent_at_pts =
nullptr)
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)
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");
368template <
int FIELD_DIM = 3>
374 boost::shared_ptr<MatrixDouble> tangents_at_pts =
nullptr)
393 boost::shared_ptr<VectorDouble> det_jac_ptr,
const FieldSpace space =
L2)
407 auto &base_fun = data.
getN(base);
408 auto &diff_base_fun = data.
getDiffN(base);
412 const auto nb_int_pts = base_fun.size1();
414 if (nb_int_pts != det_vec.size())
416 "Number of integration pts in detJacPtr does not mush "
417 "number of integration pts in base function");
419 auto get_row = [](
auto &
m,
auto gg) {
420 return ublas::matrix_row<decltype(m)>(
m, gg);
423 for (
auto gg = 0; gg != nb_int_pts; ++gg)
424 get_row(base_fun) /= det_vec[gg];
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];
434 if (this->getFEDim() == 3) {
454 }
else if (this->getFEDim() == 2) {
471 }
else if (this->getFEDim() == 1) {
503template <
int FE_DIM,
int PROBLEM_DIM,
int SPACE_DIM>
struct AddHOOps;
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);
518 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
519 std::vector<FieldSpace> spaces, std::string geom_field_name =
"");
524 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
525 std::vector<FieldSpace> space, std::string geom_field_name =
"");
531 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
532 std::vector<FieldSpace> space, std::string geom_field_name =
"");
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);
545template <
int FIELD_DIM>
553 if (type == MBVERTEX)
554 getCoordsAtGaussPts().clear();
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) {
562 for (; bb != nb_dofs; ++bb) {
563 t_coords(
i) += t_base * t_dof(
i);
567 for (; bb < nb_base_functions; ++bb)
575template <
int FIELD_DIM>
586 &
m(0, 0), &
m(0, 1), &
m(0, 2));
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);
600 tangent1_at_gauss_pts.clear();
601 tangent2_at_gauss_pts.clear();
608 if (2 * data.
getN().size2() != data.
getDiffN().size2()) {
610 "expected two direcatives in local element coordinates");
614 "expected that number of dofs is multiplicative of field "
621 for (; nn != nb_dofs; nn++) {
627 "Data inconsistency for base %s",
635 const int nb_base_functions = data.
getN().size2();
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) {
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);
651 for (; bb != nb_base_functions; ++bb) {
664 if (moab::CN::Dimension(type) == 2) {
666 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
667 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
668 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
670 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
671 normal_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
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) {
690template <
int FIELD_DIM>
698 return *tangentsAtPts;
700 return getTangentAtGaussPts();
703 auto &tangent = get_tangent();
704 int nb_gauss_pts = getGaussPts().size2();
706 if (type == MBVERTEX) {
707 tangent.resize(nb_gauss_pts, 3,
false);
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());
722 "Inconsistent number of dofs and filed dimension");
724 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
727 for (; bb < size; ++bb) {
728 tangent_at_gauss_pts(
i) += field_data(
i) * (*diff_base_function);
730 ++diff_base_function;
732 for (; bb != nb_base_functions; ++bb) {
733 ++diff_base_function;
735 ++tangent_at_gauss_pts;
755 bool hcurl,
bool hdiv,
bool l2) {
756 std::vector<FieldSpace> spaces;
758 spaces.push_back(
H1);
760 spaces.push_back(
HCURL);
762 spaces.push_back(
HDIV);
764 spaces.push_back(
L2);
780 bool hcurl,
bool hdiv) {
781 std::vector<FieldSpace> spaces;
783 spaces.push_back(
HCURL);
785 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.
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)