v0.13.2
Loading...
Searching...
No Matches
hdiv_divergence_operator.cpp
Go to the documentation of this file.
1/**
2 * \file hdiv_divergence_operator.cpp
3 * \example hdiv_divergence_operator.cpp
4 *
5 * Using PipelineManager interface calculate the divergence of base functions,
6 * and integral of flux on the boundary. Since the h-div space is used, volume
7 * integral and boundary integral should give the same result.
8 */
9
10
11
12#include <MoFEM.hpp>
13
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
20
21 double &dIv;
22 OpVolDivergence(double &div)
24 "HDIV", UserDataOperator::OPROW),
25 dIv(div) {}
26
27 MoFEMErrorCode doWork(int side, EntityType type,
29};
30
33
34 double &dIv;
35 OpFacesFluxes(double &div)
37 "HDIV", UserDataOperator::OPROW),
38 dIv(div) {}
39
40 MoFEMErrorCode doWork(int side, EntityType type,
42};
43
44int main(int argc, char *argv[]) {
45
46 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
47
48 try {
49
50 enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
51
52 const char *list[] = {"ainsworth", "demkowicz"};
53
54 PetscBool flg;
55 PetscInt choice_value = AINSWORTH;
56 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
57 &choice_value, &flg);
58 if (flg != PETSC_TRUE) {
59 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
60 }
61
62 PetscBool ho_geometry = PETSC_FALSE;
63 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ho_geometry", &ho_geometry,
64 PETSC_NULL);
65
66 PetscInt ho_choice_value = AINSWORTH;
67 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-ho_base", list, LASTOP,
68 &ho_choice_value, &flg);
69
70 DMType dm_name = "DMMOFEM";
71 CHKERR DMRegister_MoFEM(dm_name);
72
73 // Create MoAB
74 moab::Core mb_instance; ///< database
75 moab::Interface &moab = mb_instance; ///< interface
76
77 // Create MoFEM
78 MoFEM::Core core(moab);
79 MoFEM::Interface &m_field = core;
80
81 Simple *simple_interface = m_field.getInterface<Simple>();
82 PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
83 CHKERR simple_interface->getOptions();
84 CHKERR simple_interface->loadFile("");
85
86 // fields
87 switch (choice_value) {
88 case AINSWORTH:
89 CHKERR simple_interface->addDomainField("HDIV", HDIV,
91 CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
93 break;
94 case DEMKOWICZ:
95 CHKERR simple_interface->addDomainField("HDIV", HDIV,
97 CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
99 break;
100 }
101
102 if (ho_geometry == PETSC_TRUE) {
103 switch (ho_choice_value) {
104 case AINSWORTH:
105 CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
107 case DEMKOWICZ:
108 CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
110 }
111 }
112
113 constexpr int order = 3;
114 CHKERR simple_interface->setFieldOrder("HDIV", order);
115 if (ho_geometry == PETSC_TRUE)
116 CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 2);
117 CHKERR simple_interface->setUp();
118
119 /// This has no real effect, folling line are only for atom test purpose
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();
126
127 auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
128 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
129 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
130
131 double divergence_vol = 0;
132 double divergence_skin = 0;
133
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(
138 new OpCalculateHOJac<3>(jac_ptr));
139 pipeline_mng->getOpDomainRhsPipeline().push_back(
140 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
141 pipeline_mng->getOpDomainRhsPipeline().push_back(
142 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
143 pipeline_mng->getOpDomainRhsPipeline().push_back(
144 new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
145 pipeline_mng->getOpDomainRhsPipeline().push_back(
146 new OpSetHOWeights(det_ptr));
147 pipeline_mng->getOpDomainRhsPipeline().push_back(
148 new OpVolDivergence(divergence_vol));
149
150 if (m_field.check_field("MESH_NODE_POSITIONS"))
151 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
152 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
153 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
155 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
156 new OpFacesFluxes(divergence_skin));
157
158 // project geometry form 10 node tets on higher order approx. functions
159 if (ho_geometry == PETSC_TRUE) {
160 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
161 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
162 }
163 CHKERR pipeline_mng->loopFiniteElements();
164
165 MOFEM_LOG("WORLD", Sev::inform)
166 << "divergence_vol " << std::setprecision(12) << divergence_vol;
167 MOFEM_LOG("WORLD", Sev::inform)
168 << "divergence_skin " << std::setprecision(12) << divergence_skin;
169
170 constexpr double eps = 1e-8;
171 const double error = divergence_skin - divergence_vol;
172 if (fabs(error) > eps)
173 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
174 "invalid surface flux or divergence or both error = %3.4e",
175 error);
176 }
178
180}
181
186
187 if (CN::Dimension(type) < 2)
189
190 if (data.getFieldData().size() == 0)
192
193 int nb_gauss_pts = data.getDiffN().size1();
194 int nb_dofs = data.getFieldData().size();
195
196 VectorDouble div_vec;
197 div_vec.resize(nb_dofs, 0);
198
200 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
201
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;
206 ++t_base_diff_hdiv;
207 }
208 }
209
211}
212
216
217 if (CN::Dimension(type) != 2)
219
220 int nb_gauss_pts = data.getN().size1();
221 int nb_dofs = data.getFieldData().size();
222
223 for (int gg = 0; gg < nb_gauss_pts; gg++) {
224 for (int dd = 0; dd < nb_dofs; dd++) {
225
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) {
231 w *= 0.5;
232 }
233
234 dIv += (n0 * data.getVectorN<3>(gg)(dd, 0) +
235 n1 * data.getVectorN<3>(gg)(dd, 1) +
236 n2 * data.getVectorN<3>(gg)(dd, 2)) *
237 w;
238 }
239 }
240
242}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
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_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()
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: 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.
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
static char help[]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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)
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 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 loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
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)