v0.15.0
Loading...
Searching...
No Matches
ContactOps::Monitor Struct Reference

#include "tutorials/adv-1/src/PostProcContact.hpp"

Inheritance diagram for ContactOps::Monitor:
[legend]
Collaboration diagram for ContactOps::Monitor:
[legend]

Public Member Functions

 Monitor (SmartPetscObj< DM > &dm, double scale, boost::shared_ptr< GenericElementInterface > mfront_interface=nullptr, bool is_axisymmetric=false, bool is_large_strain=false)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode operator() ()
 
MoFEMErrorCode postProcess ()
 
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 getErrorNorm (int normType)
 

Public Attributes

MoFEM::ScalarFun analyticalWavy2DPressure
 [Analytical function]
 
MoFEM::ScalarFun analyticalHertzPressureAxisymmetric
 
MoFEM::ScalarFun analyticalHertzPressurePlaneStrain
 
MoFEM::ScalarFun analyticalHertzPressurePlaneStress
 
MoFEM::VectorFun< SPACE_DIManalyticalHertzDisplacement3D
 

Private Types

enum  NORMS { TRACTION_NORM_L2 = 0 , MAG_TRACTION_NORM_L2 , TRACTION_Y_NORM_L2 , LAST_NORM }
 

Private Attributes

SmartPetscObj< DM > dM
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
boost::shared_ptr< moab::Core > postProcMesh = boost::make_shared<moab::Core>()
 
boost::shared_ptr< PostProcEleDomainpostProcDomainFe
 
boost::shared_ptr< PostProcEleBdypostProcBdyFe
 
boost::shared_ptr< BoundaryEleintegrateTraction
 
boost::shared_ptr< BoundaryEleintegrateArea
 
SmartPetscObj< Vec > normsVec
 
moab::Core mbVertexPostproc
 
moab::Interface & moabVertex
 
double lastTime
 
double deltaTime
 
int sTEP
 
bool isLargeStrain
 
boost::shared_ptr< GenericElementInterfacemfrontInterface
 

Detailed Description

Definition at line 28 of file PostProcContact.hpp.

Member Enumeration Documentation

◆ NORMS

Enumerator
TRACTION_NORM_L2 
MAG_TRACTION_NORM_L2 
TRACTION_Y_NORM_L2 
LAST_NORM 

Definition at line 873 of file PostProcContact.hpp.

Constructor & Destructor Documentation

◆ Monitor()

ContactOps::Monitor::Monitor ( SmartPetscObj< DM > & dm,
double scale,
boost::shared_ptr< GenericElementInterface > mfront_interface = nullptr,
bool is_axisymmetric = false,
bool is_large_strain = false )
inline

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.

  • For HookeOps::CommonData, maps Cauchy stress and small strain.
  • For HenckyOps::CommonData, maps first Piola-Kirchhoff stress and Hencky strain.
    Returns
    A commod data pointer to the configured post-processing hencky_or_hooke_common_data_ptr.

Definition at line 30 of file PostProcContact.hpp.

