ep.cpp

Implementation of mix-element for Large strains

Todo:

Implementation of plasticity

/**
 * \file ep.cpp
 * \example 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/>. */

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

using namespace MoFEM;
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);

  try {

    // Get mesh file
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", 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;
    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();
    // CHKERR ep.addMaterial_HMHHStVenantKirchhoff(1, LAMBDA(1, 0), MU(1, 0), 0);
    CHKERR ep.addMaterial_HMHMooneyRivlin(1, MU(1, 0) / 2., 0, LAMBDA(1, 0), 0);

    CHKERR ep.setElasticElementOps(1);
    CHKERR ep.setElasticElementToTs(ep.dmElastic);

    SmartPetscObj<Vec> x_elastic, f_elastic;
    CHKERR DMCreateGlobalVector_MoFEM(ep.dmElastic, x_elastic);
    f_elastic = smartVectorDuplicate(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 ep.setUpTSElastic(ts_elastic, m_elastic, f_elastic);
    CHKERR TSSetTime(ts_elastic, time);
    CHKERR ep.solveElastic(ts_elastic, x_elastic);
  }
  CATCH_ERRORS;

  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
}