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

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

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

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

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

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:466
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
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:1265
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:595
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
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
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:413
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:50
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:66
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:702
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:578
MoFEM::SmartPetscObj< IS >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:768
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:839
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