v0.10.0
OpPostProcElastic.hpp
Go to the documentation of this file.
1 /**
2  * \file OpPostProcElastic.hpp
3  * \example OpPostProcElastic.hpp
4  *
5  * Postprocessing
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 namespace Tutorial {
24 
25 //! [Class definition]
26 template <int DIM> struct OpPostProcElastic : public DomainEleOp {
27  OpPostProcElastic(const std::string field_name,
28  moab::Interface &post_proc_mesh,
29  std::vector<EntityHandle> &map_gauss_pts,
30  boost::shared_ptr<MatrixDouble> m_strain_ptr,
31  boost::shared_ptr<MatrixDouble> m_stress_ptr);
32  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
33 
34 private:
36  std::vector<EntityHandle> &mapGaussPts;
37  boost::shared_ptr<MatrixDouble> mStrainPtr;
38  boost::shared_ptr<MatrixDouble> mStressPtr;
39 };
40 //! [Class definition]
41 
42 //! [Postprocessing constructor]
43 template <int DIM>
45  const std::string field_name, moab::Interface &post_proc_mesh,
46  std::vector<EntityHandle> &map_gauss_pts,
47  boost::shared_ptr<MatrixDouble> m_strain_ptr,
48  boost::shared_ptr<MatrixDouble> m_stress_ptr)
49  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
50  mapGaussPts(map_gauss_pts), mStrainPtr(m_strain_ptr),
51  mStressPtr(m_stress_ptr) {
52  // Opetor is only executed for vertices
53  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
54 }
55 //! [Postprocessing constructor]
56 
57 //! [Postprocessing]
58 template <int DIM>
60  EntData &data) {
62 
63  auto get_tag = [&](const std::string name) {
64  std::array<double, 9> def;
65  std::fill(def.begin(), def.end(), 0);
66  Tag th;
67  CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
68  MB_TAG_CREAT | MB_TAG_SPARSE,
69  def.data());
70  return th;
71  };
72 
73  MatrixDouble3by3 mat(3, 3);
74 
75  auto set_matrix_symm = [&](auto &t) -> MatrixDouble3by3 & {
76  mat.clear();
77  for (size_t r = 0; r != DIM; ++r)
78  for (size_t c = 0; c != DIM; ++c)
79  mat(r, c) = t(r, c);
80  return mat;
81  };
82 
83  auto set_plain_stress_strain = [&](auto &mat, auto &t) -> MatrixDouble3by3 & {
84  mat(2, 2) = -poisson_ratio * (t(0, 0) + t(1, 1));
85  return mat;
86  };
87 
88  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
89  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
90  &*mat.data().begin());
91  };
92 
93  auto th_strain = get_tag("STRAIN");
94  auto th_stress = get_tag("STRESS");
95 
96  size_t nb_gauss_pts = data.getN().size1();
97  auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(mStrainPtr));
98  auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(mStressPtr));
99 
100  switch (DIM) {
101  case 2:
102  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
103  CHKERR set_tag(
104  th_strain, gg,
105  set_plain_stress_strain(set_matrix_symm(t_strain), t_stress));
106  CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
107  ++t_strain;
108  ++t_stress;
109  }
110  break;
111  case 3:
112  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
113  CHKERR set_tag(th_strain, gg, set_matrix_symm(t_strain));
114  CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
115  ++t_strain;
116  ++t_stress;
117  }
118  break;
119  default:
120  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
121  "Not implemeneted dimension");
122  }
123 
125 }
126 //! [Postprocessing]
127 
128 } // namespace Tutorial
poisson_ratio
constexpr double poisson_ratio
Definition: plastic.cpp:101
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
Tutorial::OpPostProcElastic::mStressPtr
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: OpPostProcElastic.hpp:38
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
Tutorial
Definition: OpPostProcElastic.hpp:23
Tutorial::OpPostProcElastic::OpPostProcElastic
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]
Definition: OpPostProcElastic.hpp:44
convert.type
type
Definition: convert.py:66
Tutorial::OpPostProcElastic::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: OpPostProcElastic.hpp:36
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
Tutorial::OpPostProcElastic
[Class definition]
Definition: OpPostProcElastic.hpp:26
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1288
DomainEleOp
ElectroPhysiology::c
const double c
Definition: ElecPhysOperators.hpp:31
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:60
Tutorial::OpPostProcElastic::postProcMesh
moab::Interface & postProcMesh
Definition: OpPostProcElastic.hpp:35
Tutorial::OpPostProcElastic::mStrainPtr
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: OpPostProcElastic.hpp:37
Tutorial::OpPostProcElastic::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing constructor]
Definition: OpPostProcElastic.hpp:59
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33