v0.9.1
second_moment_of_inertia.cpp
Go to the documentation of this file.
1 /**
2  * \file second_moment_of_inertia.cpp
3  * \example second_moment_of_inertia.cpp
4  *
5  * Using Basic interface calculate the divergence of base functions, and
6  * integral of flux on the boundary. Since the h-div space is used, volume
7  * integral and boundary integral should give the same result.
8  */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
33 
34 struct Example {
35 
36  Example(MoFEM::Interface &m_field) : mField(m_field) {}
37 
38  MoFEMErrorCode runProblem();
39 
40 private:
41  MoFEM::Interface &mField;
42 
43  MoFEMErrorCode setUP();
44  MoFEMErrorCode createCommonData();
45  MoFEMErrorCode bC();
46  MoFEMErrorCode OPs();
47  MoFEMErrorCode integrateElements();
48  MoFEMErrorCode postProcess();
49  MoFEMErrorCode checkResults();
50 
51  //! [Common data]
52  struct CommonData {
53  boost::shared_ptr<VectorDouble> rhoAtIntegrationPts;
54  enum VecElements {
55  ZERO = 0,
65  LAST_ELEMENT
66  };
68  };
69  boost::shared_ptr<CommonData> commonDataPtr;
70  //! [Common data]
71 
72  //! [Operators]
73  struct OpZero : public OpElement {
74  OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
75  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
76  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
77 
78  private:
79  boost::shared_ptr<CommonData> commonDataPtr;
80  };
81 
82  struct OpFirst : public OpElement {
83  OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
84  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
85  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
86 
87  private:
88  boost::shared_ptr<CommonData> commonDataPtr;
89  };
90 
91  struct OpSecond : public OpElement {
92  OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
93  : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {}
94  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
95 
96  private:
97  boost::shared_ptr<CommonData> commonDataPtr;
98  };
99  //! [Operators]
100 };
101 
102 //! [Run all]
105  CHKERR setUP();
106  CHKERR createCommonData();
107  CHKERR bC();
108  CHKERR OPs();
109  CHKERR integrateElements();
110  CHKERR postProcess();
111  CHKERR checkResults();
113 }
114 //! [Run all]
115 
116 //! [Set up problem]
119  Simple *simple = mField.getInterface<Simple>();
120  CHKERR simple->getOptions();
121  CHKERR simple->loadFile("");
122  // Add field
123  CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE,
124  1);
125  constexpr int order = 1;
126  CHKERR simple->setFieldOrder("rho", order);
127  CHKERR simple->setUp();
129 }
130 //! [Set up problem]
131 
132 //! [Create common data]
135  commonDataPtr = boost::make_shared<CommonData>();
136  commonDataPtr->petscVec = createSmartVectorMPI(
137  mField.get_comm(),
138  (!mField.get_comm_rank()) ? CommonData::LAST_ELEMENT : 0,
139  CommonData::LAST_ELEMENT);
140  commonDataPtr->rhoAtIntegrationPts = boost::make_shared<VectorDouble>();
142 }
143 //! [Create common data]
144 
145 //! [Distributions mass
148  auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
149  double *ycoord, double *zcoord) {
151  field_data[0] = 1;
153  };
154  FieldBlas *field_blas;
155  CHKERR mField.getInterface(field_blas);
156  CHKERR field_blas->setVertexDofs(set_density, "rho");
158 }
159 //! [Distributions mass]
160 
161 //! [Push operators to pipeline]
164  Basic *basic = mField.getInterface<Basic>();
165  basic->getOpDomainRhsPipeline().push_back(
167  commonDataPtr->rhoAtIntegrationPts));
168  basic->getOpDomainRhsPipeline().push_back(
169  new OpZero(commonDataPtr));
170  basic->getOpDomainRhsPipeline().push_back(
171  new OpFirst(commonDataPtr));
172  basic->getOpDomainRhsPipeline().push_back(
173  new OpSecond(commonDataPtr));
174  auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
175  CHKERR basic->setDomainRhsIntegrationRule(integration_rule);
177 }
178 //! [Push operators to pipeline]
179 
180 //! [Do calculations]
183  Basic *basic = mField.getInterface<Basic>();
184  CHKERR VecZeroEntries(commonDataPtr->petscVec);
185  CHKERR basic->loopFiniteElements();
186  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
187  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
189 }
190 //! [Do calculations]
191 
192 //! [Print results]
195  const double *array;
196  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
197  //! [Print results]
198  if (mField.get_comm_rank() == 0) {
199  PetscPrintf(PETSC_COMM_SELF, "Mass %6.4e\n", array[CommonData::ZERO]);
200  PetscPrintf(PETSC_COMM_SELF,
201  "First moment of inertia [ %6.4e, %6.4e, %6.4e ] \n",
202  array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
203  array[CommonData::FIRST_Z]);
204  PetscPrintf(PETSC_COMM_SELF,
205  "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
206  "%6.4e ]\n",
207  array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
208  array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
209  array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
210  }
211  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
213 }
214 //! [Print results]
215 
216 //! [Test example]
219  const double *array;
220  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
221 
222  PetscBool test = PETSC_FALSE;
223  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
224  if (mField.get_comm_rank() == 0 && test) {
225  constexpr double eps = 1e-8;
226  constexpr double expected_area = 1.;
227  constexpr double expected_first_moment = 0.;
228  constexpr double expected_second_moment = 1. / 12.;
229  if (std::abs(array[CommonData::ZERO] - expected_area) > eps)
230  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
231  "Wrong area %6.4e !+ %6.4e", expected_area,
232  array[CommonData::ZERO]);
233  for (auto i :
234  {CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z}) {
235  if (std::abs(array[i] - expected_first_moment) > eps)
236  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
237  "Wrong first moment %6.4e !+ %6.4e", expected_first_moment,
238  array[i]);
239  }
240  for (auto i : {CommonData::SECOND_XX, CommonData::SECOND_YY,
241  CommonData::SECOND_ZZ}) {
242  if (std::abs(array[i] - expected_second_moment) > eps)
243  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
244  "Wrong second moment %6.4e !+ %6.4e", expected_second_moment,
245  array[i]);
246  }
247  }
248  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
249 
251 }
252 //! [Test example]
253 
254 int main(int argc, char *argv[]) {
255 
256  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
257 
258  try {
259 
260  //! [Register MoFEM discrete manager in PETSc]
261  DMType dm_name = "DMMOFEM";
262  CHKERR DMRegister_MoFEM(dm_name);
263  //! [Register MoFEM discrete manager in PETSc
264 
265  //! [Create MoAB]
266  moab::Core mb_instance; ///< mesh database
267  moab::Interface &moab = mb_instance; ///< mesh database interface
268  //! [Create MoAB]
269 
270  //! [Create MoFEM]
271  MoFEM::Core core(moab); ///< finite element database
272  MoFEM::Interface &m_field = core; ///< finite element database insterface
273  //! [Create MoFEM]
274 
275  Example ex(m_field);
276  CHKERR ex.runProblem();
277 
278  }
279  CATCH_ERRORS;
280 
282 }
283 
284 //! [FirstOp]
285 MoFEMErrorCode Example::OpZero::doWork(int side, EntityType type, EntData &data) {
287  if (type == MBVERTEX) {
288  const int nb_integration_pts = getGaussPts().size2();
289  auto t_w = getFTensor0IntegrationWeight();
290  auto t_rho = getFTensor0FromVec(*(commonDataPtr->rhoAtIntegrationPts));
292  const double volume = getVolume();
293  double element_local_value = 0;
294  for (int gg = 0; gg != nb_integration_pts; ++gg) {
295  element_local_value += t_w * t_rho * volume;
296  ++t_w;
297  ++t_rho;
298  }
299  const int index = CommonData::ZERO;
300  CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
301  ADD_VALUES);
302  }
304 }
305 //! [FirstOps]
306 
307 MoFEMErrorCode Example::OpFirst::doWork(int side, EntityType type, EntData &data) {
309  if (type == MBVERTEX) {
310  const int nb_integration_pts = getGaussPts().size2();
311  auto t_w = getFTensor0IntegrationWeight();
312  auto t_rho = getFTensor0FromVec(*(commonDataPtr->rhoAtIntegrationPts));
313  auto t_coords = getFTensor1CoordsAtGaussPts();
314  const double volume = getVolume();
315 
316  VectorDouble3 element_local_value(3);
319  &element_local_value[0], &element_local_value[1],
320  &element_local_value[2]);
321 
322  t_s(i) = 0;
323  for (int gg = 0; gg != nb_integration_pts; ++gg) {
324  t_s(i) += t_w * t_rho * volume * t_coords(i);
325  ++t_w;
326  ++t_rho;
327  ++t_coords;
328  }
329  constexpr std::array<int, 3> indices = {
330  CommonData::FIRST_X, CommonData::FIRST_Y, CommonData::FIRST_Z};
331  CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(),
332  &element_local_value[0], ADD_VALUES);
333  }
335 }
336 
337 //! [SecondOp]
338 MoFEMErrorCode Example::OpSecond::doWork(int side, EntityType type, EntData &data) {
340  if (type == MBVERTEX) {
341  const int nb_integration_pts = getGaussPts().size2();
342  auto t_w = getFTensor0IntegrationWeight();
343  auto t_rho = getFTensor0FromVec(*(commonDataPtr->rhoAtIntegrationPts));
344  auto t_coords = getFTensor1CoordsAtGaussPts();
345  const double volume = getVolume();
346 
347  VectorDouble6 element_local_value(6);
351  &element_local_value[0], &element_local_value[1],
352  &element_local_value[2], &element_local_value[3],
353  &element_local_value[4], &element_local_value[5]);
354 
355  for (int gg = 0; gg != nb_integration_pts; ++gg) {
356  t_I(i, j) += (t_w * t_rho * volume) * (t_coords(i) ^ t_coords(j));
357  ++t_w;
358  ++t_rho;
359  ++t_coords;
360  }
361 
362  constexpr std::array<int, 6> indices = {
363  CommonData::SECOND_XX, CommonData::SECOND_XY, CommonData::SECOND_XZ,
364  CommonData::SECOND_YY, CommonData::SECOND_YZ, CommonData::SECOND_ZZ};
365  CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
366  &element_local_value[0], ADD_VALUES);
367  }
369 }
370 //! [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
Deprecated interface functions.
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< VectorDouble > rhoAtIntegrationPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[SecondOp]
static char help[]
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: Basic.hpp:466
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: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
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:90
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]
Basic algebra on fields.
Definition: FieldBlas.hpp:34
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]
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
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)
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
DataForcesAndSourcesCore::EntData EntData
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
FTensor::Index< 'i', 3 > i
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