v0.14.0
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 125 of file contact.cpp.

Constructor & Destructor Documentation

◆ Contact()

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

Definition at line 127 of file contact.cpp.

127 : mField(m_field) {}

Member Function Documentation

◆ bC()

MoFEMErrorCode Contact::bC ( )
private

[Create common data]

[Boundary condition]

Definition at line 435 of file contact.cpp.

435  {
437  auto bc_mng = mField.getInterface<BcManager>();
438  auto simple = mField.getInterface<Simple>();
439 
440  for (auto f : {"U", "SIGMA"}) {
441  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
442  "REMOVE_X", f, 0, 0);
443  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
444  "REMOVE_Y", f, 1, 1);
445  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
446  "REMOVE_Z", f, 2, 2);
447  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
448  "REMOVE_ALL", f, 0, 3);
449  }
450 
451  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
452  "SIGMA", 0, 0, false, true);
453  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
454  "SIGMA", 1, 1, false, true);
455  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
456  "SIGMA", 2, 2, false, true);
457  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
458  "SIGMA", 0, 3, false, true);
459  CHKERR bc_mng->removeBlockDOFsOnEntities(
460  simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
461 
462  // Note remove has to be always before push. Then node marking will be
463  // corrupted.
464  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
465  simple->getProblemName(), "U");
466 
468 }

◆ checkResults()

MoFEMErrorCode Contact::checkResults ( )
private

[Solve]

[Check]

Definition at line 838 of file contact.cpp.

838  {
840  if (atom_test && !mField.get_comm_rank()) {
841  const double *t_ptr;
842  CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
843  double hertz_force;
844  double fem_force;
845  double analytical_active_area = 1.0;
846  double norm = 1e-5;
847  double tol_force = 1e-3;
848  double tol_norm = 7.5; // change when analytical functions are updated
849  double tol_area = 3e-2;
850  double fem_active_area = t_ptr[3];
851 
852  switch (atom_test) {
853  case 1: // plane stress
854  hertz_force = 3.927;
855  fem_force = t_ptr[1];
856  break;
857 
858  case 2: // plane strain
859  hertz_force = 4.675;
860  fem_force = t_ptr[1];
861  norm = monitorPtr->getErrorNorm(1);
862  break;
863 
864  case 3: // Hertz 3D
865  hertz_force = 3.968;
866  tol_force = 2e-3;
867  fem_force = t_ptr[2];
868  analytical_active_area = M_PI / 4;
869  tol_area = 0.2;
870  break;
871 
872  case 4: // axisymmetric
873  tol_force = 5e-3;
874  tol_area = 0.2;
875  // analytical_active_area = M_PI;
876 
877  case 5: // axisymmetric
878  hertz_force = 15.873;
879  tol_force = 5e-3;
880  fem_force = t_ptr[1];
881  norm = monitorPtr->getErrorNorm(1);
882  analytical_active_area = M_PI;
883  break;
884 
885  case 6: // wavy 2d
886  hertz_force = 0.374;
887  fem_force = t_ptr[1];
888  break;
889 
890  case 7: // wavy 3d
891  hertz_force = 0.5289;
892  fem_force = t_ptr[2];
893  break;
894 
895  default:
896  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
897  "atom test %d does not exist", atom_test);
898  }
899  if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
900  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
901  "atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
902  atom_test, fem_force, hertz_force);
903  }
904  if (norm > tol_norm) {
905  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
906  "atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
907  atom_test, norm, tol_norm);
908  }
909  if (fabs(fem_active_area - analytical_active_area) > tol_area) {
910  SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
911  "atom test %d failed: AREA computed %3.4e but should be %3.4e",
912  atom_test, fem_active_area, analytical_active_area);
913  }
914  CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
915  }
916 
918 
920 }

◆ createCommonData()

MoFEMErrorCode Contact::createCommonData ( )
private

[Set up problem]

[Create common data]

Definition at line 314 of file contact.cpp.

