Collaboration diagram for Contact:

[legend]

Classes
struct	ScaledTimeScale

Public Member Functions
	Contact (MoFEM::Interface &m_field)

MoFEMErrorCode	runProblem ()
	[Run problem]

Private Member Functions
MoFEMErrorCode	setupProblem ()
	[Run problem]

MoFEMErrorCode	createCommonData ()
	[Set up problem]

MoFEMErrorCode	bC ()
	[Create common data]

MoFEMErrorCode	OPs ()
	[Boundary condition]

MoFEMErrorCode	tsSolve ()

MoFEMErrorCode	checkResults ()
	[Solve]

Private Attributes
MoFEM::Interface &	mField

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

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

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

boost::shared_ptr< GenericElementInterface >	mfrontInterface

boost::shared_ptr< Monitor >	monitorPtr

Detailed Description

Definition at line 126 of file contact.cpp.

Constructor & Destructor Documentation

◆ Contact()

Contact::Contact ( MoFEM::Interface & m_field )

inline

Definition at line 128 of file contact.cpp.

128: mField(m_field) {}

Contact::mField

MoFEM::Interface & mField

Definition contact.cpp:133

Member Function Documentation

◆ bC()

MoFEMErrorCode Contact::bC ( )

private

[Create common data]

[Boundary condition]

Definition at line 459 of file contact.cpp.

                           {
  MoFEMFunctionBegin;
  auto bc_mng = mField.getInterface<BcManager>();
  auto simple = mField.getInterface<Simple>();
 
  for (auto f : {"U", "SIGMA"}) {
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_X", f, 0, 0);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_Y", f, 1, 1);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_Z", f, 2, 2);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_ALL", f, 0, 3);
  }
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
                                           "SIGMA", 0, 0, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
                                           "SIGMA", 1, 1, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
                                           "SIGMA", 2, 2, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
                                           "SIGMA", 0, 3, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(
      simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
 
  // Note remove has to be always before push. Then node marking will be
  // corrupted.
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}

◆ checkResults()

MoFEMErrorCode Contact::checkResults ( )

private

[Solve]

[Check]

Definition at line 872 of file contact.cpp.

                                     {
  MoFEMFunctionBegin;
  if (atom_test && !mField.get_comm_rank()) {
    const double *t_ptr;
    CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
    double hertz_force;
    double fem_force;
    double analytical_active_area = 1.0;
    double norm = 1e-5;
    double tol_force = 1e-3;
    double tol_norm = 7.5; // change when analytical functions are updated
    double tol_area = 3e-2;
    double fem_active_area = t_ptr[3];
 
    switch (atom_test) {
    case 1: // plane stress
      hertz_force = 3.927;
      fem_force = t_ptr[1];
      break;
 
    case 2: // plane strain
      hertz_force = 4.675;
      fem_force = t_ptr[1];
      norm = monitorPtr->getErrorNorm(1);
      break;
 
    case 3: // Hertz 3D
      hertz_force = 3.968;
      tol_force = 2e-3;
      fem_force = t_ptr[2];
      analytical_active_area = M_PI / 4;
      tol_area = 0.2;
      break;
 
    case 4: // axisymmetric
      tol_force = 5e-3;
      tol_area = 0.2;
      // analytical_active_area = M_PI;
 
    case 5: // axisymmetric
      hertz_force = 15.873;
      tol_force = 5e-3;
      fem_force = t_ptr[1];
      norm = monitorPtr->getErrorNorm(1);
      analytical_active_area = M_PI;
      break;
 
    case 6: // wavy 2d
      hertz_force = 0.374;
      fem_force = t_ptr[1];
      break;
 
    case 7: // wavy 3d
      hertz_force = 0.5289;
      fem_force = t_ptr[2];
      break;
 
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
               "atom test %d does not exist", atom_test);
    }
    if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
               atom_test, fem_force, hertz_force);
    }
    if (norm > tol_norm) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
               atom_test, norm, tol_norm);
    }
    if (fabs(fem_active_area - analytical_active_area) > tol_area) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "atom test %d failed: AREA computed %3.4e but should be %3.4e",
               atom_test, fem_active_area, analytical_active_area);
    }
    CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
  }
 
  ContactOps::CommonData::totalTraction.reset();
 
  MoFEMFunctionReturn(0);
}

◆ createCommonData()

MoFEMErrorCode Contact::createCommonData ( )

private

