v0.14.0
Loading...
Searching...
No Matches
simple_interface.cpp
Go to the documentation of this file.
1/**
2 * \file simple_interface.cpp
3 * \ingroup mofem_simple_interface
4 * \example simple_interface.cpp
5 *
6 * Calculate volume by integrating volume elements and using divergence theorem
7 * by integrating surface elements.
8 *
9 */
10
11
12
13#include <MoFEM.hpp>
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
18struct OpVolume : public VolumeElementForcesAndSourcesCore::UserDataOperator {
19 Vec vOl;
20 OpVolume(const std::string &field_name, Vec vol)
22 vOl(vol) {}
25
27 if (type != MBVERTEX)
29 const int nb_int_pts = getGaussPts().size2();
30 // cerr << nb_int_pts << endl;
31 auto t_w = getFTensor0IntegrationWeight();
32 double v = getMeasure();
33 double vol = 0;
34 for (int gg = 0; gg != nb_int_pts; gg++) {
35 vol += t_w * v;
36 ++t_w;
37 }
38 CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
40 }
41 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
42 EntityType col_type,
46 // PetscPrintf(PETSC_COMM_WORLD,"domain: calculate matrix\n");
48 }
49};
50
51struct OpFace : public FaceElementForcesAndSourcesCore::UserDataOperator {
52 Vec vOl;
53 OpFace(const std::string &field_name, Vec vol)
55 vOl(vol) {}
58
60 if (type != MBVERTEX)
62 const int nb_int_pts = getGaussPts().size2();
63 auto t_normal = getFTensor1NormalsAtGaussPts();
65 auto t_coords = getFTensor1CoordsAtGaussPts();
66 FTensor::Index<'i', 3> i;
67 double vol = 0;
68 for (int gg = 0; gg != nb_int_pts; gg++) {
69 vol += (t_coords(i) * t_normal(i)) * t_w;
70 ++t_normal;
71 ++t_w;
72 ++t_coords;
73 }
74 vol /= 6;
75 CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
77 }
78 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
79 EntityType col_type,
83 // PetscPrintf(PETSC_COMM_WORLD,"boundary: calculate matrix\n");
85 }
86};
87
88struct VolRule {
89 int operator()(int, int, int) const { return 2; }
90};
91struct FaceRule {
92 int operator()(int, int, int) const { return 4; }
93};
94
95int main(int argc, char *argv[]) {
96
97 // initialize petsc
98 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
99
100 try {
101
102 // Create MoAB database
103 moab::Core moab_core;
104 moab::Interface &moab = moab_core;
105
106 // Create MoFEM database and link it to MoAB
107 MoFEM::Core mofem_core(moab);
108 MoFEM::Interface &m_field = mofem_core;
109
110 // Register DM Manager
111 DMType dm_name = "DMMOFEM";
112 CHKERR DMRegister_MoFEM(dm_name);
113
114 // Simple interface
115 Simple *simple_interface;
116 CHKERR m_field.getInterface(simple_interface);
117 {
118 // get options from command line
119 CHKERR simple_interface->getOptions();
120 // load mesh file
121 CHKERR simple_interface->loadFile();
122 // add fields
123 CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
125 CHKERR simple_interface->addBoundaryField("MESH_NODE_POSITIONS", H1,
127 CHKERR simple_interface->addSkeletonField("MESH_NODE_POSITIONS", H1,
129
130 // set fields order
131 CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 1);
132 // setup problem
133 CHKERR simple_interface->setUp();
134 // Project mesh coordinate on mesh
135 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
136 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
137 // get dm
138 auto dm = simple_interface->getDM();
139 // create elements
140 auto domain_fe =
141 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
142 auto boundary_fe =
143 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
144
145 // set integration rule
146 domain_fe->getRuleHook = VolRule();
147 boundary_fe->getRuleHook = FaceRule();
148 // create distributed vector to accumulate values from processors.
149 int ghosts[] = {0};
150
151 auto vol = createGhostVector(
152 PETSC_COMM_WORLD, m_field.get_comm_rank() == 0 ? 1 : 0, 1,
153 m_field.get_comm_rank() == 0 ? 0 : 1, ghosts);
154 auto surf_vol = vectorDuplicate(vol);
155
156 // set operator to the volume element
157 auto material_grad_mat = boost::make_shared<MatrixDouble>();
158 auto material_det_vec = boost::make_shared<VectorDouble>();
159
160 domain_fe->getOpPtrVector().push_back(
161 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
162 material_grad_mat));
163 domain_fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
164 material_grad_mat, material_det_vec, nullptr));
165 domain_fe->getOpPtrVector().push_back(
166 new OpSetHOWeights(material_det_vec));
167 domain_fe->getOpPtrVector().push_back(
168 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
169 domain_fe->getOpPtrVector().push_back(
170 new OpVolume("MESH_NODE_POSITIONS", vol));
171 // set operator to the face element
172 boundary_fe->getOpPtrVector().push_back(
173 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
174 boundary_fe->getOpPtrVector().push_back(
175 new OpFace("MESH_NODE_POSITIONS", surf_vol));
176 // make integration in volume (here real calculations starts)
177 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
178 domain_fe);
179 // make integration on boundary
180 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
181 boundary_fe);
182
183 auto skeleton_fe = boost::make_shared<FEMethod>();
184 auto A = createDMMatrix(dm);
185 auto B = createDMMatrix(dm);
186 auto f = createDMVector(dm);
187 auto x = createDMVector(dm);
188 auto x_t = createDMVector(dm);
189 auto x_tt = createDMVector(dm);
190 skeleton_fe->f = f;
191 skeleton_fe->A = A;
192 skeleton_fe->B = B;
193 skeleton_fe->x = x;
194 skeleton_fe->x_t = x_t;
195 skeleton_fe->x_tt = x_tt;
196
197 skeleton_fe->preProcessHook = [&]() {
199 if (f != skeleton_fe->ts_F)
200 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
201 if (A != skeleton_fe->ts_A)
202 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
203 if (B != skeleton_fe->ts_B)
204 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
205 if (x != skeleton_fe->ts_u)
206 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
207 if (x_t != skeleton_fe->ts_u_t)
208 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
209 if (x_tt != skeleton_fe->ts_u_tt)
210 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
212 };
213
214 skeleton_fe->postProcessHook = []() { return 0; };
215 skeleton_fe->operatorHook = []() { return 0; };
216
217 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
218 skeleton_fe);
219
220 // assemble volumes from processors and accumulate on processor of rank 0
221 CHKERR VecAssemblyBegin(vol);
222 CHKERR VecAssemblyEnd(vol);
223 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
224 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
225 CHKERR VecAssemblyBegin(surf_vol);
226 CHKERR VecAssemblyEnd(surf_vol);
227 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
228 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
229 if (m_field.get_comm_rank() == 0) {
230 double *a_vol;
231 CHKERR VecGetArray(vol, &a_vol);
232 double *a_surf_vol;
233 CHKERR VecGetArray(surf_vol, &a_surf_vol);
234 cout << "Volume = " << a_vol[0] << endl;
235 cout << "Surf Volume = " << a_surf_vol[0] << endl;
236 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
237 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Should be zero");
238 }
239 CHKERR VecRestoreArray(vol, &a_vol);
240 CHKERR VecRestoreArray(vol, &a_surf_vol);
241 }
242
243 }
244 }
246
247 // finish work cleaning memory, getting statistics, etc.
249
250 return 0;
251}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#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.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:988
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 double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
constexpr AssemblyType A
constexpr auto field_name
static char help[]
Set integration rule to boundary elements.
int operator()(int, int, int) const
virtual int get_comm_rank() const =0
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)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Set inverse jacobian to base functions.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition Simple.hpp:355
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 getDM(DM *dm)
Get DM.
Definition Simple.cpp:667
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 addSkeletonField(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 skeleton.
Definition Simple.cpp:300
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:341
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
OpFace(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpVolume(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Set integration rule.
int operator()(int, int, int) const