v0.13.1
OpPostProcElastic.hpp
Go to the documentation of this file.
1/**
2 * \file OpPostProcElastic.hpp
3 * \example OpPostProcElastic.hpp
4 *
5 * Postprocessing
6 *
7 */
8
9namespace Tutorial {
10
11//! [Class definition]
12template <int DIM> struct OpPostProcElastic : public DomainEleOp {
13 OpPostProcElastic(const std::string field_name,
14 moab::Interface &post_proc_mesh,
15 std::vector<EntityHandle> &map_gauss_pts,
16 boost::shared_ptr<MatrixDouble> m_strain_ptr,
17 boost::shared_ptr<MatrixDouble> m_stress_ptr);
19
20private:
22 std::vector<EntityHandle> &mapGaussPts;
23 boost::shared_ptr<MatrixDouble> mStrainPtr;
24 boost::shared_ptr<MatrixDouble> mStressPtr;
25};
26//! [Class definition]
27
28//! [Postprocessing constructor]
29template <int DIM>
31 const std::string field_name, moab::Interface &post_proc_mesh,
32 std::vector<EntityHandle> &map_gauss_pts,
33 boost::shared_ptr<MatrixDouble> m_strain_ptr,
34 boost::shared_ptr<MatrixDouble> m_stress_ptr)
35 : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
36 mapGaussPts(map_gauss_pts), mStrainPtr(m_strain_ptr),
37 mStressPtr(m_stress_ptr) {
38 // Opetor is only executed for vertices
39 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
40}
41//! [Postprocessing constructor]
42
43//! [Postprocessing]
44template <int DIM>
46 EntData &data) {
48
49 auto get_tag = [&](const std::string name) {
50 std::array<double, 9> def;
51 std::fill(def.begin(), def.end(), 0);
52 Tag th;
53 CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
54 MB_TAG_CREAT | MB_TAG_SPARSE,
55 def.data());
56 return th;
57 };
58
59 MatrixDouble3by3 mat(3, 3);
60
61 auto set_float_precision = [](const double x) {
62 if (std::abs(x) < std::numeric_limits<float>::epsilon())
63 return 0.;
64 else
65 return x;
66 };
67
68 auto set_matrix_symm = [&](auto &t) -> MatrixDouble3by3 & {
69 mat.clear();
70 for (size_t r = 0; r != DIM; ++r)
71 for (size_t c = 0; c != DIM; ++c)
72 mat(r, c) = t(r, c);
73 return mat;
74 };
75
76 auto set_plain_stress_strain = [&](auto &mat, auto &t) -> MatrixDouble3by3 & {
77 mat(2, 2) = -poisson_ratio * (t(0, 0) + t(1, 1));
78 return mat;
79 };
80
81 auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
82 for(auto &v : mat.data())
83 v = set_float_precision(v);
84 return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
85 &*mat.data().begin());
86 };
87
88 auto th_strain = get_tag("STRAIN");
89 auto th_stress = get_tag("STRESS");
90
91 size_t nb_gauss_pts = data.getN().size1();
92 auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(mStrainPtr));
93 auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(mStressPtr));
94
95 switch (DIM) {
96 case 2:
97 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
98 CHKERR set_tag(
99 th_strain, gg,
100 set_plain_stress_strain(set_matrix_symm(t_strain), t_stress));
101 CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
102 ++t_strain;
103 ++t_stress;
104 }
105 break;
106 case 3:
107 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
108 CHKERR set_tag(th_strain, gg, set_matrix_symm(t_strain));
109 CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
110 ++t_strain;
111 ++t_stress;
112 }
113 break;
114 default:
115 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
116 "Not implemeneted dimension");
117 }
118
120}
121//! [Postprocessing]
122
123} // namespace Tutorial
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_NOT_IMPLEMENTED
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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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....
[Class definition]
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