ep.cpp

Go to the documentation of this file.

/**
 * \file ep.cpp
 * \example mofem/users_modules/eshelbian_plasticity/ep.cpp
 *
 * \brief Implementation of mix-element for Large strains
 *
 * \todo Implementation of plasticity
 *
 */
 
/* This file is part of MoFEM.
 * 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/>. */
 
constexpr int adolc_tag = 1;
 
#define SINGULARITY
 
#include <MoFEM.hpp>
using namespace MoFEM;
using namespace boost::multi_index;
using namespace boost::multiprecision;
using namespace boost::numeric;
#include <cmath>
 
#include <cholesky.hpp>
#ifdef ENABLE_PYTHON_BINDING 
  #include <boost/python.hpp>
  #include <boost/python/def.hpp>
  #include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
 
#include <MatOps.hpp>
#include <MatElastic.hpp>
#include <EshelbianPlasticity.hpp>
using namespace EshelbianPlasticity;
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  // initialize petsc
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
  LogManager::setLog("EP");
  MOFEM_LOG_TAG("EP", "ep");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
  LogManager::setLog("EPSELF");
  MOFEM_LOG_TAG("EPSELF", "ep");
 
  #ifdef ENABLE_PYTHON_BINDING
  Py_Initialize();
  np::initialize();
  MOFEM_LOG("EP", Sev::inform) << "Python initialised";
#else
  MOFEM_LOG("EP", Sev::inform) << "Python NOT initialised";