35
36 MoFEM::Interface *m_field_ptr;
37 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
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)
43 : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
44 mPtr(m_ptr), scale(s) {}
45 MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
46 *mPtr *= 1./scale;
47 return 0;
48 }
49
50 private:
51 boost::shared_ptr<MatrixDouble> mPtr;
52 double scale;
53 };
54
55 auto push_domain_ops = [&](auto &pip) {
57 pip, {H1, HDIV}, "GEOMETRY")),
58 "Apply base transform");
59 // Define variant type that holds either HookeOps or HenckyOps
60 // CommonData.
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
68 if (!isLargeStrain) { // use HookeOps
69 auto hooke_common_data_ptr =
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 { // use HenckyOps
77 auto hencky_common_data_ptr =
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 // evaluate traction
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> & {
101 if constexpr (SPACE_DIM == 3) {
102 auto op_loop_side = new OpLoopSide<SideEle>(
103 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
104 boost::make_shared<
105 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
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 // Get common data pointer (either Hooke or Hencky)
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 * @brief Adds operators for evaluating stress and strain based on the material model: i.e.
130 * Visit the variant to handle either HookeOps::CommonData or HenckyOps::CommonData. i.e.
131 * - For `HookeOps::CommonData`, maps Cauchy stress and small strain.
132 * - For `HenckyOps::CommonData`, maps first Piola-Kirchhoff stress and Hencky strain.
133 * @return A commod data pointer to the configured post-processing hencky_or_hooke_common_data_ptr.
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>>) {
140 pip.push_back(new OpPPMap(
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}}, //Displacement gradient
146 {{"STRESS", common_data_ptr->getMatCauchyStress()},
147 {"STRAIN", common_data_ptr->getMatStrain()}}));
148 }
149 // Handle HenckyOps::CommonData
150 else if constexpr (std::is_same_v<
151 std::decay_t<decltype(common_data_ptr)>,
152 boost::shared_ptr<
154
155 pip.push_back(new OpPPMap(
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
169 if (SPACE_DIM == 3) {
170
172 pip, {HDIV}, "GEOMETRY")),
173 "Apply transform");
174 auto common_data_ptr = push_bdy_ops(pip);
175
176 pip.push_back(
177
178 new OpPPMap(
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 // create OP which run element on side
215 auto op_loop_side = new OpLoopSide<SideEle>(
216 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
217 boost::make_shared<
218 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
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
230 new OpPPMap(
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 // We have to integrate on curved face geometry, thus integration weight have to adjusted.
262 integrate_traction->getOpPtrVector().push_back(
263 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
264 integrate_traction->getRuleHook = [](int, int, int approx_order) {
265 return 2 * approx_order + geom_order - 1;
266 };
267
270 integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
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 // We have to integrate on curved face geometry, thus integration weight
284 // have to adjusted.
285 integrate_area->getOpPtrVector().push_back(
286 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
287 integrate_area->getRuleHook = [](int, int, int approx_order) {
288 return 2 * approx_order + geom_order - 1;
289 };
290 Range contact_range;
291 for (auto m :
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;
296 CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
297 meshset, SPACE_DIM - 1, contact_meshset_range, true);
298
299 CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
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(),
308 SPACE_DIM);
310 op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
311
314 integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
315 is_axisymmetric, contact_range_ptr)),
316 "push operators to calculate area");
317
318 return integrate_area;
319 };
320
321 postProcDomainFe = get_post_proc_domain_fe();
322 if constexpr (SPACE_DIM == 2)
323 postProcBdyFe = get_post_proc_bdy_fe();
324
325 integrateTraction = get_integrate_traction();
326 integrateArea = get_integrate_area();
327
329 m_field_ptr->get_comm(),
330 (m_field_ptr->get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
331 }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
PetscBool is_large_strain
Definition contact.cpp:94
PetscBool is_axisymmetric
Definition contact.cpp:93
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:410
MoFEMErrorCode opFactoryCalculateArea(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, OpLoopSide< SideEle > *op_loop_side, std::string sigma, std::string u, bool is_axisymmetric=false, boost::shared_ptr< Range > contact_range_ptr=nullptr)
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
Definition HookeOps.hpp:208
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition plastic.cpp:123
int geom_order
Order if fixed.
Definition plastic.cpp:141
static constexpr int approx_order
FTensor::Index< 'm', 3 > m
moab::Interface & moabVertex
boost::shared_ptr< BoundaryEle > integrateTraction
SmartPetscObj< Vec > normsVec
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
boost::shared_ptr< GenericElementInterface > mfrontInterface
SmartPetscObj< DM > dM
boost::shared_ptr< moab::Core > postProcMesh
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
boost::shared_ptr< BoundaryEle > integrateArea
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

Member Function Documentation

◆ getErrorNorm()

MoFEMErrorCode ContactOps::Monitor::getErrorNorm ( int normType)
inline

Definition at line 851 of file PostProcContact.hpp.

851 {
852 const double *norm;
853 CHKERR VecGetArrayRead(normsVec, &norm);
854 double norm_val = std::sqrt(norm[normType]);
855 CHKERR VecRestoreArrayRead(normsVec, &norm);
856 return norm_val;
857 }

◆ operator()()

MoFEMErrorCode ContactOps::Monitor::operator() ( )
inline

Definition at line 334 of file PostProcContact.hpp.

334{ return 0; }

◆ postProcess()

MoFEMErrorCode ContactOps::Monitor::postProcess ( )
inline

Definition at line 336 of file PostProcContact.hpp.

336 {
338 MoFEM::Interface *m_field_ptr;
339 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
340
341 auto post_proc = [&]() {
343
344 if (!mfrontInterface) {
345 auto post_proc_begin =
346 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
348 auto post_proc_end =
349 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
351
353 post_proc_begin->getFEMethod());
354 if (!postProcBdyFe) {
355 postProcDomainFe->copyTs(*this); // this here is a Monitor
357 } else {
358 postProcDomainFe->copyTs(*this); // this here is a Monitor
359 postProcBdyFe->copyTs(*this);
362 }
364 post_proc_end->getFEMethod());
365
366 CHKERR post_proc_end->writeFile(
367 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
368 } else {
369 CHKERR mfrontInterface->postProcessElement(
370 ts_step, dM,
371 m_field_ptr->getInterface<Simple>()->getDomainFEName());
372 }
373
375 };
376
377 auto calculate_force = [&] {
379 CHKERR VecZeroEntries(CommonData::totalTraction);
381 CHKERR VecAssemblyBegin(CommonData::totalTraction);
382 CHKERR VecAssemblyEnd(CommonData::totalTraction);
384 };
385
386 auto calculate_area = [&] {
388 integrateArea->copyTs(*this);
390 CHKERR VecAssemblyBegin(CommonData::totalTraction);
391 CHKERR VecAssemblyEnd(CommonData::totalTraction);
393 };
394
395 auto calculate_reactions = [&]() {
397
398 auto res = createDMVector(dM);
399
400 auto assemble_domain = [&]() {
402 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
403 auto &pip = fe_rhs->getOpPtrVector();
404 fe_rhs->f = res;
405
406 auto integration_rule = [](int, int, int approx_order) {
407 return 2 * approx_order + geom_order - 1;
408 };
409 fe_rhs->getRuleHook = integration_rule;
410
412 "GEOMETRY");
413 CHKERR
415 pip, "SIGMA", "U", is_axisymmetric);
416
417 if (!mfrontInterface) {
418 CHKERR
420 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
421 } else {
422 CHKERR mfrontInterface->opFactoryDomainRhs(pip);
423 }
424 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
425
426 CHKERR VecAssemblyBegin(res);
427 CHKERR VecAssemblyEnd(res);
428 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
429 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
430
432 };
433
434 auto assemble_boundary = [&]() {
436 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
437 auto &pip = fe_rhs->getOpPtrVector();
438 fe_rhs->f = res;
439
440 auto integration_rule = [](int, int, int approx_order) {
441 return 2 * approx_order + geom_order - 1;
442 };
443 fe_rhs->getRuleHook = integration_rule;
444
446 "GEOMETRY");
447 // We have to integrate on curved face geometry, thus integration weight
448 // have to adjusted.
449 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
450
451 auto u_disp = boost::make_shared<MatrixDouble>();
452 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
453 pip.push_back(
454 new OpSpringRhs("U", u_disp, [this](double, double, double) {
455 return spring_stiffness;
456 }));
457
458 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
459
461 };
462
463 CHKERR assemble_domain();
464 CHKERR assemble_boundary();
465
466 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
467 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
468 m_field_ptr]() {
470 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
471 *m_field_ptr, fe_post_proc_ptr, res)();
473 };
474 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
475 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
476
478 };
479
480 auto print_max_min = [&](auto &tuple, const std::string msg) {
482 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
483 INSERT_VALUES, SCATTER_FORWARD);
484 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
485 INSERT_VALUES, SCATTER_FORWARD);
486 double max, min;
487 CHKERR VecMax(std::get<0>(tuple), PETSC_NULLPTR, &max);
488 CHKERR VecMin(std::get<0>(tuple), PETSC_NULLPTR, &min);
489 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %6.4e min %6.4e max %6.4e",
490 msg.c_str(), ts_t, min, max);
492 };
493
494 auto print_force_and_area = [&]() {
496 MoFEM::Interface *m_field_ptr;
497 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
498 if (!m_field_ptr->get_comm_rank()) {
499 const double *t_ptr;
500 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
501 MOFEM_LOG_C("CONTACT", Sev::inform,
502 "Contact force: time %6.3e Fx: %6.6e Fy: %6.6e Fz: %6.6e",
503 ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
504 MOFEM_LOG_C("CONTACT", Sev::inform,
505 "Contact area: time %6.3e Active: %6.6e Potential: %6.6e",
506 ts_t, t_ptr[3], t_ptr[4]);
507 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
508 }
510 };
511
512 if (mfrontInterface) {
513 CHKERR mfrontInterface->updateElementVariables(
514 dM, m_field_ptr->getInterface<Simple>()->getDomainFEName());
515 }
516
517 auto calculate_error = [&](MoFEM::ScalarFun &fun) {
519 struct OpCalcTractions : public BoundaryEleOp {
520 OpCalcTractions(boost::shared_ptr<MatrixDouble> m_ptr,
521 boost::shared_ptr<VectorDouble> p_ptr,
522 boost::shared_ptr<VectorDouble> mag_ptr,
523 boost::shared_ptr<VectorDouble> traction_y_ptr,
524 boost::shared_ptr<MatrixDouble> t_ptr,
525 boost::shared_ptr<MatrixDouble> grad_sdf_ptr)
526 : BoundaryEleOp(NOSPACE, OPSPACE), mPtr(m_ptr), pPtr(p_ptr),
527 magPtr(mag_ptr), tyPtr(traction_y_ptr), tPtr(t_ptr),
528 gradSDFPtr(grad_sdf_ptr) {}
529 MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
532 mPtr->resize(SPACE_DIM, pPtr->size());
533 mPtr->clear();
534 magPtr->resize(pPtr->size());
535 magPtr->clear();
536 tyPtr->resize(pPtr->size());
537 tyPtr->clear();
538
539 auto t_traction = getFTensor1FromMat<SPACE_DIM>(*mPtr);
540 auto t_contact_traction = getFTensor1FromMat<SPACE_DIM>(*tPtr);
541 auto t_p = getFTensor0FromVec(*pPtr);
542 int nb_gauss_pts = pPtr->size();
543 auto t_normal = getFTensor1FromMat<SPACE_DIM>(*gradSDFPtr);
544 auto t_normal_at_gauss = getFTensor1NormalsAtGaussPts();
545 auto t_mag = getFTensor0FromVec(*magPtr);
546 auto t_ty = getFTensor0FromVec(*tyPtr);
547
548 for (int gg = 0; gg != nb_gauss_pts; gg++) {
550 t_traction(i) = t_p * (-(t_normal(i) / t_normal.l2()));
551 t_mag = t_contact_traction.l2();
552 t_ty = t_contact_traction(1);
553
554 ++t_normal;
555 ++t_traction;
556 ++t_p;
557 ++t_mag;
558 ++t_contact_traction;
559 ++t_ty;
560 ++t_normal_at_gauss;
561 }
563 }
564
565 private:
566 boost::shared_ptr<MatrixDouble> mPtr;
567 boost::shared_ptr<VectorDouble> pPtr;
568 boost::shared_ptr<VectorDouble> magPtr;
569 boost::shared_ptr<MatrixDouble> tPtr;
570 boost::shared_ptr<VectorDouble> tyPtr;
571 boost::shared_ptr<MatrixDouble> gradSDFPtr;
572 };
573
574 auto post_proc_norm_fe = boost::make_shared<BoundaryEle>(*m_field_ptr);
575 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
576 Range contact_range;
577 for (auto m :
578 m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
579 std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
580 auto meshset = m->getMeshset();
581 Range contact_meshset_range;
582 CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
583 meshset, SPACE_DIM - 1, contact_meshset_range, true);
584
585 CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
586 contact_meshset_range);
587 contact_range.merge(contact_meshset_range);
588 }
589
592 post_proc_norm_fe->getOpPtrVector(), {HDIV}, "GEOMETRY")),
593 "Apply transform");
594 // We have to integrate on curved face geometry, thus integration weight
595 // have to adjusted.
596 post_proc_norm_fe->getOpPtrVector().push_back(
597 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
598 post_proc_norm_fe->getRuleHook = [](int, int, int approx_order) {
599 return 2 * approx_order + geom_order - 1;
600 };
601
602 post_proc_norm_fe->getOpPtrVector().push_back(
603 new OpCalculateVectorFieldValues<SPACE_DIM>(
604 "U", common_data_ptr->contactDispPtr()));
605 post_proc_norm_fe->getOpPtrVector().push_back(
606 new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
607 "SIGMA", common_data_ptr->contactTractionPtr()));
608 using C = ContactIntegrators<BoundaryEleOp>;
609 post_proc_norm_fe->getOpPtrVector().push_back(
610 new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
611 common_data_ptr));
612
613 auto analytical_traction_ptr = boost::make_shared<MatrixDouble>();
614 auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
615 auto mag_traction_ptr = boost::make_shared<VectorDouble>();
616 auto traction_y_ptr = boost::make_shared<VectorDouble>();
617 auto contact_range_ptr = boost::make_shared<Range>(contact_range);
618
619 post_proc_norm_fe->getOpPtrVector().push_back(
620 new OpGetTensor0fromFunc(analytical_pressure_ptr, fun));
621
622 post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcTractions(
623 analytical_traction_ptr, analytical_pressure_ptr, mag_traction_ptr,
624 traction_y_ptr, common_data_ptr->contactTractionPtr(),
625 common_data_ptr->gradSdfPtr()));
626
627 post_proc_norm_fe->getOpPtrVector().push_back(
628 new OpCalcNormL2Tensor1<SPACE_DIM>(
629 common_data_ptr->contactTractionPtr(), normsVec, TRACTION_NORM_L2,
630 analytical_traction_ptr, contact_range_ptr));
631
632 // calculate magnitude of traction
633
634 post_proc_norm_fe->getOpPtrVector().push_back(new OpCalcNormL2Tensor0(
635 mag_traction_ptr, normsVec, MAG_TRACTION_NORM_L2,
636 analytical_pressure_ptr, contact_range_ptr));
637
638 post_proc_norm_fe->getOpPtrVector().push_back(
639 new OpCalcNormL2Tensor0(traction_y_ptr, normsVec, TRACTION_Y_NORM_L2,
640 analytical_pressure_ptr, contact_range_ptr));
641
642 CHKERR VecZeroEntries(normsVec);
643 post_proc_norm_fe->copyTs(*this); // set time as is in Monitor
644 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", post_proc_norm_fe);
645 CHKERR VecAssemblyBegin(normsVec);
646 CHKERR VecAssemblyEnd(normsVec);
647
648 MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
649 if (m_field_ptr->get_comm_rank() == 0) {
650 const double *norms;
651 CHKERR VecGetArrayRead(normsVec, &norms);
652 MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
653 << "norm_traction: " << std::scientific
654 << std::sqrt(norms[TRACTION_NORM_L2]);
655 MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
656 << "norm_mag_traction: " << std::scientific
657 << std::sqrt(norms[MAG_TRACTION_NORM_L2]);
658 MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
659 << "norm_traction_y: " << std::scientific
660 << std::sqrt(norms[TRACTION_Y_NORM_L2]);
661 CHKERR VecRestoreArrayRead(normsVec, &norms);
662 }
664 };
665
666 int se = 1;
667 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-save_every", &se, PETSC_NULLPTR);
668
669 if (!(ts_step % se)) {
670 MOFEM_LOG("CONTACT", Sev::inform)
671 << "Write file at time " << ts_t << " write step " << sTEP;
672 CHKERR post_proc();
673 }
674 CHKERR calculate_force();
675 CHKERR calculate_area();
676
677 CHKERR calculate_reactions();
678
679 if (atom_test && sTEP) {
680 switch (atom_test) {
681 case 1:
683 break;
684 case 2:
686 break;
687 case 5:
689 break;
690 case 6:
691 CHKERR calculate_error(analyticalWavy2DPressure);
692 break;
693 default:
694 break;
695 }
696 }
697
698 CHKERR print_max_min(uXScatter, "Ux");
699 CHKERR print_max_min(uYScatter, "Uy");
700 if (SPACE_DIM == 3)
701 CHKERR print_max_min(uZScatter, "Uz");
702 CHKERR print_force_and_area();
703 ++sTEP;
704
706 }
#define MOFEM_LOG_C(channel, severity, format,...)
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition contact.cpp:74
double spring_stiffness
Definition contact.cpp:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
auto fun
Function to approximate.
auto integration_rule
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
int atom_test
Atom test.
Definition plastic.cpp:121
static SmartPetscObj< Vec > totalTraction
MoFEM::ScalarFun analyticalWavy2DPressure
[Analytical function]
MoFEM::ScalarFun analyticalHertzPressurePlaneStress
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
MoFEM::ScalarFun analyticalHertzPressureAxisymmetric
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEM::ScalarFun analyticalHertzPressurePlaneStrain

