Adds operators for evaluating stress and strain based on the material model: i.e. Visit the variant to handle either HookeOps::CommonData or HenckyOps::CommonData. i.e.
35
38
39 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
40
41 struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
42 OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
44 mPtr(m_ptr),
scale(s) {}
45 MoFEMErrorCode doWork(
int, EntityType, EntitiesFieldData::EntData &) {
47 return 0;
48 }
49
50 private:
51 boost::shared_ptr<MatrixDouble> mPtr;
53 };
54
55 auto push_domain_ops = [&](auto &pip) {
57 pip, {
H1,
HDIV},
"GEOMETRY")),
58 "Apply base transform");
59
60
61 using CommonDataVariant =
62 std::variant<boost::shared_ptr<HookeOps::CommonData>,
63 boost::shared_ptr<HenckyOps::CommonData>>;
64 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
65 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
66 "SIGMA", contact_stress_ptr));
67
69 auto hooke_common_data_ptr =
70 HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
71 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
72
73 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
74 return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
75 contact_stress_ptr);
76 } else {
77 auto hencky_common_data_ptr =
78 HenckyOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
79 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
80 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
81 return std::make_tuple(CommonDataVariant(hencky_common_data_ptr),
82 contact_stress_ptr);
83 }
84 };
85 auto push_bdy_ops = [&](auto &pip) {
86
87 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
88 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
89 "U", common_data_ptr->contactDispPtr()));
90 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
91 "SIGMA", common_data_ptr->contactTractionPtr()));
92 pip.push_back(
new OpScale(common_data_ptr->contactTractionPtr(),
scale));
93 using C = ContactIntegrators<BoundaryEleOp>;
94 pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
95 common_data_ptr));
96 return common_data_ptr;
97 };
98
99 auto get_domain_pip = [&](auto &pip)
100 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
102 auto op_loop_side = new OpLoopSide<SideEle>(
103 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
104 boost::make_shared<
106 pip.push_back(op_loop_side);
107 return op_loop_side->getOpPtrVector();
108 } else {
109 return pip;
110 }
111 };
112
113
114
115 auto get_post_proc_domain_fe = [&]() {
116 auto post_proc_fe =
117 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
118 auto &pip = post_proc_fe->getOpPtrVector();
119
120 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
121 push_domain_ops(get_domain_pip(pip));
122
123 auto u_ptr = boost::make_shared<MatrixDouble>();
124 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
125 auto X_ptr = boost::make_shared<MatrixDouble>();
126 pip.push_back(
127 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
128
129
130
131
132
133
134
135 std::visit(
136 [&](auto &common_data_ptr) {
137 if constexpr (std::is_same_v<
138 std::decay_t<decltype(common_data_ptr)>,
139 boost::shared_ptr<HookeOps::CommonData>>) {
141 post_proc_fe->getPostProcMesh(),
142 post_proc_fe->getMapGaussPts(), {},
143 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
144 {{"SIGMA", contact_stress_ptr},
145 {"G", common_data_ptr->matGradPtr}},
146 {{"STRESS", common_data_ptr->getMatCauchyStress()},
147 {"STRAIN", common_data_ptr->getMatStrain()}}));
148 }
149
150 else if constexpr (std::is_same_v<
151 std::decay_t<decltype(common_data_ptr)>,
152 boost::shared_ptr<
154
156 post_proc_fe->getPostProcMesh(),
157 post_proc_fe->getMapGaussPts(), {},
158 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
159 {{"SIGMA", contact_stress_ptr},
160 {"G", common_data_ptr->matGradPtr},
161 {"PK1", common_data_ptr->getMatFirstPiolaStress()}
162
163 },
164 {{"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
165 }
166 },
167 hencky_or_hooke_common_data_ptr);
168
170
172 pip, {
HDIV},
"GEOMETRY")),
173 "Apply transform");
174 auto common_data_ptr = push_bdy_ops(pip);
175
176 pip.push_back(
177
179
180 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
181
182 {{"SDF", common_data_ptr->sdfPtr()},
183 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
184
185 {
186
187 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
188 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
189
190 },
191
192 {},
193
194 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
195
196 )
197
198 );
199 }
200
201 return post_proc_fe;
202 };
203
204 auto get_post_proc_bdy_fe = [&]() {
205 auto post_proc_fe =
206 boost::make_shared<PostProcEleBdy>(*m_field_ptr,
postProcMesh);
207 auto &pip = post_proc_fe->getOpPtrVector();
208
210 pip, {
HDIV},
"GEOMETRY")),
211 "Apply transform");
212 auto common_data_ptr = push_bdy_ops(pip);
213
214
215 auto op_loop_side = new OpLoopSide<SideEle>(
216 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
217 boost::make_shared<
219 pip.push_back(op_loop_side);
220
221 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
222 push_domain_ops(op_loop_side->getOpPtrVector());
223
224 auto X_ptr = boost::make_shared<MatrixDouble>();
225 pip.push_back(
226 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
227
228 pip.push_back(
229
231
232 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
233
234 {{"SDF", common_data_ptr->sdfPtr()},
235 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
236
237 {{"U", common_data_ptr->contactDispPtr()},
238 {"GEOMETRY", X_ptr},
239 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
240 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
241
242 },
243
244 {},
245
246 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
247
248 )
249
250 );
251
252 return post_proc_fe;
253 };
254
255 auto get_integrate_traction = [&]() {
256 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
259 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
260 "Apply transform");
261
262 integrate_traction->getOpPtrVector().push_back(
263 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
266 };
267
269 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
271 "push operators to calculate traction");
272
273 return integrate_traction;
274 };
275
276 auto get_integrate_area = [&]() {
277 auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
278
281 integrate_area->getOpPtrVector(), {HDIV}, "GEOMETRY")),
282 "Apply transform");
283
284
285 integrate_area->getOpPtrVector().push_back(
286 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
289 };
292 m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
293 std::regex((boost::format(
"%s(.*)") %
"CONTACT").str()))) {
294 auto meshset =
m->getMeshset();
295 Range contact_meshset_range;
297 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
298
300 contact_meshset_range);
301 contact_range.merge(contact_meshset_range);
302 }
303
304 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
305
306 auto op_loop_side = new OpLoopSide<SideEle>(
307 *m_field_ptr, m_field_ptr->
getInterface<Simple>()->getDomainFEName(),
310 op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
311
313 (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
314 integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
316 "push operators to calculate area");
317
318 return integrate_area;
319 };
320
324
327
331 }