314  {
316 
317  PetscBool use_mfront = PETSC_FALSE;
318  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_mfront", &use_mfront,
319  PETSC_NULL);
320  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_axisymmetric",
321  &is_axisymmetric, PETSC_NULL);
322  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test,
323  PETSC_NULL);
324 
325  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
326  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus", &young_modulus,
327  PETSC_NULL);
328  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio", &poisson_ratio,
329  PETSC_NULL);
330  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
331  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn_contact, PETSC_NULL);
332  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
333  &spring_stiffness, PETSC_NULL);
334  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-vis_spring_stiffness",
335  &vis_spring_stiffness, PETSC_NULL);
336  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping", &alpha_damping,
337  PETSC_NULL);
338 
339  if (!mfrontInterface) {
340  MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
341  MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
342  } else {
343  MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
344  }
345 
346  MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
347  MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
348  MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
349  MOFEM_LOG("CONTACT", Sev::inform)
350  << "Vis spring_stiffness " << vis_spring_stiffness;
351 
352  MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
353 
354  PetscBool use_scale = PETSC_FALSE;
355  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_scale", &use_scale,
356  PETSC_NULL);
357  if (use_scale) {
358  scale /= young_modulus;
359  }
360 
361  MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
362 
363  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_quasi_static",
364  &is_quasi_static, PETSC_NULL);
365  MOFEM_LOG("CONTACT", Sev::inform)
366  << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
367 
368 #ifdef PYTHON_SDF
369  char sdf_file_name[255];
370  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
371  sdf_file_name, 255, PETSC_NULL);
372 
373  sdfPythonPtr = boost::make_shared<SDFPython>();
374  CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
375  sdfPythonWeakPtr = sdfPythonPtr;
376 #endif
377 
378  if (is_axisymmetric) {
379  if (SPACE_DIM == 3) {
380  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
381  "Use executable contact_2d with axisymmetric model");
382  } else {
383  if (!use_mfront) {
384  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
385  "Axisymmetric model is only available with MFront (set "
386  "use_mfront to 1)");
387  } else {
388  MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
389  }
390  }
391  } else {
392  if (SPACE_DIM == 2) {
393  MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
394  }
395  }
396 
397  if (use_mfront) {
398 #ifndef WITH_MODULE_MFRONT_INTERFACE
399  SETERRQ(
400  PETSC_COMM_SELF, MOFEM_NOT_FOUND,
401  "MFrontInterface module was not found while use_mfront was set to 1");
402 #else
403  if (SPACE_DIM == 3) {
405  boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
406  mField, "U", "GEOMETRY", true, is_quasi_static);
407  } else if (SPACE_DIM == 2) {
408  if (is_axisymmetric) {
410  boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
411  mField, "U", "GEOMETRY", true, is_quasi_static);
412  } else {
413  mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
414  mField, "U", "GEOMETRY", true, is_quasi_static);
415  }
416  }
417 #endif
418  CHKERR mfrontInterface->getCommandLineParameters();
419  }
420 
422  auto dm = simple->getDM();
423  monitorPtr =
424  boost::make_shared<Monitor>(dm, scale, mfrontInterface, is_axisymmetric);
425 
426  if (use_mfront) {
427  mfrontInterface->setMonitorPtr(monitorPtr);
428  }
429 
431 }

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

