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.
36
39
40 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
41
42 struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
43 OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
45 mPtr(m_ptr),
scale(s) {}
46 MoFEMErrorCode doWork(
int, EntityType, EntitiesFieldData::EntData &) {
48 return 0;
49 }
50
51 private:
52 boost::shared_ptr<MatrixDouble> mPtr;
54 };
55
56 auto push_domain_ops = [&](auto &pip) {
58 pip, {
H1,
HDIV},
"GEOMETRY")),
59 "Apply base transform");
60
61
62 using CommonDataVariant =
63 std::variant<boost::shared_ptr<HookeOps::CommonData>,
64 boost::shared_ptr<HenckyOps::CommonData>>;
65 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
66 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
67 "SIGMA", contact_stress_ptr));
68
70 auto hooke_common_data_ptr =
71 HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
72 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
73
74 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
75 return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
76 contact_stress_ptr);
77 } else {
78 auto hencky_common_data_ptr =
79 HenckyOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
80 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform,
scale);
81 pip.push_back(
new OpScale(contact_stress_ptr,
scale));
82 return std::make_tuple(CommonDataVariant(hencky_common_data_ptr),
83 contact_stress_ptr);
84 }
85 };
86 auto push_bdy_ops = [&](auto &pip) {
87
88 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
89 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
90 "U", common_data_ptr->contactDispPtr()));
91 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
92 "SIGMA", common_data_ptr->contactTractionPtr()));
93 pip.push_back(
new OpScale(common_data_ptr->contactTractionPtr(),
scale));
94 using C = ContactIntegrators<BoundaryEleOp>;
95 pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
96 common_data_ptr));
97 return common_data_ptr;
98 };
99
100 auto get_domain_pip = [&](auto &pip)
101 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
103 auto op_loop_side = new OpLoopSide<SideEle>(
104 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
105 boost::make_shared<
106 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
107 pip.push_back(op_loop_side);
108 return op_loop_side->getOpPtrVector();
109 } else {
110 return pip;
111 }
112 };
113
114 auto get_post_proc_domain_fe = [&]() {
115 auto post_proc_fe =
116 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
117 auto &pip = post_proc_fe->getOpPtrVector();
118
119 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
120 push_domain_ops(get_domain_pip(pip));
121
122 auto u_ptr = boost::make_shared<MatrixDouble>();
123 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
124 auto X_ptr = boost::make_shared<MatrixDouble>();
125 pip.push_back(
126 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
127
128
129
130
131
132
133
134
135
136
137 std::visit(
138 [&](auto &common_data_ptr) {
139 if constexpr (std::is_same_v<
140 std::decay_t<decltype(common_data_ptr)>,
141 boost::shared_ptr<HookeOps::CommonData>>) {
143 post_proc_fe->getPostProcMesh(),
144 post_proc_fe->getMapGaussPts(), {},
145 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
146 {{"SIGMA", contact_stress_ptr},
147 {"G", common_data_ptr->matGradPtr}},
148 {{"STRESS", common_data_ptr->getMatCauchyStress()},
149 {"STRAIN", common_data_ptr->getMatStrain()}}));
150 }
151
152 else if constexpr (std::is_same_v<
153 std::decay_t<decltype(common_data_ptr)>,
154 boost::shared_ptr<HenckyOps::CommonData>>) {
155
157 post_proc_fe->getPostProcMesh(),
158 post_proc_fe->getMapGaussPts(), {},
159 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
160 {{"SIGMA", contact_stress_ptr},
161 {"G", common_data_ptr->matGradPtr},
162 {"PK1", common_data_ptr->getMatFirstPiolaStress()}
163
164 },
165 {{"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
166 }
167 },
168 hencky_or_hooke_common_data_ptr);
169
171
173 pip, {
HDIV},
"GEOMETRY")),
174 "Apply transform");
175 auto common_data_ptr = push_bdy_ops(pip);
176
177 pip.push_back(
178
180
181 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
182
183 {{"SDF", common_data_ptr->sdfPtr()},
184 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
185
186 {
187
188 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
189 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
190
191 },
192
193 {},
194
195 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
196
197 )
198
199 );
200 }
201
202 return post_proc_fe;
203 };
204
205 auto get_post_proc_bdy_fe = [&]() {
206 auto post_proc_fe =
207 boost::make_shared<PostProcEleBdy>(*m_field_ptr,
postProcMesh);
208 auto &pip = post_proc_fe->getOpPtrVector();
209
211 pip, {
HDIV},
"GEOMETRY")),
212 "Apply transform");
213 auto common_data_ptr = push_bdy_ops(pip);
214
215
216 auto op_loop_side = new OpLoopSide<SideEle>(
217 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
218 boost::make_shared<
219 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
220 pip.push_back(op_loop_side);
221
222 auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
223 push_domain_ops(op_loop_side->getOpPtrVector());
224
225 auto X_ptr = boost::make_shared<MatrixDouble>();
226 pip.push_back(
227 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
228
229 pip.push_back(
230
232
233 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
234
235 {{"SDF", common_data_ptr->sdfPtr()},
236 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
237
238 {{"U", common_data_ptr->contactDispPtr()},
239 {"GEOMETRY", X_ptr},
240 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
241 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
242
243 },
244
245 {},
246
247 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
248
249 )
250
251 );
252
253 return post_proc_fe;
254 };
255
256 auto get_integrate_traction = [&]() {
257 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
259 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
260 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
261 "Apply transform");
262
263
264 integrate_traction->getOpPtrVector().push_back(
265 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
268 };
269
271 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
273 "push operators to calculate traction");
274
275 return integrate_traction;
276 };
277
278 auto get_integrate_area = [&]() {
279 auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
280
282 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
283 integrate_area->getOpPtrVector(), {HDIV}, "GEOMETRY")),
284 "Apply transform");
285
286
287 integrate_area->getOpPtrVector().push_back(
288 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
291 };
294 m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
295 std::regex((boost::format(
"%s(.*)") %
"CONTACT").str()))) {
296 auto meshset =
m->getMeshset();
297 Range contact_meshset_range;
299 meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
300
302 contact_meshset_range);
303 contact_range.merge(contact_meshset_range);
304 }
305
306 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
307
308 auto op_loop_side = new OpLoopSide<SideEle>(
309 *m_field_ptr, m_field_ptr->
getInterface<Simple>()->getDomainFEName(),
311 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
312 op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
313
315 (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
316 integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
318 "push operators to calculate area");
319
320 return integrate_area;
321 };
322
326
329 }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ HDIV
field with continuous normal traction
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int geom_order
Order if fixed.
static constexpr int approx_order
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.