v0.10.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  * 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 #ifndef __POSTPROCSTRESSES_HPP__
24 #define __POSTPROCSTRESSES_HPP__
25 
26 #ifndef WITH_ADOL_C
27 #error "MoFEM need to be compiled with ADOL-C"
28 #endif
29 
32 
34  std::vector<EntityHandle> &mapGaussPts;
35 
38  const bool fieldDisp;
40  const double maxVal;
41 
43  std::vector<EntityHandle> &map_gauss_pts,
44  const std::string field_name,
47  const bool field_disp = false,
48  const bool replace_nonanumber_by_max_value = false,
49  const double max_val = 1e16)
52  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), dAta(data),
53  commonData(common_data), fieldDisp(field_disp),
54  replaceNonANumberByMaxValue(replace_nonanumber_by_max_value),
55  maxVal(max_val) {}
56 
58 
59  MoFEMErrorCode doWork(int side, EntityType type,
62 
63  if (type != MBVERTEX)
65  if (data.getIndices().size() == 0)
67  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
68  dAta.tEts.end()) {
70  }
71 
72  const auto &dof_ptr = data.getFieldDofs()[0];
73 
74  int id = dAta.iD;
75 
76  Tag th_id;
77  int def_block_id = -1;
78  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
79  MB_TAG_CREAT | MB_TAG_SPARSE,
80  &def_block_id);
81  Range::iterator tit = commonData.tEts.begin();
82  for (; tit != commonData.tEts.end(); tit++) {
83  CHKERR postProcMesh.tag_set_data(th_id, &*tit, 1, &id);
84  }
85 
86  string tag_name_piola1 = dof_ptr->getName() + "_PIOLA1_STRESS";
87  string tag_name_energy = dof_ptr->getName() + "_ENERGY_DENSITY";
88 
89  int tag_length = 9;
90  double def_VAL[tag_length];
91  bzero(def_VAL, tag_length * sizeof(double));
92  Tag th_piola1, th_energy;
93  CHKERR postProcMesh.tag_get_handle(tag_name_piola1.c_str(), tag_length,
94  MB_TYPE_DOUBLE, th_piola1,
95  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
96  CHKERR postProcMesh.tag_get_handle(tag_name_energy.c_str(), 1,
97  MB_TYPE_DOUBLE, th_energy,
98  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
99 
100  int nb_gauss_pts = data.getN().size1();
101  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
102  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
103  "Nb. of integration points is not equal to number points on "
104  "post-processing mesh");
105  }
106  if (commonData.gradMap[rowFieldName].size() != (unsigned int)nb_gauss_pts) {
107  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
108  "Gradient of field not found, filed <%s> not found",
109  rowFieldName.c_str());
110  }
111 
112  MatrixDouble3by3 H, invH;
113  double detH;
114 
116  dAta.materialDoublePtr->opPtr = this;
117  CHKERR dAta.materialDoublePtr->getDataOnPostProcessor(commonData.fieldMap,
119 
122 
123  MatrixDouble3by3 maxP(3, 3);
124  maxP.clear();
125 
126  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
127 
128  dAta.materialDoublePtr->gG = gg;
129  dAta.materialDoublePtr->F.resize(3, 3);
130  noalias(dAta.materialDoublePtr->F) =
132  if (fieldDisp) {
133  for (int dd = 0; dd != 3; dd++) {
134  dAta.materialDoublePtr->F(dd, dd) += 1;
135  }
136  }
137  if (commonData.gradMap["MESH_NODE_POSITIONS"].size() ==
138  (unsigned int)nb_gauss_pts) {
139  H.resize(3, 3);
140  invH.resize(3, 3);
141  noalias(H) = (commonData.gradMap["MESH_NODE_POSITIONS"])[gg];
142  CHKERR dAta.materialDoublePtr->dEterminant(H, detH);
143  CHKERR dAta.materialDoublePtr->iNvert(detH, H, invH);
144  noalias(dAta.materialDoublePtr->F) =
145  prod(dAta.materialDoublePtr->F, invH);
146  }
147 
148  int nb_active_variables = 9;
149  CHKERR dAta.materialDoublePtr->setUserActiveVariables(
150  nb_active_variables);
151  CHKERR dAta.materialDoublePtr->calculateP_PiolaKirchhoffI(
153  CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
155  CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
156  &dAta.materialDoublePtr->P(0, 0));
157  CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
158  &dAta.materialDoublePtr->eNergy);
159  }
160 
162  MatrixDouble3by3 P(3, 3);
163  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
164  double val_energy;
165  CHKERR postProcMesh.tag_get_data(th_energy, &mapGaussPts[gg], 1,
166  &val_energy);
167  if (!std::isnormal(val_energy)) {
168  CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
169  &maxVal);
170  CHKERR postProcMesh.tag_get_data(th_piola1, &mapGaussPts[gg], 1,
171  &P(0, 0));
172  for (unsigned int r = 0; r != P.size1(); ++r) {
173  for (unsigned int c = 0; c != P.size2(); ++c) {
174  if (!std::isnormal(P(r, c)))
175  P(r, c) = copysign(maxVal, P(r, c));
176  }
177  }
178  CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
179  &P(0, 0));
180  }
181  }
182  }
183 
185  }
186 };
187 
188 /// \deprecated Use PostProcStress
190 
191 #endif //__POSTPROCSTRESSES_HPP__
const bool fieldDisp
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr
NonlinearElasticElement::BlockData & dAta
std::map< std::string, std::vector< MatrixDouble > > gradMap
Range tEts
const double maxVal
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
moab::Interface & postProcMesh
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
constexpr double H
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
const bool replaceNonANumberByMaxValue
NonlinearElasticElement::CommonData nonLinearElementCommonData
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
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)
Range tEts
constrains elements in block set
std::map< std::string, std::vector< VectorDouble > > fieldMap
DEPRECATED typedef PostProcStress PostPorcStress
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
common data used by volume elements
PostProcCommonOnRefMesh::CommonDataForVolume & commonData
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< EntityHandle > & mapGaussPts
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
#define DEPRECATED
Definition: definitions.h:29
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
data for calculation heat conductivity and heat capacity elements
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...