[Set up problem]

[Create common data]

Definition at line 326 of file contact.cpp.

                                         {
  MoFEMFunctionBegin;
 
  PetscBool use_mfront = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_mfront", &use_mfront,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_axisymmetric",
                             &is_axisymmetric, PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                            PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus", &young_modulus,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio", &poisson_ratio,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn", &cn_contact, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-spring_stiffness",
                               &spring_stiffness, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-vis_spring_stiffness",
                               &vis_spring_stiffness, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping", &alpha_damping,
                               PETSC_NULLPTR);
 
  if (!mfrontInterface) {
    MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
    MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
  } else {
    MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
  }
 
  MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
  MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
  MOFEM_LOG("CONTACT", Sev::inform) << "Spring stiffness " << spring_stiffness;
  MOFEM_LOG("CONTACT", Sev::inform)
      << "Vis spring_stiffness " << vis_spring_stiffness;
 
  MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
 
  PetscBool use_scale = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_scale", &use_scale,
                             PETSC_NULLPTR);
  if (use_scale) {
    scale /= young_modulus;
  }
 
  MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_quasi_static",
                             &is_quasi_static, PETSC_NULLPTR);
  MOFEM_LOG("CONTACT", Sev::inform)
      << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
 
#ifdef ENABLE_PYTHON_BINDING
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
  char sdf_file_name[255] = "sdf.py";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
                               sdf_file_name, 255, PETSC_NULLPTR);
 
  if (file_exists(sdf_file_name)) {
    MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
    sdfPythonPtr = boost::make_shared<SDFPython>();
    CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
    sdfPythonWeakPtr = sdfPythonPtr;
  } else {
    MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
  }
#endif
 
  if (is_axisymmetric) {
    if (SPACE_DIM == 3) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Use executable contact_2d with axisymmetric model");
    } else {
      if (!use_mfront) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
                "Axisymmetric model is only available with MFront (set "
                "use_mfront to 1)");
      } else {
        MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
      }
    }
  } else {
    if (SPACE_DIM == 2) {
      MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
    }
  }
 
  if (use_mfront) {
#ifndef WITH_MODULE_MFRONT_INTERFACE
    SETERRQ(
        PETSC_COMM_SELF, MOFEM_NOT_FOUND,
        "MFrontInterface module was not found while use_mfront was set to 1");
#else
    if (SPACE_DIM == 3) {
      mfrontInterface =
          boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
              mField, "U", "GEOMETRY", true, is_quasi_static);
    } else if (SPACE_DIM == 2) {
      if (is_axisymmetric) {
        mfrontInterface =
            boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
                mField, "U", "GEOMETRY", true, is_quasi_static);
      } else {
        mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
            mField, "U", "GEOMETRY", true, is_quasi_static);
      }
    }
#endif
    CHKERR mfrontInterface->getCommandLineParameters();
  }
 
  Simple *simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
  monitorPtr =
      boost::make_shared<Monitor>(dm, scale, mfrontInterface, is_axisymmetric, is_large_strain);
 
  if (use_mfront) {
    mfrontInterface->setMonitorPtr(monitorPtr);
  }
 
  MoFEMFunctionReturn(0);
}

◆ OPs()

MoFEMErrorCode Contact::OPs ( )

private

[Boundary condition]

[Push operators to pip]

[Only used for dynamics]

[Operators used for contact]

