crack_mesh_cut.cpp

Go to the documentation of this file.

/* 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/>. */
 
#ifdef WITH_TETGEN
 
#include <tetgen.h>
#ifdef REAL
#undef REAL
#endif
 
#endif
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <NeoHookean.hpp>
#include <Hooke.hpp>
 
#include <CrackFrontElement.hpp>
#include <ComplexConstArea.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
#include <VolumeLengthQuality.hpp>
#include <CrackPropagation.hpp>
#include <CPSolvers.hpp>
#include <CPMeshCut.hpp>
#include <AnalyticalFun.hpp>
 
static char help[] = "Calculate crack release energy and crack propagation"
                     "\n\n";
 
using namespace FractureMechanics;
 
bool CrackPropagation::debug = false;
 
bool CrackPropagation::parallelReadAndBroadcast =
    false; ///< That is unstable development, for
           ///< some meshses (propably generated by
           ///< cubit) this is not working. Error can
           ///< be attributed to bug in MOAB.
 
int main(int argc, char *argv[]) {
 
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  try {
 
    PetscBool flg = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-debug", &flg, PETSC_NULL);
    if (flg == PETSC_TRUE)
      CrackPropagation::debug = true;
 
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
                                 255, &flg);
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    // Read mesh file
    if (flg) {
      const char *option;
      option = "";
      CHKERR moab.load_file(mesh_file_name, 0, option);
    }
 
    // Create mofem database
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    auto add_last_twenty = [&]() {
      MoFEMFunctionBegin;
      BitRefLevel mask;
      for (int ll = 20; ll != BITREFLEVEL_SIZE; ++ll)
        mask.set(ll);
      CHKERR m_field.getInterface<BitRefManager>()
          ->addToDatabaseBitRefLevelByType(MBTET, mask, BitRefLevel().set());
      MoFEMFunctionReturn(0);
    };
 
    auto add_first_four = [&]() {
      MoFEMFunctionBegin;
      BitRefLevel mask;
      for (int ll = 0; ll != 5; ++ll)
        mask.set(ll);
      CHKERR m_field.getInterface<BitRefManager>()
          ->addToDatabaseBitRefLevelByType(MBTET, mask, BitRefLevel().set());
      MoFEMFunctionReturn(0);
    };
 
    CHKERR add_last_twenty();
    CHKERR add_first_four();
 
    if (CrackPropagation::debug)
      CHKERR m_field.getInterface<BitRefManager>()
          ->writeEntitiesAllBitLevelsByType(BitRefLevel().set(), MBTET,
                                            "all_start.vtk", "VTK", "");
 
    // Create data structure for crack propagation
    CrackPropagation cp(m_field);
    CHKERR cp.getOptions();
 
    if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
      Range ents;
      CHKERR m_field.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
          ents);
      if (!ents.empty()) {
        EntityHandle meshset;
        CHKERR moab.create_meshset(MESHSET_SET, meshset);
        CHKERR moab.add_entities(meshset, ents);
        CHKERR moab.write_file("ents_not_in_database_to_delete.vtk", "VTK", "",
                               &meshset, 1);
        CHKERR moab.delete_entities(&meshset, 1);
 
        Range meshsets;
        CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets, false);
        for (auto m : meshsets)
          CHKERR moab.remove_entities(m, ents);
        CHKERR moab.delete_entities(ents);
      }
      CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
          cp.crackAccelerationFactor, "cutting_surface.vtk", QUIET,
          CrackPropagation::debug);
 
      CHKERR cp.getInterface<CPMeshCut>()->refineMesh(nullptr, false, QUIET,
                                                      CrackPropagation::debug);
    } else {
 
      if (!cp.getInterface<CPMeshCut>()->getCutSurfMeshName().empty())
        CHKERR cp.getInterface<CPMeshCut>()->setCutSurfaceFromFile();
      CHKERR cp.getInterface<CPMeshCut>()->copySurface("cutting_surface.vtk");
      Range vol;
      CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
      CHKERR cp.getInterface<CPMeshCut>()->refineMesh(&vol, true, QUIET,
                                                      CrackPropagation::debug);
    }
 
    CHKERR cp.getInterface<CPMeshCut>()->cutRefineAndSplit(
        QUIET, CrackPropagation::debug);
 
    CHKERR moab.write_mesh("cut_mesh_out.h5m");
 
    CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
        cp.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTET,
        "out_level0.vtk", "VTK", "");
    CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
        cp.mapBitLevel["spatial_domain"], BitRefLevel().set(), MBTET,
        "out_level1_tets.vtk", "VTK", "");
    CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
        cp.mapBitLevel["spatial_domain"], BitRefLevel().set(), MBPRISM,
        "out_level1_prisms.vtk", "VTK", "");
 
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}