small_strain_plasticity_problem_paraboloidal.cpp

Go to the documentation of this file.

/** \file small_strain_plasticity.cpp
 * \ingroup small_strain_plasticity
 * \brief Small strain plasticity example
 *
 */
 
/* 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 <MoFEM.hpp>
using namespace MoFEM;
 
#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
 
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
 
#include <adolc/adolc.h>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
 
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <SurfacePressure.hpp>
#include <NodalForce.hpp>
#include <EdgeForce.hpp>
#include <DirichletBC.hpp>
 
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
 
#include <string>
 
using namespace boost::numeric;
 
static char help[] = "...\n"
                     "\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, "-my_file", mesh_file_name, 255,
                                 &flg);
    if (flg != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
    }
 
    PetscInt order;
    CHKERR PetscOptionsGetInt(PETSC_NULL, "-my_order", &order, &flg);
    if (flg != PETSC_TRUE) {
      order = 2;
    }
    PetscInt bubble_order;
    CHKERR PetscOptionsGetInt(PETSC_NULL, "-my_bubble_order", &bubble_order,
                              &flg);
    if (flg != PETSC_TRUE) {
      bubble_order = order;
    }
 
    // use this if your mesh is partitioned and you run code on parts,
    // you can solve very big problems
    PetscBool is_partitioned = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "-my_is_partitioned",
                               &is_partitioned, &flg);
 
    // Read mesh to MOAB
    if (is_partitioned == PETSC_TRUE) {
      const char *option;
      option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
               "PARALLEL_PARTITION;";
      rval = moab.load_file(mesh_file_name, 0, option);
      CHKERRQ_MOAB(rval);
    } else {
      const char *option;
      option = "";
      rval = moab.load_file(mesh_file_name, 0, option);
      CHKERRQ_MOAB(rval);
    }
 
    SmallStrainParaboloidalPlasticity cp;
    {
      cp.tAgs.resize(3);
      cp.tAgs[0] = 0;
      cp.tAgs[1] = 1;
      cp.tAgs[2] = 2;
      cp.tOl = 1e-12;
 
      double young_modulus = 1;
      double poisson_ratio = 0.25;
      cp.sIgma_yt = 1;
      cp.sIgma_yc = 1;
 
      cp.Ht = 0.1;
      cp.Hc = 0.1;
 
      cp.nup = 0.3;
 
      {
        CHKERR PetscOptionsGetReal(0, "-my_young_modulus", &young_modulus, 0);
        CHKERR PetscOptionsGetReal(0, "-my_poisson_ratio", &poisson_ratio, 0);
        cp.mu = MU(young_modulus, poisson_ratio);
        cp.lambda = LAMBDA(young_modulus, poisson_ratio);
        CHKERR PetscOptionsGetReal(0, "-my_sigma_yt", &cp.sIgma_yt, 0);
        CHKERR PetscOptionsGetReal(0, "-my_sigma_yc", &cp.sIgma_yc, 0);
 
        CHKERR PetscOptionsGetReal(0, "-my_Ht", &cp.Ht, 0);
        CHKERR PetscOptionsGetReal(0, "-my_Hc", &cp.Hc, 0);
 
        CHKERR PetscOptionsGetReal(0, "-my_nt", &cp.nt, 0);
        CHKERR PetscOptionsGetReal(0, "-my_nc", &cp.nc, 0);
 
        CHKERR PetscOptionsGetReal(0, "-my_nup", &cp.nup, 0);
      }
 
      PetscPrintf(PETSC_COMM_WORLD, "young_modulus = %4.2e \n", young_modulus);
      PetscPrintf(PETSC_COMM_WORLD, "poisson_ratio = %4.2e \n", poisson_ratio);
 
      PetscPrintf(PETSC_COMM_WORLD, "sIgma_yt = %4.2e \n", cp.sIgma_yt);
      PetscPrintf(PETSC_COMM_WORLD, "sIgma_yc = %4.2e \n", cp.sIgma_yt);
 
      PetscPrintf(PETSC_COMM_WORLD, "nup = %4.2e \n", cp.nup);
 
      cp.sTrain.resize(6, false);
      cp.sTrain.clear();
      cp.plasticStrain.resize(6, false);
      cp.plasticStrain.clear();
      cp.internalVariables.resize(2, false);
      cp.internalVariables.clear();
      cp.createMatAVecR();
      cp.snesCreate();
      // CHKERR SNESSetFromOptions(cp.sNes);
    }
 
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    // ref meshset ref level 0
    CHKERR m_field.seed_ref_level_3D(0, BitRefLevel().set(0));
    vector<BitRefLevel> bit_levels;
    bit_levels.push_back(BitRefLevel().set(0));
 
    {
 
      BitRefLevel problem_bit_level = bit_levels.back();
      Range CubitSideSets_meshsets;
      CHKERR m_field.get_cubit_meshsets(SIDESET, CubitSideSets_meshsets);
 
      // Fields
      CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
      CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1,
                               AINSWORTH_LEGENDRE_BASE, 3);
      // add entitities (by tets) to the field
      CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
      CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
 
      // set app. order
      CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", bubble_order);
      CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
      CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
      CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
 
      CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS",
                                     order > 1 ? 2 : 1);
      CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS",
                                     order > 1 ? 2 : 1);
      CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS",
                                     order > 1 ? 2 : 1);
      CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
 
      // build field
      CHKERR m_field.build_fields();
      {
        // 10 node tets
        Projection10NodeCoordsOnField ent_method_material(
            m_field, "MESH_NODE_POSITIONS");
        CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material, 0);
      }
 
      // FE
      CHKERR m_field.add_finite_element("PLASTIC");
      // Define rows/cols and element data
      CHKERR m_field.modify_finite_element_add_field_row("PLASTIC",
                                                         "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_col("PLASTIC",
                                                         "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_data("PLASTIC",
                                                          "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_data(
          "PLASTIC", "MESH_NODE_POSITIONS");
      CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTET, "PLASTIC");
 
      // Add Neumann forces
      CHKERR MetaNeummanForces::addNeumannBCElements(m_field, "DISPLACEMENT");
      CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
      CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
 
      // build finite elements
      CHKERR m_field.build_finite_elements();
      // build adjacencies
      CHKERR m_field.build_adjacencies(problem_bit_level);
 
      DMType dm_name = "PLASTIC_PROB";
      CHKERR DMRegister_MoFEM(dm_name);
      // craete dm instance
      DM dm;
      CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
      CHKERR DMSetType(dm, dm_name);
      {
        // set dm datastruture which created mofem datastructures
        CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, problem_bit_level);
        CHKERR DMSetFromOptions(dm);
        CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
        // add elements to dm
        CHKERR DMMoFEMAddElement(dm, "PLASTIC");
        CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
        CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
        CHKERR DMSetUp(dm);
      }
 
      // create matrices
      Vec F, D;
      Mat Aij;
      {
        CHKERR DMCreateGlobalVector_MoFEM(dm, &D);
        CHKERR VecDuplicate(D, &F);
        CHKERR DMCreateMatrix_MoFEM(dm, &Aij);
        CHKERR VecZeroEntries(D);
        CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
        CHKERR MatZeroEntries(Aij);
      }
 
      // assemble Aij and F
      DirichletDisplacementBc dirichlet_bc(m_field, "DISPLACEMENT", Aij, D, F);
      dirichlet_bc.methodsOp.push_back(
          new TimeForceScale("-my_displacements_history", false));
 
      SmallStrainPlasticity small_strain_plasticity(m_field);
      {
        PetscBool bbar = PETSC_TRUE;
        CHKERR PetscOptionsGetBool(0, "-my_bbar", &bbar, 0);
        small_strain_plasticity.commonData.bBar = bbar;
      }
      {
        small_strain_plasticity.feRhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        small_strain_plasticity.feRhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpCalculateStress(
                m_field, "DISPLACEMENT", small_strain_plasticity.commonData, cp,
                false));
        small_strain_plasticity.feRhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpAssembleRhs(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        small_strain_plasticity.feLhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        small_strain_plasticity.feLhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpCalculateStress(
                m_field, "DISPLACEMENT", small_strain_plasticity.commonData, cp,
                false));
        small_strain_plasticity.feLhs.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpAssembleLhs(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpCalculateStress(
                m_field, "DISPLACEMENT", small_strain_plasticity.commonData, cp,
                false));
        small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
            new SmallStrainPlasticity::OpUpdate(
                "DISPLACEMENT", small_strain_plasticity.commonData));
        CHKERR small_strain_plasticity.createInternalVariablesTag();
      }
 
      // Setting finite element method for applying tractions
      boost::ptr_map<string, NeummanForcesSurface> neumann_forces;
      boost::ptr_map<string, NodalForce> nodal_forces;
      boost::ptr_map<string, EdgeForce> edge_forces;
      TimeForceScale time_force_scale("-my_load_history", false);
      {
        // forces and pressures on surface
        CHKERR MetaNeummanForces::setMomentumFluxOperators(
            m_field, neumann_forces, PETSC_NULL, "DISPLACEMENT");
        // noadl forces
        CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, PETSC_NULL,
                                             "DISPLACEMENT");
        // edge forces
        CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, PETSC_NULL,
                                            "DISPLACEMENT");
        for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
                 neumann_forces.begin();
             mit != neumann_forces.end(); mit++) {
          mit->second->methodsOp.push_back(
              new TimeForceScale("-my_load_history", false));
        }
        for (boost::ptr_map<string, NodalForce>::iterator mit =
                 nodal_forces.begin();
             mit != nodal_forces.end(); mit++) {
          mit->second->methodsOp.push_back(
              new TimeForceScale("-my_load_history", false));
        }
        for (boost::ptr_map<string, EdgeForce>::iterator mit =
                 edge_forces.begin();
             mit != edge_forces.end(); mit++) {
          mit->second->methodsOp.push_back(
              new TimeForceScale("-my_load_history", false));
        }
      }
 
      // Adding elements to DMSnes
      {
        // Rhs
        CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, &dirichlet_bc,
                                      NULL);
        {
          boost::ptr_map<string, NeummanForcesSurface>::iterator fit;
          fit = neumann_forces.begin();
          for (; fit != neumann_forces.end(); fit++) {
            CHKERR DMMoFEMSNESSetFunction(
                dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
          }
          fit = nodal_forces.begin();
          for (; fit != nodal_forces.end(); fit++) {
            CHKERR DMMoFEMSNESSetFunction(
                dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
          }
          fit = edge_forces.begin();
          for (; fit != edge_forces.end(); fit++) {
            CHKERR DMMoFEMSNESSetFunction(
                dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
          }
        }
        CHKERR DMMoFEMSNESSetFunction(dm, "PLASTIC",
                                      &small_strain_plasticity.feRhs,
                                      PETSC_NULL, PETSC_NULL);
        CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, NULL,
                                      &dirichlet_bc);
 
        // Lhs
        CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, &dirichlet_bc,
                                      NULL);
        CHKERR DMMoFEMSNESSetJacobian(
            dm, "PLASTIC", &small_strain_plasticity.feLhs, NULL, NULL);
        CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, NULL,
                                      &dirichlet_bc);
      }
 
      // Create Time Stepping solver
      SNES snes;
      SnesCtx *snes_ctx;
      {
        CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
        // CHKERR SNESSetDM(snes,dm);
        CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
        CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
        CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
        CHKERR SNESSetFromOptions(snes);
      }
 
      PostProcVolumeOnRefinedMesh post_proc(m_field);
      {
        CHKERR post_proc.generateReferenceElementMesh();
        CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
        CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
        CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
      }
 
      {
        double final_time = 1, delta_time = 0.1;
        CHKERR PetscOptionsGetReal(0, "-my_final_time", &final_time, 0);
        CHKERR PetscOptionsGetReal(0, "-my_delta_time", &delta_time, 0);
        double delta_time0 = delta_time;
 
        // CHKERR MatView(Aij,PETSC_VIEWER_DRAW_SELF);
        // std::string wait;
        // std::cin >> wait;
 
        Vec D0;
        CHKERR VecDuplicate(D, &D0);
 
        NumeredDofEntity_multiIndex_uid_view_ordered load_path_dofs_view;
        {
          Range node_set;
          for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field, "LoadPath", cit)) {
            EntityHandle meshset = cit->getMeshset();
            Range nodes;
            rval = moab.get_entities_by_type(meshset, MBVERTEX, nodes, true);
            CHKERRQ_MOAB(rval);
            node_set.merge(nodes);
          }
          PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
                      node_set.size());
          const Problem *problem_ptr;
          CHKERR m_field.get_problem("PLASTIC_PROB", &problem_ptr);
          boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
              problem_ptr->getNumeredDofsRows();
          Range::iterator nit = node_set.begin();
          for (; nit != node_set.end(); nit++) {
            NumeredDofEntityByEnt::iterator dit, hi_dit;
            dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
            hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
            for (; dit != hi_dit; dit++) {
              load_path_dofs_view.insert(*dit);
            }
          }
        }
 
        int step = 0;
        double t = 0;
        SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
        for (; t < final_time; step++) {
 
          t += delta_time;
          PetscPrintf(PETSC_COMM_WORLD, "Step %d Time %6.4g final time %3.2g\n",
                      step, t, final_time);
 
          // set time
          {
            dirichlet_bc.ts_t = t;
            {
              boost::ptr_map<string, NeummanForcesSurface>::iterator fit;
              fit = neumann_forces.begin();
              for (; fit != neumann_forces.end(); fit++) {
                fit->second->getLoopFe().ts_t = t;
              }
              fit = nodal_forces.begin();
              for (; fit != nodal_forces.end(); fit++) {
                fit->second->getLoopFe().ts_t = t;
              }
              fit = edge_forces.begin();
              for (; fit != edge_forces.end(); fit++) {
                fit->second->getLoopFe().ts_t = t;
              }
            }
          }
 
          // solve problem at step
          {
            // dirichlet_bc.snes_B = Aij;
            // small_strain_plasticity.feLhs.snes_B = Aij;
            CHKERR VecAssemblyBegin(D);
            CHKERR VecAssemblyEnd(D);
            CHKERR VecCopy(D, D0);
            if (step == 0 || reason < 0) {
              CHKERR SNESSetLagJacobian(snes, -2);
            } else {
              CHKERR SNESSetLagJacobian(snes, -1);
            }
            CHKERR SNESSolve(snes, PETSC_NULL, D);
 
            int its;
            CHKERR SNESGetIterationNumber(snes, &its);
 
            // CHKERR MatView(Aij,PETSC_VIEWER_DRAW_SELF);
            // std::string wait;
            // std::cin >> wait;
 
            CHKERR PetscPrintf(PETSC_COMM_WORLD,
                               "number of Newton iterations = %D\n", its);
 
            CHKERR SNESGetConvergedReason(snes, &reason);
 
            if (reason < 0) {
              t -= delta_time;
              delta_time *= 0.1;
              CHKERR VecCopy(D0, D);
            } else {
              const int its_d = 7;
              const double gamma = 0.5;
              const double max_reudction = 1;
              const double min_reduction = 1e-1;
              double reduction;
              reduction = pow((double)its_d / (double)(its + 1), gamma);
              if (delta_time >= max_reudction * delta_time0 && reduction > 1) {
                reduction = 1;
              } else if (delta_time <= min_reduction * delta_time0 &&
                         reduction < 1) {
                reduction = 1;
              }
 
              CHKERR PetscPrintf(
                  PETSC_COMM_WORLD,
                  "reduction delta_time = %6.4e delta_time = %6.4e\n",
                  reduction, delta_time);
              delta_time *= reduction;
              if (reduction > 1 && delta_time < min_reduction * delta_time0) {
                delta_time = min_reduction * delta_time0;
              }
 
              CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES,
                                              SCATTER_REVERSE);
              CHKERR DMoFEMLoopFiniteElements(
                  dm, "PLASTIC", &small_strain_plasticity.feUpdate);
 
              {
                double scale;
                CHKERR time_force_scale.getForceScale(t, scale);
                NumeredDofEntity_multiIndex_uid_view_ordered::iterator it,
                    hi_it;
                it = load_path_dofs_view.begin();
                hi_it = load_path_dofs_view.end();
                for (; it != hi_it; it++) {
                  PetscPrintf(PETSC_COMM_WORLD,
                              "load_path %s [ %d ] %6.4e %6.4e %6.4e\n",
                              (*it)->getName().c_str(), (*it)->getDofCoeffIdx(),
                              (*it)->getFieldData(), t, scale);
                }
              }
 
              // Save data on mesh
              {
                CHKERR DMoFEMLoopFiniteElements(dm, "PLASTIC", &post_proc);
                string out_file_name;
                {
                  std::ostringstream stm;
                  stm << "out_" << step << ".h5m";
                  out_file_name = stm.str();
                }
                CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
                                   out_file_name.c_str());
                rval = post_proc.postProcMesh.write_file(
                    out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
                CHKERRQ_MOAB(rval);
              }
            }
          }
        }
 
        CHKERR VecDestroy(&D0);
      }
 
      {
        CHKERR MatDestroy(&Aij);
        CHKERR VecDestroy(&F);
        CHKERR VecDestroy(&D);
      }
    }
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
  return 0;
}