v0.13.0
integration.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson1_moment_of_inertia.cpp
3  * \example lesson1_moment_of_inertia.cpp
4  *
5  * \brief Calculate zero, first and second moments of inertia.
6  *
7  * Example intended to show how to write user data operator and integrate
8  * scalar field.
9  *
10  */
11 
12 /* This file is part of MoFEM.
13  * MoFEM is free software: you can redistribute it and/or modify it under
14  * the terms of the GNU Lesser General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21  * License for more details.
22  *
23  * You should have received a copy of the GNU Lesser General Public
24  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #include <MoFEM.hpp>
27 
28 using namespace MoFEM;
29 
30 static char help[] = "...\n\n";
31 
35 
36 //! [Example]
37 struct Example {
38 
39  Example(MoFEM::Interface &m_field) : mField(m_field) {}
40 
42 
43 private:
44  MoFEM::Interface &mField;
45 
46  MoFEMErrorCode setUp();
48  MoFEMErrorCode setFieldValues();
49  MoFEMErrorCode pushOperators();
50  MoFEMErrorCode integrateElements();
53 
54  struct CommonData;
55 
56  boost::shared_ptr<CommonData> commonDataPtr;
57 
58  struct OpZero;
59 
60  struct OpFirst;
61 
62  struct OpSecond;
63 };
64 //! [Example]
65 
66 //! [Common data]
68  : public boost::enable_shared_from_this<Example::CommonData> {
69 
70  VectorDouble rhoAtIntegrationPts; ///< Storing density at integration points
71 
72  inline boost::shared_ptr<VectorDouble> getRhoAtIntegrationPtsPtr() {
73  return boost::shared_ptr<VectorDouble>(shared_from_this(),
74  &rhoAtIntegrationPts);
75  }
76 
77  /**
78  * @brief Vector to indicate indices for storing zero, first and second
79  * moments of inertia.
80  *
81  */
82  enum VecElements {
83  ZERO = 0,
93  LAST_ELEMENT
94  };
95 
97  petscVec; ///< Smart pointer which stores PETSc distributed vector
98 };
99 //! [Common data]
100 
101 struct Example::OpZero : public OpElement {
102  OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
103  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
104  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
105  }
106  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
107 
108 private:
109  boost::shared_ptr<CommonData> commonDataPtr;
110 };
111 
112 struct Example::OpFirst : public OpElement {
113  OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
114  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
115  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
116  }
117  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
118 
119 private:
120  boost::shared_ptr<CommonData> commonDataPtr;
121 };
122 
123 //! [Operator]
124 struct Example::OpSecond : public OpElement {
125  OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
126  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
127  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
128  }
129  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
130 
131 private:
132  boost::shared_ptr<CommonData> commonDataPtr;
133 };
134 //! [Operator]
135 
136 //! [Run all]
139  CHKERR setUp();
140  CHKERR createCommonData();
141  CHKERR setFieldValues();
142  CHKERR pushOperators();
143  CHKERR integrateElements();
144  CHKERR postProcess();
145  CHKERR checkResults();
147 }
148 //! [Run all]
149 
150 //! [Set up problem]
153  Simple *simple = mField.getInterface<Simple>();
154  CHKERR simple->getOptions();
155  CHKERR simple->loadFile();
156  // Add field
157  CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
158  constexpr int order = 1;
159  CHKERR simple->setFieldOrder("rho", order);
160  CHKERR simple->setUp();
162 }
163 //! [Set up problem]
164 
165 //! [Create common data]
168  commonDataPtr = boost::make_shared<CommonData>();
169 
170  int local_size;
171  if (mField.get_comm_rank() == 0) // get_comm_rank() gets processor number
172  // processor 0
173  local_size = CommonData::LAST_ELEMENT; // last element gives size of vector
174  else
175  // other processors (e.g. 1, 2, 3, etc.)
176  local_size = 0; // local size of vector is zero on other processors
177 
178  commonDataPtr->petscVec = createSmartVectorMPI(mField.get_comm(), local_size,
179  CommonData::LAST_ELEMENT);
180 
182 }
183 //! [Create common data]
184 
185 //! [Set density distribution]
188  auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
189  double *ycoord, double *zcoord) {
191  field_data[0] = 1;
193  };
194  FieldBlas *field_blas;
195  CHKERR mField.getInterface(field_blas);
196  CHKERR field_blas->setVertexDofs(set_density, "rho");
198 }
199 //! [Set density distribution]
200 
201 //! [Push operators to pipeline]
204  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
205 
206  // Push an operator which calculates values of density at integration points
207  pipeline_mng->getOpDomainRhsPipeline().push_back(
209  "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
210 
211  // Push an operator to pipeline to calculate zero moment of inertia (mass)
212  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
213 
214  // Push an operator to the pipeline to calculate first moment of inertaia
215  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
216 
217  // Push an operator to the pipeline to calculate second moment of inertaia
218  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
219 
220  // Set integration rule. Integration rule is equal to the polynomial order of
221  // the density field plus 2, since under the integral of the second moment of
222  // inertia term x*x is present
223  auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
224 
225  // Add integration rule to the element
228 }
229 //! [Push operators to pipeline]
230 
231 //! [Integrate]
234  // Zero global vector
235  CHKERR VecZeroEntries(commonDataPtr->petscVec);
236 
237  // Integrate elements by executing operators in the pipeline
238  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
239  CHKERR pipeline_mng->loopFiniteElements();
240 
241  // Assemble MPI vector
242  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
243  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
245 }
246 //! [Integrate]
247 
248 //! [Print results]
251  const double *array;
252  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
253  if (mField.get_comm_rank() == 0) {
254  MOFEM_LOG_C("SELF", Sev::inform, "Mass %6.4e", array[CommonData::ZERO]);
255  MOFEM_LOG_C("SELF", Sev::inform,
256  "First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
257  array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
258  array[CommonData::FIRST_Z]);
259  MOFEM_LOG_C("SELF", Sev::inform,
260  "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
261  "%6.4e ]",
262  array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
263  array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
264  array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
265  }
266  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
268 }
269 //! [Print results]
270 
271 //! [Test example]
274  const double *array;
275  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
276 
277  PetscBool test = PETSC_FALSE;
278  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
279  if (mField.get_comm_rank() == 0 && test) {
280  constexpr double eps = 1e-8;
281  constexpr double expected_volume = 1.;
282  constexpr double expected_first_moment = 0.;
283  constexpr double expected_second_moment = 1. / 12.;
284  if (std::abs(array[CommonData::ZERO] - expected_volume) > eps)
285  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
286  "Wrong area %6.4e != %6.4e", expected_volume,
287  array[CommonData::ZERO]);
288  for (auto i :
289  {CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z}) {
290  if (std::abs(array[i] - expected_first_moment) > eps)
291  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
292  "Wrong first moment %6.4e != %6.4e", expected_first_moment,
293  array[i]);
294  }
295  for (auto i : {CommonData::SECOND_XX, CommonData::SECOND_YY,
296  CommonData::SECOND_ZZ}) {
297  if (std::abs(array[i] - expected_second_moment) > eps)
298  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
299  "Wrong second moment %6.4e != %6.4e", expected_second_moment,
300  array[i]);
301  }
302  }
303  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
304 
306 }
307 //! [Test example]
308 
309 //! [main]
310 int main(int argc, char *argv[]) {
311 
312  // Initialisation of MoFEM/PETSc and MOAB data structures
313  const char param_file[] = "param_file.petsc";
314  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
315 
316  // Error handling
317  try {
318 
319  //! [Register MoFEM discrete manager in PETSc]
320  DMType dm_name = "DMMOFEM";
321  CHKERR DMRegister_MoFEM(dm_name);
322  //! [Register MoFEM discrete manager in PETSc]
323 
324  //! [Create MoAB]
325  moab::Core mb_instance; ///< mesh database
326  moab::Interface &moab = mb_instance; ///< mesh database interface
327  //! [Create MoAB]
328 
329  //! [Create MoFEM]
330  MoFEM::Core core(moab); ///< finite element database
331  MoFEM::Interface &m_field = core; ///< finite element database interface
332  //! [Create MoFEM]
333 
334  Example ex(m_field);
335  CHKERR ex.runProblem();
336  }
337  CATCH_ERRORS;
338 
340 }
341 //! [main]
342 
343 //! [ZeroOp]
345  EntData &data) {
347  if (type == MBVERTEX) {
348 
349  const int nb_integration_pts =
350  getGaussPts().size2(); // Number of integration points
351  auto t_w = getFTensor0IntegrationWeight(); // Integration weights
352  auto t_rho = getFTensor0FromVec(
353  commonDataPtr->rhoAtIntegrationPts); // Density at integration weights
354  const double volume = getMeasure();
355 
356  // Integrate area of the element
357  double element_local_value = 0;
358  for (int gg = 0; gg != nb_integration_pts; ++gg) {
359  element_local_value += t_w * t_rho * volume;
360  ++t_w;
361  ++t_rho;
362  }
363 
364  // Assemble area of the element into global PETSc vector
365  const int index = CommonData::ZERO;
366  CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
367  ADD_VALUES);
368  }
370 }
371 //! [ZeroOp]
372 
373 //! [FirstOp]
375  EntData &data) {
377  const int nb_integration_pts = getGaussPts().size2();
378  auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
379  auto t_rho = getFTensor0FromVec(
380  commonDataPtr->rhoAtIntegrationPts); ///< Density at integration points
381  auto t_x =
382  getFTensor1CoordsAtGaussPts(); ///< Coordinates at integration points
383  const double volume = getMeasure(); ///< Get Volume of element
384 
386 
388  t_s; ///< First moment of inertia is tensor of rank 1, i.e. vector.
389  t_s(i) = 0; // Zero entries
390 
391  // Integrate
392  for (int gg = 0; gg != nb_integration_pts; ++gg) {
393 
394  t_s(i) += t_w * t_rho * volume * t_x(i);
395 
396  ++t_w; // move weight to next integration pts
397  ++t_rho; // move density
398  ++t_x; // move coordinate
399  }
400 
401  // Set array of indices
402  constexpr std::array<int, 3> indices = {
403  CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z};
404 
405  // Assemble first moment of inertia
406  CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(), &t_s(0),
407  ADD_VALUES);
409 }
410 //! [FirstOp]
411 
412 //! [SecondOp]
414  EntData &data) {
416 
417  const int nb_integration_pts = getGaussPts().size2();
418  auto t_w = getFTensor0IntegrationWeight();
419  auto t_rho = getFTensor0FromVec(commonDataPtr->rhoAtIntegrationPts);
420  auto t_x = getFTensor1CoordsAtGaussPts();
421  const double volume = getMeasure();
422 
423  // Create storage for a symmetric tensor
424  std::array<double, 6> element_local_value;
425  std::fill(element_local_value.begin(), element_local_value.end(), 0.0);
426 
427  // Create symmetric tensor using pointers to the storage
429  &element_local_value[0], &element_local_value[1], &element_local_value[2],
430  &element_local_value[3], &element_local_value[4],
431  &element_local_value[5]);
432 
435 
436  // Integrate
437  for (int gg = 0; gg != nb_integration_pts; ++gg) {
438 
439  // Symbol "^" indicates outer product which yields a symmetric tensor
440  t_I(i, j) += (t_w * t_rho * volume) * (t_x(i) ^ t_x(j));
441 
442  ++t_w; // move weight pointer to the next integration point
443  ++t_rho; // move density pointer
444  ++t_x; // move coordinate pointer
445  }
446 
447  // Set array of indices
448  constexpr std::array<int, 6> indices = {
449  CommonData::SECOND_XX, CommonData::SECOND_XY, CommonData::SECOND_XZ,
450  CommonData::SECOND_YY, CommonData::SECOND_YZ, CommonData::SECOND_ZZ};
451 
452  // Assemble second moment of inertia
453  CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
454  &element_local_value[0], ADD_VALUES);
455 
457 }
458 //! [SecondOp]
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static const double eps
#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
#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 integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FTensor::Index< 'i', SPACE_DIM > i
int main(int argc, char *argv[])
[Test example]
EntitiesFieldData::EntData EntData
Definition: integration.cpp:34
static char help[]
Definition: integration.cpp:30
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:126
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto createSmartVectorMPI
Create MPI Vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< VectorDouble > getRhoAtIntegrationPtsPtr()
Definition: integration.cpp:72
VecElements
Vector to indicate indices for storing zero, first and second moments of inertia.
Definition: integration.cpp:82
SmartPetscObj< Vec > petscVec
Smart pointer which stores PETSc distributed vector.
Definition: integration.cpp:97
VectorDouble rhoAtIntegrationPts
Storing density at integration points.
Definition: integration.cpp:70
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[ZeroOp]
OpFirst(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOp]
boost::shared_ptr< CommonData > commonDataPtr
OpSecond(boost::shared_ptr< CommonData > &common_data_ptr)
[Common data]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[main]
boost::shared_ptr< CommonData > commonDataPtr
OpZero(boost::shared_ptr< CommonData > &common_data_ptr)
[Example]
MoFEMErrorCode pushOperators()
[Set density distribution]
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:54
MoFEMErrorCode setUp()
[Run all]
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: integration.cpp:39
MoFEMErrorCode runProblem()
MoFEMErrorCode postProcess()
MoFEMErrorCode setFieldValues()
[Create common data]
MoFEMErrorCode integrateElements()
[Push operators to pipeline]
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)
Basic algebra on fields.
Definition: FieldBlas.hpp:31
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:333
Get value at integration points for scalar field.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.