v0.13.2
Loading...
Searching...
No Matches
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#include <MoFEM.hpp>
13
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
21
22//! [Example]
23struct Example {
24
25 Example(MoFEM::Interface &m_field) : mField(m_field) {}
26
28
29private:
31
39
40 struct CommonData;
41
42 boost::shared_ptr<CommonData> commonDataPtr;
43
44 struct OpZero;
45
46 struct OpFirst;
47
48 struct OpSecond;
49};
50//! [Example]
51
52//! [Common data]
54 : public boost::enable_shared_from_this<Example::CommonData> {
55
56 VectorDouble rhoAtIntegrationPts; ///< Storing density at integration points
57
58 inline boost::shared_ptr<VectorDouble> getRhoAtIntegrationPtsPtr() {
59 return boost::shared_ptr<VectorDouble>(shared_from_this(),
61 }
62
63 /**
64 * @brief Vector to indicate indices for storing zero, first and second
65 * moments of inertia.
66 *
67 */
69 ZERO = 0,
80 };
81
83 petscVec; ///< Smart pointer which stores PETSc distributed vector
84};
85//! [Common data]
86
87struct Example::OpZero : public OpElement {
88 OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
89 : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
90 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
91 }
92 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
93
94private:
95 boost::shared_ptr<CommonData> commonDataPtr;
96};
97
98struct Example::OpFirst : public OpElement {
99 OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
100 : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
101 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
102 }
103 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
104
105private:
106 boost::shared_ptr<CommonData> commonDataPtr;
107};
108
109//! [Operator]
111 OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
112 : OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
113 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
114 }
115 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
116
117private:
118 boost::shared_ptr<CommonData> commonDataPtr;
119};
120//! [Operator]
121
122//! [Run all]
125 CHKERR setUp();
133}
134//! [Run all]
135
136//! [Set up problem]
140 CHKERR simple->getOptions();
141 CHKERR simple->loadFile();
142 // Add field
143 CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
144 constexpr int order = 1;
145 CHKERR simple->setFieldOrder("rho", order);
146 CHKERR simple->setUp();
148}
149//! [Set up problem]
150
151//! [Create common data]
154 commonDataPtr = boost::make_shared<CommonData>();
155
156 int local_size;
157 if (mField.get_comm_rank() == 0) // get_comm_rank() gets processor number
158 // processor 0
159 local_size = CommonData::LAST_ELEMENT; // last element gives size of vector
160 else
161 // other processors (e.g. 1, 2, 3, etc.)
162 local_size = 0; // local size of vector is zero on other processors
163
164 commonDataPtr->petscVec = createSmartVectorMPI(mField.get_comm(), local_size,
166
168}
169//! [Create common data]
170
171//! [Set density distribution]
174 auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
175 double *ycoord, double *zcoord) {
177 field_data[0] = 1;
179 };
180 FieldBlas *field_blas;
181 CHKERR mField.getInterface(field_blas);
182 CHKERR field_blas->setVertexDofs(set_density, "rho");
184}
185//! [Set density distribution]
186
187//! [Push operators to pipeline]
191
192 // Push an operator which calculates values of density at integration points
193 pipeline_mng->getOpDomainRhsPipeline().push_back(
195 "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
196
197 // Push an operator to pipeline to calculate zero moment of inertia (mass)
198 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
199
200 // Push an operator to the pipeline to calculate first moment of inertaia
201 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
202
203 // Push an operator to the pipeline to calculate second moment of inertaia
204 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
205
206 // Set integration rule. Integration rule is equal to the polynomial order of
207 // the density field plus 2, since under the integral of the second moment of
208 // inertia term x*x is present
209 auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
210
211 // Add integration rule to the element
214}
215//! [Push operators to pipeline]
216
217//! [Integrate]
220 // Zero global vector
221 CHKERR VecZeroEntries(commonDataPtr->petscVec);
222
223 // Integrate elements by executing operators in the pipeline
225 CHKERR pipeline_mng->loopFiniteElements();
226
227 // Assemble MPI vector
228 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
229 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
231}
232//! [Integrate]
233
234//! [Print results]
237 const double *array;
238 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
239 if (mField.get_comm_rank() == 0) {
240 MOFEM_LOG_C("SELF", Sev::inform, "Mass %6.4e", array[CommonData::ZERO]);
241 MOFEM_LOG_C("SELF", Sev::inform,
242 "First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
244 array[CommonData::FIRST_Z]);
245 MOFEM_LOG_C("SELF", Sev::inform,
246 "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
247 "%6.4e ]",
251 }
252 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
254}
255//! [Print results]
256
257//! [Test example]
260 const double *array;
261 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
262
263 PetscBool test = PETSC_FALSE;
264 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
265 if (mField.get_comm_rank() == 0 && test) {
266 constexpr double eps = 1e-8;
267 constexpr double expected_volume = 1.;
268 constexpr double expected_first_moment = 0.;
269 constexpr double expected_second_moment = 1. / 12.;
270 if (std::abs(array[CommonData::ZERO] - expected_volume) > eps)
271 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
272 "Wrong area %6.4e != %6.4e", expected_volume,
273 array[CommonData::ZERO]);
274 for (auto i :
276 if (std::abs(array[i] - expected_first_moment) > eps)
277 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
278 "Wrong first moment %6.4e != %6.4e", expected_first_moment,
279 array[i]);
280 }
283 if (std::abs(array[i] - expected_second_moment) > eps)
284 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
285 "Wrong second moment %6.4e != %6.4e", expected_second_moment,
286 array[i]);
287 }
288 }
289 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
290
292}
293//! [Test example]
294
295//! [main]
296int main(int argc, char *argv[]) {
297
298 // Initialisation of MoFEM/PETSc and MOAB data structures
299 const char param_file[] = "param_file.petsc";
301
302 // Error handling
303 try {
304
305 //! [Register MoFEM discrete manager in PETSc]
306 DMType dm_name = "DMMOFEM";
307 CHKERR DMRegister_MoFEM(dm_name);
308 //! [Register MoFEM discrete manager in PETSc]
309
310 //! [Create MoAB]
311 moab::Core mb_instance; ///< mesh database
312 moab::Interface &moab = mb_instance; ///< mesh database interface
313 //! [Create MoAB]
314
315 //! [Create MoFEM]
316 MoFEM::Core core(moab); ///< finite element database
317 MoFEM::Interface &m_field = core; ///< finite element database interface
318 //! [Create MoFEM]
319
320 Example ex(m_field);
321 CHKERR ex.runProblem();
322 }
324
326}
327//! [main]
328
329//! [ZeroOp]
331 EntData &data) {
333 if (type == MBVERTEX) {
334
335 const int nb_integration_pts =
336 getGaussPts().size2(); // Number of integration points
337 auto t_w = getFTensor0IntegrationWeight(); // Integration weights
338 auto t_rho = getFTensor0FromVec(
339 commonDataPtr->rhoAtIntegrationPts); // Density at integration weights
340 const double volume = getMeasure();
341
342 // Integrate area of the element
343 double element_local_value = 0;
344 for (int gg = 0; gg != nb_integration_pts; ++gg) {
345 element_local_value += t_w * t_rho * volume;
346 ++t_w;
347 ++t_rho;
348 }
349
350 // Assemble area of the element into global PETSc vector
351 const int index = CommonData::ZERO;
352 CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
353 ADD_VALUES);
354 }
356}
357//! [ZeroOp]
358
359//! [FirstOp]
361 EntData &data) {
363 const int nb_integration_pts = getGaussPts().size2();
364 auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
365 auto t_rho = getFTensor0FromVec(
366 commonDataPtr->rhoAtIntegrationPts); ///< Density at integration points
367 auto t_x =
368 getFTensor1CoordsAtGaussPts(); ///< Coordinates at integration points
369 const double volume = getMeasure(); ///< Get Volume of element
370
372
374 t_s; ///< First moment of inertia is tensor of rank 1, i.e. vector.
375 t_s(i) = 0; // Zero entries
376
377 // Integrate
378 for (int gg = 0; gg != nb_integration_pts; ++gg) {
379
380 t_s(i) += t_w * t_rho * volume * t_x(i);
381
382 ++t_w; // move weight to next integration pts
383 ++t_rho; // move density
384 ++t_x; // move coordinate
385 }
386
387 // Set array of indices
388 constexpr std::array<int, 3> indices = {
390
391 // Assemble first moment of inertia
392 CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(), &t_s(0),
393 ADD_VALUES);
395}
396//! [FirstOp]
397
398//! [SecondOp]
400 EntData &data) {
402
403 const int nb_integration_pts = getGaussPts().size2();
404 auto t_w = getFTensor0IntegrationWeight();
405 auto t_rho = getFTensor0FromVec(commonDataPtr->rhoAtIntegrationPts);
406 auto t_x = getFTensor1CoordsAtGaussPts();
407 const double volume = getMeasure();
408
409 // Create storage for a symmetric tensor
410 std::array<double, 6> element_local_value;
411 std::fill(element_local_value.begin(), element_local_value.end(), 0.0);
412
413 // Create symmetric tensor using pointers to the storage
415 &element_local_value[0], &element_local_value[1], &element_local_value[2],
416 &element_local_value[3], &element_local_value[4],
417 &element_local_value[5]);
418
421
422 // Integrate
423 for (int gg = 0; gg != nb_integration_pts; ++gg) {
424
425 // Symbol "^" indicates outer product which yields a symmetric tensor
426 t_I(i, j) += (t_w * t_rho * volume) * (t_x(i) ^ t_x(j));
427
428 ++t_w; // move weight pointer to the next integration point
429 ++t_rho; // move density pointer
430 ++t_x; // move coordinate pointer
431 }
432
433 // Set array of indices
434 constexpr std::array<int, 6> indices = {
437
438 // Assemble second moment of inertia
439 CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
440 &element_local_value[0], ADD_VALUES);
441
443}
444//! [SecondOp]
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
static char help[]
Definition: integration.cpp:16
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
boost::shared_ptr< VectorDouble > getRhoAtIntegrationPtsPtr()
Definition: integration.cpp:58
VecElements
Vector to indicate indices for storing zero, first and second moments of inertia.
Definition: integration.cpp:68
SmartPetscObj< Vec > petscVec
Smart pointer which stores PETSc distributed vector.
Definition: integration.cpp:83
VectorDouble rhoAtIntegrationPts
Storing density at integration points.
Definition: integration.cpp:56
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[ZeroOp]
OpFirst(boost::shared_ptr< CommonData > &common_data_ptr)
Definition: integration.cpp:99
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[FirstOp]
boost::shared_ptr< CommonData > commonDataPtr
OpSecond(boost::shared_ptr< CommonData > &common_data_ptr)
[Common data]
Definition: integration.cpp:87
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[main]
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:95
OpZero(boost::shared_ptr< CommonData > &common_data_ptr)
Definition: integration.cpp:88
[Example]
Definition: plastic.cpp:141
MoFEMErrorCode pushOperators()
[Set density distribution]
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:150
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:42
MoFEMErrorCode setUp()
[Run all]
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: integration.cpp:25
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:148
MoFEMErrorCode postProcess()
MoFEMErrorCode setFieldValues()
[Create common data]
MoFEMErrorCode integrateElements()
[Push operators to pipeline]
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:320
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get value at integration points for scalar field.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.