v0.15.0
Loading...
Searching...
No Matches
mofem/tutorials/fun-1/integration.cpp

Calculate zero, first and second moments of inertia.

Calculate zero, first and second moments of inertia.Example intended to show how to write user data operator and integrate scalar field.

/**
* \file lesson1_moment_of_inertia.cpp
* \example mofem/tutorials/fun-1/integration.cpp
*
* \brief Calculate zero, first and second moments of inertia.
*
* Example intended to show how to write user data operator and integrate
* scalar field.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
//! [Example]
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
struct CommonData;
boost::shared_ptr<CommonData> commonDataPtr;
struct OpZero;
struct OpFirst;
struct OpSecond;
};
//! [Example]
//! [Common data]
: public boost::enable_shared_from_this<Example::CommonData> {
VectorDouble rhoAtIntegrationPts; ///< Storing density at integration points
inline boost::shared_ptr<VectorDouble> getRhoAtIntegrationPtsPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(),
}
/**
* @brief Vector to indicate indices for storing zero, first and second
* moments of inertia.
*
*/
enum VecElements {
ZERO = 0,
};
petscVec; ///< Smart pointer which stores PETSc distributed vector
};
//! [Common data]
struct Example::OpZero : public OpElement {
OpZero(boost::shared_ptr<CommonData> &common_data_ptr)
: OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct Example::OpFirst : public OpElement {
OpFirst(boost::shared_ptr<CommonData> &common_data_ptr)
: OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
//! [Operator]
struct Example::OpSecond : public OpElement {
OpSecond(boost::shared_ptr<CommonData> &common_data_ptr)
: OpElement("rho", OPROW), commonDataPtr(common_data_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
//! [Operator]
//! [Run all]
}
//! [Run all]
//! [Set up problem]
CHKERR simple->loadFile();
// Add field
CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
constexpr int order = 1;
CHKERR simple->setFieldOrder("rho", order);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Create common data]
commonDataPtr = boost::make_shared<CommonData>();
int local_size;
if (mField.get_comm_rank() == 0) // get_comm_rank() gets processor number
// processor 0
local_size = CommonData::LAST_ELEMENT; // last element gives size of vector
else
// other processors (e.g. 1, 2, 3, etc.)
local_size = 0; // local size of vector is zero on other processors
commonDataPtr->petscVec = createVectorMPI(mField.get_comm(), local_size,
}
//! [Create common data]
//! [Set density distribution]
auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
double *ycoord, double *zcoord) {
field_data[0] = 1;
};
FieldBlas *field_blas;
CHKERR mField.getInterface(field_blas);
CHKERR field_blas->setVertexDofs(set_density, "rho");
}
//! [Set density distribution]
//! [Push operators to pipeline]
// Push an operator which calculates values of density at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
"rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
// Push an operator to pipeline to calculate zero moment of inertia (mass)
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
// Push an operator to the pipeline to calculate first moment of inertaia
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
// Push an operator to the pipeline to calculate second moment of inertaia
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
// Set integration rule. Integration rule is equal to the polynomial order of
// the density field plus 2, since under the integral of the second moment of
// inertia term x*x is present
auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
// Add integration rule to the element
}
//! [Push operators to pipeline]
//! [Integrate]
// Zero global vector
CHKERR VecZeroEntries(commonDataPtr->petscVec);
// Integrate elements by executing operators in the pipeline
CHKERR pipeline_mng->loopFiniteElements();
// Assemble MPI vector
CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
}
//! [Integrate]
//! [Print results]
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
if (mField.get_comm_rank() == 0) {
MOFEM_LOG_C("SELF", Sev::inform, "Mass %6.4e", array[CommonData::ZERO]);
MOFEM_LOG_C("SELF", Sev::inform,
"First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
MOFEM_LOG_C("SELF", Sev::inform,
"Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
"%6.4e ]",
}
CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
}
//! [Print results]
//! [Test example]
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
PetscBool test = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test, PETSC_NULLPTR);
if (mField.get_comm_rank() == 0 && test) {
constexpr double eps = 1e-8;
constexpr double expected_volume = 1.;
constexpr double expected_first_moment = 0.;
constexpr double expected_second_moment = 1. / 12.;
if (std::abs(array[CommonData::ZERO] - expected_volume) > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong area %6.4e != %6.4e", expected_volume,
for (auto i :
if (std::abs(array[i] - expected_first_moment) > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong first moment %6.4e != %6.4e", expected_first_moment,
array[i]);
}
if (std::abs(array[i] - expected_second_moment) > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong second moment %6.4e != %6.4e", expected_second_moment,
array[i]);
}
}
CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
}
//! [Test example]
//! [main]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Error handling
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc]
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
Example ex(m_field);
CHKERR ex.runProblem();
}
}
//! [main]
//! [ZeroOp]
MoFEMErrorCode Example::OpZero::doWork(int side, EntityType type,
EntData &data) {
if (type == MBVERTEX) {
const int nb_integration_pts =
getGaussPts().size2(); // Number of integration points
auto t_w = getFTensor0IntegrationWeight(); // Integration weights
auto t_rho = getFTensor0FromVec(
commonDataPtr->rhoAtIntegrationPts); // Density at integration weights
const double volume = getMeasure();
// Integrate area of the element
double element_local_value = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
element_local_value += t_w * t_rho * volume;
++t_w;
++t_rho;
}
// Assemble area of the element into global PETSc vector
const int index = CommonData::ZERO;
CHKERR VecSetValue(commonDataPtr->petscVec, index, element_local_value,
ADD_VALUES);
}
}
//! [ZeroOp]
//! [FirstOp]
MoFEMErrorCode Example::OpFirst::doWork(int side, EntityType type,
EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
auto t_rho = getFTensor0FromVec(
commonDataPtr->rhoAtIntegrationPts); ///< Density at integration points
auto t_x =
getFTensor1CoordsAtGaussPts(); ///< Coordinates at integration points
const double volume = getMeasure(); ///< Get Volume of element
FTensor::Index<'i', 3> i;
t_s; ///< First moment of inertia is tensor of rank 1, i.e. vector.
t_s(i) = 0; // Zero entries
// Integrate
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_s(i) += t_w * t_rho * volume * t_x(i);
++t_w; // move weight to next integration pts
++t_rho; // move density
++t_x; // move coordinate
}
// Set array of indices
constexpr std::array<int, 3> indices = {
// Assemble first moment of inertia
CHKERR VecSetValues(commonDataPtr->petscVec, 3, indices.data(), &t_s(0),
ADD_VALUES);
}
//! [FirstOp]
//! [SecondOp]
MoFEMErrorCode Example::OpSecond::doWork(int side, EntityType type,
EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_rho = getFTensor0FromVec(commonDataPtr->rhoAtIntegrationPts);
auto t_x = getFTensor1CoordsAtGaussPts();
const double volume = getMeasure();
// Create storage for a symmetric tensor
std::array<double, 6> element_local_value;
std::fill(element_local_value.begin(), element_local_value.end(), 0.0);
// Create symmetric tensor using pointers to the storage
&element_local_value[0], &element_local_value[1], &element_local_value[2],
&element_local_value[3], &element_local_value[4],
&element_local_value[5]);
FTensor::Index<'i', 3> i;
FTensor::Index<'j', 3> j;
// Integrate
for (int gg = 0; gg != nb_integration_pts; ++gg) {
// Symbol "^" indicates outer product which yields a symmetric tensor
t_I(i, j) += (t_w * t_rho * volume) * (t_x(i) ^ t_x(j));
++t_w; // move weight pointer to the next integration point
++t_rho; // move density pointer
++t_x; // move coordinate pointer
}
// Set array of indices
constexpr std::array<int, 6> indices = {
// Assemble second moment of inertia
CHKERR VecSetValues(commonDataPtr->petscVec, 6, indices.data(),
&element_local_value[0], ADD_VALUES);
}
//! [SecondOp]
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
static char help[]
int main()
static const double eps
#define CATCH_ERRORS
Catch errors.
@ 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()
@ 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 ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
UBlasVector< double > VectorDouble
Definition Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< VectorDouble > getRhoAtIntegrationPtsPtr()
VecElements
Vector to indicate indices for storing zero, first and second moments of inertia.
SmartPetscObj< Vec > petscVec
Smart pointer which stores PETSc distributed vector.
VectorDouble rhoAtIntegrationPts
Storing density at integration points.
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]
Definition plastic.cpp:216
MoFEMErrorCode pushOperators()
[Set density distribution]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode setUp()
[Run all]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Set up problem]
Definition plastic.cpp:480
MoFEMErrorCode runProblem()
[Run problem]
Definition plastic.cpp:256
MoFEM::Interface & mField
Definition plastic.cpp:226
MoFEMErrorCode postProcess()
[Integrate]
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:118
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.
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
Specialization for double precision scalar field values calculation.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.