v0.13.1
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
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__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
DEPRECATED typedef PostProcStress PostPorcStress
EntitiesFieldData::EntData EntData
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define DEPRECATED
Definition: definitions.h:17
#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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1210
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
const double r
rate factor
double H
Definition: plastic.cpp:100
constexpr auto field_name
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPROW
operator doWork function is executed on FE rows
data for calculation heat conductivity and heat capacity elements
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr
Range tEts
constrains elements in block set
common data used by volume elements
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Range tEts
std::map< std::string, std::vector< VectorDouble > > fieldMap
std::map< std::string, std::vector< MatrixDouble > > gradMap
PostProcCommonOnRefMesh::CommonDataForVolume & commonData
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)
NonlinearElasticElement::BlockData & dAta
const bool printCauchy
const bool replaceNonANumberByMaxValue
std::vector< EntityHandle > & mapGaussPts
const double maxVal
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
const bool fieldDisp
moab::Interface & postProcMesh
NonlinearElasticElement::CommonData nonLinearElementCommonData