472  {
474  auto simple = mField.getInterface<Simple>();
475  auto *pip_mng = mField.getInterface<PipelineManager>();
476  auto time_scale = boost::make_shared<ScaledTimeScale>();
477  auto body_force_time_scale =
478  boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
479 
480  auto integration_rule_vol = [](int, int, int approx_order) {
481  return 2 * approx_order + geom_order - 1;
482  };
483  auto integration_rule_boundary = [](int, int, int approx_order) {
484  return 2 * approx_order + geom_order - 1;
485  };
486 
487  auto add_domain_base_ops = [&](auto &pip) {
490  "GEOMETRY");
492  };
493 
494  auto hencky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
495  hencky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
496  hencky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
497 
498  auto add_domain_ops_lhs = [&](auto &pip) {
500 
501  //! [Only used for dynamics]
504  //! [Only used for dynamics]
505  if (is_quasi_static == PETSC_FALSE) {
506 
507  auto *pip_mng = mField.getInterface<PipelineManager>();
508  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
509 
510  auto get_inertia_and_mass_damping =
511  [this, fe_domain_lhs](const double, const double, const double) {
512  return (rho * scale) * fe_domain_lhs->ts_aa +
513  (alpha_damping * scale) * fe_domain_lhs->ts_a;
514  };
515  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
516  } else {
517 
518  auto *pip_mng = mField.getInterface<PipelineManager>();
519  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
520 
521  auto get_inertia_and_mass_damping =
522  [this, fe_domain_lhs](const double, const double, const double) {
523  return (alpha_damping * scale) * fe_domain_lhs->ts_a;
524  };
525  pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
526  }
527 
528  if (!mfrontInterface) {
529  CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
530  mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
531  } else {
532  CHKERR mfrontInterface->opFactoryDomainLhs(pip);
533  }
534 
536  };
537 
538  auto add_domain_ops_rhs = [&](auto &pip) {
540 
542  pip, mField, "U", {body_force_time_scale}, Sev::inform);
543 
544  //! [Only used for dynamics]
546  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
547  //! [Only used for dynamics]
548 
549  // only in case of dynamics
550  if (is_quasi_static == PETSC_FALSE) {
551  auto mat_acceleration = boost::make_shared<MatrixDouble>();
553  "U", mat_acceleration));
554  pip.push_back(
555  new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
556  return rho * scale;
557  }));
558  }
559 
560  // only in case of viscosity
561  if (alpha_damping > 0) {
562  auto mat_velocity = boost::make_shared<MatrixDouble>();
563  pip.push_back(
564  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
565  pip.push_back(
566  new OpInertiaForce("U", mat_velocity, [](double, double, double) {
567  return alpha_damping * scale;
568  }));
569  }
570 
571  if (!mfrontInterface) {
572  CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
573  mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
574  } else {
575  CHKERR mfrontInterface->opFactoryDomainRhs(pip);
576  }
577 
578  CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
579  pip, "SIGMA", "U", is_axisymmetric);
580 
582  };
583 
584  auto add_boundary_base_ops = [&](auto &pip) {
587  "GEOMETRY");
588  // We have to integrate on curved face geometry, thus integration weight
589  // have to adjusted.
590  pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
592  };
593 
594  auto add_boundary_ops_lhs = [&](auto &pip) {
596 
597  //! [Operators used for contact]
600  //! [Operators used for contact]
601 
602  // Add Natural BCs to LHS
604  pip, mField, "U", Sev::inform);
605 
606  if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
607 
608  auto *pip_mng = mField.getInterface<PipelineManager>();
609  auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
610 
611  pip.push_back(new OpSpringLhs(
612  "U", "U",
613 
614  [this, fe_boundary_lhs](double, double, double) {
615  return spring_stiffness * scale +
616  (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
617  }
618 
619  ));
620  }
621 
622  CHKERR
623  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
624  pip, "SIGMA", "U", is_axisymmetric);
626  DomainEle>(
627  mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
628  integration_rule_vol, is_axisymmetric);
629 
631  };
632 
633  auto add_boundary_ops_rhs = [&](auto &pip) {
635 
636  //! [Operators used for contact]
638  AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
639  //! [Operators used for contact]
640 
641  // Add Natural BCs to RHS
643  pip, mField, "U", {time_scale}, Sev::inform);
644 
645  if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
646  auto u_disp = boost::make_shared<MatrixDouble>();
647  auto dot_u_disp = boost::make_shared<MatrixDouble>();
648  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
649  pip.push_back(
650  new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
651  pip.push_back(
652  new OpSpringRhs("U", u_disp, [this](double, double, double) {
653  return spring_stiffness * scale;
654  }));
655  pip.push_back(
656  new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
657  return vis_spring_stiffness * scale;
658  }));
659  }
660 
661  CHKERR
662  ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
663  pip, "SIGMA", "U", is_axisymmetric);
664 
666  };
667 
668  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
669  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
670  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
671  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
672 
673  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
674  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
675  CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
676  CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
677 
678  if (mfrontInterface) {
679  CHKERR mfrontInterface->setUpdateElementVariablesOperators();
680  }
681 
682  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
683  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
684  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
685  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
686 
688 }

◆ runProblem()

MoFEMErrorCode Contact::runProblem ( )

[Run problem]

Definition at line 161 of file contact.cpp.

161  {
165  CHKERR bC();
166  CHKERR OPs();
167  CHKERR tsSolve();
170 }

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

