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
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
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();
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->meshPositionsFieldName = "none";
161 domain_fe->getOpPtrVector().push_back(
162 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
163 material_grad_mat));
164 domain_fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
165 material_grad_mat, material_det_vec, nullptr));
166 domain_fe->getOpPtrVector().push_back(
167 new OpSetHOWeights(material_det_vec));
168 domain_fe->getOpPtrVector().push_back(
169 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
170 domain_fe->getOpPtrVector().push_back(
171 new OpVolume("MESH_NODE_POSITIONS", vol));
172 // set operator to the face element
173 boundary_fe->getOpPtrVector().push_back(
174 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
175 boundary_fe->getOpPtrVector().push_back(
176 new OpFace("MESH_NODE_POSITIONS", surf_vol));
177 // make integration in volume (here real calculations starts)
178 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
179 domain_fe);
180 // make integration on boundary
181 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
182 boundary_fe);
183
184 auto skeleton_fe = boost::make_shared<FEMethod>();
185 auto A = createDMMatrix(dm);
186 auto B = createDMMatrix(dm);
187 auto f = createDMVector(dm);
188 auto x = createDMVector(dm);
189 auto x_t = createDMVector(dm);
190 auto x_tt = createDMVector(dm);
191 skeleton_fe->f = f;
192 skeleton_fe->A = A;
193 skeleton_fe->B = B;
194 skeleton_fe->x = x;
195 skeleton_fe->x_t = x_t;
196 skeleton_fe->x_tt = x_tt;
197
198 skeleton_fe->preProcessHook = [&]() {
200 if (f != skeleton_fe->ts_F)
201 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
202 if (A != skeleton_fe->ts_A)
203 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
204 if (B != skeleton_fe->ts_B)
205 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
206 if (x != skeleton_fe->ts_u)
207 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
208 if (x_t != skeleton_fe->ts_u_t)
209 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
210 if (x_tt != skeleton_fe->ts_u_tt)
211 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
213 };
214
215 skeleton_fe->postProcessHook = []() { return 0; };
216 skeleton_fe->operatorHook = []() { return 0; };
217
218 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
219 skeleton_fe);
220
221 // assemble volumes from processors and accumulate on processor of rank 0
222 CHKERR VecAssemblyBegin(vol);
223 CHKERR VecAssemblyEnd(vol);
224 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
225 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
226 CHKERR VecAssemblyBegin(surf_vol);
227 CHKERR VecAssemblyEnd(surf_vol);
228 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
229 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
230 if (m_field.get_comm_rank() == 0) {
231 double *a_vol;
232 CHKERR VecGetArray(vol, &a_vol);
233 double *a_surf_vol;
234 CHKERR VecGetArray(surf_vol, &a_surf_vol);
235 cout << "Volume = " << a_vol[0] << endl;
236 cout << "Surf Volume = " << a_surf_vol[0] << endl;
237 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
238 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Should be zero");
239 }
240 CHKERR VecRestoreArray(vol, &a_vol);
241 CHKERR VecRestoreArray(vol, &a_surf_vol);
242 }
243
244 }
245 }
247
248 // finish work cleaning memory, getting statistics, etc.
250
251 return 0;
252}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ 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()
Definition: definitions.h:447
@ 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 ...
Definition: definitions.h:346
@ 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()
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
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.
Definition: Exceptions.hpp:56
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:334
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:341
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 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: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:327
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