v0.14.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 
9 
10 /**
11  * \brief Operator post-procesing stresses for Hook isotropic material
12 
13  * Example how to use it
14 
15  \code
16  PostProcVolumeOnRefinedMesh post_proc(m_field);
17  {
18  CHKERR post_proc.generateReferenceElementMesh();
19  CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
20  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
21  CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
22  //add postprocessing for stresses
23  post_proc.getOpPtrVector().push_back(
24  new PostProcHookStress(
25  m_field,
26  post_proc.postProcMesh,
27  post_proc.mapGaussPts,
28  "DISPLACEMENT",
29  post_proc.commonData,
30  &elastic.setOfBlocks
31  )
32  );
33  CHKERR DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc);
34  CHKERR post_proc.writeFile("out.h5m");
35  }
36 
37  \endcode
38 
39  */
42 
45  std::vector<EntityHandle> &mapGaussPts;
47 
48 #ifdef __NONLINEAR_ELASTIC_HPP
49  /// Material block data, ket is block id
50  const std::map<int, NonlinearElasticElement::BlockData>
51  *setOfBlocksMaterialDataPtr;
52 #endif //__NONLINEAR_ELASTIC_HPP
53 
55 
56  /**
57  * Constructor
58  */
60  std::vector<EntityHandle> &map_gauss_pts,
61  const std::string field_name,
63 #ifdef __NONLINEAR_ELASTIC_HPP
64  const std::map<int, NonlinearElasticElement::BlockData>
65  *set_of_block_data_ptr = NULL,
66 #endif
67  const bool is_field_disp = true)
68  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
70  mField(m_field), postProcMesh(post_proc_mesh),
71  mapGaussPts(map_gauss_pts),
72 #ifdef __NONLINEAR_ELASTIC_HPP
73  setOfBlocksMaterialDataPtr(set_of_block_data_ptr),
74 #endif //__NONLINEAR_ELASTIC_HPP
75  commonData(common_data), isFieldDisp(is_field_disp) {
76  }
77 
78  /**
79  * \brief get material parameter
80 
81  * Material parameters are read form BlockSet, however if block data are
82  present,
83  * use data how are set for elastic element operators.
84 
85  * @param _lambda elastic material constant
86  * @param _mu elastic material constant
87  * @param _block_id block id
88  * @return error code
89 
90  */
91  MoFEMErrorCode getMatParameters(double *_lambda, double *_mu,
92  int *_block_id) {
94 
95  *_lambda = 1;
96  *_mu = 1;
97 
100  mField, BLOCKSET | MAT_ELASTICSET, it)) {
101  Mat_Elastic mydata;
102  CHKERR it->getAttributeDataStructure(mydata);
103 
104  Range meshsets;
105  CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBENTITYSET,
106  meshsets, false);
107  meshsets.insert(it->meshset);
108  for (Range::iterator mit = meshsets.begin(); mit != meshsets.end();
109  mit++) {
110  if (mField.get_moab().contains_entities(*mit, &ent, 1)) {
111  *_lambda = LAMBDA(mydata.data.Young, mydata.data.Poisson);
112  *_mu = MU(mydata.data.Young, mydata.data.Poisson);
113  *_block_id = it->getMeshsetId();
114 #ifdef __NONLINEAR_ELASTIC_HPP
115  if (setOfBlocksMaterialDataPtr) {
116  *_lambda =
117  LAMBDA(setOfBlocksMaterialDataPtr->at(*_block_id).E,
118  setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
119  *_mu = MU(setOfBlocksMaterialDataPtr->at(*_block_id).E,
120  setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
121  }
122 #endif //__NONLINEAR_ELASTIC_HPP
124  }
125  }
126  }
127 
128  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
129  "Element is not in elastic block, however you run linear elastic "
130  "analysis with that element\n"
131  "top tip: check if you update block sets after mesh refinements or "
132  "interface insertion");
133 
135  }
136 
137  /**
138  * \brief Here real work is done
139  */
140  MoFEMErrorCode doWork(int side, EntityType type,
143 
144  if (type != MBVERTEX)
146  if (data.getFieldData().size() == 0)
148 
149  int id;
150  double lambda, mu;
151  CHKERR getMatParameters(&lambda, &mu, &id);
152 
153  MatrixDouble D_lambda, D_mu, D;
154  D_lambda.resize(6, 6);
155  D_lambda.clear();
156  for (int rr = 0; rr < 3; rr++) {
157  for (int cc = 0; cc < 3; cc++) {
158  D_lambda(rr, cc) = 1;
159  }
160  }
161  D_mu.resize(6, 6);
162  D_mu.clear();
163  for (int rr = 0; rr < 6; rr++) {
164  D_mu(rr, rr) = rr < 3 ? 2 : 1;
165  }
166  D = lambda * D_lambda + mu * D_mu;
167 
168  int tag_length = 9;
169  double def_VAL[tag_length];
170  bzero(def_VAL, tag_length * sizeof(double));
171  Tag th_stress;
172  CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
173  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
174 
175  Tag th_id;
176  int def_block_id = -1;
177  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
178  MB_TAG_CREAT | MB_TAG_SPARSE,
179  &def_block_id);
180  Range::iterator tit = commonData.tEts.begin();
181  for (; tit != commonData.tEts.end(); tit++)
182  CHKERR postProcMesh.tag_set_data(th_id, &*tit, 1, &id);
183 
184  VectorDouble strain;
185  VectorDouble stress;
186  MatrixDouble Stress;
187 
188  int nb_gauss_pts = data.getN().size1();
189  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts)
190  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
191 
192  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
193 
194  strain.resize(6);
195  strain[0] = (commonData.gradMap[rowFieldName][gg])(0, 0);
196  strain[1] = (commonData.gradMap[rowFieldName][gg])(1, 1);
197  strain[2] = (commonData.gradMap[rowFieldName][gg])(2, 2);
198  strain[3] = (commonData.gradMap[rowFieldName][gg])(0, 1) +
199  (commonData.gradMap[rowFieldName][gg])(1, 0);
200  strain[4] = (commonData.gradMap[rowFieldName][gg])(1, 2) +
201  (commonData.gradMap[rowFieldName][gg])(2, 1);
202  strain[5] = (commonData.gradMap[rowFieldName][gg])(0, 2) +
203  (commonData.gradMap[rowFieldName][gg])(2, 0);
204 
205  if (!isFieldDisp) {
206  strain[0] -= 1.0;
207  strain[1] -= 1.0;
208  strain[2] -= 1.0;
209  }
210 
211  stress.resize(6);
212  noalias(stress) = prod(D, strain);
213 
214  Stress.resize(3, 3);
215  Stress(0, 0) = stress[0];
216  Stress(1, 1) = stress[1];
217  Stress(2, 2) = stress[2];
218  Stress(0, 1) = Stress(1, 0) = stress[3];
219  Stress(1, 2) = Stress(2, 1) = stress[4];
220  Stress(2, 0) = Stress(0, 2) = stress[5];
221 
222  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
223  &Stress(0, 0));
224  }
225 
227  }
228 };
229 
230 /// \deprecated Class name with spelling mistake
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
EntityHandle
PostProcHookStress::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: PostProcHookStresses.hpp:45
PostProcCommonOnRefMesh::CommonDataForVolume
Definition: PostProcOnRefMesh.hpp:37
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PostProcHookStress::isFieldDisp
bool isFieldDisp
Definition: PostProcHookStresses.hpp:46
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
PostProcHookStress::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcHookStresses.hpp:44
PostProcHookStress::getMatParameters
MoFEMErrorCode getMatParameters(double *_lambda, double *_mu, int *_block_id)
get material parameter
Definition: PostProcHookStresses.hpp:91
MoFEM::ForcesAndSourcesCore::UserDataOperator::rowFieldName
std::string rowFieldName
Definition: ForcesAndSourcesCore.hpp:577
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::ForcesAndSourcesCore::UserDataOperator::ForcesAndSourcesCore
friend class ForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:990
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PostProcHookStress::PostProcHookStress
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)
Definition: PostProcHookStresses.hpp:59
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
convert.type
type
Definition: convert.py:64
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
PostProcHookStress::commonData
PostProcVolumeOnRefinedMesh::CommonData & commonData
Definition: PostProcHookStresses.hpp:54
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
PostProcCommonOnRefMesh::CommonDataForVolume::tEts
Range tEts
Definition: PostProcOnRefMesh.hpp:38
PostProcHookStress::mField
MoFEM::Interface & mField
Definition: PostProcHookStresses.hpp:43
Range
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
PostProcHookStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Here real work is done.
Definition: PostProcHookStresses.hpp:140
PostProcCommonOnRefMesh::CommonData::gradMap
std::map< std::string, std::vector< MatrixDouble > > gradMap
Definition: PostProcOnRefMesh.hpp:34
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MU
#define MU(E, NU)
Definition: fem_tools.h:23
PostPorcHookStress
DEPRECATED typedef PostProcHookStress PostPorcHookStress
Definition: PostProcHookStresses.hpp:231
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567