* \file OpPostProcElastic.hpp
* \example OpPostProcElastic.hpp
* Postprocessing
namespace Tutorial {
//! [Class definition]
template <int DIM> struct OpPostProcElastic : public DomainEleOp {
OpPostProcElastic(const std::string field_name,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
boost::shared_ptr<MatrixDouble> m_strain_ptr,
boost::shared_ptr<MatrixDouble> m_stress_ptr);
std::vector<EntityHandle> &mapGaussPts;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
//! [Class definition]
//! [Postprocessing constructor]
template <int DIM>
const std::string field_name, moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
boost::shared_ptr<MatrixDouble> m_strain_ptr,
boost::shared_ptr<MatrixDouble> m_stress_ptr)
: DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts), mStrainPtr(m_strain_ptr),
mStressPtr(m_stress_ptr) {
// Opetor is only executed for vertices
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
//! [Postprocessing constructor]
//! [Postprocessing]
template <int DIM>
EntData &data) {
auto get_tag = [&](const std::string name) {
std::array<double, 9> def;
std::fill(def.begin(), def.end(), 0);
Tag th;
CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
return th;
MatrixDouble3by3 mat(3, 3);
auto set_float_precision = [](const double x) {
if (std::abs(x) < std::numeric_limits<float>::epsilon())
return 0.;
return x;
auto set_matrix_symm = [&](auto &t) -> MatrixDouble3by3 & {
for (size_t r = 0; r != DIM; ++r)
for (size_t c = 0; c != DIM; ++c)
mat(r, c) = t(r, c);
return mat;
auto set_plain_stress_strain = [&](auto &mat, auto &t) -> MatrixDouble3by3 & {
mat(2, 2) = -poisson_ratio * (t(0, 0) + t(1, 1));
return mat;
auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
for(auto &v : mat.data())
v = set_float_precision(v);
return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
auto th_strain = get_tag("STRAIN");
auto th_stress = get_tag("STRESS");
size_t nb_gauss_pts = data.getN().size1();
auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(mStrainPtr));
auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(mStressPtr));
switch (DIM) {
case 2:
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
CHKERR set_tag(
th_strain, gg,
set_plain_stress_strain(set_matrix_symm(t_strain), t_stress));
CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
case 3:
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
CHKERR set_tag(th_strain, gg, set_matrix_symm(t_strain));
CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
"Not implemeneted dimension");
//! [Postprocessing]
} // namespace Tutorial
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
Definition: definitions.h:32
#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
double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double r
rate factor
constexpr double t
plate stiffness
Definition: plate.cpp:62
constexpr auto field_name
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
operator doWork function is executed on FE rows
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing constructor]
boost::shared_ptr< MatrixDouble > mStressPtr
moab::Interface & postProcMesh
boost::shared_ptr< MatrixDouble > mStrainPtr
OpPostProcElastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< MatrixDouble > m_strain_ptr, boost::shared_ptr< MatrixDouble > m_stress_ptr)
[Class definition]
std::vector< EntityHandle > & mapGaussPts