mortar_contact.cpp

Implementation of mortar contact between surfaces with non-matching meshes

/** \file mortar_contact.cpp
 * \example mortar_contact.cpp
 *
 * Implementation of mortar contact between surfaces with non-matching meshes
 *
 **/
 
/* MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
 */
 
#include <Mortar.hpp>
#include <MortarCtestFunctions.hpp>
 
// #include <boost/program_options.hpp>
// namespace po = boost::program_options;
#include <ElasticMaterials.hpp>
#include <NonlinearElasticElementInterface.hpp>
 
static char help[] = "\n";
 
double MortarContactProblem::LoadScale::lAmbda = 1;
int main(int argc, char *argv[]) {
 
  const string default_options = "-ksp_type fgmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_type mumps \n"
                                 "-mat_mumps_icntl_13 1 \n"
                                 "-mat_mumps_icntl_14 800 \n"
                                 "-mat_mumps_icntl_20 0 \n"
                                 "-mat_mumps_icntl_24 1 \n"
                                 "-snes_type newtonls \n"
                                 "-snes_linesearch_type basic \n"
                                 "-snes_divergence_tolerance 0 \n"
                                 "-snes_max_it 50 \n"
                                 "-snes_atol 1e-8 \n"
                                 "-snes_rtol 1e-10 \n"
                                 "-snes_monitor \n"
                                 "-ksp_monitor \n"
                                 "-snes_converged_reason \n"
                                 "-order 1 \n"
                                 "-order_lambda 1 \n"
                                 "-order_contact 2 \n"
                                 "-ho_levels_num 1 \n"
                                 "-step_num 1 \n"
                                 "-cn_value 1. \n"
                                 "-alm_flag 0 \n";
 
  string param_file = "param_file.petsc";
  if (!static_cast<bool>(ifstream(param_file))) {
    std::ofstream file(param_file.c_str(), std::ios::ate);
    if (file.is_open()) {
      file << default_options;
      file.close();
    }
  }
 
  // Initialize MoFEM
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
 
  // Create mesh database
  moab::Core mb_instance;              // create database
  moab::Interface &moab = mb_instance; // create interface to database
 
  try {
    PetscBool flg_file;
 
    char mesh_file_name[255];
    PetscInt order = 1;
    PetscInt nb_sub_steps = 1;
    PetscBool is_partitioned = PETSC_FALSE;
    PetscInt test_num = 0;
    PetscBool wave_surf_flag = PETSC_FALSE;
    PetscInt wave_dim = 2;
    PetscBool is_displacement_field = PETSC_FALSE;
    PetscInt wave_surf_block_id = 1;
    PetscReal wave_length = 1.0;
    PetscReal wave_ampl = 0.01;
    PetscReal mesh_height = 1.0;
 
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
 
    CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_file);
 
    CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
                           "", 1, &order, PETSC_NULL);
 
    CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
                           &nb_sub_steps, PETSC_NULL);
 
    CHKERR PetscOptionsBool("-is_partitioned",
                            "set if mesh is partitioned (this result that each "
                            "process keeps only part of the mesh",
                            "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
 
    CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
                           PETSC_NULL);
 
    CHKERR PetscOptionsBool("-is_displacement_field", "use displacement field",
                            "", PETSC_FALSE, &is_displacement_field,
                            PETSC_NULL);
    CHKERR PetscOptionsBool("-wave_surf",
                            "if set true, make one of the surfaces wavy", "",
                            PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
    CHKERR PetscOptionsInt("-wave_surf_block_id",
                           "make wavy surface of the block with this id", "",
                           wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
    CHKERR PetscOptionsInt("-wave_dim", "dimension (2 or 3)", "", wave_dim,
                           &wave_dim, PETSC_NULL);
    CHKERR PetscOptionsReal("-wave_length", "profile wavelength", "",
                            wave_length, &wave_length, PETSC_NULL);
    CHKERR PetscOptionsReal("-wave_ampl", "profile amplitude", "", wave_ampl,
                            &wave_ampl, PETSC_NULL);
    CHKERR PetscOptionsReal("-mesh_height", "vertical dimension of the mesh ",
                            "", mesh_height, &mesh_height, PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
 
    // Check if mesh file was provided
    if (flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
    }
    // FIXME: really???
    if (is_partitioned == PETSC_TRUE)
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Partitioned mesh is not supported");
 
    // CHKERR moab.load_file(mesh_file_name, 0, option);
 
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    auto simple = m_field.getInterface<Simple>();
    simple->loadFile("", mesh_file_name);
 
    MeshsetsManager *mmanager_ptr;
    CHKERR m_field.getInterface(mmanager_ptr);
    CHKERR mmanager_ptr->printDisplacementSet();
    CHKERR mmanager_ptr->printForceSet();
    // print block sets with materials
    CHKERR mmanager_ptr->printMaterialsSet();
 
    string primary_field_name;
 
    if (!is_displacement_field)
      primary_field_name = "SPATIAL_POSITION";
    else
      primary_field_name = "U";
    // create interfaces (could be based on cmake file ?)
 
    MortarContactInterface m_contact(m_field, primary_field_name,
                                     "MESH_NODE_POSITIONS",
                                     is_displacement_field);
    NonlinearElasticElementInterface m_elastic(m_field, primary_field_name,
                                               "MESH_NODE_POSITIONS",
                                               is_displacement_field);
 
    BasicBoundaryConditionsInterface m_boundary(
        m_field, primary_field_name, "MESH_NODE_POSITIONS", "CONTACT_PROB",
        "ELASTIC", is_displacement_field, true,
        &MortarContactProblem::LoadScale::lAmbda, is_partitioned);
 
    CHKERR m_boundary.getCommandLineParameters();
    CHKERR m_contact.getCommandLineParameters();
    CHKERR m_elastic.getCommandLineParameters();
 
    CHKERR m_boundary.addElementFields();
    CHKERR m_contact.addElementFields();
    CHKERR m_elastic.addElementFields();
 
    // build field
    CHKERR m_field.build_fields();
 
    CHKERR m_boundary.createElements();
    CHKERR m_contact.createElements();
    CHKERR m_elastic.createElements();
 
    // Projection on "x" field
    if (!is_displacement_field) {
      Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
      CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
    }
    // Projection on "X" field
    {
      Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
      CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
    }
 
    // build finite elemnts
    CHKERR m_field.build_finite_elements();
 
    auto bit_ref = m_contact.getBitRefLevel();
    // build adjacencies
    CHKERR m_field.build_adjacencies(bit_ref);
 
    // define problems
    CHKERR m_field.add_problem("CONTACT_PROB");
 
    // set refinement level for problem
    CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB", bit_ref);
 
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    SmartPetscObj<DM> dm;
    dm = createSmartDM(m_field.get_comm(), dm_name);
 
    // create dm instance
    CHKERR DMSetType(dm, dm_name);
 
    // set dm datastruture which created mofem datastructures
    CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_ref);
    CHKERR DMSetFromOptions(dm);
    CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
    // add elements to dm
 
    CHKERR m_boundary.setOperators();
    CHKERR m_contact.setOperators();
    CHKERR m_elastic.setOperators();
 
    CHKERR m_boundary.addElementsToDM(dm);
    CHKERR m_contact.addElementsToDM(dm);
    CHKERR m_elastic.addElementsToDM(dm);
 
    CHKERR DMSetUp(dm);
 
    CHKERR m_boundary.setupSolverFunctionSNES();
    CHKERR m_contact.setupSolverFunctionSNES();
    CHKERR m_elastic.setupSolverFunctionSNES();
 
    CHKERR m_boundary.setupSolverJacobianSNES();
    CHKERR m_contact.setupSolverJacobianSNES();
    CHKERR m_elastic.setupSolverJacobianSNES();
 
// #ifndef NDEBUG
//     cout << "<<< checkMPIAIJWithArraysMatrixFillIn <<< " << endl;
//     CHKERR m_field.getInterface<MatrixManager>()
//         ->checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
//             "CONTACT_PROB", -2, -2, 0);
// #endif // NDEBUG
 
    auto D = smartCreateDMVector(dm);
    if (test_num) {
      char testing_options[] = "-ksp_type fgmres "
                               "-pc_type lu "
                               "-pc_factor_mat_solver_type mumps "
                               "-snes_type newtonls "
                               "-snes_linesearch_type basic "
                               "-snes_max_it 20 "
                               "-snes_atol 1e-8 "
                               "-snes_rtol 1e-8 ";
      CHKERR PetscOptionsInsertString(NULL, testing_options);
    }
 
    auto snes = MoFEM::createSNES(m_field.get_comm());
    CHKERR SNESSetDM(snes, dm);
    CHKERR SNESSetFromOptions(snes);
 
    for (int ss = 0; ss != nb_sub_steps; ++ss) {
 
      MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
      MOFEM_LOG_C("WORLD", Sev::inform, "Load lambda %6.4e",
                  MortarContactProblem::LoadScale::lAmbda);
 
      CHKERR SNESSolve(snes, PETSC_NULL, D);
 
      // save on mesh
      CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
      CHKERR m_elastic.postProcessElement(ss);
 
      auto common_data_post = m_contact.commonDataPostProc;
      std::ofstream ofs((std::string("test_simple_contact") + ".txt").c_str());
 
      m_contact.getPostProcSimpleContactPipeline().push_back(
          new MortarContactProblem::OpMakeTestTextFile(
              m_field, primary_field_name, common_data_post, ofs));
 
      CHKERR m_contact.postProcessElement(ss);
      CHKERR m_boundary.postProcessElement(ss);
 
      double energy = m_elastic.elasticElementPtr->getLoopFeEnergy().eNergy;
 
      ofs << "Elastic energy: " << energy << endl;
      ofs.flush();
      ofs.close();
 
      if (test_num) {
 
        auto &nb_gauss_pts = m_contact.nbGaussPts;
        auto &contact_area = m_contact.contactArea;
 
        int expected_nb_gauss_pts;
        double expected_energy, expected_contact_area;
 
        MortarCtestFunctions::check_tests(ss, test_num, expected_energy,
                                          expected_contact_area,
                                          expected_nb_gauss_pts);
 
        constexpr double eps = 1e-6;
        if (m_field.get_comm_rank() == 0) {
          if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong number of active gauss pts %d != %d",
                     expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
          }
          if (std::abs(contact_area[0] - expected_contact_area) > eps) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong active contact area %6.4e != %6.4e",
                     expected_contact_area, contact_area[0]);
          }
        }
 
        if (std::abs(energy - expected_energy) > eps)
          SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                   "Wrong energy %6.4e != %6.4e", expected_energy, energy);
      }
    }
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc
  MoFEM::Core::Finalize();
 
  return 0;
}