v0.14.0
PostProcStresses.hpp
Go to the documentation of this file.
1 /** \file PostProcStresses.hpp
2  * \brief Post-processing stresses for non-linear analysis
3  * \ingroup nonlinear_elastic_elem
4  *
5  * Implementation of method for post-processing stresses.
6  */
7 
8 
9 
10 #ifndef __POSTPROCSTRESSES_HPP__
11 #define __POSTPROCSTRESSES_HPP__
12 
13 #ifndef WITH_ADOL_C
14 #error "MoFEM need to be compiled with ADOL-C"
15 #endif
16 
19 
21  std::vector<EntityHandle> &mapGaussPts;
22 
25  const bool fieldDisp;
27  const double maxVal;
28  const bool printCauchy;
29 
31  std::vector<EntityHandle> &map_gauss_pts,
32  const std::string field_name,
35  const bool field_disp = false,
36  const bool replace_nonanumber_by_max_value = false,
37  const double max_val = 1e16,
38  const bool print_cauchy_stress = false)
39  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
41  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), dAta(data),
42  commonData(common_data), fieldDisp(field_disp),
43  replaceNonANumberByMaxValue(replace_nonanumber_by_max_value),
44  maxVal(max_val), printCauchy(print_cauchy_stress) {}
45 
47 
48  MoFEMErrorCode doWork(int side, EntityType type,
51 
52  if (type != MBVERTEX)
54  if (data.getIndices().size() == 0)
56  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
57  dAta.tEts.end()) {
59  }
60 
61  const auto &dof_ptr = data.getFieldDofs()[0];
62 
63  int id = dAta.iD;
64 
65  Tag th_id;
66  int def_block_id = -1;
67  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
68  MB_TAG_CREAT | MB_TAG_SPARSE,
69  &def_block_id);
70  Range::iterator tit = commonData.tEts.begin();
71  for (; tit != commonData.tEts.end(); tit++) {
72  CHKERR postProcMesh.tag_set_data(th_id, &*tit, 1, &id);
73  }
74 
75  string tag_name_piola1 = dof_ptr->getName() + "_PIOLA1_STRESS";
76  string tag_name_energy = dof_ptr->getName() + "_ENERGY_DENSITY";
77 
78  int tag_length = 9;
79  double def_VAL[tag_length];
80  bzero(def_VAL, tag_length * sizeof(double));
81  Tag th_piola1, th_energy, th_cauchy;
82  CHKERR postProcMesh.tag_get_handle(tag_name_piola1.c_str(), tag_length,
83  MB_TYPE_DOUBLE, th_piola1,
84  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
85  CHKERR postProcMesh.tag_get_handle(tag_name_energy.c_str(), 1,
86  MB_TYPE_DOUBLE, th_energy,
87  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
88 
89  if (printCauchy) {
90  string tag_name_cauchy = "MED_" + dof_ptr->getName() + "_CAUCHY_STRESS";
91  CHKERR postProcMesh.tag_get_handle(tag_name_cauchy.c_str(), tag_length,
92  MB_TYPE_DOUBLE, th_cauchy,
93  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
94  }
95 
96  int nb_gauss_pts = data.getN().size1();
97  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
98  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
99  "Nb. of integration points is not equal to number points on "
100  "post-processing mesh");
101  }
102  if (commonData.gradMap[rowFieldName].size() != (unsigned int)nb_gauss_pts) {
103  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
104  "Gradient of field not found, filed <%s> not found",
105  rowFieldName.c_str());
106  }
107 
108  MatrixDouble3by3 H, invH;
109  double detH;
110 
112  dAta.materialDoublePtr->opPtr = this;
113  CHKERR dAta.materialDoublePtr->getDataOnPostProcessor(commonData.fieldMap,
115 
118 
119  MatrixDouble3by3 maxP(3, 3);
120  maxP.clear();
121 
122  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
123 
124  dAta.materialDoublePtr->gG = gg;
125  dAta.materialDoublePtr->F.resize(3, 3);
126  noalias(dAta.materialDoublePtr->F) =
128  if (fieldDisp) {
129  for (int dd = 0; dd != 3; dd++) {
130  dAta.materialDoublePtr->F(dd, dd) += 1;
131  }
132  }
133  if (commonData.gradMap["MESH_NODE_POSITIONS"].size() ==
134  (unsigned int)nb_gauss_pts) {
135  H.resize(3, 3);
136  invH.resize(3, 3);
137  noalias(H) = (commonData.gradMap["MESH_NODE_POSITIONS"])[gg];
138  detH = determinantTensor3by3(H);
139  CHKERR invertTensor3by3(H, detH, invH);
140  noalias(dAta.materialDoublePtr->F) =
141  prod(dAta.materialDoublePtr->F, invH);
142  }
143 
144  int nb_active_variables = 9;
145  CHKERR dAta.materialDoublePtr->setUserActiveVariables(
146  nb_active_variables);
147  CHKERR dAta.materialDoublePtr->calculateP_PiolaKirchhoffI(
149  CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
151  CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
152  &dAta.materialDoublePtr->P(0, 0));
153  CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
154  &dAta.materialDoublePtr->eNergy);
155  if (printCauchy) {
156  dAta.materialDoublePtr->sigmaCauchy.resize(3, 3);
157  CHKERR dAta.materialDoublePtr->calculateCauchyStress(
159  CHKERR postProcMesh.tag_set_data(
160  th_cauchy, &mapGaussPts[gg], 1,
161  &dAta.materialDoublePtr->sigmaCauchy(0, 0));
162  }
163  }
164 
166  MatrixDouble3by3 P(3, 3);
167  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
168  double val_energy;
169  CHKERR postProcMesh.tag_get_data(th_energy, &mapGaussPts[gg], 1,
170  &val_energy);
171  if (!std::isnormal(val_energy)) {
172  CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
173  &maxVal);
174  CHKERR postProcMesh.tag_get_data(th_piola1, &mapGaussPts[gg], 1,
175  &P(0, 0));
176  for (unsigned int r = 0; r != P.size1(); ++r) {
177  for (unsigned int c = 0; c != P.size2(); ++c) {
178  if (!std::isnormal(P(r, c)))
179  P(r, c) = copysign(maxVal, P(r, c));
180  }
181  }
182  CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
183  &P(0, 0));
184  }
185  }
186  }
187 
189  }
190 };
191 
192 /// \deprecated Use PostProcStress
194 
195 #endif //__POSTPROCSTRESSES_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
PostProcCommonOnRefMesh::CommonData::fieldMap
std::map< std::string, std::vector< VectorDouble > > fieldMap
Definition: PostProcOnRefMesh.hpp:33
PostProcCommonOnRefMesh::CommonDataForVolume
Definition: PostProcOnRefMesh.hpp:37
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PostProcStress::commonData
PostProcCommonOnRefMesh::CommonDataForVolume & commonData
Definition: PostProcStresses.hpp:24
PostProcStress::PostProcStress
PostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, NonlinearElasticElement::BlockData &data, PostProcCommonOnRefMesh::CommonDataForVolume &common_data, const bool field_disp=false, const bool replace_nonanumber_by_max_value=false, const double max_val=1e16, const bool print_cauchy_stress=false)
Definition: PostProcStresses.hpp:30
MoFEM::ForcesAndSourcesCore::UserDataOperator::rowFieldName
std::string rowFieldName
Definition: ForcesAndSourcesCore.hpp:577
sdf.r
int r
Definition: sdf.py:8
PostProcStress::replaceNonANumberByMaxValue
const bool replaceNonANumberByMaxValue
Definition: PostProcStresses.hpp:26
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1566
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
PostProcStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: PostProcStresses.hpp:48
MoFEM::ForcesAndSourcesCore::UserDataOperator::ForcesAndSourcesCore
friend class ForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:989
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
PostProcStress::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: PostProcStresses.hpp:21
PostProcStress::dAta
NonlinearElasticElement::BlockData & dAta
Definition: PostProcStresses.hpp:23
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
PostProcStress
Definition: PostProcStresses.hpp:17
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
convert.type
type
Definition: convert.py:64
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
PostProcStress::nonLinearElementCommonData
NonlinearElasticElement::CommonData nonLinearElementCommonData
Definition: PostProcStresses.hpp:46
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:999
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
PostProcStress::fieldDisp
const bool fieldDisp
Definition: PostProcStresses.hpp:25
PostProcCommonOnRefMesh::CommonDataForVolume::tEts
Range tEts
Definition: PostProcOnRefMesh.hpp:38
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
PostProcStress::printCauchy
const bool printCauchy
Definition: PostProcStresses.hpp:28
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
NonlinearElasticElement::BlockData::iD
int iD
Definition: HookeElement.hpp:33
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
NonlinearElasticElement::BlockData::tEts
Range tEts
constrains elements in block set
Definition: HookeElement.hpp:36
PostProcStress::maxVal
const double maxVal
Definition: PostProcStresses.hpp:27
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
NonlinearElasticElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition: NonLinearElasticElement.hpp:107
PostProcCommonOnRefMesh::CommonData::gradMap
std::map< std::string, std::vector< MatrixDouble > > gradMap
Definition: PostProcOnRefMesh.hpp:34
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
NonlinearElasticElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition: NonLinearElasticElement.hpp:108
NonlinearElasticElement::BlockData::materialDoublePtr
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr
Definition: NonLinearElasticElement.hpp:94
H
double H
Hardening.
Definition: plastic.cpp:124
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
PostProcStress::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcStresses.hpp:20
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
NonlinearElasticElement::CommonData
common data used by volume elements
Definition: NonLinearElasticElement.hpp:105