v0.15.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
Contact Struct Reference
Collaboration diagram for Contact:
[legend]

Classes

struct  ScaledTimeScale
 

Public Member Functions

 Contact (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run problem]
 

Private Member Functions

MoFEMErrorCode setupProblem ()
 [Run problem]
 
MoFEMErrorCode createCommonData ()
 [Set up problem]
 
MoFEMErrorCode bC ()
 [Create common data]
 
MoFEMErrorCode OPs ()
 [Boundary condition]
 
MoFEMErrorCode tsSolve ()
 
MoFEMErrorCode checkResults ()
 [Solve]
 

Private Attributes

MoFEM::InterfacemField
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
boost::shared_ptr< MFrontInterfacemfrontInterfacePtr
 
boost::shared_ptr< MonitormonitorPtr
 

Detailed Description

Definition at line 120 of file contact.cpp.

Constructor & Destructor Documentation

◆ Contact()

Contact::Contact ( MoFEM::Interface m_field)
inline

Definition at line 122 of file contact.cpp.

123 : mField(m_field), mfrontInterfacePtr(nullptr) {}
boost::shared_ptr< MFrontInterface > mfrontInterfacePtr
Definition contact.cpp:141
MoFEM::Interface & mField
Definition contact.cpp:128

Member Function Documentation

◆ bC()

MoFEMErrorCode Contact::bC ( )
private

[Create common data]

[Boundary condition]

Definition at line 443 of file contact.cpp.

443 {
445 auto bc_mng = mField.getInterface<BcManager>();
447
448 for (auto f : {"U", "SIGMA"}) {
449 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
450 "REMOVE_X", f, 0, 0);
451 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
452 "REMOVE_Y", f, 1, 1);
453 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
454 "REMOVE_Z", f, 2, 2);
455 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
456 "REMOVE_ALL", f, 0, 3);
457 }
458
459 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
460 "SIGMA", 0, 0, false, true);
461 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
462 "SIGMA", 1, 1, false, true);
463 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
464 "SIGMA", 2, 2, false, true);
465 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
466 "SIGMA", 0, 3, false, true);
467 CHKERR bc_mng->removeBlockDOFsOnEntities(
468 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
469
470 // Note remove has to be always before push. Then node marking will be
471 // corrupted.
472 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
473 simple->getProblemName(), "U");
474
476}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
Boundary condition manager for finite element problem setup.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

◆ checkResults()

MoFEMErrorCode Contact::checkResults ( )
private

[Solve]

[Check]

Definition at line 860 of file contact.cpp.

860 {
862 if (atom_test && !mField.get_comm_rank()) {
863 const double *t_ptr;
864 CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
865 double hertz_force;
866 double fem_force;
867 double analytical_active_area = 1.0;
868 double norm = 1e-5;
869 double tol_force = 1e-3;
870 double tol_norm = 7.5; // change when analytical functions are updated
871 double tol_area = 3e-2;
872 double fem_active_area = t_ptr[3];
873
874 switch (atom_test) {
875 case 1: // plane stress
876 hertz_force = 3.927;
877 fem_force = t_ptr[1];
878 break;
879
880 case 2: // plane strain
881 hertz_force = 4.675;
882 fem_force = t_ptr[1];
883 norm = monitorPtr->getErrorNorm(1);
884 break;
885
886 case 3: // Hertz 3D
887 hertz_force = 3.968;
888 tol_force = 2e-3;
889 fem_force = t_ptr[2];
890 analytical_active_area = M_PI / 4;
891 tol_area = 0.2;
892 break;
893
894 case 4: // axisymmetric
895 tol_force = 5e-3;
896 tol_area = 0.2;
897 // analytical_active_area = M_PI;
898
899 case 5: // axisymmetric
900 hertz_force = 15.873;
901 tol_force = 5e-3;
902 fem_force = t_ptr[1];
903 norm = monitorPtr->getErrorNorm(1);
904 analytical_active_area = M_PI;
905 break;
906
907 case 6: // wavy 2d
908 hertz_force = 0.374;
909 fem_force = t_ptr[1];
910 break;
911
912 case 7: // wavy 3d
913 hertz_force = 0.5289;
914 fem_force = t_ptr[2];
915 break;
916
917 default:
918 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
919 "atom test %d does not exist", atom_test);
920 }
921 if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
922 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
923 "atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
924 atom_test, fem_force, hertz_force);
925 }
926 if (norm > tol_norm) {
927 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
928 "atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
929 atom_test, norm, tol_norm);
930 }
931 if (fabs(fem_active_area - analytical_active_area) > tol_area) {
932 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
933 "atom test %d failed: AREA computed %3.4e but should be %3.4e",
934 atom_test, fem_active_area, analytical_active_area);
935 }
936 CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
937 }
938
940
942}
int atom_test
Definition contact.cpp:94
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_INVALID_DATA
Definition definitions.h:36
static SmartPetscObj< Vec > totalTraction
boost::shared_ptr< Monitor > monitorPtr
Definition contact.cpp:142
virtual int get_comm_rank() const =0