174  {
177 
178  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
179  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-contact_order", &contact_order,
180  PETSC_NULL);
181  sigma_order = std::max(order, contact_order) - 1;
182  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-sigma_order", &sigma_order,
183  PETSC_NULL);
184  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
185  PETSC_NULL);
186 
187  MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
188  MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
189  MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
190  MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
191 
192  // Select base
193  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
194  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
195  PetscInt choice_base_value = AINSWORTH;
196  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
197  LASBASETOPT, &choice_base_value, PETSC_NULL);
198 
200  switch (choice_base_value) {
201  case AINSWORTH:
203  MOFEM_LOG("CONTACT", Sev::inform)
204  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
205  break;
206  case DEMKOWICZ:
207  base = DEMKOWICZ_JACOBI_BASE;
208  MOFEM_LOG("CONTACT", Sev::inform)
209  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
210  break;
211  default:
212  base = LASTBASE;
213  break;
214  }
215 
216  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
217  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
218  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
219  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
220  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
221  SPACE_DIM);
222  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
223  SPACE_DIM);
224  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
225 
226  CHKERR simple->setFieldOrder("U", order);
227  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
228 
229  auto get_skin = [&]() {
230  Range body_ents;
231  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
232  Skinner skin(&mField.get_moab());
233  Range skin_ents;
234  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
235  return skin_ents;
236  };
237 
238  auto filter_blocks = [&](auto skin) {
239  bool is_contact_block = false;
240  Range contact_range;
241  for (auto m :
243 
244  (boost::format("%s(.*)") % "CONTACT").str()
245 
246  ))
247 
248  ) {
249  is_contact_block =
250  true; ///< blocs interation is collective, so that is set irrespective
251  ///< if there are entities in given rank or not in the block
252  MOFEM_LOG("CONTACT", Sev::inform)
253  << "Find contact block set: " << m->getName();
254  auto meshset = m->getMeshset();
255  Range contact_meshset_range;
256  CHKERR mField.get_moab().get_entities_by_dimension(
257  meshset, SPACE_DIM - 1, contact_meshset_range, true);
258 
259  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
260  contact_meshset_range);
261  contact_range.merge(contact_meshset_range);
262  }
263  if (is_contact_block) {
264  MOFEM_LOG("SYNC", Sev::inform)
265  << "Nb entities in contact surface: " << contact_range.size();
267  skin = intersect(skin, contact_range);
268  }
269  return skin;
270  };
271 
272  auto filter_true_skin = [&](auto skin) {
273  Range boundary_ents;
274  ParallelComm *pcomm =
275  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
276  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
277  PSTATUS_NOT, -1, &boundary_ents);
278  return boundary_ents;
279  };
280 
281  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
282  CHKERR simple->setFieldOrder("SIGMA", 0);
283  CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
284 
285  if (contact_order > order) {
286  Range ho_ents;
287 
288  CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
289  moab::Interface::UNION);
290 
291  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
292  CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
293  CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
294  }
295 
296  CHKERR simple->setUp();
297 
298  auto project_ho_geometry = [&]() {
299  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
300  return mField.loop_dofs("GEOMETRY", ent_method);
301  };
302 
303  PetscBool project_geometry = PETSC_TRUE;
304  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
305  &project_geometry, PETSC_NULL);
306  if (project_geometry) {
307  CHKERR project_ho_geometry();
308  }
309 
311 } //! [Set up problem]

◆ tsSolve()

MoFEMErrorCode Contact::tsSolve ( )
private

Definition at line 701 of file contact.cpp.

701  {
703 
706  ISManager *is_manager = mField.getInterface<ISManager>();
707 
708  auto set_section_monitor = [&](auto solver) {
710  SNES snes;
711  CHKERR TSGetSNES(solver, &snes);
712  PetscViewerAndFormat *vf;
713  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
714  PETSC_VIEWER_DEFAULT, &vf);
715  CHKERR SNESMonitorSet(
716  snes,
717  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
718  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
720  };
721 
722  auto scatter_create = [&](auto D, auto coeff) {
724  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
725  ROW, "U", coeff, coeff, is);
726  int loc_size;
727  CHKERR ISGetLocalSize(is, &loc_size);
728  Vec v;
729  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
730  VecScatter scatter;
731  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
732  return std::make_tuple(SmartPetscObj<Vec>(v),
733  SmartPetscObj<VecScatter>(scatter));
734  };
735 
736  auto set_time_monitor = [&](auto dm, auto solver) {
738  monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
739  boost::shared_ptr<ForcesAndSourcesCore> null;
740  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
741  monitorPtr, null, null);
743  };
744 
745  auto set_essential_bc = [&]() {
747  // This is low level pushing finite elements (pipelines) to solver
748  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
749  auto pre_proc_ptr = boost::make_shared<FEMethod>();
750  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
751  auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
752 
753  // Add boundary condition scaling
754  auto time_scale = boost::make_shared<TimeScale>();
755 
756  auto get_bc_hook_rhs = [&]() {
758  {time_scale}, false);
759  return hook;
760  };
761  pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
762 
763  auto get_post_proc_hook_rhs = [&]() {
765  mField, post_proc_rhs_ptr, 1.);
766  };
767  auto get_post_proc_hook_lhs = [&]() {
769  mField, post_proc_lhs_ptr, 1.);
770  };
771  post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
772 
773  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
774  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
775  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
776  post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
777  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
779  };
780 
781  // Set up Schur preconditioner
782  auto set_schur_pc = [&](auto solver) {
783  boost::shared_ptr<SetUpSchur> schur_ptr;
784  if (AT == AssemblyType::BLOCK_SCHUR) {
785  // Set up Schur preconditioner
787  CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
788  }
789  return schur_ptr;
790  };
791 
792  auto dm = simple->getDM();
793  auto D = createDMVector(dm);
794 
796 
797  uXScatter = scatter_create(D, 0);
798  uYScatter = scatter_create(D, 1);
799  if (SPACE_DIM == 3)
800  uZScatter = scatter_create(D, 2);
801 
802  // Add extra finite elements to SNES solver pipelines to resolve essential
803  // boundary conditions
804  CHKERR set_essential_bc();
805 
806  if (is_quasi_static == PETSC_TRUE) {
807  auto solver = pip_mng->createTSIM();
808  CHKERR TSSetFromOptions(solver);
809  auto schur_pc_ptr = set_schur_pc(solver);
810 
811  auto D = createDMVector(dm);
812  CHKERR set_section_monitor(solver);
813  CHKERR set_time_monitor(dm, solver);
814  CHKERR TSSetSolution(solver, D);
815  CHKERR TSSetUp(solver);
816  CHKERR TSSolve(solver, NULL);
817  } else {
818  auto solver = pip_mng->createTSIM2();
819  CHKERR TSSetFromOptions(solver);
820  auto schur_pc_ptr = set_schur_pc(solver);
821 
822  auto dm = simple->getDM();
823  auto D = createDMVector(dm);
824  auto DD = vectorDuplicate(D);
825 
826  CHKERR set_section_monitor(solver);
827  CHKERR set_time_monitor(dm, solver);
828  CHKERR TS2SetSolution(solver, D, DD);
829  CHKERR TSSetUp(solver);
830  CHKERR TSSolve(solver, NULL);
831  }
832 
834 }