◆ preProcess()

MoFEMErrorCode ContactOps::Monitor::preProcess ( )
inline

Definition at line 333 of file PostProcContact.hpp.

333{ return 0; }

◆ setScatterVectors()

MoFEMErrorCode ContactOps::Monitor::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 )
inline

Definition at line 840 of file PostProcContact.hpp.

843 {
845 uXScatter = ux_scatter;
846 uYScatter = uy_scatter;
847 uZScatter = uz_scatter;
849 }

Member Data Documentation

◆ analyticalHertzDisplacement3D

MoFEM::VectorFun<SPACE_DIM> ContactOps::Monitor::analyticalHertzDisplacement3D

Definition at line 783 of file PostProcContact.hpp.

785 {
786 // update to atom test values
787 double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
788 // Radius
789 double R = 100.;
790 // Contact area radius
791 double a = 1;
792 // max pressure
793 double p_0 = (2. * E_star * a) / (M_PI * R);
794 // current radius
795 double r = std::sqrt((x * x) + (y * y));
797 std::vector<double> v_u;
798
799 double u_z = 0.;
800 double u_r = 0.;
801 // outside contact zone
802 if (r > a) {
803 u_z = (1. - std::pow(poisson_ratio, 2.)) / young_modulus *
804 ((p_0) / (2. * a)) *
805 ((2. * std::pow(a, 2.) - std::pow(r, 2.)) * asin(a / r) +
806 std::pow(r, 2.) * (a / r) *
807 std::pow(1 - (std::pow(a, 2.) / std::pow(r, 2.)), 2.));
808 u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
809 (3. * young_modulus) * ((std::pow(a, 2) / r)) * p_0;
810
811 if (SPACE_DIM == 2)
812 v_u = {u_r, u_z};
813 else
814 v_u = {u_r, u_z, u_r};
815
816 for (int i = 0; i < SPACE_DIM; ++i)
817 u(i) = v_u[i];
818
819 return u;
820 }
821
822 // In contact zone
823 u_z = ((1. - std::pow(poisson_ratio, 2.)) / young_modulus) *
824 ((M_PI * p_0) / 4. * a) * (2. * std::pow(a, 2.) - std::pow(r, 2.));
825 u_r = -((1. - 2. * poisson_ratio) * (1. + poisson_ratio)) /
826 (3. * young_modulus) * ((std::pow(a, 2.) / r)) * p_0 *
827 (1 - std::pow(1 - (std::pow(r, 2.) / std::pow(a, 2.)), 1.5));
828
829 if (SPACE_DIM == 2)
830 v_u = {u_r, u_z};
831 else
832 v_u = {u_r, u_z, u_r};
833
834 for (int i = 0; i < SPACE_DIM; ++i)
835 u(i) = v_u[i];
836
837 return u;
838 };
constexpr double a
@ R
int r
Definition sdf.py:8
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126

