v0.14.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
24#ifdef PYTHON_SDF
25#include <boost/python.hpp>
26#include <boost/python/def.hpp>
27#include <boost/python/numpy.hpp>
28namespace bp = boost::python;
29namespace np = boost::python::numpy;
30#endif
31
32using namespace MoFEM;
33
34constexpr AssemblyType AT =
35 (SCHUR_ASSEMBLE) ? AssemblyType::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. */
64
65//! [Specialisation for assembly]
66
67// Assemble to A matrix, by default, however, some terms are assembled only to
68// preconditioning.
69
70template <>
74 const EntitiesFieldData::EntData &row_data,
75 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
76 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
77 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
78 };
79
80template <>
84 const EntitiesFieldData::EntData &row_data,
85 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
86 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
87 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
88 };
89
90/**
91 * @brief Element used to specialise assembly
92 *
93 */
95 using BoundaryEleOp::BoundaryEleOp;
96};
97
98/**
99 * @brief Specialise assembly for Stabilised matrix
100 *
101 * @tparam
102 */
103template <>
107 const EntitiesFieldData::EntData &row_data,
108 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
109 return MatSetValues<AssemblyTypeSelector<SCHUR>>(
110 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
111 };
112//! [Specialisation for assembly]
113
115
116//! [Operators used for contact]
120 IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
121//! [Operators used for contact]
122
123PetscBool is_quasi_static = PETSC_TRUE;
124
125int order = 2;
127double young_modulus = 100;
128double poisson_ratio = 0.25;
129double rho = 0.0;
130double spring_stiffness = 0.0;
132double alpha_damping = 0;
133
134double scale = 1.;
135
136PetscBool use_mfront = PETSC_FALSE;
137PetscBool is_axisymmetric = PETSC_FALSE;
138
139int atom_test = 0;
140
141namespace ContactOps {
142double cn_contact = 0.1;
143}; // namespace ContactOps
144
145//#define HENCKY_SMALL_STRAIN
146
147#include <HenckyOps.hpp>
148using namespace HenckyOps;
149#include <ContactOps.hpp>
150#include <PostProcContact.hpp>
151#include <ContactNaturalBC.hpp>
152
153#ifdef WITH_MODULE_MFRONT_INTERFACE
154#include <MFrontMoFEMInterface.hpp>
155#endif
156
166
167using namespace ContactOps;
168
169struct Contact {
170
171 Contact(MoFEM::Interface &m_field) : mField(m_field) {
172 mfrontInterface = nullptr;
173 }
174
176
177private:
179
186
187 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
188 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
189 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
190
191 boost::shared_ptr<GenericElementInterface> mfrontInterface;
192 boost::shared_ptr<Monitor> monitorPtr;
193
194#ifdef PYTHON_SDF
195 boost::shared_ptr<SDFPython> sdfPythonPtr;
196#endif
197
200 double getScale(const double time) {
201 return scale * MoFEM::TimeScale::getScale(time);
202 };
203 };
204};
205
206//! [Run problem]
211 CHKERR bC();
212 CHKERR OPs();
213 CHKERR tsSolve();
216}
217//! [Run problem]
218
219//! [Set up problem]
223
224 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
225 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
226 PETSC_NULL);
227 MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
228 MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
229
230 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_mfront", &use_mfront,
231 PETSC_NULL);
232 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_axisymmetric",
233 &is_axisymmetric, PETSC_NULL);
234 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test, PETSC_NULL);
235
236 // Select base
237 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
238 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
239 PetscInt choice_base_value = AINSWORTH;
240 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
241 LASBASETOPT, &choice_base_value, PETSC_NULL);
242
244 switch (choice_base_value) {
245 case AINSWORTH:
247 MOFEM_LOG("CONTACT", Sev::inform)
248 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
249 break;
250 case DEMKOWICZ:
252 MOFEM_LOG("CONTACT", Sev::inform)
253 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
254 break;
255 default:
256 base = LASTBASE;
257 break;
258 }
259
260 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
261 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
262 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
263 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
264 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
265 SPACE_DIM);
266 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
267 SPACE_DIM);
268 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
269
270 CHKERR simple->setFieldOrder("U", order);
271 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
272
273 auto get_skin = [&]() {
274 Range body_ents;
275 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
276 Skinner skin(&mField.get_moab());
277 Range skin_ents;
278 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
279 return skin_ents;
280 };
281
282 auto filter_blocks = [&](auto skin) {
283 bool is_contact_block = false;
284 Range contact_range;
285 for (auto m :
287
288 (boost::format("%s(.*)") % "CONTACT").str()
289
290 ))
291
292 ) {
293 is_contact_block =
294 true; ///< bloks interation is collectibe, so that is set irrespective
295 ///< if there are enerities in given rank or not in the block
296 MOFEM_LOG("CONTACT", Sev::inform)
297 << "Find contact block set: " << m->getName();
298 auto meshset = m->getMeshset();
299 Range contact_meshset_range;
300 CHKERR mField.get_moab().get_entities_by_dimension(
301 meshset, SPACE_DIM - 1, contact_meshset_range, true);
302
303 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
304 contact_meshset_range);
305 contact_range.merge(contact_meshset_range);
306 }
307 if (is_contact_block) {
308 MOFEM_LOG("SYNC", Sev::inform)
309 << "Nb entities in contact surface: " << contact_range.size();
311 skin = intersect(skin, contact_range);
312 }
313 return skin;
314 };
315
316 auto filter_true_skin = [&](auto skin) {
317 Range boundary_ents;
318 ParallelComm *pcomm =
319 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
320 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
321 PSTATUS_NOT, -1, &boundary_ents);
322 return boundary_ents;
323 };
324
325 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
326 CHKERR simple->setFieldOrder("SIGMA", 0);
327 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
328
329 if (is_axisymmetric) {
330 if (SPACE_DIM == 3) {
331 SETERRQ(PETSC_COMM_SELF, 1,
332 "Use executable contact_2d with axisymmetric model");
333 } else {
334 if (!use_mfront) {
335 SETERRQ(PETSC_COMM_SELF, 1,
336 "Axisymmetric model is only available with MFront (set "
337 "use_mfront to 1)");
338 } else {
339 MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
340 }
341 }
342 } else {
343 if (SPACE_DIM == 2) {
344 MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
345 }
346 }
347
348 if (!use_mfront) {
349 CHKERR simple->setUp();
350 } else {
351#ifndef WITH_MODULE_MFRONT_INTERFACE
352 SETERRQ(
353 PETSC_COMM_SELF, 1,
354 "MFrontInterface module was not found while use_mfront was set to 1");
355#else
356 if (SCHUR_ASSEMBLE) {
357 SETERRQ(PETSC_COMM_SELF, 1,
358 "MFrontInterface module is not compatible with Schur assembly");
359 }
360 if (SPACE_DIM == 3) {
362 boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
363 mField, "U", "GEOMETRY", true, is_quasi_static);
364 } else if (SPACE_DIM == 2) {
365 if (is_axisymmetric) {
367 boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
368 mField, "U", "GEOMETRY", true, is_quasi_static);
369 } else {
370 mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
371 mField, "U", "GEOMETRY", true, is_quasi_static);
372 }
373 }
374#endif
375
376 CHKERR mfrontInterface->getCommandLineParameters();
377 CHKERR mfrontInterface->addElementFields();
378 CHKERR mfrontInterface->createElements();
379
380 CHKERR simple->defineFiniteElements();
381 CHKERR simple->defineProblem(PETSC_TRUE);
382 CHKERR simple->buildFields();
383 CHKERR simple->buildFiniteElements();
384
386 CHKERR mfrontInterface->addElementsToDM(simple->getDM());
387
388 CHKERR simple->buildProblem();
389 }
390
391 auto dm = simple->getDM();
392 monitorPtr = boost::make_shared<Monitor>(dm, use_mfront, mfrontInterface,
394 if (use_mfront) {
395 mfrontInterface->setMonitorPtr(monitorPtr);
396 }
397
398 auto project_ho_geometry = [&]() {
399 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
400 return mField.loop_dofs("GEOMETRY", ent_method);
401 };
402
403 PetscBool project_geometry = PETSC_TRUE;
404 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
405 &project_geometry, PETSC_NULL);
406 if (project_geometry){
407 CHKERR project_ho_geometry();
408 }
409
411} //! [Set up problem]
412
413//! [Create common data]
416
417 auto get_options = [&]() {
419 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
420 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
421 &young_modulus, PETSC_NULL);
422 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
423 &poisson_ratio, PETSC_NULL);
424 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
425 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn_contact,
426 PETSC_NULL);
427 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
428 &spring_stiffness, PETSC_NULL);
429 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-vis_spring_stiffness",
430 &vis_spring_stiffness, PETSC_NULL);
431 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
432 &alpha_damping, PETSC_NULL);
433
434 if (!use_mfront) {
435 MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
436 MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
437 } else {
438 MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
439 }
440
441 MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
442 MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
443 MOFEM_LOG("CONTACT", Sev::inform)
444 << "Spring stiffness " << spring_stiffness;
445 MOFEM_LOG("CONTACT", Sev::inform)
446 << "Vis spring_stiffness " << vis_spring_stiffness;
447
448 MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
449
450 PetscBool use_scale = PETSC_FALSE;
451 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_scale", &use_scale,
452 PETSC_NULL);
453 if (use_scale) {
455 }
456
457 MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
458
459 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_quasi_static",
460 &is_quasi_static, PETSC_NULL);
461 MOFEM_LOG("CONTACT", Sev::inform)
462 << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
463
465 };
466
467 CHKERR get_options();
468
469#ifdef PYTHON_SDF
470 char sdf_file_name[255];
471 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
472 sdf_file_name, 255, PETSC_NULL);
473
474 sdfPythonPtr = boost::make_shared<SDFPython>();
475 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
476 sdfPythonWeakPtr = sdfPythonPtr;
477#endif
478
480}
481//! [Create common data]
482
483//! [Boundary condition]
486 auto bc_mng = mField.getInterface<BcManager>();
488
489 for (auto f : {"U", "SIGMA"}) {
490 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
491 "REMOVE_X", f, 0, 0);
492 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
493 "REMOVE_Y", f, 1, 1);
494 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
495 "REMOVE_Z", f, 2, 2);
496 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
497 "REMOVE_ALL", f, 0, 3);
498 }
499
500 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
501 "SIGMA", 0, 0, false, true);
502 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
503 "SIGMA", 1, 1, false, true);
504 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
505 "SIGMA", 2, 2, false, true);
506 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
507 "SIGMA", 0, 3, false, true);
508 CHKERR bc_mng->removeBlockDOFsOnEntities(
509 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
510
511 if (is_axisymmetric) {
512 CHKERR bc_mng->removeBlockDOFsOnEntities(
513 simple->getProblemName(), "REMOVE_X", "SIGMA", 0, 3, false, true);
514 }
515
516 // Note remove has to be always before push. Then node marking will be
517 // corrupted.
518 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
519 simple->getProblemName(), "U");
520
522}
523//! [Boundary condition]
524
525//! [Push operators to pip]
529 auto *pip_mng = mField.getInterface<PipelineManager>();
530 auto bc_mng = mField.getInterface<BcManager>();
531 auto time_scale = boost::make_shared<ScaledTimeScale>();
532 auto body_force_time_scale =
533 boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
534
535 auto integration_rule_vol = [](int, int, int approx_order) {
536 return 2 * approx_order + geom_order - 1;
537 };
538 auto integration_rule_boundary = [](int, int, int approx_order) {
539 return 2 * approx_order + geom_order - 1;
540 };
541
542 auto add_domain_base_ops = [&](auto &pip) {
545 "GEOMETRY");
547 };
548
549 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
550 henky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
551 henky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
552
553 auto add_domain_ops_lhs = [&](auto &pip) {
555
556 //! [Only used for dynamics]
559 //! [Only used for dynamics]
560 if (is_quasi_static == PETSC_FALSE) {
561
562 auto *pip_mng = mField.getInterface<PipelineManager>();
563 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
564
565 auto get_inertia_and_mass_damping =
566 [this, fe_domain_lhs](const double, const double, const double) {
567 return (rho * scale) * fe_domain_lhs->ts_aa +
568 (alpha_damping * scale) * fe_domain_lhs->ts_a;
569 };
570 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
571 } else {
572
573 auto *pip_mng = mField.getInterface<PipelineManager>();
574 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
575
576 auto get_inertia_and_mass_damping =
577 [this, fe_domain_lhs](const double, const double, const double) {
578 return (alpha_damping * scale) * fe_domain_lhs->ts_a;
579 };
580 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
581
582 }
583
584 if (!use_mfront) {
585 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
586 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
587 }
588
590 };
591
592 auto add_domain_ops_rhs = [&](auto &pip) {
594
596 pip, mField, "U", {body_force_time_scale}, Sev::inform);
597
598 //! [Only used for dynamics]
601 //! [Only used for dynamics]
602
603 // only in case of dynamics
604 if (is_quasi_static == PETSC_FALSE) {
605 auto mat_acceleration = boost::make_shared<MatrixDouble>();
607 "U", mat_acceleration));
608 pip.push_back(
609 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
610 return rho * scale;
611 }));
612
613 }
614
615 // only in case of viscosity
616 if (alpha_damping > 0) {
617 auto mat_velocity = boost::make_shared<MatrixDouble>();
618 pip.push_back(
619 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
620 pip.push_back(
621 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
622 return alpha_damping * scale;
623 }));
624 }
625
626 if (!use_mfront) {
627 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
628 mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
629 }
630
631 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
632 pip, "SIGMA", "U", is_axisymmetric);
633
635 };
636
637 auto add_boundary_base_ops = [&](auto &pip) {
640 "GEOMETRY");
641 // We have to integrate on curved face geometry, thus integration weight
642 // have to adjusted.
643 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
645 };
646
647 auto add_boundary_ops_lhs = [&](auto &pip) {
649
650 //! [Operators used for contact]
653 //! [Operators used for contact]
654
655 // Add Natural BCs to LHS
657 pip, mField, "U", Sev::inform);
658
659 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
660
661 auto *pip_mng = mField.getInterface<PipelineManager>();
662 auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
663
664 pip.push_back(new OpSpringLhs(
665 "U", "U",
666
667 [this, fe_boundary_lhs](double, double, double) {
668 return spring_stiffness * scale +
669 (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
670 }
671
672 ));
673 }
674
675 CHKERR
676 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
677 pip, "SIGMA", "U", is_axisymmetric);
679 DomainEle>(
680 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
681 integration_rule_vol, is_axisymmetric);
682
684 };
685
686 auto add_boundary_ops_rhs = [&](auto &pip) {
688
689 //! [Operators used for contact]
692 //! [Operators used for contact]
693
694 // Add Natural BCs to RHS
696 pip, mField, "U", {time_scale}, Sev::inform);
697
698 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
699 auto u_disp = boost::make_shared<MatrixDouble>();
700 auto dot_u_disp = boost::make_shared<MatrixDouble>();
701 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
702 pip.push_back(
703 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
704 pip.push_back(
705 new OpSpringRhs("U", u_disp, [this](double, double, double) {
706 return spring_stiffness * scale;
707 }));
708 pip.push_back(
709 new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
711 }));
712 }
713
714 CHKERR
715 ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
716 pip, "SIGMA", "U", is_axisymmetric);
717
719 };
720
721 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
722 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
723 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
724 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
725
726 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
727 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
728 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
729 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
730
731 if (use_mfront) {
732 auto t_type = GenericElementInterface::IM2;
733 if (is_quasi_static == PETSC_TRUE)
735
736 mfrontInterface->setOperators();
737 mfrontInterface->setupSolverFunctionTS(t_type);
738 mfrontInterface->setupSolverJacobianTS(t_type);
739 }
740
741 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
742 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
743 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
744 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
745
747}
748//! [Push operators to pip]
749
750//! [Solve]
751struct SetUpSchur {
752 static boost::shared_ptr<SetUpSchur>
755
756protected:
757 SetUpSchur() = default;
758};
759
762
765 ISManager *is_manager = mField.getInterface<ISManager>();
766
767 auto set_section_monitor = [&](auto solver) {
769 SNES snes;
770 CHKERR TSGetSNES(solver, &snes);
771 PetscViewerAndFormat *vf;
772 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
773 PETSC_VIEWER_DEFAULT, &vf);
774 CHKERR SNESMonitorSet(
775 snes,
776 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
777 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
779 };
780
781 auto scatter_create = [&](auto D, auto coeff) {
783 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
784 ROW, "U", coeff, coeff, is);
785 int loc_size;
786 CHKERR ISGetLocalSize(is, &loc_size);
787 Vec v;
788 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
789 VecScatter scatter;
790 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
791 return std::make_tuple(SmartPetscObj<Vec>(v),
793 };
794
795 auto set_time_monitor = [&](auto dm, auto solver) {
797 monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
798 boost::shared_ptr<ForcesAndSourcesCore> null;
799 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
800 monitorPtr, null, null);
802 };
803
804 auto set_essential_bc = [&]() {
806 // This is low level pushing finite elements (pipelines) to solver
807 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
808 auto pre_proc_ptr = boost::make_shared<FEMethod>();
809 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
810 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
811
812 // Add boundary condition scaling
813 auto time_scale = boost::make_shared<TimeScale>();
814
815 auto get_bc_hook_rhs = [&]() {
817 {time_scale}, false);
818 return hook;
819 };
820 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
821
822 auto get_post_proc_hook_rhs = [&]() {
824 mField, post_proc_rhs_ptr, 1.);
825 };
826 auto get_post_proc_hook_lhs = [&]() {
828 mField, post_proc_lhs_ptr, 1.);
829 };
830 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
831
832 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
833 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
834 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
835 if (AT != AssemblyType::SCHUR) {
836 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
837 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
838 }
840 };
841
842 auto set_schur_pc = [&](auto solver) {
843 boost::shared_ptr<SetUpSchur> schur_ptr;
844 if (AT == AssemblyType::SCHUR) {
846 CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
847 }
848 return schur_ptr;
849 };
850
851 auto dm = simple->getDM();
852 auto D = createDMVector(dm);
853
855
856 uXScatter = scatter_create(D, 0);
857 uYScatter = scatter_create(D, 1);
858 if (SPACE_DIM == 3)
859 uZScatter = scatter_create(D, 2);
860
861 // Add extra finite elements to SNES solver pipelines to resolve essential
862 // boundary conditions
863 CHKERR set_essential_bc();
864
865 if (is_quasi_static == PETSC_TRUE) {
866 auto solver = pip_mng->createTSIM();
867 CHKERR TSSetFromOptions(solver);
868
869 auto D = createDMVector(dm);
870 auto schur_pc_ptr = set_schur_pc(solver);
871
872 CHKERR set_section_monitor(solver);
873 CHKERR set_time_monitor(dm, solver);
874 CHKERR TSSetSolution(solver, D);
875 CHKERR TSSetUp(solver);
876 CHKERR TSSolve(solver, NULL);
877 } else {
878 auto solver = pip_mng->createTSIM2();
879 CHKERR TSSetFromOptions(solver);
880
881 auto dm = simple->getDM();
882 auto D = createDMVector(dm);
883 auto DD = vectorDuplicate(D);
884 auto schur_pc_ptr = set_schur_pc(solver);
885
886 CHKERR set_section_monitor(solver);
887 CHKERR set_time_monitor(dm, solver);
888 CHKERR TS2SetSolution(solver, D, DD);
889 CHKERR TSSetUp(solver);
890 CHKERR TSSolve(solver, NULL);
891 }
892
894
896}
897//! [Solve]
898
899//! [Check]
901//! [Check]
902
903static char help[] = "...\n\n";
904
905int main(int argc, char *argv[]) {
906
907#ifdef PYTHON_SDF
908 Py_Initialize();
909 np::initialize();
910#endif
911
912 // Initialisation of MoFEM/PETSc and MOAB data structures
913 const char param_file[] = "param_file.petsc";
915
916 // Add logging channel for CONTACT
917 auto core_log = logging::core::get();
918 core_log->add_sink(
920 LogManager::setLog("CONTACT");
921 MOFEM_LOG_TAG("CONTACT", "Indent");
922
923 try {
924
925 //! [Register MoFEM discrete manager in PETSc]
926 DMType dm_name = "DMMOFEM";
927 CHKERR DMRegister_MoFEM(dm_name);
928 //! [Register MoFEM discrete manager in PETSc
929
930 //! [Create MoAB]
931 moab::Core mb_instance; ///< mesh database
932 moab::Interface &moab = mb_instance; ///< mesh database interface
933 //! [Create MoAB]
934
935 //! [Create MoFEM]
936 MoFEM::Core core(moab); ///< finite element database
937 MoFEM::Interface &m_field = core; ///< finite element database interface
938 //! [Create MoFEM]
939
940 //! [Load mesh]
941 Simple *simple = m_field.getInterface<Simple>();
942 CHKERR simple->getOptions();
943 CHKERR simple->loadFile("");
944 //! [Load mesh]
945
946 //! [CONTACT]
947 Contact ex(m_field);
948 CHKERR ex.runProblem();
949 //! [CONTACT]
950 }
952
954
955#ifdef PYTHON_SDF
956 if (Py_FinalizeEx() < 0) {
957 exit(120);
958 }
959#endif
960
961 return 0;
962}
963
964struct SetUpSchurImpl : public SetUpSchur {
965
967
968 virtual ~SetUpSchurImpl() {
969 A.reset();
970 P.reset();
971 S.reset();
972 }
973
975
976private:
979 MoFEMErrorCode setPC(PC pc);
980
982
986
988
990};
991
995 auto pip = mField.getInterface<PipelineManager>();
996 auto dm = simple->getDM();
997
998 SNES snes;
999 CHKERR TSGetSNES(solver, &snes);
1000 KSP ksp;
1001 CHKERR SNESGetKSP(snes, &ksp);
1002 CHKERR KSPSetFromOptions(ksp);
1003
1004 PC pc;
1005 CHKERR KSPGetPC(ksp, &pc);
1006
1007 PetscBool is_pcfs = PETSC_FALSE;
1008 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1009 if (is_pcfs) {
1010
1011 MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
1012
1013 if (A || P || S) {
1016 "It is expected that Schur matrix is not allocated. This is "
1017 "possible only if PC is set up twice");
1018 }
1019
1020 A = createDMMatrix(dm);
1021 P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1022 subDM = createSubDM();
1024 CHKERR MatSetBlockSize(S, SPACE_DIM);
1025
1026 auto ts_ctx_ptr = getDMTsCtx(dm);
1027 CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1028
1030 CHKERR setPC(pc);
1031
1032 } else {
1033 MOFEM_LOG("CONTACT", Sev::inform) << "No Schur pc";
1034
1035 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1036 pip->getOpBoundaryLhsPipeline().push_back(
1037 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1038 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1039 pip->getOpDomainLhsPipeline().push_back(
1040 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1041
1042 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1043 post_proc_schur_lhs_ptr->postProcessHook = [this,
1044 post_proc_schur_lhs_ptr]() {
1047 mField, post_proc_schur_lhs_ptr, 1.)();
1049 };
1050 auto ts_ctx_ptr = getDMTsCtx(dm);
1051 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1052 }
1054}
1055
1057 auto simple = mField.getInterface<Simple>();
1058 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1059 auto set_up = [&]() {
1061 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
1062 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1063 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1064 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
1065 CHKERR DMSetUp(sub_dm);
1067 };
1068 CHK_THROW_MESSAGE(set_up(), "sey up dm");
1069 return sub_dm;
1070}
1071
1074
1075 double eps_stab = 1e-4;
1076 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1077 PETSC_NULL);
1078
1079 using B =
1081 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1082
1083 auto pip = mField.getInterface<PipelineManager>();
1084 // Boundary
1085 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
1086 auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
1087
1088 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1089 pip->getOpBoundaryLhsPipeline().push_back(
1090 new OpMassStab("SIGMA", "SIGMA",
1091 [eps_stab](double, double, double) { return eps_stab; }));
1092 pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1093 {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1094
1095 // Domain
1096 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1097 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1098 {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1099
1100 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1101 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1102
1103 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1105 CHKERR MatZeroEntries(A);
1106 CHKERR MatZeroEntries(P);
1107 CHKERR MatZeroEntries(S);
1108 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1110 };
1111
1112 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1113 post_proc_schur_lhs_ptr]() {
1115 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1116
1117 *post_proc_schur_lhs_ptr->matAssembleSwitch = false;
1118
1119 auto print_mat_norm = [this](auto a, std::string prefix) {
1121 double nrm;
1122 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1123 MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1125 };
1126
1127 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1128 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1130 mField, post_proc_schur_lhs_ptr, 1, A)();
1131
1132 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1133 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1134 CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1135
1136 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1137 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1138
1140 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1141
1142#ifndef NDEBUG
1143 CHKERR print_mat_norm(A, "A");
1144 CHKERR print_mat_norm(P, "P");
1145 CHKERR print_mat_norm(S, "S");
1146#endif // NDEBUG
1147
1148 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1150 };
1151
1152 auto simple = mField.getInterface<Simple>();
1153 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1154 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1155 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1156
1158}
1159
1162 auto simple = mField.getInterface<Simple>();
1164 mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1165 simple->getProblemName(), ROW, "SIGMA", 0, SPACE_DIM, is);
1166 CHKERR PCFieldSplitSetIS(pc, NULL, is);
1167 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1169}
1170
1171boost::shared_ptr<SetUpSchur>
1173 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1174}
std::string param_file
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
double young_modulus
Definition: contact.cpp:127
constexpr AssemblyType AT
Definition: contact.cpp:34
constexpr IntegrationType IT
Definition: contact.cpp:37
static char help[]
[Check]
Definition: contact.cpp:903
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
double spring_stiffness
Definition: contact.cpp:130
double rho
Definition: contact.cpp:129
int atom_test
Definition: contact.cpp:139
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:14
PetscBool is_quasi_static
[Operators used for contact]
Definition: contact.cpp:123
double alpha_damping
Definition: contact.cpp:132
PetscBool use_mfront
Definition: contact.cpp:136
constexpr int SPACE_DIM
Definition: contact.cpp:53
double poisson_ratio
Definition: contact.cpp:128
double vis_spring_stiffness
Definition: contact.cpp:131
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: contact.cpp:60
double scale
Definition: contact.cpp:134
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:118
int order
Definition: contact.cpp:125
PetscBool is_axisymmetric
Definition: contact.cpp:137
int geom_order
Definition: contact.cpp:126
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:120
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition: contact.cpp:114
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
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.
Definition: definitions.h:595
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
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
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:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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)
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
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.
Definition: Exceptions.hpp:56
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:1042
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:165
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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.
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061
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)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
static constexpr int approx_order
Element used to specialise assembly.
Definition: contact.cpp:94
double getScale(const double time)
Get scaling at a given time.
Definition: contact.cpp:200
boost::shared_ptr< Monitor > monitorPtr
Definition: contact.cpp:192
MoFEMErrorCode tsSolve()
Definition: contact.cpp:760
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:188
MoFEMErrorCode runProblem()
[Run problem]
Definition: contact.cpp:207
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:220
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:526
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:414
MoFEMErrorCode checkResults()
[Solve]
Definition: contact.cpp:900
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: contact.cpp:191
MoFEM::Interface & mField
Definition: contact.cpp:178
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:484
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:187
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:189
Contact(MoFEM::Interface &m_field)
Definition: contact.cpp:171
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() 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:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
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.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
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)
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
Definition: plastic.cpp:762
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1482
virtual MoFEMErrorCode setUp(SmartPetscObj< TS > solver)=0
SmartPetscObj< Mat > A
Definition: contact.cpp:983
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: contact.cpp:966
SmartPetscObj< DM > createSubDM()
Definition: contact.cpp:1056
MoFEMErrorCode setEntities()
Definition: elastic.cpp:942
virtual ~SetUpSchurImpl()
Definition: contact.cpp:968
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1160
MoFEMErrorCode setOperator()
Definition: contact.cpp:1072
MoFEM::Interface & mField