16static char help[] =
"...\n\n";
44int main(
int argc,
char *argv[]) {
50 enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
52 const char *list[] = {
"ainsworth",
"demkowicz"};
55 PetscInt choice_value = AINSWORTH;
58 if (flg != PETSC_TRUE) {
62 PetscBool ho_geometry = PETSC_FALSE;
66 PetscInt ho_choice_value = AINSWORTH;
68 &ho_choice_value, &flg);
70 DMType dm_name =
"DMMOFEM";
74 moab::Core mb_instance;
75 moab::Interface &moab = mb_instance;
87 switch (choice_value) {
102 if (ho_geometry == PETSC_TRUE) {
103 switch (ho_choice_value) {
113 constexpr int order = 3;
115 if (ho_geometry == PETSC_TRUE)
120 pipeline_mng->getDomainLhsFE().reset();
121 pipeline_mng->getDomainRhsFE().reset();
122 pipeline_mng->getBoundaryLhsFE().reset();
123 pipeline_mng->getBoundaryRhsFE().reset();
124 pipeline_mng->getSkeletonLhsFE().reset();
125 pipeline_mng->getSkeletonRhsFE().reset();
131 double divergence_vol = 0;
132 double divergence_skin = 0;
134 auto jac_ptr = boost::make_shared<MatrixDouble>();
135 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
136 auto det_ptr = boost::make_shared<VectorDouble>();
137 pipeline_mng->getOpDomainRhsPipeline().push_back(
139 pipeline_mng->getOpDomainRhsPipeline().push_back(
141 pipeline_mng->getOpDomainRhsPipeline().push_back(
143 pipeline_mng->getOpDomainRhsPipeline().push_back(
145 pipeline_mng->getOpDomainRhsPipeline().push_back(
147 pipeline_mng->getOpDomainRhsPipeline().push_back(
151 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
153 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
155 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
159 if (ho_geometry == PETSC_TRUE) {
163 CHKERR pipeline_mng->loopFiniteElements();
166 <<
"divergence_vol " << std::setprecision(12) << divergence_vol;
168 <<
"divergence_skin " << std::setprecision(12) << divergence_skin;
170 constexpr double eps = 1e-8;
171 const double error = divergence_skin - divergence_vol;
172 if (fabs(error) >
eps)
174 "invalid surface flux or divergence or both error = %3.4e",
187 if (CN::Dimension(type) < 2)
193 int nb_gauss_pts = data.
getDiffN().size1();
197 div_vec.resize(nb_dofs, 0);
202 for (
size_t gg = 0; gg < nb_gauss_pts; gg++) {
203 for (
size_t dd = 0; dd != div_vec.size(); dd++) {
204 double w = getGaussPts()(3, gg) * getVolume();
205 dIv += t_base_diff_hdiv(
i,
i) * w;
217 if (CN::Dimension(type) != 2)
220 int nb_gauss_pts = data.
getN().size1();
223 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
224 for (
int dd = 0; dd < nb_dofs; dd++) {
226 double w = getGaussPts()(2, gg);
227 const double n0 = getNormalsAtGaussPts(gg)[0];
228 const double n1 = getNormalsAtGaussPts(gg)[1];
229 const double n2 = getNormalsAtGaussPts(gg)[2];
230 if (getFEType() == MBTRI) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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_ATOM_TEST_INVALID
#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 ...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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)
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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
default operator for TRI element
Calculate normals at Gauss points of triangle element.
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.
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.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
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.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpFacesFluxes(double &div)
OpVolDivergence(double &div)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)