v0.9.2
lesson1_moment_of_inertia.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 
41  MoFEMErrorCode runProblem();
42 
43 private:
45 
46  MoFEMErrorCode setUp();
47  MoFEMErrorCode createCommonData();
48  MoFEMErrorCode setFieldValues();
49  MoFEMErrorCode pushOperators();
50  MoFEMErrorCode integrateElements();
51  MoFEMErrorCode postProcess();
52  MoFEMErrorCode checkResults();
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(new OpCalculateScalarFieldValues(
208  "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
209 
210  // Push an operator to pipeline to calculate zero moment of inertia (mass)
211  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
212 
213  // Push an operator to the pipeline to calculate first moment of inertaia
214  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
215 
216  // Push an operator to the pipeline to calculate second moment of inertaia
217  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
218 
219  // Set integration rule. Integration rule is equal to the polynomial order of
220  // the density field plus 2, since under the integral of the second moment of
221  // inertia term x*x is present
222  auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
223 
224  // Add integration rule to the element
225  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
227 }
228 //! [Push operators to pipeline]
229 
230 //! [Integrate]
233  // Zero global vector
234  CHKERR VecZeroEntries(commonDataPtr->petscVec);
235 
236  // Integrate elements by executing operators in the pipeline
237  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
238  CHKERR pipeline_mng->loopFiniteElements();
239 
240  // Assemble MPI vector
241  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
242  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
244 }
245 //! [Integrate]
246 
247 //! [Print results]
250  const double *array;
251  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
252  if (mField.get_comm_rank() == 0) {
253  PetscPrintf(PETSC_COMM_SELF, "Mass %6.4e\n", array[CommonData::ZERO]);
254  PetscPrintf(PETSC_COMM_SELF,
255  "First moment of inertia [ %6.4e, %6.4e, %6.4e ] \n",
256  array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
257  array[CommonData::FIRST_Z]);
258  PetscPrintf(PETSC_COMM_SELF,
259  "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
260  "%6.4e ]\n",
261  array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
262  array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
263  array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
264  }
265  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
267 }
268 //! [Print results]
269 
270 //! [Test example]
273  const double *array;
274  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
275 
276  PetscBool test = PETSC_FALSE;
277  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
278  if (mField.get_comm_rank() == 0 && test) {
279  constexpr double eps = 1e-8;
280  constexpr double expected_volume = 1.;
281  constexpr double expected_first_moment = 0.;
282  constexpr double expected_second_moment = 1. / 12.;
283  if (std::abs(array[CommonData::ZERO] - expected_volume) > eps)
284  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
285  "Wrong area %6.4e != %6.4e", expected_volume,
286  array[CommonData::ZERO]);
287  for (auto i :
288  {CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z}) {
289  if (std::abs(array[i] - expected_first_moment) > eps)
290  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
291  "Wrong first moment %6.4e != %6.4e", expected_first_moment,
292  array[i]);
293  }
294  for (auto i : {CommonData::SECOND_XX, CommonData::SECOND_YY,
295  CommonData::SECOND_ZZ}) {
296  if (std::abs(array[i] - expected_second_moment) > eps)
297  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
298  "Wrong second moment %6.4e != %6.4e", expected_second_moment,
299  array[i]);
300  }
301  }
302  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
303 
305 }
306 //! [Test example]
307 
308 //! [main]
309 int main(int argc, char *argv[]) {
310 
311  // Initialisation of MoFEM/PETSc and MoAB data structures
312  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
313 
314  // Error handling
315  try {
316 
317  //! [Register MoFEM discrete manager in PETSc]
318  DMType dm_name = "DMMOFEM";
319  CHKERR DMRegister_MoFEM(dm_name);
320  //! [Register MoFEM discrete manager in PETSc]
321 
322  //! [Create MoAB]
323  moab::Core mb_instance; ///< mesh database
324  moab::Interface &moab = mb_instance; ///< mesh database interface
325  //! [Create MoAB]
326 
327  //! [Create MoFEM]
328  MoFEM::Core core(moab); ///< finite element database
329  MoFEM::Interface &m_field = core; ///< finite element database interface
330  //! [Create MoFEM]
331 
332  Example ex(m_field);
333  CHKERR ex.runProblem();
334  }
335  CATCH_ERRORS;
336 
338 }
339 //! [main]
340 
341 //! [ZeroOp]
343  EntData &data) {
345  if (type == MBVERTEX) {
346 
347  const int nb_integration_pts =
348  getGaussPts().size2(); // Number of integration points
349  auto t_w = getFTensor0IntegrationWeight(); // Integration weights
350  auto t_rho = getFTensor0FromVec(
351  commonDataPtr->rhoAtIntegrationPts); // Density at integration weights
352  const double volume = getMeasure();
353 
354  // Integrate area of the element
355  double element_local_value = 0;
356  for (int gg = 0; gg != nb_integration_pts; ++gg) {
357  element_local_value += t_w * t_rho * volume;
358  ++t_w;
359  ++t_rho;
360  }
361 
362  // Assemble area of the element into global PETSc vector
363  const int index = CommonData::ZERO;
364  CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
365  ADD_VALUES);
366  }
368 }
369 //! [ZeroOp]
370 
371 //! [FirstOp]
373  EntData &data) {
375  const int nb_integration_pts = getGaussPts().size2();
376  auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
377  auto t_rho = getFTensor0FromVec(
378  commonDataPtr->rhoAtIntegrationPts); ///< Density at integration points
379  auto t_x =
380  getFTensor1CoordsAtGaussPts(); ///< Coordinates at integration points
381  const double volume = getMeasure(); ///< Get Volume of element
382 
384 
386  t_s; ///< First moment of inertia is tensor of rank 1, i.e. vector.
387  t_s(i) = 0; // Zero entries
388 
389  // Integrate
390  for (int gg = 0; gg != nb_integration_pts; ++gg) {
391 
392  t_s(i) += t_w * t_rho * volume * t_x(i);
393 
394  ++t_w; // move weight to next integration pts
395  ++t_rho; // move density
396  ++t_x; // move coordinate
397  }
398 
399  // Set array of indices
400  constexpr std::array<int, 3> indices = {
401  CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z};
402 
403  // Assemble first moment of inertia
404  CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(), &t_s(0),
405  ADD_VALUES);
407 }
408 //! [FirstOp]
409 
410 //! [SecondOp]
412  EntData &data) {
414 
415  const int nb_integration_pts = getGaussPts().size2();
416  auto t_w = getFTensor0IntegrationWeight();
417  auto t_rho = getFTensor0FromVec(commonDataPtr->rhoAtIntegrationPts);
418  auto t_x = getFTensor1CoordsAtGaussPts();
419  const double volume = getMeasure();
420 
421  // Create storage for a symmetric tensor
422  std::array<double, 6> element_local_value;
423  std::fill(element_local_value.begin(), element_local_value.end(), 0.0);
424 
425  // Create symmetric tensor using pointers to the storage
427  &element_local_value[0], &element_local_value[1], &element_local_value[2],
428  &element_local_value[3], &element_local_value[4],
429  &element_local_value[5]);
430 
433 
434  // Integrate
435  for (int gg = 0; gg != nb_integration_pts; ++gg) {
436 
437  // Symbol "^" indicates outer product which yields a symmetric tensor
438  t_I(i, j) += (t_w * t_rho * volume) * (t_x(i) ^ t_x(j));
439 
440  ++t_w; // move weight pointer to the next integration point
441  ++t_rho; // move density pointer
442  ++t_x; // move coordinate pointer
443  }
444 
445  // Set array of indices
446  constexpr std::array<int, 6> indices = {
447  CommonData::SECOND_XX, CommonData::SECOND_XY, CommonData::SECOND_XZ,
448  CommonData::SECOND_YY, CommonData::SECOND_YZ, CommonData::SECOND_ZZ};
449 
450  // Assemble second moment of inertia
451  CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
452  &element_local_value[0], ADD_VALUES);
453 
455 }
456 //! [SecondOp]
OpZero(boost::shared_ptr< CommonData > &common_data_ptr)
Volume finite element baseUser is implementing own operator at Gauss point level, by class derived fr...
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Deprecated interface functions.
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< VectorDouble > getRhoAtIntegrationPtsPtr()
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOp]
MoFEM::Interface & mField
int main(int argc, char *argv[])
[Test example]
MoFEMErrorCode integrateElements()
[Push operators to pipeline]
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
VectorDouble rhoAtIntegrationPts
Storing density at integration points.
boost::shared_ptr< CommonData > commonDataPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Example(MoFEM::Interface &m_field)
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
boost::shared_ptr< CommonData > commonDataPtr
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
OpSecond(boost::shared_ptr< CommonData > &common_data_ptr)
auto createSmartVectorMPI
Create MPI Vector.
Definition: AuxPETSc.hpp:224
MoFEMErrorCode runProblem()
[Operator]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[main]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode checkResults()
[Print results]
Basic algebra on fields.
Definition: FieldBlas.hpp:34
MoFEMErrorCode pushOperators()
[Set density distribution]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
MoFEMErrorCode setFieldValues()
[Create common data]
DataForcesAndSourcesCore::EntData EntData
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
PipelineManager interface.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
static char help[]
constexpr int order
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[ZeroOp]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< CommonData > commonDataPtr
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
VecElements
Vector to indicate indices for storing zero, first and second moments of inertia.
#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
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
OpFirst(boost::shared_ptr< CommonData > &common_data_ptr)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Get value at integration points for scalar field.
SmartPetscObj< Vec > petscVec
Smart pointer which stores PETSc distributed vector.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:69
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user functionExample:
Definition: FieldBlas.cpp:172
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
MoFEMErrorCode setUp()
[Run all]