Definition at line 496 of file contact.cpp.

                            {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto *pip_mng = mField.getInterface<PipelineManager>();
  auto time_scale = boost::make_shared<ScaledTimeScale>();
  auto body_force_time_scale =
      boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
 
  auto integration_rule_vol = [](int, int, int approx_order) {
    return 2 * approx_order + geom_order - 1;
  };
  auto integration_rule_boundary = [](int, int, int approx_order) {
    return 2 * approx_order + geom_order - 1;
  };
 
  auto add_domain_base_ops = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                          "GEOMETRY");
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Only used for dynamics]
    using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
        GAUSS>::OpMass<1, SPACE_DIM>;
    //! [Only used for dynamics]
    if (is_quasi_static == PETSC_FALSE) {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_domain_lhs = pip_mng->getDomainLhsFE();
 
      auto get_inertia_and_mass_damping =
          [this, fe_domain_lhs](const double, const double, const double) {
            return (rho * scale) * fe_domain_lhs->ts_aa +
                   (alpha_damping * scale) * fe_domain_lhs->ts_a;
          };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
    } else {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_domain_lhs = pip_mng->getDomainLhsFE();
 
      auto get_inertia_and_mass_damping =
          [this, fe_domain_lhs](const double, const double, const double) {
            return (alpha_damping * scale) * fe_domain_lhs->ts_a;
          };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
    }
 
    if (!mfrontInterface) {
      if (!is_large_strain) {
        CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
            mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
      } else {
        CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
            mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
      }
    } else {
      CHKERR mfrontInterface->opFactoryDomainLhs(pip);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
        pip, mField, "U", {body_force_time_scale}, Sev::inform);
 
    //! [Only used for dynamics]
    using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
    //! [Only used for dynamics]
 
    // only in case of dynamics
    if (is_quasi_static == PETSC_FALSE) {
      auto mat_acceleration = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
          "U", mat_acceleration));
      pip.push_back(
          new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
            return rho * scale;
          }));
    }
 
    // only in case of viscosity
    if (alpha_damping > 0) {
      auto mat_velocity = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
      pip.push_back(
          new OpInertiaForce("U", mat_velocity, [](double, double, double) {
            return alpha_damping * scale;
          }));
    }
 
    if (!mfrontInterface) {
      if (!is_large_strain) {
        CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
            mField, pip, "U", "MAT_ELASTIC", Sev::inform, 1, scale);
      } else {
        CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
            mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
      }
    } else {
      CHKERR mfrontInterface->opFactoryDomainRhs(pip);
    }
 
    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_base_ops = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
                                                              "GEOMETRY");
    // We have to integrate on curved face geometry, thus integration weight
    // have to adjusted.
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Operators used for contact]
    using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
        AT>::BiLinearForm<IT>::OpMass<1, SPACE_DIM>;
    //! [Operators used for contact]
 
    // Add Natural BCs to LHS
    CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
        pip, mField, "U", Sev::inform);
 
    if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
 
      pip.push_back(new OpSpringLhs(
          "U", "U",
 
          [this, fe_boundary_lhs](double, double, double) {
            return spring_stiffness * scale +
                   (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
          }
 
          ));
    }
 
    CHKERR
    ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
    CHKERR ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, GAUSS,
                                                    DomainEle>(
        mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
        integration_rule_vol, is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Operators used for contact]
    using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
        AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
    //! [Operators used for contact]
 
    // Add Natural BCs to RHS
    CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
        pip, mField, "U", {time_scale}, Sev::inform);
 
    if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
      auto u_disp = boost::make_shared<MatrixDouble>();
      auto dot_u_disp = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
      pip.push_back(
          new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
      pip.push_back(
          new OpSpringRhs("U", u_disp, [this](double, double, double) {
            return spring_stiffness * scale;
          }));
      pip.push_back(
          new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
            return vis_spring_stiffness * scale;
          }));
    }
 
    CHKERR
    ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
 
  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
  CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
 
  if (mfrontInterface) {
    CHKERR mfrontInterface->setUpdateElementVariablesOperators();
  }
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
 
  MoFEMFunctionReturn(0);
}

◆ runProblem()

MoFEMErrorCode Contact::runProblem ( )

[Run problem]

Definition at line 162 of file contact.cpp.

                                   {
  MoFEMFunctionBegin;
  CHKERR setupProblem();
  CHKERR createCommonData();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

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

                                     {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strain", &is_large_strain,
                             PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain", &is_plane_strain,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-contact_order", &contact_order,
                            PETSC_NULLPTR);
  sigma_order = std::max(order, contact_order) - 1;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-sigma_order", &sigma_order,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
                            PETSC_NULLPTR);
  if (!is_large_strain) {
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Selected material model: Hooke (small strain)";
  } else {
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Selected deafult material model: Hencky (finite strain)";
  }
  MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
  MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
  MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
  MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
 
  // Select base
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                SPACE_DIM);
  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                  SPACE_DIM);
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_blocks = [&](auto skin) {
    bool is_contact_block = false;
    Range contact_range;
    for (auto m :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % "CONTACT").str()
 
                 ))
 
    ) {
      is_contact_block =
          true; ///< blocs interation is collective, so that is set irrespective
                ///< if there are entities in given rank or not in the block
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Find contact block set:  " << m->getName();
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR mField.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    if (is_contact_block) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Nb entities in contact surface: " << contact_range.size();
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
      skin = intersect(skin, contact_range);
    }
    return skin;
  };
 
  auto filter_true_skin = [&](auto skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
  CHKERR simple->setFieldOrder("SIGMA", 0);
  CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
 
  if (contact_order > order) {
    Range ho_ents;
 
    CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
                                             moab::Interface::UNION);
 
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
    CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
  }
 
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
 
  PetscBool project_geometry = PETSC_TRUE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
                             &project_geometry, PETSC_NULLPTR);
  if (project_geometry) {
    CHKERR project_ho_geometry();
  }
 
  MoFEMFunctionReturn(0);
} //! [Set up problem]

