v0.13.1
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/* This file is part of MoFEM.
12 * MoFEM is free software: you can redistribute it and/or modify it under
13 * the terms of the GNU Lesser General Public License as published by the
14 * Free Software Foundation, either version 3 of the License, or (at your
15 * option) any later version.
16 *
17 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20 * License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24
25#include <MoFEM.hpp>
26using namespace MoFEM;
27
28static char help[] = "...\n\n";
29
32 OpVolume(const std::string &field_name, Vec vol)
34 vOl(vol) {}
37
39 if (type != MBVERTEX)
41 const int nb_int_pts = getGaussPts().size2();
42 // cerr << nb_int_pts << endl;
43 auto t_w = getFTensor0IntegrationWeight();
44 double v = getMeasure();
45 double vol = 0;
46 for (int gg = 0; gg != nb_int_pts; gg++) {
47 vol += t_w * v;
48 ++t_w;
49 }
50 CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
52 }
53 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
54 EntityType col_type,
58 // PetscPrintf(PETSC_COMM_WORLD,"domain: calculate matrix\n");
60 }
61};
62
65 OpFace(const std::string &field_name, Vec vol)
67 vOl(vol) {}
70
72 if (type != MBVERTEX)
74 const int nb_int_pts = getGaussPts().size2();
75 auto t_normal = getFTensor1NormalsAtGaussPts();
77 auto t_coords = getFTensor1CoordsAtGaussPts();
79 double vol = 0;
80 for (int gg = 0; gg != nb_int_pts; gg++) {
81 vol += (t_coords(i) * t_normal(i)) * t_w;
82 ++t_normal;
83 ++t_w;
84 ++t_coords;
85 }
86 vol /= 6;
87 CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
89 }
90 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
91 EntityType col_type,
95 // PetscPrintf(PETSC_COMM_WORLD,"boundary: calculate matrix\n");
97 }
98};
99
100struct VolRule {
101 int operator()(int, int, int) const { return 2; }
102};
103struct FaceRule {
104 int operator()(int, int, int) const { return 4; }
105};
106
107int main(int argc, char *argv[]) {
108
109 // initialize petsc
110 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
111
112 try {
113
114 // Create MoAB database
115 moab::Core moab_core;
116 moab::Interface &moab = moab_core;
117
118 // Create MoFEM database and link it to MoAB
119 MoFEM::Core mofem_core(moab);
120 MoFEM::Interface &m_field = mofem_core;
121
122 // Register DM Manager
123 DMType dm_name = "DMMOFEM";
124 CHKERR DMRegister_MoFEM(dm_name);
125
126 // Simple interface
127 Simple *simple_interface;
128 CHKERR m_field.getInterface(simple_interface);
129 {
130 // get options from command line
131 CHKERR simple_interface->getOptions();
132 // load mesh file
133 CHKERR simple_interface->loadFile();
134 // add fields
135 CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
137 CHKERR simple_interface->addBoundaryField("MESH_NODE_POSITIONS", H1,
139 CHKERR simple_interface->addSkeletonField("MESH_NODE_POSITIONS", H1,
141
142 // set fields order
143 CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 1);
144 // setup problem
145 CHKERR simple_interface->setUp();
146 // Project mesh coordinate on mesh
147 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
148 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
149 // get dm
150 auto dm = simple_interface->getDM();
151 // create elements
152 auto domain_fe =
153 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
154 auto boundary_fe =
155 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
156
157 // set integration rule
158 domain_fe->getRuleHook = VolRule();
159 boundary_fe->getRuleHook = FaceRule();
160 // create distributed vector to accumulate values from processors.
161 int ghosts[] = {0};
162
163 auto vol = createSmartGhostVector(
164 PETSC_COMM_WORLD, m_field.get_comm_rank() == 0 ? 1 : 0, 1,
165 m_field.get_comm_rank() == 0 ? 0 : 1, ghosts);
166 auto surf_vol = smartVectorDuplicate(vol);
167
168 // set operator to the volume element
169 auto material_grad_mat = boost::make_shared<MatrixDouble>();
170 auto material_det_vec = boost::make_shared<VectorDouble>();
171
172 domain_fe->meshPositionsFieldName = "none";
173 domain_fe->getOpPtrVector().push_back(
174 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
175 material_grad_mat));
176 domain_fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
177 material_grad_mat, material_det_vec, nullptr));
178 domain_fe->getOpPtrVector().push_back(
179 new OpSetHOWeights(material_det_vec));
180 domain_fe->getOpPtrVector().push_back(
181 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
182 domain_fe->getOpPtrVector().push_back(
183 new OpVolume("MESH_NODE_POSITIONS", vol));
184 // set operator to the face element
185 boundary_fe->getOpPtrVector().push_back(
186 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
187 boundary_fe->getOpPtrVector().push_back(
188 new OpFace("MESH_NODE_POSITIONS", surf_vol));
189 // make integration in volume (here real calculations starts)
190 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
191 domain_fe);
192 // make integration on boundary
193 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
194 boundary_fe);
195
196 auto skeleton_fe = boost::make_shared<FEMethod>();
197 auto A = smartCreateDMMatrix(dm);
198 auto B = smartCreateDMMatrix(dm);
199 auto f = smartCreateDMVector(dm);
200 auto x = smartCreateDMVector(dm);
201 auto x_t = smartCreateDMVector(dm);
202 auto x_tt = smartCreateDMVector(dm);
203 skeleton_fe->f = f;
204 skeleton_fe->A = A;
205 skeleton_fe->B = B;
206 skeleton_fe->x = x;
207 skeleton_fe->x_t = x_t;
208 skeleton_fe->x_tt = x_tt;
209
210 skeleton_fe->preProcessHook = [&]() {
212 if (f != skeleton_fe->ts_F)
213 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
214 if (A != skeleton_fe->ts_A)
215 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
216 if (B != skeleton_fe->ts_B)
217 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
218 if (x != skeleton_fe->ts_u)
219 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
220 if (x_t != skeleton_fe->ts_u_t)
221 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
222 if (x_tt != skeleton_fe->ts_u_tt)
223 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong ptr");
225 };
226
227 skeleton_fe->postProcessHook = []() { return 0; };
228 skeleton_fe->operatorHook = []() { return 0; };
229
230 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
231 skeleton_fe);
232
233 // assemble volumes from processors and accumulate on processor of rank 0
234 CHKERR VecAssemblyBegin(vol);
235 CHKERR VecAssemblyEnd(vol);
236 CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
237 CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
238 CHKERR VecAssemblyBegin(surf_vol);
239 CHKERR VecAssemblyEnd(surf_vol);
240 CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
241 CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
242 if (m_field.get_comm_rank() == 0) {
243 double *a_vol;
244 CHKERR VecGetArray(vol, &a_vol);
245 double *a_surf_vol;
246 CHKERR VecGetArray(surf_vol, &a_surf_vol);
247 cout << "Volume = " << a_vol[0] << endl;
248 cout << "Surf Volume = " << a_surf_vol[0] << endl;
249 if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
250 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Should be zero");
251 }
252 CHKERR VecRestoreArray(vol, &a_vol);
253 CHKERR VecRestoreArray(vol, &a_surf_vol);
254 }
255
256 }
257 }
259
260 // finish work cleaning memory, getting statistics, etc.
262
263 return 0;
264}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:965
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
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
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double A
constexpr auto field_name
int main(int argc, char *argv[])
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:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
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:35
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:328
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:335
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:378
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:294
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:810
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:308
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:593
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:396
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:414
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:753
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:321
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