Member Data Documentation

◆ mField

MoFEM::Interface& Contact::mField
private

Definition at line 132 of file contact.cpp.

◆ mfrontInterface

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

Definition at line 145 of file contact.cpp.

◆ monitorPtr

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

Definition at line 146 of file contact.cpp.

◆ uXScatter

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

Definition at line 141 of file contact.cpp.

◆ uYScatter

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

Definition at line 142 of file contact.cpp.

◆ uZScatter

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

Definition at line 143 of file contact.cpp.


The documentation for this struct was generated from the following file:
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
spring_stiffness
double spring_stiffness
Definition: contact.cpp:87
MoFEM::CoreInterface::loop_dofs
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.
H1
@ H1
continuous field
Definition: definitions.h:85
Contact::monitorPtr
boost::shared_ptr< Monitor > monitorPtr
Definition: contact.cpp:146
Contact::bC
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:435
LASTBASE
@ LASTBASE
Definition: definitions.h:69
order
int order
Definition: contact.cpp:80
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:93
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
filter_true_skin
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Definition: EshelbianPlasticity.cpp:142
Contact::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:141
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
rho
double rho
Definition: contact.cpp:86
Contact::mField
MoFEM::Interface & mField
Definition: contact.cpp:132
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
ContactOps::opFactoryBoundaryToDomainLhs
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)
Definition: ContactOps.hpp:1263
SPACE_DIM
constexpr int SPACE_DIM
Definition: contact.cpp:55
scale
double scale
Definition: contact.cpp:91
ROW
@ ROW
Definition: definitions.h:136
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:62
atom_test
int atom_test
Definition: contact.cpp:97
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
vis_spring_stiffness
double vis_spring_stiffness
Definition: contact.cpp:88
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:29
OpSpringRhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:75
alpha_damping
double alpha_damping
Definition: contact.cpp:89
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
double
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
poisson_ratio
double poisson_ratio
Definition: contact.cpp:85
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
Contact::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:143
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
contact_order
int contact_order
Definition: contact.cpp:81
get_skin
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
Definition: EshelbianPlasticity.cpp:163
Contact::mfrontInterface
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: contact.cpp:145
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
Contact::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:314
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:31
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
Contact::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:472
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
is_quasi_static
PetscBool is_quasi_static
[Operators used for contact]
Definition: contact.cpp:78
Contact::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:174
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
young_modulus
double young_modulus
Definition: contact.cpp:84
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
approx_order
int approx_order
Definition: test_broken_space.cpp:54
MoFEM::DMMoFEMTSSetMonitor
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:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
Contact::tsSolve
MoFEMErrorCode tsSolve()
Definition: contact.cpp:701
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
geom_order
int geom_order
Definition: contact.cpp:83
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:589
MoFEM::SmartPetscObj< IS >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:802
OpSpringLhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:73
sigma_order
int sigma_order
Definition: contact.cpp:82
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
AT
constexpr AssemblyType AT
Definition: contact.cpp:36
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Contact::checkResults
MoFEMErrorCode checkResults()
[Solve]
Definition: contact.cpp:838
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition: contact.cpp:69
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Contact::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:142