v0.14.0
Loading...
Searching...
No Matches
PostProcContact.hpp
Go to the documentation of this file.
1/**
2 * \file PostProcContact.hpp
3 *
4 *
5 * @copyright Copyright (c) 2023
6 */
7
8namespace ContactOps {
9
10template <int DIM> struct PostProcEleByDim;
11
12template <> struct PostProcEleByDim<2> {
13 using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
14 using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
15 using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
16};
17
18template <> struct PostProcEleByDim<3> {
19 using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
20 using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
21 using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
22};
23
27
28struct Monitor : public FEMethod {
29
30 Monitor(SmartPetscObj<DM> &dm, bool use_mfront = false,
31 boost::shared_ptr<GenericElementInterface> mfront_interface = nullptr,
32 bool is_axisymmetric = false, int atom_test = 0)
36
37 MoFEM::Interface *m_field_ptr;
38 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
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)
44 : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
45 mPtr(m_ptr), scale(s) {}
46 MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
47 *mPtr *= 1./scale;
48 return 0;
49 }
50
51 private:
52 boost::shared_ptr<MatrixDouble> mPtr;
53 double scale;
54 };
55
56 auto push_domain_ops = [&](auto &pip) {
57 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
58 pip, {H1, HDIV}, "GEOMETRY")),
59 "Apply base transform");
60 auto henky_common_data_ptr =
61 commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
62 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
63 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
64 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
65 "SIGMA", contact_stress_ptr));
66 pip.push_back(new OpScale(contact_stress_ptr, scale));
67 return std::make_tuple(henky_common_data_ptr, contact_stress_ptr);
68 };
69
70 auto push_bdy_ops = [&](auto &pip) {
71 // evaluate traction
72 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
73 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
74 "U", common_data_ptr->contactDispPtr()));
75 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
76 "SIGMA", common_data_ptr->contactTractionPtr()));
77 pip.push_back(new OpScale(common_data_ptr->contactTractionPtr(), scale));
79 pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
80 common_data_ptr));
81 return common_data_ptr;
82 };
83
84 auto get_domain_pip = [&](auto &pip)
85 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
86 if constexpr (SPACE_DIM == 3) {
87 auto op_loop_side = new OpLoopSide<SideEle>(
88 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
89 boost::make_shared<
90 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
91 pip.push_back(op_loop_side);
92 return op_loop_side->getOpPtrVector();
93 } else {
94 return pip;
95 }
96 };
97
98 auto get_post_proc_domain_fe = [&]() {
99 auto post_proc_fe =
100 boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
101 auto &pip = post_proc_fe->getOpPtrVector();
102
103 auto [henky_common_data_ptr, contact_stress_ptr] =
104 push_domain_ops(get_domain_pip(pip));
105
106 auto u_ptr = boost::make_shared<MatrixDouble>();
107 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
108 auto X_ptr = boost::make_shared<MatrixDouble>();
109 pip.push_back(
110 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
111
112
113
114 pip.push_back(
115
116 new OpPPMap(
117
118 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
119
120 {},
121 {
122
123 {"U", u_ptr}, {"GEOMETRY", X_ptr}
124
125 },
126 {
127
128 {"SIGMA", contact_stress_ptr},
129
130 {"G", henky_common_data_ptr->matGradPtr},
131
132 {"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
133
134 },
135 {}
136
137 )
138
139 );
140
141 if (SPACE_DIM == 3) {
142
143 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
144 pip, {HDIV}, "GEOMETRY")),
145 "Apply transform");
146 auto common_data_ptr = push_bdy_ops(pip);
147
148 pip.push_back(
149
150 new OpPPMap(
151
152 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
153
154 {{"SDF", common_data_ptr->sdfPtr()},
155 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
156
157 {
158
159 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
160 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
161
162 },
163
164 {},
165
166 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
167
168 )
169
170 );
171 }
172
173 return post_proc_fe;
174 };
175
176 auto get_post_proc_bdy_fe = [&]() {
177 auto post_proc_fe =
178 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
179 auto &pip = post_proc_fe->getOpPtrVector();
180
181 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
182 pip, {HDIV}, "GEOMETRY")),
183 "Apply transform");
184 auto common_data_ptr = push_bdy_ops(pip);
185
186 // create OP which run element on side
187 auto op_loop_side = new OpLoopSide<SideEle>(
188 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
189 boost::make_shared<
190 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
191 pip.push_back(op_loop_side);
192
193 auto [henky_common_data_ptr, contact_stress_ptr] =
194 push_domain_ops(op_loop_side->getOpPtrVector());
195
196 auto X_ptr = boost::make_shared<MatrixDouble>();
197 pip.push_back(
198 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
199
200 pip.push_back(
201
202 new OpPPMap(
203
204 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
205
206 {{"SDF", common_data_ptr->sdfPtr()},
207 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
208
209 {{"U", common_data_ptr->contactDispPtr()},
210 {"GEOMETRY", X_ptr},
211 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
212 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
213
214 },
215
216 {},
217
218 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
219
220 )
221
222 );
223
224 return post_proc_fe;
225 };
226
227 auto get_integrate_traction = [&]() {
228 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
229 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
231 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
232 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
233 "Apply transform");
234 // We have to integrate on curved face geometry, thus integration weight
235 // have to adjusted.
236 integrate_traction->getOpPtrVector().push_back(
237 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
238 integrate_traction->getRuleHook = [](int, int, int approx_order) {
239 return 2 * approx_order + geom_order - 1;
240 };
241
243 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
244 integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
245 "push operators to calculate traction");
246
247 return integrate_traction;
248 };
249
250 postProcDomainFe = get_post_proc_domain_fe();
251 if constexpr (SPACE_DIM == 2)
252 postProcBdyFe = get_post_proc_bdy_fe();
253 integrateTraction = get_integrate_traction();
254 }
255
256 MoFEMErrorCode preProcess() { return 0; }
257 MoFEMErrorCode operator()() { return 0; }
258
259 MoFEMErrorCode postProcess() {
261 MoFEM::Interface *m_field_ptr;
262 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
263
264 auto post_proc = [&]() {
266
267 if (!useMFront) {
268 auto post_proc_begin =
269 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
270 postProcMesh);
271 auto post_proc_end =
272 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
273 postProcMesh);
274
275 CHKERR DMoFEMPreProcessFiniteElements(dM,
276 post_proc_begin->getFEMethod());
277 if (!postProcBdyFe) {
278 postProcDomainFe->copyTs(*this); // this here is a Monitor
279 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
280 } else {
281 postProcDomainFe->copyTs(*this); // this here is a Monitor
282 postProcBdyFe->copyTs(*this);
283 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
284 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
285 }
286 CHKERR DMoFEMPostProcessFiniteElements(dM,
287 post_proc_end->getFEMethod());
288
289 CHKERR post_proc_end->writeFile(
290 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
291 } else {
292 mfrontInterface->updateElementVariables();
293 mfrontInterface->postProcessElement(ts_step);
294 }
295
297 };
298
299 auto calculate_traction = [&] {
301 CHKERR VecZeroEntries(CommonData::totalTraction);
302 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
303 CHKERR VecAssemblyBegin(CommonData::totalTraction);
304 CHKERR VecAssemblyEnd(CommonData::totalTraction);
306 };
307
308 auto calculate_reactions = [&]() {
310
311 auto res = createDMVector(dM);
312
313 auto assemble_domain = [&]() {
315 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
316 auto &pip = fe_rhs->getOpPtrVector();
317 fe_rhs->f = res;
318
319 auto integration_rule = [](int, int, int approx_order) {
320 return 2 * approx_order + geom_order - 1;
321 };
322 fe_rhs->getRuleHook = integration_rule;
323 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
324 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
325 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
327 };
328
329 auto assemble_boundary = [&]() {
331 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
332 auto &pip = fe_rhs->getOpPtrVector();
333 fe_rhs->f = res;
334
335 auto integration_rule = [](int, int, int approx_order) {
336 return 2 * approx_order + geom_order - 1;
337 };
338 fe_rhs->getRuleHook = integration_rule;
339
340 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
341 "GEOMETRY");
342 // We have to integrate on curved face geometry, thus integration weight
343 // have to adjusted.
344 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
345
346 auto u_disp = boost::make_shared<MatrixDouble>();
347 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
348 pip.push_back(
349 new OpSpringRhs("U", u_disp, [this](double, double, double) {
350 return spring_stiffness;
351 }));
352
353 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
354
356 };
357
358 CHKERR assemble_domain();
359 CHKERR assemble_boundary();
360
361 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
362 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
363 m_field_ptr]() {
365 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
366 *m_field_ptr, fe_post_proc_ptr, res)();
368 };
369 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
370 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
371
373 };
374
375 auto print_max_min = [&](auto &tuple, const std::string msg) {
377 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
378 INSERT_VALUES, SCATTER_FORWARD);
379 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
380 INSERT_VALUES, SCATTER_FORWARD);
381 double max, min;
382 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
383 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
384 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e min %3.4e max %3.4e",
385 msg.c_str(), ts_t, min, max);
387 };
388
389 auto print_traction = [&](const std::string msg) {
391 MoFEM::Interface *m_field_ptr;
392 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
393 if (!m_field_ptr->get_comm_rank()) {
394 const double *t_ptr;
395 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
396 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
397 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
398 if (atomTest == 1 && fabs(ts_t - 1.0) < 1e-3) {
399 double hertz_tract = 158.73;
400 double tol = 4e-3;
401 if (fabs(t_ptr[1] - hertz_tract) / hertz_tract > tol) {
402 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
403 "atom test %d diverged! %3.4e != %3.4e", atom_test,
404 t_ptr[1], hertz_tract);
405 }
406 }
407 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
408 }
410 };
411
412 int se = 1;
413 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &se, PETSC_NULL);
414
415 if (!(ts_step % se)) {
416 MOFEM_LOG("CONTACT", Sev::inform)
417 << "Write file at time " << ts_t << " write step " << sTEP;
418 CHKERR post_proc();
419 }
420 CHKERR calculate_traction();
421 CHKERR calculate_reactions();
422
423 CHKERR print_max_min(uXScatter, "Ux");
424 CHKERR print_max_min(uYScatter, "Uy");
425 if (SPACE_DIM == 3)
426 CHKERR print_max_min(uZScatter, "Uz");
427 CHKERR print_traction("Contact force");
428
429 ++sTEP;
430
432 }
433
434 MoFEMErrorCode setScatterVectors(
435 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
436 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
437 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter) {
439 uXScatter = ux_scatter;
440 uYScatter = uy_scatter;
441 uZScatter = uz_scatter;
443 }
444
445private:
446 SmartPetscObj<DM> dM;
447 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
448 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
449 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
450
451 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
452
453 boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
454 boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
455
456 boost::shared_ptr<BoundaryEle> integrateTraction;
457
459 moab::Interface &moabVertex;
460
461 double lastTime;
462 double deltaTime;
463 int sTEP;
464
467 boost::shared_ptr<GenericElementInterface> mfrontInterface;
469};
470
471} // namespace ContactOps
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
constexpr int SPACE_DIM
double spring_stiffness
Definition: contact.cpp:130
int atom_test
Definition: contact.cpp:139
PetscBool use_mfront
Definition: contact.cpp:136
PetscBool is_axisymmetric
Definition: contact.cpp:137
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#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
auto integration_rule
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
double tol
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition: plastic.cpp:170
int geom_order
Order if fixed.
Definition: plastic.cpp:188
static constexpr int approx_order
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode setScatterVectors(std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter)
MoFEMErrorCode postProcess()
Monitor(SmartPetscObj< DM > &dm, bool use_mfront=false, boost::shared_ptr< GenericElementInterface > mfront_interface=nullptr, bool is_axisymmetric=false, int atom_test=0)
SmartPetscObj< DM > dM
boost::shared_ptr< GenericElementInterface > mfrontInterface
boost::shared_ptr< BoundaryEle > integrateTraction
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
MoFEMErrorCode preProcess()
moab::Interface & moabVertex
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< moab::Core > postProcMesh
MoFEMErrorCode operator()()
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomain
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleDomain
virtual int get_comm_rank() const =0
Deprecated interface functions.
Post post-proc data at points from hash maps.