◆ analyticalHertzPressureAxisymmetric

MoFEM::ScalarFun ContactOps::Monitor::analyticalHertzPressureAxisymmetric
Initial value:
= [](double x, double y,
double z) {
double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
double R = 100.;
double d = 0.01;
double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
double p_max = (3. * F) / (2. * M_PI * a * a);
double r = std::sqrt((x * x) + (y * y));
if (r > a) {
return 0.;
}
return p_max * std::sqrt(1 - ((r * r) / (a * a)));
}
@ F

Definition at line 723 of file PostProcContact.hpp.

724 {
725 // update to atom test values
726 double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
727 // Radius
728 double R = 100.;
729 // Indentation
730 double d = 0.01;
731 // Force
732 double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
733 // Contact area radius
734 double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
735 // Maximum pressure
736 double p_max = (3. * F) / (2. * M_PI * a * a);
737
738 double r = std::sqrt((x * x) + (y * y));
739
740 if (r > a) {
741 return 0.;
742 }
743 // Pressure = p_max * sqrt(1 - (r^2 / a^2))
744 return p_max * std::sqrt(1 - ((r * r) / (a * a)));
745 };

◆ analyticalHertzPressurePlaneStrain

MoFEM::ScalarFun ContactOps::Monitor::analyticalHertzPressurePlaneStrain
Initial value:
= [](double x, double y,
double z) {
double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
double R = 100.;
double a = 1;
double r = std::sqrt((x * x) + (y * y));
if (r > a) {
return 0.;
}
return E_star / (2. * R) * std::sqrt(a * a - r * r);
}

