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 =
"",
520 boost::shared_ptr<MatrixDouble> jac =
nullptr,
521 boost::shared_ptr<VectorDouble> det =
nullptr,
522 boost::shared_ptr<MatrixDouble> inv_jac =
nullptr);
528 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
529 std::vector<FieldSpace> spaces, std::string geom_field_name =
"");
534 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
535 std::vector<FieldSpace> space, std::string geom_field_name =
"");
541 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
542 std::vector<FieldSpace> space, std::string geom_field_name =
"");
548 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
549 std::vector<FieldSpace> space, std::string geom_field_name =
"",
550 boost::shared_ptr<MatrixDouble> jac =
nullptr,
551 boost::shared_ptr<VectorDouble> det =
nullptr,
552 boost::shared_ptr<MatrixDouble> inv_jac =
nullptr);
555template <
int FIELD_DIM>
563 if (type == MBVERTEX)
564 getCoordsAtGaussPts().clear();
566 auto t_coords = getFTensor1CoordsAtGaussPts();
567 const auto nb_integration_pts = data.
getN().size1();
568 const auto nb_base_functions = data.
getN().size2();
569 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
572 for (; bb != nb_dofs; ++bb) {
573 t_coords(
i) += t_base * t_dof(
i);
577 for (; bb < nb_base_functions; ++bb)
585template <
int FIELD_DIM>
596 &
m(0, 0), &
m(0, 1), &
m(0, 2));
602 int nb_gauss_pts = data.
getN().size1();
603 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
604 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
605 tangent1_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
606 tangent2_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
610 tangent1_at_gauss_pts.clear();
611 tangent2_at_gauss_pts.clear();
618 if (2 * data.
getN().size2() != data.
getDiffN().size2()) {
620 "expected two direcatives in local element coordinates");
624 "expected that number of dofs is multiplicative of field "
631 for (; nn != nb_dofs; nn++) {
637 "Data inconsistency for base %s",
645 const int nb_base_functions = data.
getN().size2();
648 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
649 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
650 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
653 for (; bb != nb_dofs /
FIELD_DIM; ++bb) {
655 t_t1(
i) += t_data(
i) * t_diff_base(
N0);
656 t_t2(
i) += t_data(
i) * t_diff_base(
N1);
661 for (; bb != nb_base_functions; ++bb) {
674 if (moab::CN::Dimension(type) == 2) {
676 auto &normal_at_gauss_pts = getNormalsAtGaussPts();
677 auto &tangent1_at_gauss_pts = getTangent1AtGaussPts();
678 auto &tangent2_at_gauss_pts = getTangent2AtGaussPts();
680 const auto nb_gauss_pts = tangent1_at_gauss_pts.size1();
681 normal_at_gauss_pts.resize(nb_gauss_pts, 3,
false);
683 auto t_normal = get_ftensor1(normal_at_gauss_pts);
684 auto t_t1 = get_ftensor1(tangent1_at_gauss_pts);
685 auto t_t2 = get_ftensor1(tangent2_at_gauss_pts);
686 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
700template <
int FIELD_DIM>
708 return *tangentsAtPts;
710 return getTangentAtGaussPts();
713 auto &tangent = get_tangent();
714 int nb_gauss_pts = getGaussPts().size2();
716 if (type == MBVERTEX) {
717 tangent.resize(nb_gauss_pts, 3,
false);
723 const int nb_base_functions = data.
getN().size2();
724 double *diff_base_function = &*data.
getDiffN().data().begin();
725 auto tangent_at_gauss_pts =
726 getFTensor1FromPtr<FIELD_DIM, 3>(&*tangent.data().begin());
732 "Inconsistent number of dofs and filed dimension");
734 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
737 for (; bb < size; ++bb) {
738 tangent_at_gauss_pts(
i) += field_data(
i) * (*diff_base_function);
740 ++diff_base_function;
742 for (; bb != nb_base_functions; ++bb) {
743 ++diff_base_function;
745 ++tangent_at_gauss_pts;
765 bool hcurl,
bool hdiv,
bool l2) {
766 std::vector<FieldSpace> spaces;
768 spaces.push_back(
H1);
770 spaces.push_back(
HCURL);
772 spaces.push_back(
HDIV);
774 spaces.push_back(
L2);
790 bool hcurl,
bool hdiv) {
791 std::vector<FieldSpace> spaces;
793 spaces.push_back(
HCURL);
795 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 > >::typ 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)