v0.13.2
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,
31 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
32 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
33 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
34 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
35 uZScatter(uz_scatter), moabVertex(mbVertexPostproc), sTEP(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 auto push_domain_ops = [&](auto &pip) {
43 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
44 pip, {H1, HDIV}, "GEOMETRY")),
45 "Apply base transform");
46 auto henky_common_data_ptr =
47 commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
48 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
49 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
50 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
51 "SIGMA", contact_stress_ptr));
52 return std::make_tuple(henky_common_data_ptr, contact_stress_ptr);
53 };
54
55 auto push_bdy_ops = [&](auto &pip, int space_dim) {
56 if (space_dim == 2) {
57 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
58 pip, {HDIV}, "GEOMETRY")),
59 "Apply transform");
60 // evaluate traction
61 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
62 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
63 "SIGMA", common_data_ptr->contactTractionPtr()));
64 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
65 "U", common_data_ptr->contactDispPtr()));
66 return common_data_ptr;
67 }
68 return boost::shared_ptr<ContactOps::CommonData>();
69 };
70
71 auto get_domain_pip = [&](auto &pip)
72 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
73 if constexpr (SPACE_DIM == 3) {
74 auto op_loop_side = new OpLoopSide<SideEle>(
75 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
76 boost::make_shared<
77 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
78 pip.push_back(op_loop_side);
79 return op_loop_side->getOpPtrVector();
80 } else {
81 return pip;
82 }
83 };
84
85 auto get_post_proc_domain_fe = [&]() {
86 auto post_proc_fe =
87 boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
88 auto &pip = post_proc_fe->getOpPtrVector();
89
90 auto common_data_ptr = push_bdy_ops(pip, SPACE_DIM - 1);
91 auto [henky_common_data_ptr, contact_stress_ptr] =
92 push_domain_ops(get_domain_pip(pip));
93
94 auto u_ptr = boost::make_shared<MatrixDouble>();
95 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
96 auto X_ptr = boost::make_shared<MatrixDouble>();
97 pip.push_back(
98 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
99
100 pip.push_back(
101
102 new OpPPMap(
103
104 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
105
106 {},
107
108 {
109
110 {"U", u_ptr},
111 {"GEOMETRY", X_ptr},
112
113 // Note: post-process tractions in 3d, i.e. when mesh is
114 // post-process on skin
115 {"TRACTION_CONTACT",
116 (common_data_ptr) ? common_data_ptr->contactTractionPtr()
117 : nullptr}
118
119 },
120
121 {
122
123 {"SIGMA", contact_stress_ptr},
124
125 {"G", henky_common_data_ptr->matGradPtr},
126
127 {"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
128
129 },
130
131 {}
132
133 )
134
135 );
136
137 return post_proc_fe;
138 };
139
140 auto get_post_proc_bdy_fe = [&]() {
141 auto post_proc_fe =
142 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
143 auto &pip = post_proc_fe->getOpPtrVector();
144
145 auto common_data_ptr = push_bdy_ops(pip, SPACE_DIM);
146 // create OP which run element on side
147 auto op_loop_side = new OpLoopSide<SideEle>(
148 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
149 boost::make_shared<
150 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
151 pip.push_back(op_loop_side);
152
153 auto [henky_common_data_ptr, contact_stress_ptr] =
154 push_domain_ops(op_loop_side->getOpPtrVector());
155
156 auto X_ptr = boost::make_shared<MatrixDouble>();
157 pip.push_back(
158 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
159
160 pip.push_back(
161
162 new OpPPMap(
163
164 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
165
166 {},
167
168 {{"U", common_data_ptr->contactDispPtr()},
169 {"GEOMETRY", X_ptr},
170 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()}},
171
172 {},
173
174 {}
175
176 )
177
178 );
179
180 return post_proc_fe;
181 };
182
183 auto get_integrate_traction = [&]() {
184 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
185 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
187 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
188 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
189 "Apply transform");
190 // We have to integrate on curved face geometry, thus integration weight
191 // have to adjusted.
192 integrate_traction->getOpPtrVector().push_back(
193 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
194 integrate_traction->getRuleHook = [](int, int, int approx_order) {
195 return 2 * approx_order + geom_order - 1;
196 };
197
199 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
200 integrate_traction->getOpPtrVector(), "SIGMA")),
201 "push operators to calculate traction");
202
203 return integrate_traction;
204 };
205
206 postProcDomainFe = get_post_proc_domain_fe();
207 if constexpr (SPACE_DIM == 2)
208 postProcBdyFe = get_post_proc_bdy_fe();
209 integrateTraction = get_integrate_traction();
210 }
211
212 MoFEMErrorCode preProcess() { return 0; }
213 MoFEMErrorCode operator()() { return 0; }
214
215 MoFEMErrorCode postProcess() {
217 MoFEM::Interface *m_field_ptr;
218 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
219
220 auto post_proc = [&]() {
222
223 auto post_proc_begin =
224 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
225 postProcMesh);
226 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
227 *m_field_ptr, postProcMesh);
228
229 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
230 if (!postProcBdyFe) {
231 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
232 } else {
233 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
234 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
235 }
236 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
237
238 CHKERR post_proc_end->writeFile(
239 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
241 };
242
243 auto calculate_traction = [&] {
245 CHKERR VecZeroEntries(CommonData::totalTraction);
246 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
247 CHKERR VecAssemblyBegin(CommonData::totalTraction);
248 CHKERR VecAssemblyEnd(CommonData::totalTraction);
250 };
251
252 auto calculate_reactions = [&]() {
254
255 auto res = createDMVector(dM);
256
257 auto assemble_domain = [&]() {
259 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
260 auto &pip = fe_rhs->getOpPtrVector();
261 fe_rhs->f = res;
262
263 auto integration_rule = [](int, int, int approx_order) {
264 return 2 * approx_order + geom_order - 1;
265 };
266 fe_rhs->getRuleHook = integration_rule;
267
268 // only in case of dynamics
269 if (!is_quasi_static) {
270 auto mat_acceleration = boost::make_shared<MatrixDouble>();
271 pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
272 "U", mat_acceleration));
273 pip.push_back(
274 new OpInertiaForce("U", mat_acceleration,
275 [](double, double, double) { return rho; }));
276 }
277 if (alpha_damping > 0) {
278 auto mat_velocity = boost::make_shared<MatrixDouble>();
279 pip.push_back(new OpCalculateVectorFieldValuesDot<SPACE_DIM>(
280 "U", mat_velocity));
281 pip.push_back(
282 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
283 return alpha_damping;
284 }));
285 }
286
287 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
288 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
289 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
291 };
292
293 auto assemble_boundary = [&]() {
295 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
296 auto &pip = fe_rhs->getOpPtrVector();
297 fe_rhs->f = res;
298
299 auto integration_rule = [](int, int, int approx_order) {
300 return 2 * approx_order + geom_order - 1;
301 };
302 fe_rhs->getRuleHook = integration_rule;
303
304 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
305 "GEOMETRY");
306 // We have to integrate on curved face geometry, thus integration weight
307 // have to adjusted.
308 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
309
310 auto u_disp = boost::make_shared<MatrixDouble>();
311 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
312 pip.push_back(
313 new OpSpringRhs("U", u_disp, [this](double, double, double) {
314 return spring_stiffness;
315 }));
316
318 };
319
320 CHKERR assemble_domain();
321
322 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
323 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
324 m_field_ptr]() {
326 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
327 *m_field_ptr, fe_post_proc_ptr, res)();
329 };
330 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
331 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
332
334 };
335
336 auto print_max_min = [&](auto &tuple, const std::string msg) {
338 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
339 INSERT_VALUES, SCATTER_FORWARD);
340 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
341 INSERT_VALUES, SCATTER_FORWARD);
342 double max, min;
343 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
344 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
345 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e min %3.4e max %3.4e",
346 msg.c_str(), ts_t, min, max);
348 };
349
350 auto print_traction = [&](const std::string msg) {
352 MoFEM::Interface *m_field_ptr;
353 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
354 if (!m_field_ptr->get_comm_rank()) {
355 const double *t_ptr;
356 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
357 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
358 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
359 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
360 }
362 };
363
364 MOFEM_LOG("CONTACT", Sev::inform)
365 << "Write file at time " << ts_t << " write step " << sTEP;
366
367 CHKERR post_proc();
368 CHKERR calculate_traction();
369 CHKERR calculate_reactions();
370
371 CHKERR print_max_min(uXScatter, "Ux");
372 CHKERR print_max_min(uYScatter, "Uy");
373 if (SPACE_DIM == 3)
374 CHKERR print_max_min(uZScatter, "Uz");
375 CHKERR print_traction("Contact force");
376
377 ++sTEP;
378
380 }
381
382private:
383 SmartPetscObj<DM> dM;
384 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
385 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
386 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
387
388 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
389
390 boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
391 boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
392
393 boost::shared_ptr<BoundaryEle> integrateTraction;
394
396 moab::Interface &moabVertex;
397
398 double lastTime;
399 double deltaTime;
400 int sTEP;
401};
402
403} // namespace ContactOps
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
constexpr int SPACE_DIM
double spring_stiffness
Definition: contact.cpp:130
double rho
Definition: contact.cpp:129
double alpha_damping
Definition: contact.cpp:131
constexpr bool is_quasi_static
Definition: contact.cpp:123
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ H1
continuous field
Definition: definitions.h:85
@ 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
#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
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
[Only used for dynamics]
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int geom_order
Order if fixed.
Definition: plastic.cpp:116
static constexpr int approx_order
moab::Interface & moabVertex
MoFEMErrorCode postProcess()
boost::shared_ptr< BoundaryEle > integrateTraction
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
Monitor(SmartPetscObj< DM > &dm, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter)
SmartPetscObj< DM > dM
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< moab::Core > postProcMesh
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode preProcess()
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
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.