Definition at line 747 of file PostProcContact.hpp.

748 {
749 // update to atom test values
750 double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
751 // Radius
752 double R = 100.;
753 // Contact area radius
754 double a = 1;
755 // current radius
756 double r = std::sqrt((x * x) + (y * y));
757
758 if (r > a) {
759 return 0.;
760 }
761 // Pressure = p_max * sqrt(1 - (r^2 / a^2))
762 return E_star / (2. * R) * std::sqrt(a * a - r * r);
763 };

◆ analyticalHertzPressurePlaneStress

MoFEM::ScalarFun ContactOps::Monitor::analyticalHertzPressurePlaneStress
Initial value:
= [](double x, double y,
double z) {
double E_star = young_modulus;
double R = 100.;
double a = 1;
double r = std::sqrt((x * x) + (y * y));
if (r > a) {
return 0.;
}
return E_star / (2. * R) * std::sqrt(a * a - r * r);
}

Definition at line 764 of file PostProcContact.hpp.

765 {
766 // update to atom test values
767 double E_star = young_modulus;
768 // Radius
769 double R = 100.;
770 // Contact area radius
771 double a = 1;
772 // current radius
773 double r = std::sqrt((x * x) + (y * y));
774
775 if (r > a) {
776 return 0.;
777 }
778 // Pressure = p_max * sqrt(1 - (r^2 / a^2))
779 return E_star / (2. * R) * std::sqrt(a * a - r * r);
780 };