◆ createCommonData()

MoFEMErrorCode Contact::createCommonData ( )
private

[Set up problem]

[Create common data]

Definition at line 310 of file contact.cpp.

310 {
312
313 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strain",
314 &is_large_strain, PETSC_NULLPTR);
315 PetscBool use_mfront = PETSC_FALSE;
316 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_mfront", &use_mfront,
317 PETSC_NULLPTR);
318 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_axisymmetric",
319 &is_axisymmetric, PETSC_NULLPTR);
320 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
321 PETSC_NULLPTR);
322
323 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale,
324 PETSC_NULLPTR);
325 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
326 &young_modulus, PETSC_NULLPTR);
327 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
328 &poisson_ratio, PETSC_NULLPTR);
329 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
330 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn", &cn_contact,
331 PETSC_NULLPTR);
332 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-spring_stiffness",
333 &spring_stiffness, PETSC_NULLPTR);
334 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-vis_spring_stiffness",
335 &vis_spring_stiffness, PETSC_NULLPTR);
336 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
337 &alpha_damping, PETSC_NULLPTR);
338
339 if (!use_mfront) {
340 if (!is_large_strain) {
341 MOFEM_LOG("CONTACT", Sev::inform)
342 << "Selected material model: Hooke (small strain)";
343 } else {
344 MOFEM_LOG("CONTACT", Sev::inform)
345 << "Selected material model: Hencky (finite strain)";
346 }
347 MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
348 MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
349 } else {
350 MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
351 }
352
353 MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
354 MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
355 MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
356 MOFEM_LOG("CONTACT", Sev::inform)
357 << "Vis spring_stiffness " << vis_spring_stiffness;
358
359 MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
360
361 PetscBool use_scale = PETSC_FALSE;
362 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_scale", &use_scale,
363 PETSC_NULLPTR);
364 if (use_scale) {
366 }
367
368 MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
369
370 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_quasi_static",
371 &is_quasi_static, PETSC_NULLPTR);
372 MOFEM_LOG("CONTACT", Sev::inform)
373 << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
374
375#ifdef ENABLE_PYTHON_BINDING
376 auto file_exists = [](std::string myfile) {
377 std::ifstream file(myfile.c_str());
378 if (file) {
379 return true;
380 }
381 return false;
382 };
383 char sdf_file_name[255] = "sdf.py";
384 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
385 sdf_file_name, 255, PETSC_NULLPTR);
386
387 if (file_exists(sdf_file_name)) {
388 MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
389 sdfPythonPtr = boost::make_shared<SDFPython>();
390 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
391 sdfPythonWeakPtr = sdfPythonPtr;
392 } else {
393 MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
394 }
395#endif
396
397 if (is_axisymmetric) {
398 if (SPACE_DIM == 3) {
399 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
400 "Use executable contact_2d with axisymmetric model");
401 } else {
402 if (!use_mfront) {
403 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
404 "Axisymmetric model is only available with MFront (set "
405 "use_mfront to 1)");
406 } else {
407 MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
408 }
409 }
410 } else {
411 if (SPACE_DIM == 2) {
412 MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
413 }
414 }
415
416 if (use_mfront) {
417 if (SPACE_DIM == 3) {
419 } else if (SPACE_DIM == 2) {
420 if (is_axisymmetric) {
422 } else {
424 }
425 }
426 CHKERR mfrontInterfacePtr->getCommandLineParameters();
427 }
428
430 auto dm = simple->getDM();
431 monitorPtr = boost::make_shared<Monitor>(dm, scale, mfrontInterfacePtr,
433
434 if (mfrontInterfacePtr) {
435 mfrontInterfacePtr->setMonitor(monitorPtr);
436 }
437
439}
double young_modulus
Definition contact.cpp:82
constexpr AssemblyType AT
Definition contact.cpp:34
double spring_stiffness
Definition contact.cpp:85
PetscBool is_large_strain
Definition contact.cpp:92
double rho
Definition contact.cpp:84
PetscBool is_quasi_static
[Operators used for contact]
Definition contact.cpp:76
double alpha_damping
Definition contact.cpp:87
constexpr int SPACE_DIM
Definition contact.cpp:53
double poisson_ratio
Definition contact.cpp:83
double vis_spring_stiffness
Definition contact.cpp:86
double scale
Definition contact.cpp:89
PetscBool is_axisymmetric
Definition contact.cpp:91
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MOFEM_LOG(channel, severity)
Log.
double cn_contact
Definition contact.cpp:97
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)
boost::shared_ptr< MFrontInterface > createMFrontInterface(MoFEM::Interface &m_field, ModelHypothesis mh, AssemblyType at=AssemblyType::PETSC)
create mfront interface
@ AXISYMMETRICAL
Axisymmetrical model hypothesis.
@ PLANESTRAIN
Plane strain model hypothesis.
@ TRIDIMENSIONAL
3D model hypothesis.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800

◆ OPs()

MoFEMErrorCode Contact::OPs ( )
private

[Boundary condition]

[Push operators to pip]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

[Operators used for contact]

[Operators used for contact]

[Operators used for contact]

[Operators used for contact]

Definition at line 480 of file contact.cpp.

480 {
483 auto *pip_mng = mField.getInterface<PipelineManager>();
484 auto time_scale = boost::make_shared<ScaledTimeScale>();
485 auto body_force_time_scale =
486 boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
487
488 auto integration_rule_vol = [](int, int, int approx_order) {
489 return 2 * approx_order + geom_order - 1;
490 };
491 auto integration_rule_boundary = [](int, int, int approx_order) {
492 return 2 * approx_order + geom_order - 1;
493 };
494
495 auto add_domain_base_ops = [&](auto &pip) {
498 "GEOMETRY");
500 };
501
502 auto add_domain_ops_lhs = [&](auto &pip) {
504
505 //! [Only used for dynamics]
508 //! [Only used for dynamics]
509 if (is_quasi_static == PETSC_FALSE) {
510
511 auto *pip_mng = mField.getInterface<PipelineManager>();
512 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
513
514 auto get_inertia_and_mass_damping =
515 [this, fe_domain_lhs](const double, const double, const double) {
516 return (rho * scale) * fe_domain_lhs->ts_aa +
517 (alpha_damping * scale) * fe_domain_lhs->ts_a;
518 };
519 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
520 } else {
521
522 auto *pip_mng = mField.getInterface<PipelineManager>();
523 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
524
525 auto get_inertia_and_mass_damping =
526 [this, fe_domain_lhs](const double, const double, const double) {
527 return (alpha_damping * scale) * fe_domain_lhs->ts_a;
528 };
529 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
530 }
531
532 if (!mfrontInterfacePtr) {
533 if (!is_large_strain) {
534 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
535 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
536 } else {
537 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
538 mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
539 }
540 } else {
541 CHKERR mfrontInterfacePtr->opFactoryDomainLhs(pip, "U");
542 }
543
545 };
546
547 auto add_domain_ops_rhs = [&](auto &pip) {
549
551 pip, mField, "U", {body_force_time_scale}, Sev::inform);
552
553 //! [Only used for dynamics]
556 //! [Only used for dynamics]
557
558 // only in case of dynamics
559 if (is_quasi_static == PETSC_FALSE) {
560 auto mat_acceleration = boost::make_shared<MatrixDouble>();
562 "U", mat_acceleration));
563 pip.push_back(
564 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
565 return rho * scale;
566 }));
567 }
568
569 // only in case of viscosity
570 if (alpha_damping > 0) {
571 auto mat_velocity = boost::make_shared<MatrixDouble>();
572 pip.push_back(
573 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
574 pip.push_back(
575 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
576 return alpha_damping * scale;
577 }));
578 }
579
580 if (!mfrontInterfacePtr) {
581 if (!is_large_strain) {
582 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
583 mField, pip, "U", "MAT_ELASTIC", Sev::inform, 1, scale);
584 } else {
585 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
586 mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
587 }
588 } else {
589 CHKERR mfrontInterfacePtr->opFactoryDomainRhs(pip, "U");
590 }
591
592 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
593 pip, "SIGMA", "U", is_axisymmetric);
594
596 };
597
598 auto add_boundary_base_ops = [&](auto &pip) {
601 "GEOMETRY");
602 // We have to integrate on curved face geometry, thus integration weight
603 // have to adjusted.
604 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
606 };
607
608 auto add_boundary_ops_lhs = [&](auto &pip) {
610
611 //! [Operators used for contact]
614 //! [Operators used for contact]
615
616 // Add Natural BCs to LHS
618 pip, mField, "U", Sev::inform);
619
620 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
621
622 auto *pip_mng = mField.getInterface<PipelineManager>();
623 auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
624
625 pip.push_back(new OpSpringLhs(
626 "U", "U",
627
628 [this, fe_boundary_lhs](double, double, double) {
629 return spring_stiffness * scale +
630 (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
631 }
632
633 ));
634 }
635
636 CHKERR
637 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
638 pip, "SIGMA", "U", is_axisymmetric);
640 DomainEle>(
641 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
642 integration_rule_vol, is_axisymmetric);
643
645 };
646
647 auto add_boundary_ops_rhs = [&](auto &pip) {
649
650 //! [Operators used for contact]
653 //! [Operators used for contact]
654
655 // Add Natural BCs to RHS
657 pip, mField, "U", {time_scale}, Sev::inform);
658
659 if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
660 auto u_disp = boost::make_shared<MatrixDouble>();
661 auto dot_u_disp = boost::make_shared<MatrixDouble>();
662 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
663 pip.push_back(
664 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
665 pip.push_back(
666 new OpSpringRhs("U", u_disp, [this](double, double, double) {
667 return spring_stiffness * scale;
668 }));
669 pip.push_back(
670 new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
672 }));
673 }
674
675 CHKERR
676 ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
677 pip, "SIGMA", "U", is_axisymmetric);
678
680 };
681
682 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
683 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
684 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
685 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
686
687 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
688 CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
689 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
690 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
691
692 if (mfrontInterfacePtr) {
693 CHKERR mfrontInterfacePtr->setUpdateInternalVariablesOperators(
694 integration_rule_vol, "U");
695 CHKERR mfrontInterfacePtr->setPostProcessOperators(
696 integration_rule_vol, simple->getDomainFEName(), "U", order);
697 }
698
699 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
700 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
701 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
702 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
703
705}
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition contact.cpp:71
int order
Definition contact.cpp:78
int geom_order
Definition contact.cpp:81
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition contact.cpp:73
@ H1
continuous field
Definition definitions.h:85
@ HDIV
field with continuous normal traction
Definition definitions.h:87
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
@ GAUSS
Gaussian quadrature integration.
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)
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56
Add operators pushing bases from local to physical configuration.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
PipelineManager interface.

◆ runProblem()

MoFEMErrorCode Contact::runProblem ( )

[Run problem]

Definition at line 157 of file contact.cpp.

157 {
161 CHKERR bC();
162 CHKERR OPs();
163 CHKERR tsSolve();
166}
MoFEMErrorCode tsSolve()
Definition contact.cpp:718
MoFEMErrorCode setupProblem()
[Run problem]
Definition contact.cpp:170
MoFEMErrorCode OPs()
[Boundary condition]
Definition contact.cpp:480
MoFEMErrorCode createCommonData()
[Set up problem]
Definition contact.cpp:310
MoFEMErrorCode checkResults()
[Solve]
Definition contact.cpp:860
MoFEMErrorCode bC()
[Create common data]
Definition contact.cpp:443

◆ setupProblem()

MoFEMErrorCode Contact::setupProblem ( )
private

[Run problem]

[Set up problem]

< blocs interation is collective, so that is set irrespective if there are entities in given rank or not in the block

Definition at line 170 of file contact.cpp.

170 {
173
174 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
175 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-contact_order", &contact_order,
176 PETSC_NULLPTR);
177 sigma_order = std::max(order, contact_order) - 1;
178 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-sigma_order", &sigma_order,
179 PETSC_NULLPTR);
180 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
181 PETSC_NULLPTR);
182
183 MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
184 MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
185 MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
186 MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
187
188 // Select base
189 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
190 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
191 PetscInt choice_base_value = AINSWORTH;
192 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
193 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
194
196 switch (choice_base_value) {
197 case AINSWORTH:
199 MOFEM_LOG("CONTACT", Sev::inform)
200 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
201 break;
202 case DEMKOWICZ:
204 MOFEM_LOG("CONTACT", Sev::inform)
205 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
206 break;
207 default:
208 base = LASTBASE;
209 break;
210 }
211
212 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
213 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
214 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
215 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
216 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
217 SPACE_DIM);
218 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
219 SPACE_DIM);
220 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
221
222 CHKERR simple->setFieldOrder("U", order);
223 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
224
225 auto get_skin = [&]() {
226 Range body_ents;
227 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
228 Skinner skin(&mField.get_moab());
229 Range skin_ents;
230 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
231 return skin_ents;
232 };
233
234 auto filter_blocks = [&](auto skin) {
235 bool is_contact_block = false;
236 Range contact_range;
237 for (auto m :
238 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
239
240 (boost::format("%s(.*)") % "CONTACT").str()
241
242 ))
243
244 ) {
245 is_contact_block =
246 true; ///< blocs interation is collective, so that is set irrespective
247 ///< if there are entities in given rank or not in the block
248 MOFEM_LOG("CONTACT", Sev::inform)
249 << "Find contact block set: " << m->getName();
250 auto meshset = m->getMeshset();
251 Range contact_meshset_range;
252 CHKERR mField.get_moab().get_entities_by_dimension(
253 meshset, SPACE_DIM - 1, contact_meshset_range, true);
254
255 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
256 contact_meshset_range);
257 contact_range.merge(contact_meshset_range);
258 }
259 if (is_contact_block) {
260 MOFEM_LOG("SYNC", Sev::inform)
261 << "Nb entities in contact surface: " << contact_range.size();
263 skin = intersect(skin, contact_range);
264 }
265 return skin;
266 };
267
268 auto filter_true_skin = [&](auto skin) {
269 Range boundary_ents;
270 ParallelComm *pcomm =
271 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
272 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
273 PSTATUS_NOT, -1, &boundary_ents);
274 return boundary_ents;
275 };
276
277 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
278 CHKERR simple->setFieldOrder("SIGMA", 0);
279 CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
280
281 if (contact_order > order) {
282 Range ho_ents;
283
284 CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
285 moab::Interface::UNION);
286
287 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
288 CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
289 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
290 }
291
292 CHKERR simple->setUp();
293
294 auto project_ho_geometry = [&]() {
295 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
296 return mField.loop_dofs("GEOMETRY", ent_method);
297 };
298
299 PetscBool project_geometry = PETSC_TRUE;
300 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
301 &project_geometry, PETSC_NULLPTR);
302 if (project_geometry) {
303 CHKERR project_ho_geometry();
304 }
305
307} //! [Set up problem]
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
int sigma_order
Definition contact.cpp:80
int contact_order
Definition contact.cpp:79
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition contact.cpp:67
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 MYPCOMM_INDEX
default communicator number PCOMM
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.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
FTensor::Index< 'm', 3 > m
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Interface for managing meshsets containing materials and boundary conditions.
Projection of edge entities with one mid-node on hierarchical basis.

