v0.14.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] More...
 

Private Member Functions

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

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< GenericElementInterfacemfrontInterface
 
boost::shared_ptr< MonitormonitorPtr
 

Detailed Description

Definition at line 169 of file contact.cpp.

Constructor & Destructor Documentation

◆ Contact()

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

Definition at line 171 of file contact.cpp.

171 : mField(m_field) {
172 mfrontInterface = nullptr;
173 }
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: contact.cpp:191
MoFEM::Interface & mField
Definition: contact.cpp:178

Member Function Documentation

◆ bC()

MoFEMErrorCode Contact::bC ( )
private

[Create common data]

[Boundary condition]

Definition at line 484 of file contact.cpp.

484 {
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}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
PetscBool is_axisymmetric
Definition: contact.cpp:137
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ checkResults()

MoFEMErrorCode Contact::checkResults ( )
private

[Solve]

[Check]

Definition at line 900 of file contact.cpp.

900{ return 0; }

◆ createCommonData()

MoFEMErrorCode Contact::createCommonData ( )
private

[Set up problem]

[Create common data]

Definition at line 414 of file contact.cpp.

414 {
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}
double young_modulus
Definition: contact.cpp:127
double spring_stiffness
Definition: contact.cpp:130
double rho
Definition: contact.cpp:129
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
double poisson_ratio
Definition: contact.cpp:128
double vis_spring_stiffness
Definition: contact.cpp:131
double scale
Definition: contact.cpp:134
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)

◆ 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 526 of file contact.cpp.

526 {
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}
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr AssemblyType AT
Definition: contact.cpp:34
constexpr int SPACE_DIM
Definition: contact.cpp:53
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:118
int geom_order
Definition: contact.cpp:126
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:120
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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)
static constexpr int approx_order
Add operators pushing bases from local to physical configuration.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
PipelineManager interface.

◆ runProblem()

MoFEMErrorCode Contact::runProblem ( )

[Run problem]

Definition at line 207 of file contact.cpp.

207 {
211 CHKERR bC();
212 CHKERR OPs();
213 CHKERR tsSolve();
216}
MoFEMErrorCode tsSolve()
Definition: contact.cpp:760
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
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:484

◆ setupProblem()

MoFEMErrorCode Contact::setupProblem ( )
private

[Run problem]

[Set up problem]

< bloks interation is collectibe, so that is set irrespective if there are enerities in given rank or not in the block

Definition at line 220 of file contact.cpp.

220 {
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]
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
int atom_test
Definition: contact.cpp:139
int order
Definition: contact.cpp:125
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition: contact.cpp:114
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
Definition: definitions.h:215
FTensor::Index< 'm', SPACE_DIM > m
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
boost::shared_ptr< Monitor > monitorPtr
Definition: contact.cpp:192
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 760 of file contact.cpp.

760 {
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}
@ ROW
Definition: definitions.h:123
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
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.
Definition: Exceptions.hpp:56
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
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:188
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:187
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:189
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
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, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1482

Member Data Documentation

◆ mField

MoFEM::Interface& Contact::mField
private

Definition at line 178 of file contact.cpp.

◆ mfrontInterface

boost::shared_ptr<GenericElementInterface> Contact::mfrontInterface
private

Definition at line 191 of file contact.cpp.

◆ monitorPtr

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

Definition at line 192 of file contact.cpp.

◆ uXScatter

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

Definition at line 187 of file contact.cpp.

◆ uYScatter

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

Definition at line 188 of file contact.cpp.

◆ uZScatter

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

Definition at line 189 of file contact.cpp.


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