◆ analyticalWavy2DPressure

MoFEM::ScalarFun ContactOps::Monitor::analyticalWavy2DPressure
Initial value:
= [](double x, double y, double z) {
double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
double delta = 0.0002;
double lambda = 2;
double p_star = M_PI * E_star * delta / lambda;
return p_star + p_star * std::cos(2. * M_PI * x / lambda);
}
static double lambda
static constexpr double delta

[Analytical function]

Definition at line 709 of file PostProcContact.hpp.

709 {
710 // update to atom test values
711 double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
712 // delta
713 double delta = 0.0002;
714 // lambda
715 double lambda = 2;
716 // pressure star
717 double p_star = M_PI * E_star * delta / lambda;
718
719 // Pressure = p_bar + p_star * cos(2 * pi * x / lambda)
720 return p_star + p_star * std::cos(2. * M_PI * x / lambda);
721 };

◆ deltaTime

double ContactOps::Monitor::deltaTime
private

Definition at line 885 of file PostProcContact.hpp.

◆ dM

SmartPetscObj<DM> ContactOps::Monitor::dM
private

Definition at line 860 of file PostProcContact.hpp.

◆ integrateArea

boost::shared_ptr<BoundaryEle> ContactOps::Monitor::integrateArea
private

Definition at line 871 of file PostProcContact.hpp.