◆ tsSolve()

MoFEMErrorCode Contact::tsSolve ( )

private

Definition at line 731 of file contact.cpp.

                                {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    PetscViewerAndFormat *vf;
    CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                      PETSC_VIEWER_DEFAULT, &vf);
    CHKERR SNESMonitorSet(
        snes,
        (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
        vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
    MoFEMFunctionReturn(0);
  };
 
  auto scatter_create = [&](auto D, auto coeff) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitorPtr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    // Add boundary condition scaling
    auto time_scale = boost::make_shared<TimeScale>();
 
    auto get_bc_hook_rhs = [&]() {
      EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
                                                     {time_scale}, false);
      return hook;
    };
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
 
    auto get_post_proc_hook_rhs = [&]() {
      return EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.);
    };
    auto get_post_proc_hook_lhs = [&]() {
      return EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
 
    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
    MoFEMFunctionReturn(0);
  };
 
  // Set up Schur preconditioner
  auto set_schur_pc = [&](auto solver) {
    boost::shared_ptr<SetUpSchur> schur_ptr;
    if (AT == AssemblyType::BLOCK_SCHUR) {
      // Set up Schur preconditioner
      schur_ptr = SetUpSchur::createSetUpSchur(mField);
      CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
    }
    return schur_ptr;
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
 
  ContactOps::CommonData::createTotalTraction(mField);
 
  uXScatter = scatter_create(D, 0);
  uYScatter = scatter_create(D, 1);
  if (SPACE_DIM == 3)
    uZScatter = scatter_create(D, 2);
 
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR set_essential_bc();
 
  if (is_quasi_static == PETSC_TRUE) {
    auto solver = pip_mng->createTSIM();
    CHKERR TSSetFromOptions(solver);
 
    auto B = createDMMatrix(dm);
    CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
    auto schur_pc_ptr = set_schur_pc(solver);
 
    auto D = createDMVector(dm);
    CHKERR set_section_monitor(solver);
    CHKERR set_time_monitor(dm, solver);
    CHKERR TSSetSolution(solver, D);
    CHKERR TSSetUp(solver);
    CHKERR TSSolve(solver, NULL);
  } else {
    auto solver = pip_mng->createTSIM2();
    CHKERR TSSetFromOptions(solver);
 
    auto B = createDMMatrix(dm);
    CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
    auto schur_pc_ptr = set_schur_pc(solver);
 
    auto D = createDMVector(dm);
    auto DD = vectorDuplicate(D);
    CHKERR set_section_monitor(solver);
    CHKERR set_time_monitor(dm, solver);
    CHKERR TS2SetSolution(solver, D, DD);
    CHKERR TSSetUp(solver);
    CHKERR TSSolve(solver, NULL);
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

◆ mField

MoFEM::Interface& Contact::mField

private

Definition at line 133 of file contact.cpp.

◆ mfrontInterface

boost::shared_ptr<GenericElementInterface> Contact::mfrontInterface

private

Definition at line 146 of file contact.cpp.

◆ monitorPtr

boost::shared_ptr<Monitor> Contact::monitorPtr

private

Definition at line 147 of file contact.cpp.

◆ uXScatter

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

private

Definition at line 142 of file contact.cpp.

◆ uYScatter

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

private

Definition at line 143 of file contact.cpp.

◆ uZScatter

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

private

Definition at line 144 of file contact.cpp.

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

tutorials/adv-1/contact.cpp

Classes

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Contact()

Member Function Documentation

◆ bC()

◆ checkResults()

◆ createCommonData()

◆ OPs()

◆ runProblem()

◆ setupProblem()

◆ tsSolve()

Member Data Documentation

◆ mField

◆ mfrontInterface

◆ monitorPtr

◆ uXScatter

◆ uYScatter

◆ uZScatter