v0.9.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>
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:673
Deprecated interface functions.
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
Get the Boundary FE Name.
Definition: Simple.hpp:288
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:707
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:507
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:238
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
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:287
virtual int get_comm_rank() const =0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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:269
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:602
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:281
auto createSmartGhostVector
Create smart ghost vector.
Definition: AuxPETSc.hpp:208
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
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.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:483
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:77
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26