v0.9.1
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;
58 
59 #ifdef __NONLINEAR_ELASTIC_HPP
60  /// Material block data, ket is block id
61  const std::map<int, NonlinearElasticElement::BlockData>
62  *setOfBlocksMaterialDataPtr;
63 #endif //__NONLINEAR_ELASTIC_HPP
64 
66 
67  /**
68  * Constructor
69  */
71  std::vector<EntityHandle> &map_gauss_pts,
72  const std::string field_name,
74 #ifdef __NONLINEAR_ELASTIC_HPP
75  ,
76  const std::map<int, NonlinearElasticElement::BlockData>
77  *set_of_block_data_ptr = NULL
78 #endif
79  )
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) {
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  // const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0];
162 
163  int id;
164  double lambda, mu;
165  CHKERR getMatParameters(&lambda, &mu, &id);
166 
167  MatrixDouble D_lambda, D_mu, D;
168  D_lambda.resize(6, 6);
169  D_lambda.clear();
170  for (int rr = 0; rr < 3; rr++) {
171  for (int cc = 0; cc < 3; cc++) {
172  D_lambda(rr, cc) = 1;
173  }
174  }
175  D_mu.resize(6, 6);
176  D_mu.clear();
177  for (int rr = 0; rr < 6; rr++) {
178  D_mu(rr, rr) = rr < 3 ? 2 : 1;
179  }
180  D = lambda * D_lambda + mu * D_mu;
181 
182  int tag_length = 9;
183  double def_VAL[tag_length];
184  bzero(def_VAL, tag_length * sizeof(double));
185  Tag th_stress;
186  CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
187  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
188 
189  Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
190  CHKERR postProcMesh.tag_get_handle("S1", 3, MB_TYPE_DOUBLE,
191  th_prin_stress_vect1,
192  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
193  CHKERR postProcMesh.tag_get_handle("S2", 3, MB_TYPE_DOUBLE,
194  th_prin_stress_vect2,
195  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
196  CHKERR postProcMesh.tag_get_handle("S3", 3, MB_TYPE_DOUBLE,
197  th_prin_stress_vect3,
198  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
199 
200  Tag th_prin_stress_vals;
201  CHKERR postProcMesh.tag_get_handle("PRINCIPAL_STRESS", 3, MB_TYPE_DOUBLE,
202  th_prin_stress_vals,
203  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
204 
205  Tag th_id;
206  int def_block_id = -1;
207  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
208  MB_TAG_CREAT | MB_TAG_SPARSE,
209  &def_block_id);
210  Range::iterator tit = commonData.tEts.begin();
211  for (; tit != commonData.tEts.end(); tit++) {
212  CHKERR postProcMesh.tag_set_data(th_id, &*tit, 1, &id);
213  }
214 
215  VectorDouble strain;
216  VectorDouble stress;
217  MatrixDouble Stress;
218 
219  // Combine eigenvalues and vectors to create principal stress vector
220  MatrixDouble prin_stress_vect(3, 3);
221  VectorDouble prin_vals_vect(3);
222 
223  int nb_gauss_pts = data.getN().size1();
224  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
225  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
226  }
227  for (int gg = 0; gg < nb_gauss_pts; gg++) {
228 
229  strain.resize(6);
230  strain[0] = (commonData.gradMap[rowFieldName][gg])(0, 0);
231  strain[1] = (commonData.gradMap[rowFieldName][gg])(1, 1);
232  strain[2] = (commonData.gradMap[rowFieldName][gg])(2, 2);
233  strain[3] = (commonData.gradMap[rowFieldName][gg])(0, 1) +
234  (commonData.gradMap[rowFieldName][gg])(1, 0);
235  strain[4] = (commonData.gradMap[rowFieldName][gg])(1, 2) +
236  (commonData.gradMap[rowFieldName][gg])(2, 1);
237  strain[5] = (commonData.gradMap[rowFieldName][gg])(0, 2) +
238  (commonData.gradMap[rowFieldName][gg])(2, 0);
239 
240  stress.resize(6);
241  noalias(stress) = prod(D, strain);
242 
243  Stress.resize(3, 3);
244  Stress(0, 0) = stress[0];
245  Stress(1, 1) = stress[1];
246  Stress(2, 2) = stress[2];
247  Stress(0, 1) = Stress(1, 0) = stress[3];
248  Stress(1, 2) = Stress(2, 1) = stress[4];
249  Stress(2, 0) = Stress(0, 2) = stress[5];
250 
251  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
252  &Stress(0, 0));
253 
254  MatrixDouble eigen_vectors = Stress;
255  VectorDouble eigen_values(3);
256 
257  // LAPACK - eigenvalues and vectors. Applied twice for initial creates
258  // memory space
259  int n = 3, lda = 3, info, lwork = -1;
260  double wkopt;
261  info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
262  &(eigen_values.data()[0]), &wkopt, lwork);
263  if (info != 0)
264  SETERRQ1(PETSC_COMM_SELF, 1,
265  "is something wrong with lapack_dsyev info = %d", info);
266  lwork = (int)wkopt;
267  double work[lwork];
268  info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
269  &(eigen_values.data()[0]), work, lwork);
270  if (info != 0)
271  SETERRQ1(PETSC_COMM_SELF, 1,
272  "is something wrong with lapack_dsyev info = %d", info);
273 
274  map<double, int> eigen_sort;
275  eigen_sort[eigen_values[0]] = 0;
276  eigen_sort[eigen_values[1]] = 1;
277  eigen_sort[eigen_values[2]] = 2;
278 
279  prin_stress_vect.clear();
280  prin_vals_vect.clear();
281 
282  int ii = 0;
283  for (map<double, int>::reverse_iterator mit = eigen_sort.rbegin();
284  mit != eigen_sort.rend(); mit++) {
285  prin_vals_vect[ii] = eigen_values[mit->second];
286  for (int dd = 0; dd != 3; dd++) {
287  prin_stress_vect(ii, dd) = eigen_vectors.data()[3 * mit->second + dd];
288  }
289  ii++;
290  }
291 
292  // Tag principle stress vectors 1, 2, 3
293  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
294  1, &prin_stress_vect(0, 0));
295  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
296  1, &prin_stress_vect(1, 0));
297  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
298  1, &prin_stress_vect(2, 0));
299  CHKERR postProcMesh.tag_set_data(th_prin_stress_vals, &mapGaussPts[gg], 1,
300  &prin_vals_vect[0]);
301  }
302 
304  }
305 };
306 
307 /// \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
Range tEts
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)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
std::vector< EntityHandle > & mapGaussPts
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
MoFEMErrorCode getMatParameters(double *_lambda, double *_mu, int *_block_id)
get material parameter
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Operator post-procesing stresses for Hook isotropic material.
MoFEM::Interface & mField
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
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
#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:601
#define MU(E, NU)
Definition: fem_tools.h:33
block name is "MAT_ELASTIC"
Definition: definitions.h:227
#define DEPRECATED
Definition: definitions.h:29
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72