v0.13.1
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";
double &dIv;
OpVolDivergence(double &div)
"HDIV", UserDataOperator::OPROW),
dIv(div) {}
};
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_IMPOSIBLE_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);
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 const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ 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()
Definition: definitions.h:447
@ 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 ...
Definition: definitions.h:346
@ MOFEM_IMPOSIBLE_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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
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
int main(int argc, char *argv[])
static char help[]
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.
Definition: Exceptions.hpp:56
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
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 geometr.
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:26
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:430
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
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:392
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
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)