mofem_manager.cpp

Go to the documentation of this file.

/** \file mofem_manager.cpp
 *
 */
/* MoFEM is free software: you can redistribute it and/or modify it under the
 * terms of the MIT License.
 */
 
#include <BasicFiniteElements.hpp>
using namespace MoFEM;
 
#include <ElasticMaterials.hpp>
#include <NonlinearElasticElementInterface.hpp>
#include <BasicBoundaryConditionsInterface.hpp>
 
#include <SurfacePressureComplexForLazy.hpp>
 
#ifdef WITH_MODULE_MORTAR_CONTACT
#include <Mortar.hpp>
#endif
 
#ifdef WITH_MODULE_MFRONT_INTERFACE
#include <MFrontMoFEMInterface.hpp>
#endif
 
#ifdef WITH_MODULE_HDIV_CONTACT
using namespace FTensor;
#include <AnalyticalSurfaces.hpp>
#include <RigidBodies.hpp>
#include <BasicFeTools.hpp>
#include <ContactOperators.hpp>
using namespace OpContactTools;
#include <HdivContactInterface.hpp>
 
#endif
 
using DomainEle = VolumeElementForcesAndSourcesCore;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = FaceElementForcesAndSourcesCore;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using PostProcEle = PostProcVolumeOnRefinedMesh;
using PostProcSkinEle = PostProcFaceOnRefinedMesh;
 
constexpr size_t SPACE_DIM = 3;
 
static char help[] = "...\n\n";
 
#include <MoFEMManager.hpp>
using namespace MoFEMManager;
 
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"
                                 "-ksp_atol 1e-10 \n"
                                 "-ksp_rtol 1e-10 \n"
                                 "-snes_monitor \n"
                                 "-snes_max_it 100 \n"
                                 "-snes_linesearch_type bt \n"
                                 "-snes_linesearch_max_it 3 \n"
                                 "-snes_atol 1e-8 \n"
                                 "-snes_rtol 1e-8 \n"
                                 "-ts_monitor \n"
                                 "-ts_alpha_radius 1 \n"
                                 "-ts_monitor \n"
                                 "-mat_mumps_icntl_20 0 \n"
                                 "-mat_mumps_icntl_14 800 \n"
                                 "-mat_mumps_icntl_24 1 \n"
                                 "-mat_mumps_icntl_13 1 \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();
    }
  }
 
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "MManager"));
  LogManager::setLog("MManager");
  MOFEM_LOG_TAG("MManager", "module_manager");
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    MFManager md(m_field);
    CHKERR md.getParams();
 
    Simple *simple = m_field.getInterface<Simple>();
    CHKERR simple->getOptions();
 
    if (md.isPartitioned)
      CHKERR simple->loadFile("");
    else
      CHKERR simple->loadFile("", "");
    simple->getProblemName() = "MoFEM Manager module";
    simple->getDomainFEName() = "ELASTIC";
    vector<string> volume_fe_elements{simple->getDomainFEName()};
    simple->getBoundaryFEName() = "BOUNDARY_ELEMENT";
 
    // setup the modules
    boost::ptr_vector<GenericElementInterface> m_modules;
 
    // Basic Boundary Conditions module should always be first (dirichlet)
    m_modules.push_back(new BasicBoundaryConditionsInterface(
        m_field, "U", "MESH_NODE_POSITIONS", simple->getProblemName(),
        simple->getDomainFEName(), true, md.isQuasiStatic, nullptr,
        md.isPartitioned));
 
#ifdef WITH_MODULE_MORTAR_CONTACT
    m_modules.push_back(new MortarContactInterface(
        m_field, "U", "MESH_NODE_POSITIONS", true, md.isQuasiStatic));
#endif
 
    // // Nonlinear Elastic Element module
    m_modules.push_back(new NonlinearElasticElementInterface(
        m_field, "U", "MESH_NODE_POSITIONS", true, md.isQuasiStatic));
 
#ifdef WITH_MODULE_MFRONT_INTERFACE
    m_modules.push_back(new MFrontMoFEMInterface(
        m_field, "U", "MESH_NODE_POSITIONS", true, md.isQuasiStatic));
    volume_fe_elements.push_back("MFRONT_EL");
#endif
 
#ifdef WITH_MODULE_HDIV_CONTACT
    m_modules.push_back(new HdivContactInterface(
        m_field, "U", "MESH_NODE_POSITIONS", simple->getProblemName(),
        volume_fe_elements, {simple->getBoundaryFEName()}, PETSC_TRUE,
        md.isQuasiStatic, md.isPartitioned));
#endif
 
    CHKERR md.setMainField();
 
    for (auto &&mod : m_modules) {
      CHKERR mod.getCommandLineParameters();
      CHKERR mod.addElementFields();
      CHKERR mod.createElements();
    }
 
    CHKERR m_field.build_fields();
    CHKERR md.defineDM();
 
    for (auto &&mod : m_modules)
      CHKERR mod.addElementsToDM(md.dM);
 
    Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
    CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
    CHKERR m_field.build_finite_elements();
    CHKERR m_field.build_adjacencies(simple->getBitRefLevel());
 
    CHKERR DMSetUp(md.dM);
 
    md.monitorPtr->preProcessHook = []() { return 0; };
    md.monitorPtr->operatorHook = []() { return 0; };
    md.monitorPtr->postProcessHook = [&]() {
      MoFEMFunctionBeginHot;
      auto ts_time = md.monitorPtr->ts_t;
      auto ts_step = md.monitorPtr->ts_step;
      auto &ts_u = md.monitorPtr->ts_u;
 
      for (auto &&mod : m_modules) {
        CHKERR mod.updateElementVariables();
        if (ts_step % md.saveEveryNthStep == 0)
          CHKERR mod.postProcessElement(ts_step);
      }
 
      if (md.doPrintMaxMin == PETSC_TRUE) {
        CHKERR md.printMaxMin(md.vScatter[0], "Ux", ts_u, ts_time);
        CHKERR md.printMaxMin(md.vScatter[1], "Uy", ts_u, ts_time);
        CHKERR md.printMaxMin(md.vScatter[2], "Uz", ts_u, ts_time);
      }
 
      MoFEMFunctionReturnHot(0);
    };
 
    auto t_type = md.getTSType();
    for (auto &&mod : m_modules) {
      CHKERR mod.setMonitorPtr(md.monitorPtr);
      CHKERR mod.setOperators();
      CHKERR mod.setupSolverFunctionTS(t_type);
      CHKERR mod.setupSolverJacobianTS(t_type);
    }
 
    CHKERR md.tsSetup();
    CHKERR md.tsSolve();
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}