v0.9.1
basic_moment_of_inertia.cpp
Go to the documentation of this file.
1 /**
2  * \file basic_moment_of_inertia.cpp
3  * \example basic_moment_of_inertia.cpp
4  *
5  * \brief Calculate mass and second moment of inertia.
6  *
7  * Example intend 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 
38 
39 struct Example {
40 
41  Example(MoFEM::Interface &m_field) : mField(m_field) {}
42 
43  MoFEMErrorCode runProblem();
44 
45 private:
46  MoFEM::Interface &mField;
47 
48  MoFEMErrorCode setUP();
49  MoFEMErrorCode createCommonData();
50  MoFEMErrorCode bC();
51  MoFEMErrorCode OPs();
52  MoFEMErrorCode integrateElements();
53  MoFEMErrorCode postProcess();
54  MoFEMErrorCode checkResults();
55 
56  //! [Common data]
57  struct CommonData {
58 
59  boost::shared_ptr<VectorDouble>
60  rhoAtIntegrationPts; ///< Storing density at integration point
61 
62  /**
63  * @brief Vector to indicate indices for storing, zero, first and second
64  * moment.
65  *
66  */
67  enum VecElements {
68  ZERO = 0,
78  LAST_ELEMENT
79  };
80 
82  petscVec; ///< Smart pinter which stores PETSc distributed vector
83  };
84  boost::shared_ptr<CommonData> commonDataPtr;
85 
86  //! [Common data]
87 
88  //! [Operators]
89  struct OpZero : public OpElement {
90  OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
91  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
92  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
93 
94  private:
95  boost::shared_ptr<CommonData> commonDataPtr;
96  };
97 
98  struct OpFirst : public OpElement {
99  OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
100  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
101  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
102 
103  private:
104  boost::shared_ptr<CommonData> commonDataPtr;
105  };
106 
107  struct OpSecond : public OpElement {
108  OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
109  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
110  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
111 
112  private:
113  boost::shared_ptr<CommonData> commonDataPtr;
114  };
115  //! [Operators]
116 };
117 
118 //! [Run all]
121  CHKERR setUP();
122  CHKERR createCommonData();
123  CHKERR bC();
124  CHKERR OPs();
125  CHKERR integrateElements();
126  CHKERR postProcess();
127  CHKERR checkResults();
129 }
130 //! [Run all]
131 
132 //! [Set up problem]
135  Simple *simple = mField.getInterface<Simple>();
136  CHKERR simple->getOptions();
137  CHKERR simple->loadFile("");
138  // Add field
139  CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
140  constexpr int order = 1;
141  CHKERR simple->setFieldOrder("rho", order);
142  CHKERR simple->setUp();
144 }
145 //! [Set up problem]
146 
147 //! [Create common data]
150  commonDataPtr = boost::make_shared<CommonData>();
151  commonDataPtr->petscVec = createSmartVectorMPI(
152  mField.get_comm(),
153  (!mField.get_comm_rank()) ? CommonData::LAST_ELEMENT : 0,
154  CommonData::LAST_ELEMENT);
155  commonDataPtr->rhoAtIntegrationPts = boost::make_shared<VectorDouble>();
157 }
158 //! [Create common data]
159 
160 //! [Distributions mass
163  auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
164  double *ycoord, double *zcoord) {
166  field_data[0] = 1;
168  };
169  FieldBlas *field_blas;
170  CHKERR mField.getInterface(field_blas);
171  CHKERR field_blas->setVertexDofs(set_density, "rho");
173 }
174 //! [Distributions mass]
175 
176 //! [Push operators to pipeline]
179  Basic *basic = mField.getInterface<Basic>();
180 
181  // Push operator which calculate values of densities at integration points
183  "rho", commonDataPtr->rhoAtIntegrationPts));
184 
185  // Push operator to pipeline to calculate zero moment of inertia, that is mass
186  // and when density is one everywere it is area
187  basic->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
188 
189  // Push operator to pipeline to calculate first moment of inertaia
190  basic->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
191 
192  // Push operator to pipeline to calculate second moment of inertaia
193  basic->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
194 
195  // Set integration rule. Integration rule is equal to polynomial order of
196  // density plus the 2, since to calculate second moment of inertia term x*x is
197  // present
198  auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
199 
200  // Add integration rule to element
201  CHKERR basic->setDomainRhsIntegrationRule(integration_rule);
203 }
204 //! [Push operators to pipeline]
205 
206 //! [Do calculations]
209 
210  Basic *basic = mField.getInterface<Basic>();
211  // Zero global vector
212  CHKERR VecZeroEntries(commonDataPtr->petscVec);
213 
214  // Integrate elements by executing operators in the pipeline
215  CHKERR basic->loopFiniteElements();
216 
217  // Assemble vector
218  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
219  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
221 }
222 //! [Do calculations]
223 
224 //! [Print results]
227  const double *array;
228  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
229  //! [Print results]
230  if (mField.get_comm_rank() == 0) {
231  PetscPrintf(PETSC_COMM_SELF, "Mass %6.4e\n", array[CommonData::ZERO]);
232  PetscPrintf(PETSC_COMM_SELF,
233  "First moment of inertia [ %6.4e, %6.4e, %6.4e ] \n",
234  array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
235  array[CommonData::FIRST_Z]);
236  PetscPrintf(PETSC_COMM_SELF,
237  "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
238  "%6.4e ]\n",
239  array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
240  array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
241  array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
242  }
243  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
245 }
246 //! [Print results]
247 
248 //! [Test example]
251  const double *array;
252  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
253 
254  PetscBool test = PETSC_FALSE;
255  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
256  if (mField.get_comm_rank() == 0 && test) {
257  constexpr double eps = 1e-8;
258  constexpr double expected_area = 1.;
259  constexpr double expected_first_moment = 0.;
260  constexpr double expected_second_moment = 1. / 12.;
261  if (std::abs(array[CommonData::ZERO] - expected_area) > eps)
262  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
263  "Wrong area %6.4e !+ %6.4e", expected_area,
264  array[CommonData::ZERO]);
265  for (auto i :
266  {CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z}) {
267  if (std::abs(array[i] - expected_first_moment) > eps)
268  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
269  "Wrong first moment %6.4e !+ %6.4e", expected_first_moment,
270  array[i]);
271  }
272  for (auto i : {CommonData::SECOND_XX, CommonData::SECOND_YY,
273  CommonData::SECOND_ZZ}) {
274  if (std::abs(array[i] - expected_second_moment) > eps)
275  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
276  "Wrong second moment %6.4e !+ %6.4e", expected_second_moment,
277  array[i]);
278  }
279  }
280  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
281 
283 }
284 //! [Test example]
285 
286 int main(int argc, char *argv[]) {
287 
288  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
289 
290  try {
291 
292  //! [Register MoFEM discrete manager in PETSc]
293  DMType dm_name = "DMMOFEM";
294  CHKERR DMRegister_MoFEM(dm_name);
295  //! [Register MoFEM discrete manager in PETSc
296 
297  //! [Create MoAB]
298  moab::Core mb_instance; ///< mesh database
299  moab::Interface &moab = mb_instance; ///< mesh database interface
300  //! [Create MoAB]
301 
302  //! [Create MoFEM]
303  MoFEM::Core core(moab); ///< finite element database
304  MoFEM::Interface &m_field = core; ///< finite element database insterface
305  //! [Create MoFEM]
306 
307  Example ex(m_field);
308  CHKERR ex.runProblem();
309  }
310  CATCH_ERRORS;
311 
313 }
314 
315 //! [FirstOp]
317  EntData &data) {
319  if (type == MBVERTEX) {
320 
321  const int nb_integration_pts =
322  getGaussPts().size2(); // Number of integration points
323  auto t_w = getFTensor0IntegrationWeight(); // Integration weights
324  auto t_rho = getFTensor0FromVec(*(
325  commonDataPtr->rhoAtIntegrationPts)); // Density at integration weights
326  const double volume = getMeasure();
327 
328  // Integrate area of the element
329  double element_local_value = 0;
330  for (int gg = 0; gg != nb_integration_pts; ++gg) {
331  element_local_value += t_w * t_rho * volume;
332  ++t_w;
333  ++t_rho;
334  }
335 
336  // Assemble area of the element into global PETSc vector
337  const int index = CommonData::ZERO;
338  CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
339  ADD_VALUES);
340  }
342 }
343 //! [FirstOps]
344 
346  EntData &data) {
348  if (type == MBVERTEX) {
349  const int nb_integration_pts = getGaussPts().size2();
350  auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
351  auto t_rho = getFTensor0FromVec(*(
352  commonDataPtr->rhoAtIntegrationPts)); ///< Density at integration points
353  auto t_coords =
354  getFTensor1CoordsAtGaussPts(); ///< Coordinates at integration points
355  const double volume = getMeasure(); ///< Get Volume of element
356 
358  t_s; ///< First moment of inertia is tensor of rank 1, i.e. vector.
359  t_s(i) = 0; // Zero entries
360 
361  // Integrate
362  for (int gg = 0; gg != nb_integration_pts; ++gg) {
363 
364  t_s(i) += t_w * t_rho * volume * t_coords(i);
365 
366  ++t_w; // move weight to next integration pts
367  ++t_rho; // move density
368  ++t_coords; // move coordinate
369  }
370 
371  // Set array of indices
372  constexpr std::array<int, 3> indices = {
373  CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z};
374 
375  // Assemble first moment of inertia
376  CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(), &t_s(0),
377  ADD_VALUES);
378  }
380 }
381 
382 //! [SecondOp]
384  EntData &data) {
386  if (type == MBVERTEX) {
387 
388  const int nb_integration_pts = getGaussPts().size2();
389  auto t_w = getFTensor0IntegrationWeight();
390  auto t_rho = getFTensor0FromVec(*(commonDataPtr->rhoAtIntegrationPts));
391  auto t_coords = getFTensor1CoordsAtGaussPts();
392  const double volume = getMeasure();
393 
394  // Create storage for symmetric tensor
395  std::array<double, 6> element_local_value;
396 
397  // Crate symmetric tensor with points to the storrage
399  &element_local_value[CommonData::SECOND_XX],
400  &element_local_value[CommonData::SECOND_XY],
401  &element_local_value[CommonData::SECOND_XZ],
402  &element_local_value[CommonData::SECOND_YY],
403  &element_local_value[CommonData::SECOND_YZ],
404  &element_local_value[CommonData::SECOND_ZZ]);
405 
406  // Integate
407  for (int gg = 0; gg != nb_integration_pts; ++gg) {
408 
409  // Symbol "^" indicate multiplication which yield symmetric tensor
410  t_I(i, j) += (t_w * t_rho * volume) * (t_coords(i) ^ t_coords(j));
411 
412  ++t_w;
413  ++t_rho;
414  ++t_coords;
415  }
416 
417  // Set array of indices
418  constexpr std::array<int, 6> indices = {
419  CommonData::SECOND_XX, CommonData::SECOND_XY, CommonData::SECOND_XZ,
420  CommonData::SECOND_YY, CommonData::SECOND_YZ, CommonData::SECOND_ZZ};
421 
422  // Assemble second moment of inertia
423  CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
424  &element_local_value[0], ADD_VALUES);
425  }
427 }
428 //! [SecondOp]
FTensor::Index< 'i', 3 > i
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
static char help[]
Deprecated interface functions.
FTensor::Index< 'j', 3 > j
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< VectorDouble > rhoAtIntegrationPts
Storing density at integration point.
DataForcesAndSourcesCore::EntData EntData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[SecondOp]
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: Basic.hpp:466
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:506
boost::shared_ptr< CommonData > commonDataPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:513
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()
[Run problem]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOp]
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 OPs()
[Boundary condition]
Basic interface.
Definition: Basic.hpp:36
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode bC()
[Create common data]
int main(int argc, char *argv[])
[Test example]
Basic algebra on fields.
Definition: FieldBlas.hpp:34
constexpr int order
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:151
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:104
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: Basic.cpp:36
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode setUP()
[Set up problem]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOps]
MoFEMErrorCode postProcess()
[Solve]
boost::shared_ptr< CommonData > commonDataPtr
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
Data on single entity (This is passed as argument to DataOperator::doWork)
VecElements
Vector to indicate indices for storing, zero, first and second moment.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: Basic.hpp:311
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
continuous field
Definition: definitions.h:176
static const double eps
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 pinter which stores PETSc distributed vector.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user functionExample:
Definition: FieldBlas.cpp:175