◆ tsSolve()

MoFEMErrorCode Contact::tsSolve ( )
private

Definition at line 718 of file contact.cpp.

718 {
720
723 ISManager *is_manager = mField.getInterface<ISManager>();
724
725 auto set_section_monitor = [&](auto solver) {
727 SNES snes;
728 CHKERR TSGetSNES(solver, &snes);
729 PetscViewerAndFormat *vf;
730 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
731 PETSC_VIEWER_DEFAULT, &vf);
732 CHKERR SNESMonitorSet(
733 snes,
734 (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal,
735 void *))SNESMonitorFields,
736 vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
738 };
739
740 auto scatter_create = [&](auto D, auto coeff) {
742 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
743 ROW, "U", coeff, coeff, is);
744 int loc_size;
745 CHKERR ISGetLocalSize(is, &loc_size);
746 Vec v;
747 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
748 VecScatter scatter;
749 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
750 return std::make_tuple(SmartPetscObj<Vec>(v),
752 };
753
754 auto set_time_monitor = [&](auto dm, auto solver) {
756 monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
757 boost::shared_ptr<ForcesAndSourcesCore> null;
758 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
759 monitorPtr, null, null);
761 };
762
763 auto set_essential_bc = [&]() {
765 // This is low level pushing finite elements (pipelines) to solver
766 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
767 auto pre_proc_ptr = boost::make_shared<FEMethod>();
768 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
769 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
770
771 // Add boundary condition scaling
772 auto time_scale = boost::make_shared<TimeScale>();
773
774 auto get_bc_hook_rhs = [&]() {
776 {time_scale}, false);
777 return hook;
778 };
779 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
780
781 auto get_post_proc_hook_rhs = [&]() {
783 mField, post_proc_rhs_ptr, 1.);
784 };
785 auto get_post_proc_hook_lhs = [&]() {
787 mField, post_proc_lhs_ptr, 1.);
788 };
789 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
790
791 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
792 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
793 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
794 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
795 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
797 };
798
799 // Set up Schur preconditioner
800 auto set_schur_pc = [&](auto solver) {
801 boost::shared_ptr<SetUpSchur> schur_ptr;
802 if (AT == AssemblyType::BLOCK_SCHUR) {
803 // Set up Schur preconditioner
805 CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
806 }
807 return schur_ptr;
808 };
809
810 auto dm = simple->getDM();
811 auto D = createDMVector(dm);
812
814
815 uXScatter = scatter_create(D, 0);
816 uYScatter = scatter_create(D, 1);
817 if (SPACE_DIM == 3)
818 uZScatter = scatter_create(D, 2);
819
820 // Add extra finite elements to SNES solver pipelines to resolve essential
821 // boundary conditions
822 CHKERR set_essential_bc();
823
824 if (is_quasi_static == PETSC_TRUE) {
825 auto solver = pip_mng->createTSIM();
826 CHKERR TSSetFromOptions(solver);
827
828 auto B = createDMMatrix(dm);
829 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
830 auto schur_pc_ptr = set_schur_pc(solver);
831
832 auto D = createDMVector(dm);
833 CHKERR set_section_monitor(solver);
834 CHKERR set_time_monitor(dm, solver);
835 CHKERR TSSetSolution(solver, D);
836 CHKERR TSSetUp(solver);
837 CHKERR TSSolve(solver, NULL);
838 } else {
839 auto solver = pip_mng->createTSIM2();
840 CHKERR TSSetFromOptions(solver);
841
842 auto B = createDMMatrix(dm);
843 CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
844 auto schur_pc_ptr = set_schur_pc(solver);
845
846 auto D = createDMVector(dm);
847 auto DD = vectorDuplicate(D);
848 CHKERR set_section_monitor(solver);
849 CHKERR set_time_monitor(dm, solver);
850 CHKERR TS2SetSolution(solver, D, DD);
851 CHKERR TSSetUp(solver);
852 CHKERR TSSolve(solver, NULL);
853 }
854
856}
@ ROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1191
double D
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1276
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static auto createTotalTraction(MoFEM::Interface &m_field)
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition contact.cpp:138
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition contact.cpp:137
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition contact.cpp:139
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
intrusive_ptr for managing petsc objects
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)

Member Data Documentation

◆ mField

MoFEM::Interface& Contact::mField
private

Definition at line 128 of file contact.cpp.

◆ mfrontInterfacePtr

boost::shared_ptr<MFrontInterface> Contact::mfrontInterfacePtr
private

Definition at line 141 of file contact.cpp.

◆ monitorPtr

boost::shared_ptr<Monitor> Contact::monitorPtr
private

Definition at line 142 of file contact.cpp.

◆ uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Contact::uXScatter
private

Definition at line 137 of file contact.cpp.

◆ uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Contact::uYScatter
private

Definition at line 138 of file contact.cpp.

◆ uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Contact::uZScatter
private

Definition at line 139 of file contact.cpp.


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