v0.13.1
Loading...
Searching...
No Matches
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
44 moab::Interface &postProcMesh;
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 */
59 PostProcHookStress(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh,
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
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,
141 EntitiesFieldData::EntData &data) {
143
144 if (type != MBVERTEX)
146 if (data.getFieldData().size() == 0)
148
149 int id;
150 double lambda, mu;
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
DEPRECATED typedef PostProcHookStress PostPorcHookStress
#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
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ 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
#define MU(E, NU)
Definition: fem_tools.h:23
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
constexpr double lambda
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
double D
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr auto field_name
constexpr double mu
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPROW
operator doWork function is executed on FE rows
Range tEts
std::map< std::string, std::vector< MatrixDouble > > gradMap
Operator post-procesing stresses for Hook isotropic material.
moab::Interface & postProcMesh
PostProcVolumeOnRefinedMesh::CommonData & commonData
MoFEM::Interface & mField
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Here real work is done.
MoFEMErrorCode getMatParameters(double *_lambda, double *_mu, int *_block_id)
get material parameter
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)
bool isFieldDisp
std::vector< EntityHandle > & mapGaussPts