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#include <GenericElementInterface.hpp>
27
28#ifdef ENABLE_PYTHON_BINDING
29#include <boost/python.hpp>
30#include <boost/python/def.hpp>
31#include <boost/python/numpy.hpp>
32namespace bp = boost::python;
33namespace np = boost::python::numpy;
34#endif
35
36constexpr AssemblyType AT =
37 (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
38 : AssemblyType::PETSC; //< selected assembly type
40 IntegrationType::GAUSS; //< selected integration type
41
42template <int DIM> struct ElementsAndOps;
43
45 static constexpr FieldSpace CONTACT_SPACE = HCURL;
46};
47
49 static constexpr FieldSpace CONTACT_SPACE = HDIV;
50};
51
54
55constexpr int SPACE_DIM =
56 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
57
58/* The above code is defining an alias `EntData` for the type
59`EntitiesFieldData::EntData`. This is a C++ syntax for creating a new name for
60an existing type. */
63using DomainEleOp = DomainEle::UserDataOperator;
65using BoundaryEleOp = BoundaryEle::UserDataOperator;
67//! [Specialisation for assembly]
68
70
71//! [Operators used for contact]
75 IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
76//! [Operators used for contact]
77
78PetscBool is_quasi_static = PETSC_TRUE;
79
80int order = 2; //< Order of displacements in the domain
81int contact_order = 2; //< Order of displacements in contact elements
82int sigma_order = 1; //< Order of Lagrange multiplier in side elements
83int geom_order = 1;
84double young_modulus = 100;
85double poisson_ratio = 0.25;
86double rho = 0.0;
87double spring_stiffness = 0.0;
89double alpha_damping = 0;
90
91double scale = 1.;
92
93PetscBool is_axisymmetric = PETSC_FALSE; //< flag for Axisymmetric model
94PetscBool is_large_strain = PETSC_FALSE; //< flag for large strain
95PetscBool is_plane_strain = PETSC_FALSE; //< flag for plane strain
96int atom_test = 0;
97
98namespace ContactOps {
99double cn_contact = 0.1;
100}; // namespace ContactOps
101
102#include <HenckyOps.hpp>
103#include <HookeOps.hpp>
104using namespace HenckyOps;
105using namespace HookeOps;
106#include <ContactOps.hpp>
107#include <PostProcContact.hpp>
108#include <ContactNaturalBC.hpp>
109
110#ifdef WITH_MODULE_MFRONT_INTERFACE
111#include <MFrontMoFEMInterface.hpp>
112#endif
113
123
124using namespace ContactOps;
125
126struct Contact {
127
128 Contact(MoFEM::Interface &m_field) : mField(m_field) {}
129
131
132private:
134
141
142 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
143 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
144 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
145
146 boost::shared_ptr<GenericElementInterface> mfrontInterface;
147 boost::shared_ptr<Monitor> monitorPtr;
148
149#ifdef ENABLE_PYTHON_BINDING
150 boost::shared_ptr<SDFPython> sdfPythonPtr;
151#endif
152
155 double getScale(const double time) {
156 return scale * MoFEM::TimeScale::getScale(time);
157 };
158 };
159};
160
161//! [Run problem]
172//! [Run problem]
173
174//! [Set up problem]
178
179 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strain", &is_large_strain,
180 PETSC_NULLPTR);
181
182 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain", &is_plane_strain,
183 PETSC_NULLPTR);
184 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
185 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-contact_order", &contact_order,
186 PETSC_NULLPTR);
187 sigma_order = std::max(order, contact_order) - 1;
188 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-sigma_order", &sigma_order,
189 PETSC_NULLPTR);
190 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
191 PETSC_NULLPTR);
192 if (!is_large_strain) {
193 MOFEM_LOG("CONTACT", Sev::inform)
194 << "Selected material model: Hooke (small strain)";
195 } else {
196 MOFEM_LOG("CONTACT", Sev::inform)
197 << "Selected deafult material model: Hencky (finite strain)";
198 }
199 MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
200 MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
201 MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
202 MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
203
204 // Select base
205 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
206 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
207 PetscInt choice_base_value = AINSWORTH;
208 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
209 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
210
212 switch (choice_base_value) {
213 case AINSWORTH:
215 MOFEM_LOG("CONTACT", Sev::inform)
216 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
217 break;
218 case DEMKOWICZ:
220 MOFEM_LOG("CONTACT", Sev::inform)
221 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
222 break;
223 default:
224 base = LASTBASE;
225 break;
226 }
227
228 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
229 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
230 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
231 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
232 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
233 SPACE_DIM);
234 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
235 SPACE_DIM);
236 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
237
238 CHKERR simple->setFieldOrder("U", order);
239 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
240
241 auto get_skin = [&]() {
242 Range body_ents;
243 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
244 Skinner skin(&mField.get_moab());
245 Range skin_ents;
246 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
247 return skin_ents;
248 };
249
250 auto filter_blocks = [&](auto skin) {
251 bool is_contact_block = false;
252 Range contact_range;
253 for (auto m :
255
256 (boost::format("%s(.*)") % "CONTACT").str()
257
258 ))
259
260 ) {
261 is_contact_block =
262 true; ///< blocs interation is collective, so that is set irrespective
263 ///< if there are entities in given rank or not in the block
264 MOFEM_LOG("CONTACT", Sev::inform)
265 << "Find contact block set: " << m->getName();
266 auto meshset = m->getMeshset();
267 Range contact_meshset_range;
268 CHKERR mField.get_moab().get_entities_by_dimension(
269 meshset, SPACE_DIM - 1, contact_meshset_range, true);
270
271 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
272 contact_meshset_range);
273 contact_range.merge(contact_meshset_range);
274 }
275 if (is_contact_block) {
276 MOFEM_LOG("SYNC", Sev::inform)
277 << "Nb entities in contact surface: " << contact_range.size();
279 skin = intersect(skin, contact_range);
280 }
281 return skin;
282 };
283
284 auto filter_true_skin = [&](auto skin) {
285 Range boundary_ents;
286 ParallelComm *pcomm =
287 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
288 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
289 PSTATUS_NOT, -1, &boundary_ents);
290 return boundary_ents;
291 };
292
293 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
294 CHKERR simple->setFieldOrder("SIGMA", 0);
295 CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
296
297 if (contact_order > order) {
298 Range ho_ents;
299
300 CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
301 moab::Interface::UNION);
302
303 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
304 CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
305 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
306 }
307
308 CHKERR simple->setUp();
309
310 auto project_ho_geometry = [&]() {
311 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
312 return mField.loop_dofs("GEOMETRY", ent_method);
313 };
314
315 PetscBool project_geometry = PETSC_TRUE;
316 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
317 &project_geometry, PETSC_NULLPTR);
318 if (project_geometry) {
319 CHKERR project_ho_geometry();
320 }
321
323} //! [Set up problem]
324
325//! [Create common data]
328
329 PetscBool use_mfront = PETSC_FALSE;
330 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_mfront", &use_mfront,
331 PETSC_NULLPTR);
332 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_axisymmetric",
333 &is_axisymmetric, PETSC_NULLPTR);
334 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
335 PETSC_NULLPTR);
336
337 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale, PETSC_NULLPTR);
338 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus", &young_modulus,
339 PETSC_NULLPTR);
340 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio", &poisson_ratio,
341 PETSC_NULLPTR);
342 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
343 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn", &cn_contact, PETSC_NULLPTR);
344 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-spring_stiffness",
345 &spring_stiffness, PETSC_NULLPTR);
346 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-vis_spring_stiffness",
347 &vis_spring_stiffness, PETSC_NULLPTR);
348 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping", &alpha_damping,
349 PETSC_NULLPTR);
350
351 if (!mfrontInterface) {
352 MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
353 MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
354 } else {
355 MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
356 }
357
358 MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
359 MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
360 MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
361 MOFEM_LOG("CONTACT", Sev::inform)
362 << "Vis spring_stiffness " << vis_spring_stiffness;
363
364 MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
365
366 PetscBool use_scale = PETSC_FALSE;
367 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_scale", &use_scale,
368 PETSC_NULLPTR);
369 if (use_scale) {
371 }
372
373 MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
374
375 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_quasi_static",
376 &is_quasi_static, PETSC_NULLPTR);
377 MOFEM_LOG("CONTACT", Sev::inform)
378 << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
379
380#ifdef ENABLE_PYTHON_BINDING
381 auto file_exists = [](std::string myfile) {
382 std::ifstream file(myfile.c_str());
383 if (file) {
384 return true;
385 }
386 return false;
387 };
388 char sdf_file_name[255] = "sdf.py";
389 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
390 sdf_file_name, 255, PETSC_NULLPTR);
391
392 if (file_exists(sdf_file_name)) {
393 MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
394 sdfPythonPtr = boost::make_shared<SDFPython>();
395 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
396 sdfPythonWeakPtr = sdfPythonPtr;
397 } else {
398 MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
399 }
400#endif
401
402 if (is_axisymmetric) {
403 if (SPACE_DIM == 3) {
404 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
405 "Use executable contact_2d with axisymmetric model");
406 } else {
407 if (!use_mfront) {
408 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
409 "Axisymmetric model is only available with MFront (set "
410 "use_mfront to 1)");
411 } else {
412 MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
413 }
414 }
415 } else {
416 if (SPACE_DIM == 2) {
417 MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
418 }
419 }
420
421 if (use_mfront) {
422#ifndef WITH_MODULE_MFRONT_INTERFACE
423 SETERRQ(
424 PETSC_COMM_SELF, MOFEM_NOT_FOUND,
425 "MFrontInterface module was not found while use_mfront was set to 1");
426#else
427 if (SPACE_DIM == 3) {
429 boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
430 mField, "U", "GEOMETRY", true, is_quasi_static);
431 } else if (SPACE_DIM == 2) {
432 if (is_axisymmetric) {
434 boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
435 mField, "U", "GEOMETRY", true, is_quasi_static);
436 } else {
437 mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
438 mField, "U", "GEOMETRY", true, is_quasi_static);
439 }
440 }
441#endif
442 CHKERR mfrontInterface->getCommandLineParameters();
443 }
444
446 auto dm = simple->getDM();
447 monitorPtr =
448 boost::make_shared<Monitor>(dm, scale, mfrontInterface, is_axisymmetric, is_large_strain);
449
450 if (use_mfront) {
451 mfrontInterface->setMonitorPtr(monitorPtr);
452 }
453
455}
456//! [Create common data]
457
458//! [Boundary condition]
461 auto bc_mng = mField.getInterface<BcManager>();
463
464 for (auto f : {"U", "SIGMA"}) {
465 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
466 "REMOVE_X", f, 0, 0);
467 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
468 "REMOVE_Y", f, 1, 1);
469 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
470 "REMOVE_Z", f, 2, 2);
471 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
472 "REMOVE_ALL", f, 0, 3);
473 }
474
475 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
476 "SIGMA", 0, 0, false, true);
477 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
478 "SIGMA", 1, 1, false, true);
479 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
480 "SIGMA", 2, 2, false, true);
481 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
482 "SIGMA", 0, 3, false, true);
483 CHKERR bc_mng->removeBlockDOFsOnEntities(
484 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
485
486 // Note remove has to be always before push. Then node marking will be
487 // corrupted.
488 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
489 simple->getProblemName(), "U");
490
492}
493//! [Boundary condition]
494
495//! [Push operators to pip]
499 auto *pip_mng = mField.getInterface<PipelineManager>();
500 auto time_scale = boost::make_shared<ScaledTimeScale>();
501 auto body_force_time_scale =
502 boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
503
504 auto integration_rule_vol = [](int, int, int approx_order) {
505 return 2 * approx_order + geom_order - 1;
506 };
507 auto integration_rule_boundary = [](int, int, int approx_order) {
508 return 2 * approx_order + geom_order - 1;
509 };
510
511 auto add_domain_base_ops = [&](auto &pip) {
514 "GEOMETRY");
516 };
517
518 auto add_domain_ops_lhs = [&](auto &pip) {
520
521 //! [Only used for dynamics]
524 //! [Only used for dynamics]
525 if (is_quasi_static == PETSC_FALSE) {
526
527 auto *pip_mng = mField.getInterface<PipelineManager>();
528 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
529
530 auto get_inertia_and_mass_damping =
531 [this, fe_domain_lhs](const double, const double, const double) {
532 return (rho * scale) * fe_domain_lhs->ts_aa +
533 (alpha_damping * scale) * fe_domain_lhs->ts_a;
534 };
535 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
536 } else {
537
538 auto *pip_mng = mField.getInterface<PipelineManager>();
539 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
540
541 auto get_inertia_and_mass_damping =
542 [this, fe_domain_lhs](const double, const double, const double) {
543 return (alpha_damping * scale) * fe_domain_lhs->ts_a;
544 };
545 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
546 }
547
548 if (!mfrontInterface) {
549 if (!is_large_strain) {
551 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
552 } else {
554 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
555 }
556 } else {
557 CHKERR mfrontInterface->opFactoryDomainLhs(pip);
558 }
559
561 };
562
563 auto add_domain_ops_rhs = [&](auto &pip) {
565
567 pip, mField, "U", {body_force_time_scale}, Sev::inform);
568
569 //! [Only used for dynamics]
572 //! [Only used for dynamics]
573
574 // only in case of dynamics
575 if (is_quasi_static == PETSC_FALSE) {
576 auto mat_acceleration = boost::make_shared<MatrixDouble>();
578 "U", mat_acceleration));
579 pip.push_back(
580 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
581 return rho * scale;
582 }));
583 }
584
585 // only in case of viscosity
586 if (alpha_damping > 0) {
587 auto mat_velocity = boost::make_shared<MatrixDouble>();
588 pip.push_back(
589 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
590 pip.push_back(
591 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
592 return alpha_damping * scale;
593 }));
594 }
595
596 if (!mfrontInterface) {
597 if (!is_large_strain) {
599 mField, pip, "U", "MAT_ELASTIC", Sev::inform, 1, scale);
600 } else {
602 mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
603 }
604 } else {
605 CHKERR mfrontInterface->opFactoryDomainRhs(pip);
606 }
607
609 pip, "SIGMA", "U", is_axisymmetric);
610
612 };
613
614 auto add_boundary_base_ops = [&](auto &pip) {
617 "GEOMETRY");
618 // We have to integrate on curved face geometry, thus integration weight
619 // have to adjusted.
620 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
622 };
623
624 auto add_boundary_ops_lhs = [&](auto &pip) {
626
627 //! [Operators used for contact]
630 //! [Operators used for contact]
631
632 // Add Natural BCs to LHS
634 pip, mField, "U", Sev::inform);
635
636 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
637
638 auto *pip_mng = mField.getInterface<PipelineManager>();
639 auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
640
641 pip.push_back(new OpSpringLhs(
642 "U", "U",
643
644 [this, fe_boundary_lhs](double, double, double) {
645 return spring_stiffness * scale +
646 (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
647 }
648
649 ));
650 }
651
652 CHKERR
654 pip, "SIGMA", "U", is_axisymmetric);
656 DomainEle>(
657 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
658 integration_rule_vol, is_axisymmetric);
659
661 };
662
663 auto add_boundary_ops_rhs = [&](auto &pip) {
665
666 //! [Operators used for contact]
669 //! [Operators used for contact]
670
671 // Add Natural BCs to RHS
673 pip, mField, "U", {time_scale}, Sev::inform);
674
675 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
676 auto u_disp = boost::make_shared<MatrixDouble>();
677 auto dot_u_disp = boost::make_shared<MatrixDouble>();
678 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
679 pip.push_back(
680 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
681 pip.push_back(
682 new OpSpringRhs("U", u_disp, [this](double, double, double) {
683 return spring_stiffness * scale;
684 }));
685 pip.push_back(
686 new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
688 }));
689 }
690
691 CHKERR
693 pip, "SIGMA", "U", is_axisymmetric);
694
696 };
697
698 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
699 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
700 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
701 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
702
703 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
704 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
705 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
706 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
707
708 if (mfrontInterface) {
709 CHKERR mfrontInterface->setUpdateElementVariablesOperators();
710 }
711
712 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
713 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
714 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
715 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
716
718}
719//! [Push operators to pip]
720
721//! [Solve]
722struct SetUpSchur {
723 static boost::shared_ptr<SetUpSchur>
726
727protected:
728 SetUpSchur() = default;
729};
730
733
736 ISManager *is_manager = mField.getInterface<ISManager>();
737
738 auto set_section_monitor = [&](auto solver) {
740 SNES snes;
741 CHKERR TSGetSNES(solver, &snes);
742 PetscViewerAndFormat *vf;
743 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
744 PETSC_VIEWER_DEFAULT, &vf);
745 CHKERR SNESMonitorSet(
746 snes,
747 (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
748 vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
750 };
751
752 auto scatter_create = [&](auto D, auto coeff) {
754 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
755 ROW, "U", coeff, coeff, is);
756 int loc_size;
757 CHKERR ISGetLocalSize(is, &loc_size);
758 Vec v;
759 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
760 VecScatter scatter;
761 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
762 return std::make_tuple(SmartPetscObj<Vec>(v),
764 };
765
766 auto set_time_monitor = [&](auto dm, auto solver) {
768 monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
769 boost::shared_ptr<ForcesAndSourcesCore> null;
770 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
771 monitorPtr, null, null);
773 };
774
775 auto set_essential_bc = [&]() {
777 // This is low level pushing finite elements (pipelines) to solver
778 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
779 auto pre_proc_ptr = boost::make_shared<FEMethod>();
780 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
781 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
782
783 // Add boundary condition scaling
784 auto time_scale = boost::make_shared<TimeScale>();
785
786 auto get_bc_hook_rhs = [&]() {
788 {time_scale}, false);
789 return hook;
790 };
791 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
792
793 auto get_post_proc_hook_rhs = [&]() {
795 mField, post_proc_rhs_ptr, 1.);
796 };
797 auto get_post_proc_hook_lhs = [&]() {
799 mField, post_proc_lhs_ptr, 1.);
800 };
801 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
802
803 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
804 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
805 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
806 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
807 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
809 };
810
811 // Set up Schur preconditioner
812 auto set_schur_pc = [&](auto solver) {
813 boost::shared_ptr<SetUpSchur> schur_ptr;
814 if (AT == AssemblyType::BLOCK_SCHUR) {
815 // Set up Schur preconditioner
817 CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
818 }
819 return schur_ptr;
820 };
821
822 auto dm = simple->getDM();
823 auto D = createDMVector(dm);
824
826
827 uXScatter = scatter_create(D, 0);
828 uYScatter = scatter_create(D, 1);
829 if (SPACE_DIM == 3)
830 uZScatter = scatter_create(D, 2);
831
832 // Add extra finite elements to SNES solver pipelines to resolve essential
833 // boundary conditions
834 CHKERR set_essential_bc();
835
836 if (is_quasi_static == PETSC_TRUE) {
837 auto solver = pip_mng->createTSIM();
838 CHKERR TSSetFromOptions(solver);
839
840 auto B = createDMMatrix(dm);
841 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
842 auto schur_pc_ptr = set_schur_pc(solver);
843
844 auto D = createDMVector(dm);
845 CHKERR set_section_monitor(solver);
846 CHKERR set_time_monitor(dm, solver);
847 CHKERR TSSetSolution(solver, D);
848 CHKERR TSSetUp(solver);
849 CHKERR TSSolve(solver, NULL);
850 } else {
851 auto solver = pip_mng->createTSIM2();
852 CHKERR TSSetFromOptions(solver);
853
854 auto B = createDMMatrix(dm);
855 CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
856 auto schur_pc_ptr = set_schur_pc(solver);
857
858 auto D = createDMVector(dm);
859 auto DD = vectorDuplicate(D);
860 CHKERR set_section_monitor(solver);
861 CHKERR set_time_monitor(dm, solver);
862 CHKERR TS2SetSolution(solver, D, DD);
863 CHKERR TSSetUp(solver);
864 CHKERR TSSolve(solver, NULL);
865 }
866
868}
869//! [Solve]
870
871//! [Check]
874 if (atom_test && !mField.get_comm_rank()) {
875 const double *t_ptr;
876 CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
877 double hertz_force;
878 double fem_force;
879 double analytical_active_area = 1.0;
880 double norm = 1e-5;
881 double tol_force = 1e-3;
882 double tol_norm = 7.5; // change when analytical functions are updated
883 double tol_area = 3e-2;
884 double fem_active_area = t_ptr[3];
885
886 switch (atom_test) {
887 case 1: // plane stress
888 hertz_force = 3.927;
889 fem_force = t_ptr[1];
890 break;
891
892 case 2: // plane strain
893 hertz_force = 4.675;
894 fem_force = t_ptr[1];
895 norm = monitorPtr->getErrorNorm(1);
896 break;
897
898 case 3: // Hertz 3D
899 hertz_force = 3.968;
900 tol_force = 2e-3;
901 fem_force = t_ptr[2];
902 analytical_active_area = M_PI / 4;
903 tol_area = 0.2;
904 break;
905
906 case 4: // axisymmetric
907 tol_force = 5e-3;
908 tol_area = 0.2;
909 // analytical_active_area = M_PI;
910
911 case 5: // axisymmetric
912 hertz_force = 15.873;
913 tol_force = 5e-3;
914 fem_force = t_ptr[1];
915 norm = monitorPtr->getErrorNorm(1);
916 analytical_active_area = M_PI;
917 break;
918
919 case 6: // wavy 2d
920 hertz_force = 0.374;
921 fem_force = t_ptr[1];
922 break;
923
924 case 7: // wavy 3d
925 hertz_force = 0.5289;
926 fem_force = t_ptr[2];
927 break;
928
929 default:
930 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
931 "atom test %d does not exist", atom_test);
932 }
933 if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
934 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
935 "atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
936 atom_test, fem_force, hertz_force);
937 }
938 if (norm > tol_norm) {
939 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
940 "atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
941 atom_test, norm, tol_norm);
942 }
943 if (fabs(fem_active_area - analytical_active_area) > tol_area) {
944 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
945 "atom test %d failed: AREA computed %3.4e but should be %3.4e",
946 atom_test, fem_active_area, analytical_active_area);
947 }
948 CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
949 }
950
952
954}
955//! [Check]
956
957static char help[] = "...\n\n";
958
959int main(int argc, char *argv[]) {
960
961#ifdef ENABLE_PYTHON_BINDING
962 Py_Initialize();
963 np::initialize();
964#endif
965
966 // Initialisation of MoFEM/PETSc and MOAB data structures
967 const char param_file[] = "param_file.petsc";
968 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
969
970 // Add logging channel for CONTACT
971 auto core_log = logging::core::get();
972 core_log->add_sink(
974 LogManager::setLog("CONTACT");
975 MOFEM_LOG_TAG("CONTACT", "Indent");
976
977 try {
978
979 //! [Register MoFEM discrete manager in PETSc]
980 DMType dm_name = "DMMOFEM";
981 CHKERR DMRegister_MoFEM(dm_name);
982 DMType dm_name_mg = "DMMOFEM_MG";
984 //! [Register MoFEM discrete manager in PETSc
985
986 //! [Create MoAB]
987 moab::Core mb_instance; ///< mesh database
988 moab::Interface &moab = mb_instance; ///< mesh database interface
989 //! [Create MoAB]
990
991 //! [Create MoFEM]
992 MoFEM::Core core(moab); ///< finite element database
993 MoFEM::Interface &m_field = core; ///< finite element database interface
994 //! [Create MoFEM]
995
996 //! [Load mesh]
997 Simple *simple = m_field.getInterface<Simple>();
999 CHKERR simple->loadFile("");
1000 //! [Load mesh]
1001
1002 //! [CONTACT]
1003 Contact ex(m_field);
1004 CHKERR ex.runProblem();
1005 //! [CONTACT]
1006 }
1008
1010
1011#ifdef ENABLE_PYTHON_BINDING
1012 if (Py_FinalizeEx() < 0) {
1013 exit(120);
1014 }
1015#endif
1016
1017 return 0;
1018}
1019
1020struct SetUpSchurImpl : public SetUpSchur {
1021
1023
1024 virtual ~SetUpSchurImpl() {}
1025
1027
1028private:
1031 MoFEMErrorCode setPC(PC pc);
1033
1035
1039};
1040
1043 auto simple = mField.getInterface<Simple>();
1044 auto pip = mField.getInterface<PipelineManager>();
1045
1046 SNES snes;
1047 CHKERR TSGetSNES(solver, &snes);
1048 KSP ksp;
1049 CHKERR SNESGetKSP(snes, &ksp);
1050 CHKERR KSPSetFromOptions(ksp);
1051
1052 PC pc;
1053 CHKERR KSPGetPC(ksp, &pc);
1054
1055 PetscBool is_pcfs = PETSC_FALSE;
1056 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1057 if (is_pcfs) {
1058
1059 MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
1060
1061 if (S) {
1064 "It is expected that Schur matrix is not allocated. This is "
1065 "possible only if PC is set up twice");
1066 }
1067
1069
1070 // Add data to DM storage
1072 CHKERR MatSetBlockSize(S, SPACE_DIM);
1073 // CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
1074
1075 // Set DM to use shell block matrix
1076 DM solver_dm;
1077 CHKERR TSGetDM(solver, &solver_dm);
1078 CHKERR DMSetMatType(solver_dm, MATSHELL);
1079
1080 auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1081 auto A = createDMBlockMat(simple->getDM());
1082 auto P = createDMNestSchurMat(simple->getDM());
1083
1084 if (is_quasi_static == PETSC_TRUE) {
1085 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1086 Mat A, Mat B, void *ctx) {
1087 return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1088 };
1089 CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1090 } else {
1091 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1092 PetscReal a, PetscReal aa, Mat A, Mat B,
1093 void *ctx) {
1094 return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1095 };
1096 CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1097 }
1098 CHKERR KSPSetOperators(ksp, A, P);
1099
1101 CHKERR setPC(pc);
1102 CHKERR TSSetUp(solver);
1103 CHKERR KSPSetUp(ksp);
1105
1106 } else {
1107 MOFEM_LOG("CONTACT", Sev::inform) << "No Schur PC";
1108 pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1109 pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
1110 pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1111 pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
1112 }
1114}
1115
1118 auto simple = mField.getInterface<Simple>();
1119
1120 auto create_dm = [&](const char *name, const char *field_name, auto dm_type) {
1121 auto dm = createDM(mField.get_comm(), dm_type);
1122 auto create_dm_imp = [&]() {
1124 CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
1125 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
1126 CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
1129 CHKERR DMSetUp(dm);
1131 };
1133 create_dm_imp(),
1134 "Error in creating schurDM. It is possible that schurDM is "
1135 "already created");
1136 return dm;
1137 };
1138
1139 // Note: here we can make block with bubbles of "U" and "SIGMA" fields. See
1140 // vec-0 where bubbles are added.
1141
1142 schurDM = create_dm("SCHUR", "U", "DMMOFEM_MG");
1143 blockDM = create_dm("BLOCK", "SIGMA", "DMMOFEM");
1144
1145 if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1146
1147 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1148 auto block_mat_data = createBlockMatStructure(
1149 simple->getDM(),
1150
1151 {{
1152
1153 simple->getDomainFEName(),
1154
1155 {
1156
1157 {"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
1158
1159 }}}
1160
1161 );
1162
1164
1165 {schur_dm, block_dm}, block_mat_data,
1166
1167 {"SIGMA"}, {nullptr}, true
1168
1169 );
1170 };
1171
1172 auto nested_mat_data = get_nested_mat_data(schurDM, blockDM);
1173 CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1174
1175 } else {
1176 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1177 "Only BLOCK_SCHUR is implemented");
1178 }
1179
1181}
1182
1185
1186 double eps_stab = 1e-4;
1187 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
1188 PETSC_NULLPTR);
1189
1192 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1193
1194 auto simple = mField.getInterface<Simple>();
1195 auto pip = mField.getInterface<PipelineManager>();
1196
1197 auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
1198 auto ao_up = createAOMappingIS(dm_is, PETSC_NULLPTR);
1199
1200 // Boundary
1201 pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
1202 pip->getOpBoundaryLhsPipeline().push_back(
1203 new OpMassStab("SIGMA", "SIGMA",
1204 [eps_stab](double, double, double) { return eps_stab; }));
1205 pip->getOpBoundaryLhsPipeline().push_back(
1206
1207 createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
1208
1209 );
1210
1211 // Domain
1212 pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1213 pip->getOpDomainLhsPipeline().push_back(
1214
1215 createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
1216
1217 );
1218
1219 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1220 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1221
1222 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1224 CHKERR MatZeroEntries(S);
1225 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1227 };
1228
1229 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1230 post_proc_schur_lhs_ptr]() {
1232 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1233 auto print_mat_norm = [this](auto a, std::string prefix) {
1235 double nrm;
1236 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1237 MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1239 };
1240 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1241 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1243 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1244#ifndef NDEBUG
1245 CHKERR print_mat_norm(S, "S");
1246#endif // NDEBUG
1247 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1249 };
1250
1251 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1252 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1253 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1254
1256}
1257
1260 auto block_is = getDMSubData(blockDM)->getSmartRowIs();
1261 CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
1262 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1264}
1265
1268 KSP *subksp;
1269 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
1270 auto get_pc = [](auto ksp) {
1271 PC pc_raw;
1272 CHKERR KSPGetPC(ksp, &pc_raw);
1273 return SmartPetscObj<PC>(pc_raw, true); // bump reference
1274 };
1275 CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1276
1277 auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
1279 CHKERR PCSetDM(pc, dm);
1280 PetscBool same = PETSC_FALSE;
1281 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
1282 if (same) {
1284 pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
1285 CHKERR PCSetFromOptions(pc);
1286 }
1288 };
1289
1290 auto set_pc_ksp = [&](auto dm, auto pc, auto S) {
1292 PetscBool same = PETSC_FALSE;
1293 PetscObjectTypeCompare((PetscObject)pc, PCKSP, &same);
1294 if (same) {
1295 CHKERR PCSetFromOptions(pc);
1296 KSP inner_ksp;
1297 CHKERR PCKSPGetKSP(pc, &inner_ksp);
1298 CHKERR KSPSetFromOptions(inner_ksp);
1299 PC inner_pc;
1300 CHKERR KSPGetPC(inner_ksp, &inner_pc);
1301 CHKERR PCSetFromOptions(inner_pc);
1302 CHKERR set_pc_p_mg(dm, inner_pc, S);
1303 }
1305 };
1306
1307 CHKERR set_pc_ksp(schurDM, get_pc(subksp[1]), S);
1308 CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
1309
1310 CHKERR PetscFree(subksp);
1312}
1313
1314boost::shared_ptr<SetUpSchur>
1316 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1317}
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 >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
double young_modulus
Definition contact.cpp:84
int sigma_order
Definition contact.cpp:82
constexpr AssemblyType AT
Definition contact.cpp:36
constexpr IntegrationType IT
Definition contact.cpp:39
static char help[]
[Check]
Definition contact.cpp:957
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
PetscBool is_plane_strain
Definition contact.cpp:95
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition contact.cpp:74
double spring_stiffness
Definition contact.cpp:87
PetscBool is_large_strain
Definition contact.cpp:94
double rho
Definition contact.cpp:86
int atom_test
Definition contact.cpp:96
#define EXECUTABLE_DIMENSION
Definition contact.cpp:14
PetscBool is_quasi_static
[Operators used for contact]
Definition contact.cpp:78
double alpha_damping
Definition contact.cpp:89
constexpr int SPACE_DIM
Definition contact.cpp:55
double poisson_ratio
Definition contact.cpp:85
double vis_spring_stiffness
Definition contact.cpp:88
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition contact.cpp:62
double scale
Definition contact.cpp:91
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition contact.cpp:72
int contact_order
Definition contact.cpp:81
int order
Definition contact.cpp:80
PetscBool is_axisymmetric
Definition contact.cpp:93
int geom_order
Definition contact.cpp:83
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition contact.cpp:69
@ 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 nme:nme847.
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_NOT_FOUND
Definition definitions.h:33
@ 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:1102
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
@ BLOCK_PRECONDITIONER_SCHUR
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.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
double D
const double v
phase velocity of light in medium (cm/ns)
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
double cn_contact
Definition contact.cpp:99
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)
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
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)
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev)
Assemble domain LHS K factory (LHS first overload) Initializes and pushes operators to assemble the L...
Definition HookeOps.hpp:302
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev, bool is_non_linear=false)
Factory function to create and push internal force operators for the domain RHS. (RHS first overload)...
Definition HookeOps.hpp:242
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:1144
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:1160
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.
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:1086
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:1079
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:55
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:155
boost::shared_ptr< Monitor > monitorPtr
Definition contact.cpp:147
MoFEMErrorCode tsSolve()
Definition contact.cpp:731
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition contact.cpp:143
MoFEMErrorCode runProblem()
[Run problem]
Definition contact.cpp:162
MoFEMErrorCode setupProblem()
[Run problem]
Definition contact.cpp:175
MoFEMErrorCode OPs()
[Boundary condition]
Definition contact.cpp:496
MoFEMErrorCode createCommonData()
[Set up problem]
Definition contact.cpp:326
MoFEMErrorCode checkResults()
[Solve]
Definition contact.cpp:872
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition contact.cpp:146
MoFEM::Interface & mField
Definition contact.cpp:133
MoFEMErrorCode bC()
[Create common data]
Definition contact.cpp:459
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition contact.cpp:142
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition contact.cpp:144
Contact(MoFEM::Interface &m_field)
Definition contact.cpp:128
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
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.
Get values at integration pts for tensor field rank 1, i.e. vector field.
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:1116
MoFEMErrorCode setDiagonalPC(PC pc)
Definition contact.cpp:1266
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition contact.cpp:1022
SmartPetscObj< DM > schurDM
Definition contact.cpp:1037
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
SmartPetscObj< DM > blockDM
Definition contact.cpp:1038
virtual ~SetUpSchurImpl()
Definition contact.cpp:1024
MoFEMErrorCode setPC(PC pc)
Definition contact.cpp:1258
MoFEMErrorCode setOperator()
Definition contact.cpp:1183
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