v0.9.0
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>
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
31  Vec vOl;
32  OpVolume(const std::string &field_name, Vec vol)
34  vOl(vol) {}
35  MoFEMErrorCode doWork(int side, EntityType type,
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  auto t_ho_det = getFTenosr0HoMeasure();
45  double v = getMeasure();
46  double vol = 0;
47  for (int gg = 0; gg != nb_int_pts; gg++) {
48  vol += t_w * t_ho_det * v;
49  // cerr << t_ho_det << endl;
50  ++t_w;
51  ++t_ho_det;
52  }
53  CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
55  }
56  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
57  EntityType col_type,
61  // PetscPrintf(PETSC_COMM_WORLD,"domain: calculate matrix\n");
63  }
64 };
65 
67  Vec vOl;
68  OpFace(const std::string &field_name, Vec vol)
70  vOl(vol) {}
71  MoFEMErrorCode doWork(int side, EntityType type,
73 
75  if (type != MBVERTEX)
77  const int nb_int_pts = getGaussPts().size2();
78  auto t_normal = getFTensor1NormalsAtGaussPts();
79  auto t_w = getFTensor0IntegrationWeight();
80  auto t_coords = getFTensor1HoCoordsAtGaussPts();
82  double vol = 0;
83  for (int gg = 0; gg != nb_int_pts; gg++) {
84  vol += (t_coords(i) * t_normal(i)) * t_w;
85  ++t_normal;
86  ++t_w;
87  ++t_coords;
88  }
89  vol /= 6;
90  CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
92  }
93  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
94  EntityType col_type,
98  // PetscPrintf(PETSC_COMM_WORLD,"boundary: calculate matrix\n");
100  }
101 };
102 
103 struct VolRule {
104  int operator()(int, int, int) const { return 2; }
105 };
106 struct FaceRule {
107  int operator()(int, int, int) const { return 4; }
108 };
109 
110 int main(int argc, char *argv[]) {
111 
112  // initialize petsc
113  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
114 
115  try {
116 
117  // Create MoAB database
118  moab::Core moab_core;
119  moab::Interface &moab = moab_core;
120 
121  // Create MoFEM database and link it to MoAB
122  MoFEM::Core mofem_core(moab);
123  MoFEM::Interface &m_field = mofem_core;
124 
125  // Register DM Manager
126  DMType dm_name = "DMMOFEM";
127  CHKERR DMRegister_MoFEM(dm_name);
128 
129  // Simple interface
130  Simple *simple_interface;
131  CHKERR m_field.getInterface(simple_interface);
132  {
133  // get options from command line
134  CHKERR simple_interface->getOptions();
135  // load mesh file
136  CHKERR simple_interface->loadFile();
137  // add fields
138  CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
140  CHKERR simple_interface->addBoundaryField("MESH_NODE_POSITIONS", H1,
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  boost::shared_ptr<ForcesAndSourcesCore> domain_fe =
153  boost::shared_ptr<ForcesAndSourcesCore>(
154  new VolumeElementForcesAndSourcesCore(m_field));
155  boost::shared_ptr<ForcesAndSourcesCore> boundary_fe =
156  boost::shared_ptr<ForcesAndSourcesCore>(
157  new FaceElementForcesAndSourcesCore(m_field));
158  // set integration rule
159  domain_fe->getRuleHook = VolRule();
160  boundary_fe->getRuleHook = FaceRule();
161  // create distributed vector to accumulate values from processors.
162  int ghosts[] = {0};
163 
164  auto vol = createSmartGhostVector(
165  PETSC_COMM_WORLD, m_field.get_comm_rank() == 0 ? 1 : 0, 1,
166  m_field.get_comm_rank() == 0 ? 0 : 1, ghosts);
167  auto surf_vol = smartVectorDuplicate(vol);
168 
169  // set operator to the volume element
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 OpFace("MESH_NODE_POSITIONS", surf_vol));
175  // make integration in volume (here real calculations starts)
176  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
177  domain_fe);
178  // make integration on boundary
179  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
180  boundary_fe);
181  // assemble volumes from processors and accumulate on processor of rank 0
182  CHKERR VecAssemblyBegin(vol);
183  CHKERR VecAssemblyEnd(vol);
184  CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
185  CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
186  CHKERR VecAssemblyBegin(surf_vol);
187  CHKERR VecAssemblyEnd(surf_vol);
188  CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
189  CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
190  if (m_field.get_comm_rank() == 0) {
191  double *a_vol;
192  CHKERR VecGetArray(vol, &a_vol);
193  double *a_surf_vol;
194  CHKERR VecGetArray(surf_vol, &a_surf_vol);
195  cout << "Volume = " << a_vol[0] << endl;
196  cout << "Surf Volume = " << a_surf_vol[0] << endl;
197  if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
198  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Should be zero");
199  }
200  CHKERR VecRestoreArray(vol, &a_vol);
201  CHKERR VecRestoreArray(vol, &a_surf_vol);
202  }
203 
204  }
205  }
206  CATCH_ERRORS;
207 
208  // finish work cleaning memory, getting statistics, etc.
210 
211  return 0;
212 }
Set integration rule.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:507
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
int main(int argc, char *argv[])
OpFace(const std::string &field_name, Vec vol)
const std::string getBoundaryFEName() const
Definition: Simple.hpp:222
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:517
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
int operator()(int, int, int) const
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
OpVolume(const std::string &field_name, Vec vol)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Core (interface) class.
Definition: Core.hpp:50
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:209
MoFEMErrorCode loadFile(const std::string options="PARALLEL=READ_PART;" "PARALLEL_RESOLVE_SHARED_ENTS;" "PARTITION=PARALLEL_PARTITION;")
Load mesh file.
Definition: Simple.cpp:62
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
Projection of edge entities with one mid-node on hierarchical basis.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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:116
virtual int get_comm_rank() const =0
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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:101
static char help[]
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.
Set integration rule to boundary elements.
int operator()(int, int, int) const
#define CHKERR
Inline error check.
Definition: definitions.h:596
const std::string getDomainFEName() const
Definition: Simple.hpp:221
ForcesAndSourcesCore::UserDataOperator UserDataOperator
auto createSmartGhostVector
Create smart ghost vector.
Definition: AuxPETSc.hpp:194
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:48
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:256
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61