#endif
 
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
  LogManager::setLog("EPSYNC");
  MOFEM_LOG_TAG("EPSYNC", "ep");
 
  try {
 
    // Get mesh file
    PetscBool flg = PETSC_FALSE;
    char mesh_file_name[255] = "";
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-my_file", mesh_file_name,
                                 255, &flg);
    if (!flg) {
      CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-file_name",
                                   mesh_file_name, 255, &flg);
    }
    if (!flg) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Missing required mesh option: use -file_name or -my_file");
    }
    char restart_file[255] = "no_restart_file";
    PetscBool restart_flg = PETSC_TRUE;
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file, 255,
                                 &restart_flg);
    double time = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-time", &time, PETSC_NULLPTR);
 
 
    MOFEM_LOG("EP", Sev::inform) << "Mesh file: " << mesh_file_name;
    MOFEM_LOG("EP", Sev::inform) << "Restart file: " << restart_file;
    MOFEM_LOG("EP", Sev::inform) << "Time: " << time;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    DMType dm_name_mg = "DMMOFEM_MG";
    CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
 
    // Create MoAB database
    moab::Core moab_core;
    moab::Interface &moab = moab_core;
 
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
    // Read mesh to MOAB
    PetscBool fully_distributed = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
                               &fully_distributed, PETSC_NULLPTR);
    if (fully_distributed) {
      const char *option;
      if (pcomm->proc_config().proc_size() == 1)
        option = "";
      else
        option = "PARALLEL=READ_PART;"
                 "PARALLEL_RESOLVE_SHARED_ENTS;"
                 "PARTITION=PARALLEL_PARTITION";
      CHKERR moab.load_file(mesh_file_name, 0, option);
    } else {
      CHKERR CommInterface::loadFileRootProcAllRestDistributed(
          moab, mesh_file_name, SPACE_DIM);
    }
 
    // Create MoFEM database and link it to MoAB
    MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
    MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
 
    CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
 
    // Register mofem DM
    CHKERR DMRegister_MoFEM("DMMOFEM");
 
    BitRefLevel bit_level0 = BitRefLevel().set(0);
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, bit_level0);
 
    // Data stuctures
    EshelbianCore ep(m_field);
 
    auto meshset_ptr = get_temp_meshset_ptr(moab);
    CHKERR moab.add_entities(
        *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
    auto get_adj = [&](Range ents, int dim) {
      Range adj;
      CHKERR moab.get_adjacencies(ents, dim, false, adj,
                                  moab::Interface::UNION);
      return adj;
    };
    CHKERR moab.add_entities(
        *meshset_ptr,
        get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
                    .subset_by_dimension(SPACE_DIM),
                2));
    CHKERR moab.add_entities(
        *meshset_ptr,
        get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
                    .subset_by_dimension(SPACE_DIM),
                1));
    CHKERR moab.add_entities(
        *meshset_ptr,
        get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
                    .subset_by_dimension(SPACE_DIM),
                0));
 
    // create growing crack surface meshset
    CHKERR ep.createCrackSurfaceMeshset();
 
    CHKERR ep.getSpatialDispBc();
    CHKERR ep.getSpatialRotationBc();
    CHKERR ep.getSpatialTractionBc();
    CHKERR ep.getSpatialTractionFreeBc();
    CHKERR ep.getExternalStrain();
    CHKERR ep.setBlockTagsOnSkin();
 
    CHKERR ep.createExchangeVectors(Sev::inform);
 
    CHKERR ep.projectMaterialTags();
    CHKERR ep.addFields(*meshset_ptr);
    CHKERR ep.projectGeometry(*meshset_ptr, 0.0);
    CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
    CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
    CHKERR ep.addDMs();
 
    const char *list_materials[EshelbianCore::MaterialModel::LastMaterial] = {
        "stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean", "meta"};
    PetscInt choice_material = EshelbianCore::materialModel;
    CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-material",
                                list_materials,
                                EshelbianCore::MaterialModel::LastMaterial,
                                &choice_material, PETSC_NULLPTR);
    EshelbianCore::materialModel =
        static_cast<EshelbianCore::MaterialModel>(choice_material);
 
    auto material_model_names_dictionary = std::map<std::string, std::string>{
        {"grad", "stretchH1AtPts"}, {"P", "PAtPts"}, {"P_dF", "PAtPts_du"}};
    MatOps::MatElastic::useDeformationGradient =
        true; // by default use deformation gradient as common data
 
    switch (choice_material) {
    case EshelbianCore::MaterialModel::StVenantKirchhoff: {
      MOFEM_LOG("EP", Sev::inform) << "St. Venant-Kirchhoff material model";
      MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
      ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
      auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
          MatOps::ELASTICITY::STVENANTKIRCHHOFF, MatOps::MODEL_3D>(
          MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                      material_model_names_dictionary),
          MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
      CHKERR physical_equations_ptr->getOptions(&m_field);
      CHKERR physical_equations_ptr->recordTape();
      CHKERR ep.addMaterial_Core(physical_equations_ptr);
    };
    case EshelbianCore::MaterialModel::MooneyRivlin: {
      MOFEM_LOG("EP", Sev::inform) << "Mooney-Rivlin material model";
      ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
      auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
          MatOps::ELASTICITY::MOONEYRIVLINWRIGGERSEQ63, MatOps::MODEL_3D>(
          MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                      material_model_names_dictionary),
          MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
      CHKERR physical_equations_ptr->getOptions(&m_field);
      CHKERR physical_equations_ptr->recordTape();
      CHKERR ep.addMaterial_Core(physical_equations_ptr);
    }; break;
    case EshelbianCore::MaterialModel::Hencky:
      MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
      CHKERR ep.addMaterial_Hencky(5., 0.25);
      break;
    case EshelbianCore::MaterialModel::Neohookean:
      MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
      CHKERR ep.addMaterial_HMHNeohookean(adolc_tag, 1.0, 1.0);
      break;
    case EshelbianCore::MaterialModel::Meta: {
      MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
      ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
      auto physical_equations_ptr =
          MatOps::createMatOpsPhysicalEquationsPtr<MatOps::ELASTICITY::META,
                                                  MatOps::MODEL_3D>(
              MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                          material_model_names_dictionary),
              MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
      CHKERR physical_equations_ptr->getOptions(&m_field);
      CHKERR physical_equations_ptr->recordTape();
      CHKERR ep.addMaterial_Core(physical_equations_ptr);
    }; break;
    default:
      SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
      break;
    }
 
    CHKERR ep.setElasticElementOps(adolc_tag);
    CHKERR ep.setElasticElementToTs(ep.dmElastic);
 
    auto x_elastic = createDMVector(ep.dmElastic);
    auto f_elastic = vectorDuplicate(x_elastic);
    CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
 
    PetscInt start_step = 0;
    if (restart_flg) {
      PetscViewer viewer;
      CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
                                   FILE_MODE_READ, &viewer);
      CHKERR VecLoad(x_elastic, viewer);
      CHKERR PetscViewerDestroy(&viewer);
      CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
 
      std::string restart_file_str(restart_file);
      std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
      std::smatch match;
      if (std::regex_search(restart_file_str, match, restart_pattern)) {
        start_step = std::stoi(match[1]);
      } else {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Restart file name must be in the format restart_##.dat");
      }
    }
 
    auto ts_elastic = createTS(PETSC_COMM_WORLD);
    CHKERR TSSetType(ts_elastic, TSBEULER);
    CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
    TSAdapt adapt;
    CHKERR TSGetAdapt(ts_elastic, &adapt);
    CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
 
    switch (EshelbianCore::solverType) {
    case EshelbianCore::SolverType::TimeSolver:
      MOFEM_LOG("EP", Sev::inform) << "Solver type: TS";
      CHKERR TSSetTime(ts_elastic, time);
      CHKERR TSSetStepNumber(ts_elastic, start_step);
      CHKERR ep.solveElastic(ts_elastic, x_elastic); 
      break;
    case EshelbianCore::SolverType::DynamicRelaxation:
      MOFEM_LOG("EP", Sev::inform) << "Solver type: Dynamic Relaxation";
      CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
      break;
    case EshelbianCore::SolverType::Cohesive:
      MOFEM_LOG("EP", Sev::inform) << "Solver type: Cohesive";
      CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR ep.solveCohesiveCrackGrowth(ts_elastic, x_elastic, start_step,
                                         time);
      break;
    case EshelbianCore::SolverType::LoadFactor:
      MOFEM_LOG("EP", Sev::inform) << "Solver type: Load Factor";
      CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR ep.solveLoadFactor(ts_elastic, x_elastic, start_step, time);
      break;
    case EshelbianCore::SolverType::ShapeOptimization:
      MOFEM_LOG("EP", Sev::inform) << "Solver type: Shape Optimization";
      CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR ep.solveSchapeOptimisation(ts_elastic, x_elastic, start_step,
                                          time);
      break;
    case EshelbianCore::SolverType::TestTopologicalDerivative:
      MOFEM_LOG("EP", Sev::inform)
          << "Solver type: Test Topological Derivative";
      CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR ep.solveTestTopologicalDerivative(ts_elastic, x_elastic,
                                               start_step, time);
      break;
    default:
      SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Unknown solver type");
      break;
    }
 
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
#ifdef ENABLE_PYTHON_BINDING
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
#endif
}