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);

    // 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);

    Vec x_elastic, f_elastic;
    CHKERR DMCreateGlobalVector(ep.dmElastic, &x_elastic);
    CHKERR VecDuplicate(x_elastic, &f_elastic);
    CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
    Mat m_elastic;
    // CHKERR DMSetMatType(ep.dmElastic, MATAIJ);
    // CHKERR DMSetMatType(ep.dmElasticSchurStreach, MATAIJ);
    // CHKERR DMSetMatType(ep.dmElasticSchurBubble, MATAIJ);
    CHKERR DMCreateMatrix(ep.dmElastic, &m_elastic);


    TS ts_elastic;
    CHKERR TSCreate(PETSC_COMM_WORLD, &ts_elastic);
    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 ep.solveElastic(ts_elastic, x_elastic);

    CHKERR TSDestroy(&ts_elastic);
    CHKERR VecDestroy(&x_elastic);
    CHKERR VecDestroy(&f_elastic);
    CHKERR MatDestroy(&m_elastic);

  }
  CATCH_ERRORS;

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