v0.15.0
Loading...
Searching...
No Matches
contact.cpp
Go to the documentation of this file.
1/**
2 * \file contact.cpp
3 * \CONTACT contact.cpp
4 *
5 * CONTACT of contact problem
6 *
7 * @copyright Copyright (c) 2023
8 */
9
10/* The above code is a preprocessor directive in C++ that checks if the macro
11"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it is set
12to 3" */
13#ifndef EXECUTABLE_DIMENSION
14 #define EXECUTABLE_DIMENSION 3
15#endif
16
17#ifndef SCHUR_ASSEMBLE
18 #define SCHUR_ASSEMBLE 0
19#endif
20
21#include <MoFEM.hpp>
22#include <MatrixFunction.hpp>
23
24using namespace MoFEM;
25
26#ifdef ENABLE_PYTHON_BINDING
27 #include <boost/python.hpp>
28 #include <boost/python/def.hpp>
29 #include <boost/python/numpy.hpp>
30namespace bp = boost::python;
31namespace np = boost::python::numpy;
32#endif
33
34constexpr AssemblyType AT =
35 (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
36 : AssemblyType::PETSC; //< selected assembly type
38 IntegrationType::GAUSS; //< selected integration type
39
40template <int DIM> struct ElementsAndOps;
41
43 static constexpr FieldSpace CONTACT_SPACE = HCURL;
44};
45
47 static constexpr FieldSpace CONTACT_SPACE = HDIV;
48};
49
52
53constexpr int SPACE_DIM =
54 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
55
56/* The above code is defining an alias `EntData` for the type
57`EntitiesFieldData::EntData`. This is a C++ syntax for creating a new name for
58an existing type. */
61using DomainEleOp = DomainEle::UserDataOperator;
63using BoundaryEleOp = BoundaryEle::UserDataOperator;
65//! [Specialisation for assembly]
66
68
69//! [Operators used for contact]
73 IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
74//! [Operators used for contact]
75
76PetscBool is_quasi_static = PETSC_TRUE;
77
78int order = 2; //< Order of displacements in the domain
79int contact_order = 2; //< Order of displacements in contact elements
80int sigma_order = 1; //< Order of Lagrange multiplier in side elements
81int geom_order = 1;
82double young_modulus = 100;
83double poisson_ratio = 0.25;
84double rho = 0.0;
85double spring_stiffness = 0.0;
87double alpha_damping = 0;
88
89double scale = 1.;
90
91PetscBool is_axisymmetric = PETSC_FALSE; //< flag for Axisymmetric model
92PetscBool is_large_strain = PETSC_FALSE; //< flag for large strain
93
94int atom_test = 0;
95
96namespace ContactOps {
97double cn_contact = 0.1;
98}; // namespace ContactOps
99
100#include <HenckyOps.hpp>
101#include <HookeOps.hpp>
102using namespace HenckyOps;
103using namespace HookeOps;
104#include <ContactOps.hpp>
105#include <PostProcContact.hpp>
106#include <ContactNaturalBC.hpp>
107
117
118using namespace ContactOps;
119
120struct Contact {
121
123 : mField(m_field), mfrontInterfacePtr(nullptr) {}
124
126
127private:
129
136
137 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
138 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
139 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
140
141 boost::shared_ptr<MFrontInterface> mfrontInterfacePtr;
142 boost::shared_ptr<Monitor> monitorPtr;
143
144#ifdef ENABLE_PYTHON_BINDING
145 boost::shared_ptr<SDFPython> sdfPythonPtr;
146#endif
147
150 double getScale(const double time) {
151 return scale * MoFEM::TimeScale::getScale(time);
152 };
153 };
154};
155
156//! [Run problem]
167//! [Run problem]
168
169//! [Set up problem]
173
174 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
175 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-contact_order", &contact_order,
176 PETSC_NULLPTR);
177 sigma_order = std::max(order, contact_order) - 1;
178 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-sigma_order", &sigma_order,
179 PETSC_NULLPTR);
180 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
181 PETSC_NULLPTR);
182
183 MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
184 MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
185 MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
186 MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
187
188 // Select base
189 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
190 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
191 PetscInt choice_base_value = AINSWORTH;
192 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
193 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
194
196 switch (choice_base_value) {
197 case AINSWORTH:
199 MOFEM_LOG("CONTACT", Sev::inform)
200 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
201 break;
202 case DEMKOWICZ:
204 MOFEM_LOG("CONTACT", Sev::inform)
205 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
206 break;
207 default:
208 base = LASTBASE;
209 break;
210 }
211
212 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
213 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
214 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
215 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
216 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
217 SPACE_DIM);
218 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
219 SPACE_DIM);
220 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
221
222 CHKERR simple->setFieldOrder("U", order);
223 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
224
225 auto get_skin = [&]() {
226 Range body_ents;
227 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
228 Skinner skin(&mField.get_moab());
229 Range skin_ents;
230 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
231 return skin_ents;
232 };
233
234 auto filter_blocks = [&](auto skin) {
235 bool is_contact_block = false;
236 Range contact_range;
237 for (auto m :
238 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
239
240 (boost::format("%s(.*)") % "CONTACT").str()
241
242 ))
243
244 ) {
245 is_contact_block =
246 true; ///< blocs interation is collective, so that is set irrespective
247 ///< if there are entities in given rank or not in the block
248 MOFEM_LOG("CONTACT", Sev::inform)
249 << "Find contact block set: " << m->getName();
250 auto meshset = m->getMeshset();
251 Range contact_meshset_range;
252 CHKERR mField.get_moab().get_entities_by_dimension(
253 meshset, SPACE_DIM - 1, contact_meshset_range, true);
254
255 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
256 contact_meshset_range);
257 contact_range.merge(contact_meshset_range);
258 }
259 if (is_contact_block) {
260 MOFEM_LOG("SYNC", Sev::inform)
261 << "Nb entities in contact surface: " << contact_range.size();
263 skin = intersect(skin, contact_range);
264 }
265 return skin;
266 };
267
268 auto filter_true_skin = [&](auto skin) {
269 Range boundary_ents;
270 ParallelComm *pcomm =
271 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
272 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
273 PSTATUS_NOT, -1, &boundary_ents);
274 return boundary_ents;
275 };
276
277 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
278 CHKERR simple->setFieldOrder("SIGMA", 0);
279 CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
280
281 if (contact_order > order) {
282 Range ho_ents;
283
284 CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
285 moab::Interface::UNION);
286
287 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
288 CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
289 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
290 }
291
292 CHKERR simple->setUp();
293
294 auto project_ho_geometry = [&]() {
295 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
296 return mField.loop_dofs("GEOMETRY", ent_method);
297 };
298
299 PetscBool project_geometry = PETSC_TRUE;
300 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
301 &project_geometry, PETSC_NULLPTR);
302 if (project_geometry) {
303 CHKERR project_ho_geometry();
304 }
305
307} //! [Set up problem]
308
309//! [Create common data]
312
313 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strain",
314 &is_large_strain, PETSC_NULLPTR);
315 PetscBool use_mfront = PETSC_FALSE;
316 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_mfront", &use_mfront,
317 PETSC_NULLPTR);
318 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_axisymmetric",
319 &is_axisymmetric, PETSC_NULLPTR);
320 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
321 PETSC_NULLPTR);
322
323 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale,
324 PETSC_NULLPTR);
325 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
326 &young_modulus, PETSC_NULLPTR);
327 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
328 &poisson_ratio, PETSC_NULLPTR);
329 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
330 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn", &cn_contact,
331 PETSC_NULLPTR);
332 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-spring_stiffness",
333 &spring_stiffness, PETSC_NULLPTR);
334 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-vis_spring_stiffness",
335 &vis_spring_stiffness, PETSC_NULLPTR);
336 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
337 &alpha_damping, PETSC_NULLPTR);
338
339 if (!use_mfront) {
340 if (!is_large_strain) {
341 MOFEM_LOG("CONTACT", Sev::inform)
342 << "Selected material model: Hooke (small strain)";
343 } else {
344 MOFEM_LOG("CONTACT", Sev::inform)
345 << "Selected material model: Hencky (finite strain)";
346 }
347 MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
348 MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
349 } else {
350 MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
351 }
352
353 MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
354 MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
355 MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
356 MOFEM_LOG("CONTACT", Sev::inform)
357 << "Vis spring_stiffness " << vis_spring_stiffness;
358
359 MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
360
361 PetscBool use_scale = PETSC_FALSE;
362 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_scale", &use_scale,
363 PETSC_NULLPTR);
364 if (use_scale) {
366 }
367
368 MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
369
370 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_quasi_static",
371 &is_quasi_static, PETSC_NULLPTR);
372 MOFEM_LOG("CONTACT", Sev::inform)
373 << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
374
375#ifdef ENABLE_PYTHON_BINDING
376 auto file_exists = [](std::string myfile) {
377 std::ifstream file(myfile.c_str());
378 if (file) {
379 return true;
380 }
381 return false;
382 };
383 char sdf_file_name[255] = "sdf.py";
384 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
385 sdf_file_name, 255, PETSC_NULLPTR);
386
387 if (file_exists(sdf_file_name)) {
388 MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
389 sdfPythonPtr = boost::make_shared<SDFPython>();
390 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
391 sdfPythonWeakPtr = sdfPythonPtr;
392 } else {
393 MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
394 }
395#endif
396
397 if (is_axisymmetric) {
398 if (SPACE_DIM == 3) {
399 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
400 "Use executable contact_2d with axisymmetric model");
401 } else {
402 if (!use_mfront) {
403 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
404 "Axisymmetric model is only available with MFront (set "
405 "use_mfront to 1)");
406 } else {
407 MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
408 }
409 }
410 } else {
411 if (SPACE_DIM == 2) {
412 MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
413 }
414 }
415
416 if (use_mfront) {
417 if (SPACE_DIM == 3) {
419 } else if (SPACE_DIM == 2) {
420 if (is_axisymmetric) {
422 } else {
424 }
425 }
426 CHKERR mfrontInterfacePtr->getCommandLineParameters();
427 }
428
430 auto dm = simple->getDM();
431 monitorPtr = boost::make_shared<Monitor>(dm, scale, mfrontInterfacePtr,
433
434 if (mfrontInterfacePtr) {
435 mfrontInterfacePtr->setMonitor(monitorPtr);
436 }
437
439}
440//! [Create common data]
441
442//! [Boundary condition]
445 auto bc_mng = mField.getInterface<BcManager>();
447
448 for (auto f : {"U", "SIGMA"}) {
449 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
450 "REMOVE_X", f, 0, 0);
451 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
452 "REMOVE_Y", f, 1, 1);
453 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
454 "REMOVE_Z", f, 2, 2);
455 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
456 "REMOVE_ALL", f, 0, 3);
457 }
458
459 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
460 "SIGMA", 0, 0, false, true);
461 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
462 "SIGMA", 1, 1, false, true);
463 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
464 "SIGMA", 2, 2, false, true);
465 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
466 "SIGMA", 0, 3, false, true);
467 CHKERR bc_mng->removeBlockDOFsOnEntities(
468 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
469
470 // Note remove has to be always before push. Then node marking will be
471 // corrupted.
472 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
473 simple->getProblemName(), "U");
474
476}
477//! [Boundary condition]
478
479//! [Push operators to pip]
483 auto *pip_mng = mField.getInterface<PipelineManager>();
484 auto time_scale = boost::make_shared<ScaledTimeScale>();
485 auto body_force_time_scale =
486 boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
487
488 auto integration_rule_vol = [](int, int, int approx_order) {
489 return 2 * approx_order + geom_order - 1;
490 };
491 auto integration_rule_boundary = [](int, int, int approx_order) {
492 return 2 * approx_order + geom_order - 1;
493 };
494
495 auto add_domain_base_ops = [&](auto &pip) {
498 "GEOMETRY");
500 };
501
502 auto add_domain_ops_lhs = [&](auto &pip) {
504
505 //! [Only used for dynamics]
508 //! [Only used for dynamics]
509 if (is_quasi_static == PETSC_FALSE) {
510
511 auto *pip_mng = mField.getInterface<PipelineManager>();
512 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
513
514 auto get_inertia_and_mass_damping =
515 [this, fe_domain_lhs](const double, const double, const double) {
516 return (rho * scale) * fe_domain_lhs->ts_aa +
517 (alpha_damping * scale) * fe_domain_lhs->ts_a;
518 };
519 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
520 } else {
521
522 auto *pip_mng = mField.getInterface<PipelineManager>();
523 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
524
525 auto get_inertia_and_mass_damping =
526 [this, fe_domain_lhs](const double, const double, const double) {
527 return (alpha_damping * scale) * fe_domain_lhs->ts_a;
528 };
529 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
530 }
531
532 if (!mfrontInterfacePtr) {
533 if (!is_large_strain) {
534 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
535 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
536 } else {
537 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
538 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
539 }
540 } else {
541 CHKERR mfrontInterfacePtr->opFactoryDomainLhs(pip, "U");
542 }
543
545 };
546
547 auto add_domain_ops_rhs = [&](auto &pip) {
549
551 pip, mField, "U", {body_force_time_scale}, Sev::inform);
552
553 //! [Only used for dynamics]
556 //! [Only used for dynamics]
557
558 // only in case of dynamics
559 if (is_quasi_static == PETSC_FALSE) {
560 auto mat_acceleration = boost::make_shared<MatrixDouble>();
562 "U", mat_acceleration));
563 pip.push_back(
564 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
565 return rho * scale;
566 }));
567 }
568
569 // only in case of viscosity
570 if (alpha_damping > 0) {
571 auto mat_velocity = boost::make_shared<MatrixDouble>();
572 pip.push_back(
573 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
574 pip.push_back(
575 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
576 return alpha_damping * scale;
577 }));
578 }
579
580 if (!mfrontInterfacePtr) {
581 if (!is_large_strain) {
582 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
583 mField, pip, "U", "MAT_ELASTIC", Sev::inform, 1, scale);
584 } else {
585 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
586 mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
587 }
588 } else {
589 CHKERR mfrontInterfacePtr->opFactoryDomainRhs(pip, "U");
590 }
591
592 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
593 pip, "SIGMA", "U", is_axisymmetric);
594
596 };
597
598 auto add_boundary_base_ops = [&](auto &pip) {
601 "GEOMETRY");
602 // We have to integrate on curved face geometry, thus integration weight
603 // have to adjusted.
604 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
606 };
607
608 auto add_boundary_ops_lhs = [&](auto &pip) {
610
611 //! [Operators used for contact]
614 //! [Operators used for contact]
615
616 // Add Natural BCs to LHS
618 pip, mField, "U", Sev::inform);
619
620 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
621
622 auto *pip_mng = mField.getInterface<PipelineManager>();
623 auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
624
625 pip.push_back(new OpSpringLhs(
626 "U", "U",
627
628 [this, fe_boundary_lhs](double, double, double) {
629 return spring_stiffness * scale +
630 (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
631 }
632
633 ));
634 }
635
636 CHKERR
637 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
638 pip, "SIGMA", "U", is_axisymmetric);
640 DomainEle>(
641 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
642 integration_rule_vol, is_axisymmetric);
643
645 };
646
647 auto add_boundary_ops_rhs = [&](auto &pip) {
649
650 //! [Operators used for contact]
653 //! [Operators used for contact]
654
655 // Add Natural BCs to RHS
657 pip, mField, "U", {time_scale}, Sev::inform);
658
659 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
660 auto u_disp = boost::make_shared<MatrixDouble>();
661 auto dot_u_disp = boost::make_shared<MatrixDouble>();
662 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
663 pip.push_back(
664 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
665 pip.push_back(
666 new OpSpringRhs("U", u_disp, [this](double, double, double) {
667 return spring_stiffness * scale;
668 }));
669 pip.push_back(
670 new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
672 }));
673 }
674
675 CHKERR
676 ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
677 pip, "SIGMA", "U", is_axisymmetric);
678
680 };
681
682 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
683 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
684 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
685 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
686
687 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
688 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
689 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
690 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
691
692 if (mfrontInterfacePtr) {
693 CHKERR mfrontInterfacePtr->setUpdateInternalVariablesOperators(
694 integration_rule_vol, "U");
695 CHKERR mfrontInterfacePtr->setPostProcessOperators(
696 integration_rule_vol, simple->getDomainFEName(), "U", order);
697 }
698
699 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
700 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
701 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
702 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
703
705}
706//! [Push operators to pip]
707
708//! [Solve]
709struct SetUpSchur {
710 static boost::shared_ptr<SetUpSchur>
713
714protected:
715 SetUpSchur() = default;
716};
717
720
723 ISManager *is_manager = mField.getInterface<ISManager>();
724
725 auto set_section_monitor = [&](auto solver) {
727 SNES snes;
728 CHKERR TSGetSNES(solver, &snes);
729 PetscViewerAndFormat *vf;
730 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
731 PETSC_VIEWER_DEFAULT, &vf);
732 CHKERR SNESMonitorSet(
733 snes,
734 (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal,
735 void *))SNESMonitorFields,
736 vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
738 };
739
740 auto scatter_create = [&](auto D, auto coeff) {
742 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
743 ROW, "U", coeff, coeff, is);
744 int loc_size;
745 CHKERR ISGetLocalSize(is, &loc_size);
746 Vec v;
747 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
748 VecScatter scatter;
749 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
750 return std::make_tuple(SmartPetscObj<Vec>(v),
752 };
753
754 auto set_time_monitor = [&](auto dm, auto solver) {
756 monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
757 boost::shared_ptr<ForcesAndSourcesCore> null;
758 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
759 monitorPtr, null, null);
761 };
762
763 auto set_essential_bc = [&]() {
765 // This is low level pushing finite elements (pipelines) to solver
766 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
767 auto pre_proc_ptr = boost::make_shared<FEMethod>();
768 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
769 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
770
771 // Add boundary condition scaling
772 auto time_scale = boost::make_shared<TimeScale>();
773
774 auto get_bc_hook_rhs = [&]() {
776 {time_scale}, false);
777 return hook;
778 };
779 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
780
781 auto get_post_proc_hook_rhs = [&]() {
783 mField, post_proc_rhs_ptr, 1.);
784 };
785 auto get_post_proc_hook_lhs = [&]() {
787 mField, post_proc_lhs_ptr, 1.);
788 };
789 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
790
791 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
792 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
793 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
794 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
795 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
797 };
798
799 // Set up Schur preconditioner
800 auto set_schur_pc = [&](auto solver) {
801 boost::shared_ptr<SetUpSchur> schur_ptr;
802 if (AT == AssemblyType::BLOCK_SCHUR) {
803 // Set up Schur preconditioner
805 CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
806 }
807 return schur_ptr;
808 };
809
810 auto dm = simple->getDM();
811 auto D = createDMVector(dm);
812
814
815 uXScatter = scatter_create(D, 0);
816 uYScatter = scatter_create(D, 1);
817 if (SPACE_DIM == 3)
818 uZScatter = scatter_create(D, 2);
819
820 // Add extra finite elements to SNES solver pipelines to resolve essential
821 // boundary conditions
822 CHKERR set_essential_bc();
823
824 if (is_quasi_static == PETSC_TRUE) {
825 auto solver = pip_mng->createTSIM();
826 CHKERR TSSetFromOptions(solver);
827
828 auto B = createDMMatrix(dm);
829 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
830 auto schur_pc_ptr = set_schur_pc(solver);
831
832 auto D = createDMVector(dm);
833 CHKERR set_section_monitor(solver);
834 CHKERR set_time_monitor(dm, solver);
835 CHKERR TSSetSolution(solver, D);
836 CHKERR TSSetUp(solver);
837 CHKERR TSSolve(solver, NULL);
838 } else {
839 auto solver = pip_mng->createTSIM2();
840 CHKERR TSSetFromOptions(solver);
841
842 auto B = createDMMatrix(dm);
843 CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
844 auto schur_pc_ptr = set_schur_pc(solver);
845
846 auto D = createDMVector(dm);
847 auto DD = vectorDuplicate(D);
848 CHKERR set_section_monitor(solver);
849 CHKERR set_time_monitor(dm, solver);
850 CHKERR TS2SetSolution(solver, D, DD);
851 CHKERR TSSetUp(solver);
852 CHKERR TSSolve(solver, NULL);
853 }
854
856}
857//! [Solve]
858
859//! [Check]
862 if (atom_test && !mField.get_comm_rank()) {
863 const double *t_ptr;
864 CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
865 double hertz_force;
866 double fem_force;
867 double analytical_active_area = 1.0;
868 double norm = 1e-5;
869 double tol_force = 1e-3;
870 double tol_norm = 7.5; // change when analytical functions are updated
871 double tol_area = 3e-2;
872 double fem_active_area = t_ptr[3];
873
874 switch (atom_test) {
875 case 1: // plane stress
876 hertz_force = 3.927;
877 fem_force = t_ptr[1];
878 break;
879
880 case 2: // plane strain
881 hertz_force = 4.675;
882 fem_force = t_ptr[1];
883 norm = monitorPtr->getErrorNorm(1);
884 break;
885
886 case 3: // Hertz 3D
887 hertz_force = 3.968;
888 tol_force = 2e-3;
889 fem_force = t_ptr[2];
890 analytical_active_area = M_PI / 4;
891 tol_area = 0.2;
892 break;
893
894 case 4: // axisymmetric
895 tol_force = 5e-3;
896 tol_area = 0.2;
897 // analytical_active_area = M_PI;
898
899 case 5: // axisymmetric
900 hertz_force = 15.873;
901 tol_force = 5e-3;
902 fem_force = t_ptr[1];
903 norm = monitorPtr->getErrorNorm(1);
904 analytical_active_area = M_PI;
905 break;
906
907 case 6: // wavy 2d
908 hertz_force = 0.374;
909 fem_force = t_ptr[1];
910 break;
911
912 case 7: // wavy 3d
913 hertz_force = 0.5289;
914 fem_force = t_ptr[2];
915 break;
916
917 default:
918 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
919 "atom test %d does not exist", atom_test);
920 }
921 if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
922 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
923 "atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
924 atom_test, fem_force, hertz_force);
925 }
926 if (norm > tol_norm) {
927 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
928 "atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
929 atom_test, norm, tol_norm);
930 }
931 if (fabs(fem_active_area - analytical_active_area) > tol_area) {
932 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
933 "atom test %d failed: AREA computed %3.4e but should be %3.4e",
934 atom_test, fem_active_area, analytical_active_area);
935 }
936 CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
937 }
938
940
942}
943//! [Check]
944
945static char help[] = "...\n\n";
946
947int main(int argc, char *argv[]) {
948
949#ifdef ENABLE_PYTHON_BINDING
950 Py_Initialize();
951 np::initialize();
952#endif
953
954 // Initialisation of MoFEM/PETSc and MOAB data structures
955 const char param_file[] = "param_file.petsc";
956 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
957
958 // Add logging channel for CONTACT
959 auto core_log = logging::core::get();
960 core_log->add_sink(
962 LogManager::setLog("CONTACT");
963 MOFEM_LOG_TAG("CONTACT", "Indent");
964
965 try {
966
967 //! [Register MoFEM discrete manager in PETSc]
968 DMType dm_name = "DMMOFEM";
969 CHKERR DMRegister_MoFEM(dm_name);
970 DMType dm_name_mg = "DMMOFEM_MG";
972 //! [Register MoFEM discrete manager in PETSc
973
974 //! [Create MoAB]
975 moab::Core mb_instance; ///< mesh database
976 moab::Interface &moab = mb_instance; ///< mesh database interface
977 //! [Create MoAB]
978
979 //! [Create MoFEM]
980 MoFEM::Core core(moab); ///< finite element database
981 MoFEM::Interface &m_field = core; ///< finite element database interface
982 //! [Create MoFEM]
983
984 //! [Load mesh]
985 Simple *simple = m_field.getInterface<Simple>();
987 CHKERR simple->loadFile("");
988 //! [Load mesh]
989
990 //! [CONTACT]
991 Contact ex(m_field);
992 CHKERR ex.runProblem();
993 //! [CONTACT]
994 }
996
998
999#ifdef ENABLE_PYTHON_BINDING
1000 if (Py_FinalizeEx() < 0) {
1001 exit(120);
1002 }
1003#endif
1004
1005 return 0;
1006}
1007
1008struct SetUpSchurImpl : public SetUpSchur {
1009
1011
1012 virtual ~SetUpSchurImpl() {}
1013
1015
1016private:
1019 MoFEMErrorCode setPC(PC pc);
1021
1023
1027};
1028
1031 auto simple = mField.getInterface<Simple>();
1032 auto pip = mField.getInterface<PipelineManager>();
1033
1034 SNES snes;
1035 CHKERR TSGetSNES(solver, &snes);
1036 KSP ksp;
1037 CHKERR SNESGetKSP(snes, &ksp);
1038 CHKERR KSPSetFromOptions(ksp);
1039
1040 PC pc;
1041 CHKERR KSPGetPC(ksp, &pc);
1042
1043 PetscBool is_pcfs = PETSC_FALSE;
1044 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1045 if (is_pcfs) {
1046
1047 MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
1048
1049 if (S) {
1052 "It is expected that Schur matrix is not allocated. This is "
1053 "possible only if PC is set up twice");
1054 }
1055
1057
1058 // Add data to DM storage
1060 CHKERR MatSetBlockSize(S, SPACE_DIM);
1061 // CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
1062
1063 // Set DM to use shell block matrix
1064 DM solver_dm;
1065 CHKERR TSGetDM(solver, &solver_dm);
1066 CHKERR DMSetMatType(solver_dm, MATSHELL);
1067
1068 auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1069 auto A = createDMBlockMat(simple->getDM());
1070 auto P = createDMNestSchurMat(simple->getDM());
1071
1072 if (is_quasi_static == PETSC_TRUE) {
1073 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1074 Mat A, Mat B, void *ctx) {
1075 return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1076 };
1077 CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1078 } else {
1079 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1080 PetscReal a, PetscReal aa, Mat A, Mat B,
1081 void *ctx) {
1082 return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1083 };
1084 CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1085 }
1086 CHKERR KSPSetOperators(ksp, A, P);
1087
1089 CHKERR setPC(pc);
1090 CHKERR TSSetUp(solver);
1091 CHKERR KSPSetUp(ksp);
1093
1094 } else {
1095 MOFEM_LOG("CONTACT", Sev::inform) << "No Schur PC";
1096 pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1097 pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
1098 pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1099 pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
1100 }
1102}
1103
1106 auto simple = mField.getInterface<Simple>();
1107
1108 auto create_dm = [&](const char *name, const char *field_name, auto dm_type) {
1109 auto dm = createDM(mField.get_comm(), dm_type);
1110 auto create_dm_imp = [&]() {
1112 CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
1113 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
1114 CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
1117 CHKERR DMSetUp(dm);
1119 };
1121 create_dm_imp(),
1122 "Error in creating schurDM. It is possible that schurDM is "
1123 "already created");
1124 return dm;
1125 };
1126
1127 // Note: here we can make block with bubbles of "U" and "SIGMA" fields. See
1128 // vec-0 where bubbles are added.
1129
1130 schurDM = create_dm("SCHUR", "U", "DMMOFEM_MG");
1131 blockDM = create_dm("BLOCK", "SIGMA", "DMMOFEM");
1132
1133 if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1134
1135 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1136 auto block_mat_data = createBlockMatStructure(
1137 simple->getDM(),
1138
1139 {{
1140
1141 simple->getDomainFEName(),
1142
1143 {
1144
1145 {"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
1146
1147 }}}
1148
1149 );
1150
1152
1153 {schur_dm, block_dm}, block_mat_data,
1154
1155 {"SIGMA"}, {nullptr}, true
1156
1157 );
1158 };
1159
1160 auto nested_mat_data = get_nested_mat_data(schurDM, blockDM);
1161 CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1162
1163 } else {
1164 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1165 "Only BLOCK_SCHUR is implemented");
1166 }
1167
1169}
1170
1173
1174 double eps_stab = 1e-4;
1175 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
1176 PETSC_NULLPTR);
1177
1180 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1181
1182 auto simple = mField.getInterface<Simple>();
1183 auto pip = mField.getInterface<PipelineManager>();
1184
1185 auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
1186 auto ao_up = createAOMappingIS(dm_is, PETSC_NULLPTR);
1187
1188 // Boundary
1189 pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1190 pip->getOpBoundaryLhsPipeline().push_back(
1191 new OpMassStab("SIGMA", "SIGMA",
1192 [eps_stab](double, double, double) { return eps_stab; }));
1193 pip->getOpBoundaryLhsPipeline().push_back(
1194
1195 createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
1196
1197 );
1198
1199 // Domain
1200 pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1201 pip->getOpDomainLhsPipeline().push_back(
1202
1203 createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
1204
1205 );
1206
1207 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1208 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1209
1210 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1212 CHKERR MatZeroEntries(S);
1213 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1215 };
1216
1217 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1218 post_proc_schur_lhs_ptr]() {
1220 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1221 auto print_mat_norm = [this](auto a, std::string prefix) {
1223 double nrm;
1224 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1225 MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1227 };
1228 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1229 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1231 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1232#ifndef NDEBUG
1233 CHKERR print_mat_norm(S, "S");
1234#endif // NDEBUG
1235 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1237 };
1238
1239 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1240 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1241 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1242
1244}
1245
1248 auto block_is = getDMSubData(blockDM)->getSmartRowIs();
1249 CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
1250 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1252}
1253
1256 KSP *subksp;
1257 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
1258 auto get_pc = [](auto ksp) {
1259 PC pc_raw;
1260 CHKERR KSPGetPC(ksp, &pc_raw);
1261 return SmartPetscObj<PC>(pc_raw, true); // bump reference
1262 };
1263 CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1264
1265 auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
1267 CHKERR PCSetDM(pc, dm);
1268 PetscBool same = PETSC_FALSE;
1269 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
1270 if (same) {
1272 pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
1273 CHKERR PCSetFromOptions(pc);
1274 }
1276 };
1277
1278 auto set_pc_ksp = [&](auto dm, auto pc, auto S) {
1280 PetscBool same = PETSC_FALSE;
1281 PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
1282 if (same) {
1283 CHKERR PCSetFromOptions(pc);
1284 KSP inner_ksp;
1285 CHKERR PCKSPGetKSP(pc, &inner_ksp);
1286 CHKERR KSPSetFromOptions(inner_ksp);
1287 PC inner_pc;
1288 CHKERR KSPGetPC(inner_ksp, &inner_pc);
1289 CHKERR PCSetFromOptions(inner_pc);
1290 CHKERR set_pc_p_mg(dm, inner_pc, S);
1291 }
1293 };
1294
1295 CHKERR set_pc_ksp(schurDM, get_pc(subksp[1]), S);
1296 CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
1297
1298 CHKERR PetscFree(subksp);
1300}
1301
1302boost::shared_ptr<SetUpSchur>
1304 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1305}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Implementation of Hookes operator Hookes for linear elastic problems in MoFEM.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr double a
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
double young_modulus
Definition contact.cpp:82
int sigma_order
Definition contact.cpp:80
constexpr AssemblyType AT
Definition contact.cpp:34
constexpr IntegrationType IT
Definition contact.cpp:37
static char help[]
[Check]
Definition contact.cpp:945
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
double spring_stiffness
Definition contact.cpp:85
PetscBool is_large_strain
Definition contact.cpp:92
double rho
Definition contact.cpp:84
int atom_test
Definition contact.cpp:94
#define EXECUTABLE_DIMENSION
Definition contact.cpp:14
PetscBool is_quasi_static
[Operators used for contact]
Definition contact.cpp:76
double alpha_damping
Definition contact.cpp:87
constexpr int SPACE_DIM
Definition contact.cpp:53
double poisson_ratio
Definition contact.cpp:83
double vis_spring_stiffness
Definition contact.cpp:86
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition contact.cpp:60
double scale
Definition contact.cpp:89
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition contact.cpp:71
int contact_order
Definition contact.cpp:79
int order
Definition contact.cpp:78
PetscBool is_axisymmetric
Definition contact.cpp:91
int geom_order
Definition contact.cpp:81
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition contact.cpp:73
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition contact.cpp:67
@ ROW
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
FieldSpace
approximation spaces
Definition definitions.h:82
@ H1
continuous field
Definition definitions.h:85
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1191
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
@ GAUSS
Gaussian quadrature integration.
@ PETSC
Standard PETSc assembly.
@ BLOCK_PRECONDITIONER_SCHUR
Block preconditioner Schur assembly.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
double D
const double v
phase velocity of light in medium (cm/ns)
double cn_contact
Definition contact.cpp:97
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition DMMoFEM.cpp:1046
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:169
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1276
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition Schur.cpp:2585
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition Schur.cpp:2627
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition DMMoFEM.hpp:1292
boost::shared_ptr< MFrontInterface > createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh, AssemblyType at=AssemblyType::PETSC)
create mfront interface
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition TsCtx.cpp:519
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
@ AXISYMMETRICAL
Axisymmetrical model hypothesis.
@ PLANESTRAIN
Plane strain model hypothesis.
@ TRIDIMENSIONAL
3D model hypothesis.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition Schur.cpp:1082
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition Schur.cpp:2343
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition DMMoFEM.cpp:1554
auto createDMNestSchurMat(DM dm)
Definition DMMoFEM.hpp:1218
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition Schur.cpp:2580
auto createDMBlockMat(DM dm)
Definition DMMoFEM.hpp:1211
constexpr AssemblyType A
PetscBool is_quasi_static
Definition plastic.cpp:143
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61
constexpr FieldSpace CONTACT_SPACE
Definition plastic.cpp:52
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56
FTensor::Index< 'm', 3 > m
static SmartPetscObj< Vec > totalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
double getScale(const double time)
Get scaling at given time.
Definition contact.cpp:150
boost::shared_ptr< MFrontInterface > mfrontInterfacePtr
Definition contact.cpp:141
boost::shared_ptr< Monitor > monitorPtr
Definition contact.cpp:142
MoFEMErrorCode tsSolve()
Definition contact.cpp:718
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition contact.cpp:138
MoFEMErrorCode runProblem()
[Run problem]
Definition contact.cpp:157
MoFEMErrorCode setupProblem()
[Run problem]
Definition contact.cpp:170
MoFEMErrorCode OPs()
[Boundary condition]
Definition contact.cpp:480
MoFEMErrorCode createCommonData()
[Set up problem]
Definition contact.cpp:310
MoFEMErrorCode checkResults()
[Solve]
Definition contact.cpp:860
MoFEM::Interface & mField
Definition contact.cpp:128
MoFEMErrorCode bC()
[Create common data]
Definition contact.cpp:443
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition contact.cpp:137
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition contact.cpp:139
Contact(MoFEM::Interface &m_field)
Definition contact.cpp:122
Add operators pushing bases from local to physical configuration.
Boundary condition manager for finite element problem setup.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition Natural.hpp:65
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Template struct for dimension-specific finite element types.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
SmartPetscObj< Mat > S
MoFEMErrorCode createSubDM()
Definition contact.cpp:1104
MoFEMErrorCode setDiagonalPC(PC pc)
Definition contact.cpp:1254
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition contact.cpp:1010
SmartPetscObj< DM > schurDM
Definition contact.cpp:1025
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
SmartPetscObj< DM > blockDM
Definition contact.cpp:1026
virtual ~SetUpSchurImpl()
Definition contact.cpp:1012
MoFEMErrorCode setPC(PC pc)
Definition contact.cpp:1246
MoFEMErrorCode setOperator()
Definition contact.cpp:1171
MoFEM::Interface & mField
[Push operators to pipeline]
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
virtual MoFEMErrorCode setUp(SmartPetscObj< TS > solver)=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
constexpr AssemblyType AT