v0.14.0
Loading...
Searching...
No Matches
hdiv_divergence_operator.cpp

Using PipelineManager interface calculate the divergence of base functions, and integral of flux on the boundary. Since the h-div space is used, volume integral and boundary integral should give the same result.

/**
* \file hdiv_divergence_operator.cpp
* \example hdiv_divergence_operator.cpp
*
* Using PipelineManager interface calculate the divergence of base functions,
* and integral of flux on the boundary. Since the h-div space is used, volume
* integral and boundary integral should give the same result.
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
double &dIv;
OpVolDivergence(double &div)
: VolumeElementForcesAndSourcesCore::UserDataOperator(
"HDIV", UserDataOperator::OPROW),
dIv(div) {}
};
: public FaceElementForcesAndSourcesCore::UserDataOperator {
double &dIv;
OpFacesFluxes(double &div)
"HDIV", UserDataOperator::OPROW),
dIv(div) {}
};
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
const char *list[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
&choice_value, &flg);
if (flg != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
}
PetscBool ho_geometry = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ho_geometry", &ho_geometry,
PETSC_NULL);
PetscInt ho_choice_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-ho_base", list, LASTOP,
&ho_choice_value, &flg);
DMType dm_name = "DMMOFEM";
// Create MoAB
moab::Core mb_instance; ///< database
moab::Interface &moab = mb_instance; ///< interface
// Create MoFEM
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
Simple *simple_interface = m_field.getInterface<Simple>();
PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile("");
// fields
switch (choice_value) {
case AINSWORTH:
CHKERR simple_interface->addDomainField("HDIV", HDIV,
CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
break;
case DEMKOWICZ:
CHKERR simple_interface->addDomainField("HDIV", HDIV,
CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
break;
}
if (ho_geometry == PETSC_TRUE) {
switch (ho_choice_value) {
case AINSWORTH:
CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
case DEMKOWICZ:
CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
}
}
constexpr int order = 3;
CHKERR simple_interface->setFieldOrder("HDIV", order);
if (ho_geometry == PETSC_TRUE)
CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 2);
CHKERR simple_interface->setUp();
/// This has no real effect, folling line are only for atom test purpose
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getSkeletonRhsFE().reset();
auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
double divergence_vol = 0;
double divergence_skin = 0;
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJac<3>(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOWeights(det_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpVolDivergence(divergence_vol));
if (m_field.check_field("MESH_NODE_POSITIONS"))
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpFacesFluxes(divergence_skin));
// project geometry form 10 node tets on higher order approx. functions
if (ho_geometry == PETSC_TRUE) {
Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
}
CHKERR pipeline_mng->loopFiniteElements();
MOFEM_LOG("WORLD", Sev::inform)
<< "divergence_vol " << std::setprecision(12) << divergence_vol;
MOFEM_LOG("WORLD", Sev::inform)
<< "divergence_skin " << std::setprecision(12) << divergence_skin;
constexpr double eps = 1e-8;
const double error = divergence_skin - divergence_vol;
if (fabs(error) > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"invalid surface flux or divergence or both error = %3.4e",
error);
}
}
if (CN::Dimension(type) < 2)
if (data.getFieldData().size() == 0)
int nb_gauss_pts = data.getDiffN().size1();
int nb_dofs = data.getFieldData().size();
VectorDouble div_vec;
div_vec.resize(nb_dofs, 0);
FTensor::Index<'i', 3> i;
auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
for (size_t gg = 0; gg < nb_gauss_pts; gg++) {
for (size_t dd = 0; dd != div_vec.size(); dd++) {
double w = getGaussPts()(3, gg) * getVolume();
dIv += t_base_diff_hdiv(i, i) * w;
++t_base_diff_hdiv;
}
}
}
if (CN::Dimension(type) != 2)
int nb_gauss_pts = data.getN().size1();
int nb_dofs = data.getFieldData().size();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
for (int dd = 0; dd < nb_dofs; dd++) {
double w = getGaussPts()(2, gg);
const double n0 = getNormalsAtGaussPts(gg)[0];
const double n1 = getNormalsAtGaussPts(gg)[1];
const double n2 = getNormalsAtGaussPts(gg)[2];
if (getFEType() == MBTRI) {
w *= 0.5;
}
dIv += (n0 * data.getVectorN<3>(gg)(dd, 0) +
n1 * data.getVectorN<3>(gg)(dd, 1) +
n2 * data.getVectorN<3>(gg)(dd, 2)) *
w;
}
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static char help[]
int main()
static const double eps
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
@ 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 ...
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'i', SPACE_DIM > i
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
double w(const double g, const double t)
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
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
const MatrixAdaptor getVectorN(const FieldApproximationBase base, const int gg)
get Hdiv of base functions at Gauss pts
Calculate normals at Gauss points of triangle element.
transform Hdiv base fluxes from reference element to physical triangle
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode addDataField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:320
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition Simple.cpp:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:611
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)