v0.10.0
PostProcHookStresses.hpp
Go to the documentation of this file.
1 /**
2  * \file PostProcHookStress.hpp
3  * \brief Post-proc stresses for linear Hooke isotropic material
4  *
5  * \ingroup nonlinear_elastic_elem
6  */
7 
8 /* This file is part of MoFEM.
9  * MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13  *
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21 
22 /**
23  * \brief Operator post-procesing stresses for Hook isotropic material
24 
25  * Example how to use it
26 
27  \code
28  PostProcVolumeOnRefinedMesh post_proc(m_field);
29  {
30  CHKERR post_proc.generateReferenceElementMesh();
31  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
32  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
33  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
34  //add postprocessing for stresses
35  post_proc.getOpPtrVector().push_back(
36  new PostProcHookStress(
37  m_field,
38  post_proc.postProcMesh,
39  post_proc.mapGaussPts,
40  "DISPLACEMENT",
41  post_proc.commonData,
42  &elastic.setOfBlocks
43  )
44  );
45  CHKERR DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc);
46  CHKERR post_proc.writeFile("out.h5m");
47  }
48 
49  \endcode
50 
51  */
54 
57  std::vector<EntityHandle> &mapGaussPts;
59 
60 #ifdef __NONLINEAR_ELASTIC_HPP
61  /// Material block data, ket is block id
62  const std::map<int, NonlinearElasticElement::BlockData>
63  *setOfBlocksMaterialDataPtr;
64 #endif //__NONLINEAR_ELASTIC_HPP
65 
67 
68  /**
69  * Constructor
70  */
72  std::vector<EntityHandle> &map_gauss_pts,
73  const std::string field_name,
75 #ifdef __NONLINEAR_ELASTIC_HPP
76  const std::map<int, NonlinearElasticElement::BlockData>
77  *set_of_block_data_ptr = NULL,
78 #endif
79  const bool is_field_disp = true)
82  mField(m_field), postProcMesh(post_proc_mesh),
83  mapGaussPts(map_gauss_pts),
84 #ifdef __NONLINEAR_ELASTIC_HPP
85  setOfBlocksMaterialDataPtr(set_of_block_data_ptr),
86 #endif //__NONLINEAR_ELASTIC_HPP
87  commonData(common_data), isFieldDisp(is_field_disp) {
88  }
89 
90  /**
91  * \brief get material parameter
92 
93  * Material parameters are read form BlockSet, however if block data are
94  present,
95  * use data how are set for elastic element operators.
96 
97  * @param _lambda elastic material constant
98  * @param _mu elastic material constant
99  * @param _block_id block id
100  * @return error code
101 
102  */
103  MoFEMErrorCode getMatParameters(double *_lambda, double *_mu,
104  int *_block_id) {
106 
107  *_lambda = 1;
108  *_mu = 1;
109 
112  mField, BLOCKSET | MAT_ELASTICSET, it)) {
113  Mat_Elastic mydata;
114  CHKERR it->getAttributeDataStructure(mydata);
115 
116  Range meshsets;
117  CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBENTITYSET,
118  meshsets, false);
119  meshsets.insert(it->meshset);
120  for (Range::iterator mit = meshsets.begin(); mit != meshsets.end();
121  mit++) {
122  if (mField.get_moab().contains_entities(*mit, &ent, 1)) {
123  *_lambda = LAMBDA(mydata.data.Young, mydata.data.Poisson);
124  *_mu = MU(mydata.data.Young, mydata.data.Poisson);
125  *_block_id = it->getMeshsetId();
126 #ifdef __NONLINEAR_ELASTIC_HPP
127  if (setOfBlocksMaterialDataPtr) {
128  *_lambda =
129  LAMBDA(setOfBlocksMaterialDataPtr->at(*_block_id).E,
130  setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
131  *_mu = MU(setOfBlocksMaterialDataPtr->at(*_block_id).E,
132  setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
133  }
134 #endif //__NONLINEAR_ELASTIC_HPP
136  }
137  }
138  }
139 
140  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
141  "Element is not in elastic block, however you run linear elastic "
142  "analysis with that element\n"
143  "top tip: check if you update block sets after mesh refinements or "
144  "interface insertion");
145 
147  }
148 
149  /**
150  * \brief Here real work is done
151  */
152  MoFEMErrorCode doWork(int side, EntityType type,
155 
156  if (type != MBVERTEX)
158  if (data.getFieldData().size() == 0)
160 
161  int id;
162  double lambda, mu;
163  CHKERR getMatParameters(&lambda, &mu, &id);
164 
165  MatrixDouble D_lambda, D_mu, D;
166  D_lambda.resize(6, 6);
167  D_lambda.clear();
168  for (int rr = 0; rr < 3; rr++) {
169  for (int cc = 0; cc < 3; cc++) {
170  D_lambda(rr, cc) = 1;
171  }
172  }
173  D_mu.resize(6, 6);
174  D_mu.clear();
175  for (int rr = 0; rr < 6; rr++) {
176  D_mu(rr, rr) = rr < 3 ? 2 : 1;
177  }
178  D = lambda * D_lambda + mu * D_mu;
179 
180  int tag_length = 9;
181  double def_VAL[tag_length];
182  bzero(def_VAL, tag_length * sizeof(double));
183  Tag th_stress;
184  CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
185  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
186 
187  Tag th_id;
188  int def_block_id = -1;
189  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
190  MB_TAG_CREAT | MB_TAG_SPARSE,
191  &def_block_id);
192  Range::iterator tit = commonData.tEts.begin();
193  for (; tit != commonData.tEts.end(); tit++)
194  CHKERR postProcMesh.tag_set_data(th_id, &*tit, 1, &id);
195 
196  VectorDouble strain;
197  VectorDouble stress;
198  MatrixDouble Stress;
199 
200  int nb_gauss_pts = data.getN().size1();
201  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts)
202  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
203 
204  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
205 
206  strain.resize(6);
207  strain[0] = (commonData.gradMap[rowFieldName][gg])(0, 0);
208  strain[1] = (commonData.gradMap[rowFieldName][gg])(1, 1);
209  strain[2] = (commonData.gradMap[rowFieldName][gg])(2, 2);
210  strain[3] = (commonData.gradMap[rowFieldName][gg])(0, 1) +
211  (commonData.gradMap[rowFieldName][gg])(1, 0);
212  strain[4] = (commonData.gradMap[rowFieldName][gg])(1, 2) +
213  (commonData.gradMap[rowFieldName][gg])(2, 1);
214  strain[5] = (commonData.gradMap[rowFieldName][gg])(0, 2) +
215  (commonData.gradMap[rowFieldName][gg])(2, 0);
216 
217  if (!isFieldDisp) {
218  strain[0] -= 1.0;
219  strain[1] -= 1.0;
220  strain[2] -= 1.0;
221  }
222 
223  stress.resize(6);
224  noalias(stress) = prod(D, strain);
225 
226  Stress.resize(3, 3);
227  Stress(0, 0) = stress[0];
228  Stress(1, 1) = stress[1];
229  Stress(2, 2) = stress[2];
230  Stress(0, 1) = Stress(1, 0) = stress[3];
231  Stress(1, 2) = Stress(2, 1) = stress[4];
232  Stress(2, 0) = Stress(0, 2) = stress[5];
233 
234  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
235  &Stress(0, 0));
236  }
237 
239  }
240 };
241 
242 /// \deprecated Class name with spelling mistake
PostProcVolumeOnRefinedMesh::CommonData & commonData
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
Deprecated interface functions.
moab::Interface & postProcMesh
DEPRECATED typedef PostProcHookStress PostPorcHookStress
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Here real work is done.
virtual moab::Interface & get_moab()=0
std::map< std::string, std::vector< MatrixDouble > > gradMap
const double D
diffusivity
Range tEts
bool isFieldDisp
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
std::vector< EntityHandle > & mapGaussPts
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEMErrorCode getMatParameters(double *_lambda, double *_mu, int *_block_id)
get material parameter
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PostProcHookStress(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, PostProcVolumeOnRefinedMesh::CommonData &common_data, const bool is_field_disp=true)
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Operator post-procesing stresses for Hook isotropic material.
MoFEM::Interface & mField
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MU(E, NU)
Definition: fem_tools.h:33
block name is "MAT_ELASTIC"
Definition: definitions.h:228
#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
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...