|
| v0.14.0
|
Go to the documentation of this file.
6 #ifndef __HO_DATA_OPERATORS_HPP__
7 #define __HO_DATA_OPERATORS_HPP__
29 boost::shared_ptr<MatrixDouble>
jacPtr;
40 template <
int FIELD_DIM = 3>
72 template <
int DIM = 3>
74 using OpSetHOInvJacToScalarBasesImpl::OpSetHOInvJacToScalarBasesImpl;
94 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
99 "Pointer for invJacPtr not allocated");
144 template <
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");
280 template <
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");
368 template <
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) {
503 template <
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);
555 template <
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)
585 template <
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) {
700 template <
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;
763 template <
typename E>
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);
788 template <
typename E>
790 bool hcurl,
bool hdiv) {
791 std::vector<FieldSpace> spaces;
793 spaces.push_back(
HCURL);
795 spaces.push_back(
HDIV);
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Data on single entity (This is passed as argument to DataOperator::doWork)
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
Calculate tangent vector on edge form HO geometry approximation.
@ OPSPACE
operator do Work is execute on space data
OpSetHOWeights(boost::shared_ptr< VectorDouble > det_ptr)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
const static char *const ApproximationBaseNames[]
@ L2
field with C-1 continuity
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
boost::shared_ptr< MatrixDouble > tangentsAtPts
const VectorDouble & getFieldData() const
get dofs values
OpCalculateHOJacForVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Scale base functions by inverses of measure of element.
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
bool & doVertices
\deprectaed If false skip vertices
Calculate HO coordinates at gauss points.
Calculate jacobian for face element.
@ USER_BASE
user implemented approximation base
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FieldSpace
approximation spaces
boost::shared_ptr< MatrixDouble > jacPtr
implementation of Data Operators for Forces and Sources
Calculate jacobian on Hex or other volume which is not simplex.
default operator for TRI element
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.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FieldApproximationBase & getBase()
Get approximation base.
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.
OpSetHOInvJacVectorBase(const FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Calculate normals at Gauss points of triangle element.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Index< 'i', SPACE_DIM > i
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 > invJacPtr
constexpr auto field_name
Add operators pushing bases from local to physical configuration.
default operator for EDGE element
boost::shared_ptr< VectorDouble > detJacPtr
OpScaleBaseBySpaceInverseOfMeasure(boost::shared_ptr< VectorDouble > det_jac_ptr, const FieldSpace space=L2)
structure to get information form mofem into EntitiesFieldData
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.
bool & doEdges
\deprectaed If false skip edges
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coefficients.
OpGetHOTangentsOnEdge(std::string field_name, boost::shared_ptr< MatrixDouble > tangents_at_pts=nullptr)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble diffHdivInvJac
FTensor::Index< 'j', 3 > j
Set inverse jacobian to base functions.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ HCURL
field with continuous tangents
@ MOFEM_DATA_INCONSISTENCY
OpGetHONormalsOnFace(std::string field_name)
FieldApproximationBase
approximation base
FTensor::Index< 'm', 3 > m
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Modify integration weights on face to take in account higher-order geometry.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
OpCalculateHOCoords(const std::string field_name)
FTensor::Index< 'k', 3 > k
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HDIV
field with continuous normal traction
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.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
@ OPROW
operator doWork function is executed on FE rows