◆ integrateTraction

boost::shared_ptr<BoundaryEle> ContactOps::Monitor::integrateTraction
private

Definition at line 870 of file PostProcContact.hpp.

◆ isLargeStrain

bool ContactOps::Monitor::isLargeStrain
private

Definition at line 887 of file PostProcContact.hpp.

◆ lastTime

double ContactOps::Monitor::lastTime
private

Definition at line 884 of file PostProcContact.hpp.

◆ mbVertexPostproc

moab::Core ContactOps::Monitor::mbVertexPostproc
private

Definition at line 881 of file PostProcContact.hpp.

◆ mfrontInterface

boost::shared_ptr<GenericElementInterface> ContactOps::Monitor::mfrontInterface
private

Definition at line 889 of file PostProcContact.hpp.

◆ moabVertex

moab::Interface& ContactOps::Monitor::moabVertex
private

Definition at line 882 of file PostProcContact.hpp.

◆ normsVec

SmartPetscObj<Vec> ContactOps::Monitor::normsVec
private

Definition at line 879 of file PostProcContact.hpp.

◆ postProcBdyFe

boost::shared_ptr<PostProcEleBdy> ContactOps::Monitor::postProcBdyFe
private

Definition at line 868 of file PostProcContact.hpp.

◆ postProcDomainFe

boost::shared_ptr<PostProcEleDomain> ContactOps::Monitor::postProcDomainFe
private

Definition at line 867 of file PostProcContact.hpp.

◆ postProcMesh

boost::shared_ptr<moab::Core> ContactOps::Monitor::postProcMesh = boost::make_shared<moab::Core>()
private

Definition at line 865 of file PostProcContact.hpp.

◆ sTEP

int ContactOps::Monitor::sTEP
private

Definition at line 886 of file PostProcContact.hpp.

◆ uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactOps::Monitor::uXScatter
private

Definition at line 861 of file PostProcContact.hpp.

◆ uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactOps::Monitor::uYScatter
private

Definition at line 862 of file PostProcContact.hpp.

◆ uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactOps::Monitor::uZScatter
private

Definition at line 863 of file PostProcContact.hpp.


The documentation for this struct was generated from the following file: