ep.cpp File Reference

#include <cholesky.hpp>
#include <BasicFiniteElements.hpp>
#include <EshelbianPlasticity.hpp>
#include <TSEPAdapt.hpp>

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Variables
constexpr int	adolc_tag = 1
static char	help [] = "...\n\n"

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Examples

ep.cpp.

Definition at line 47 of file ep.cpp.

                                 {
 
  // 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");
 
#ifdef PYTHON_SDF
  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_TRUE;
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
                                 255, &flg);
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-file_name", mesh_file_name,
                                 255, &flg);
    char restart_file[255];
    PetscBool restart_flg = PETSC_TRUE;
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-restart", restart_file, 255,
                                 &restart_flg);
    double time = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-time", &time, PETSC_NULL);
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // 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
    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);
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
 
    // 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);
 
    CHKERR ep.getSpatialDispBc();
    CHKERR ep.getSpatialRotationBc();
    CHKERR ep.getSpatialTractionBc();
    CHKERR ep.getSpatialTractionFreeBc();
 
    CHKERR ep.addFields();
    CHKERR ep.addVolumeFiniteElement();
    CHKERR ep.addBoundaryFiniteElement();
    CHKERR ep.addDMs();
 
    enum MaterialModel { StVenantKirchhoff, MooneyRivlin, Hencky, LastMaterial };
    const char *list_materials[LastMaterial] = {"stvenant_kirchhoff",
                                                "mooney_rivlin", "hencky"};
    PetscInt choice_material = MooneyRivlin;
    CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
                                LastMaterial, &choice_material, PETSC_NULL);
 
    switch (choice_material) {
    case StVenantKirchhoff:
      MOFEM_LOG("EP", Sev::inform)
          << "StVenantKirchhoff material model";
      CHKERR ep.addMaterial_HMHHStVenantKirchhoff(adolc_tag, LAMBDA(1, 0.25),
                                                  MU(1, 0.25), 0);
      break;
    case MooneyRivlin:
      MOFEM_LOG("EP", Sev::inform) << "MooneyRivlin material model";
      CHKERR ep.addMaterial_HMHMooneyRivlin(adolc_tag, MU(1, 0) / 2., 0,
                                            LAMBDA(1, 0), 0);
      break;
    case Hencky:
      MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
      CHKERR ep.addMaterial_Hencky(5., 0.25);
      break;
    default:
      SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
      break;
    }
 
    CHKERR ep.setElasticElementOps(adolc_tag);
    CHKERR ep.setElasticElementToTs(ep.dmElastic);
 
    SmartPetscObj<Vec> x_elastic, f_elastic;
    CHKERR DMCreateGlobalVector_MoFEM(ep.dmElastic, x_elastic);
    f_elastic = vectorDuplicate(x_elastic);
    CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
    SmartPetscObj<Mat> m_elastic;
    CHKERR DMCreateMatrix_MoFEM(ep.dmElastic, m_elastic);
 
    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 DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
                                     SCATTER_REVERSE);
    }
 
    auto ts_elastic = createTS(PETSC_COMM_WORLD);
    CHKERR TSSetType(ts_elastic, TSBEULER);
    CHKERR TSAdaptRegister(TSADAPTEP, TSAdaptCreate_EP);
    TSAdapt adapt;
    CHKERR TSGetAdapt(ts_elastic, &adapt);
    CHKERR TSAdaptSetType(adapt, TSADAPTEP);
    // CHKERR TSSetFromOptions(ts_elastic);
 
    CHKERR TSSetTime(ts_elastic, time);
    CHKERR ep.solveElastic(ts_elastic, m_elastic, f_elastic, x_elastic);
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
#ifdef PYTHON_SDF
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
#endif
}

Variable Documentation

adolc_tag

constexpr int adolc_tag = 1

constexpr

Examples

ep.cpp.

Definition at line 26 of file ep.cpp.

help

char help[] = "...\n\n"

static

Examples

ep.cpp.

Definition at line 45 of file ep.cpp.