memcheck_leak.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>
 
#define BOOST_UBLAS_NDEBUG 1
 
#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);
 
  std::vector<boost::weak_ptr<RefEntity>> v_ref_entity;
  std::vector<boost::weak_ptr<RefElement>> v_ref_element;
  std::vector<boost::weak_ptr<Field>> v_field;
  boost::weak_ptr<DofEntity> arc_dof;
  boost::weak_ptr<MWLSApprox> mwls_approx;
  std::vector<boost::weak_ptr<FEMethod>> vec_fe;
  std::vector<boost::weak_ptr<NonlinearElasticElement>> v_elastic_ele_str;
 
  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
        MeshsetsManager::brodcastMeshsets = false;
    MoFEM::Core core(moab, PETSC_COMM_WORLD, VERBOSE);
    MoFEM::Interface &m_field = core;
 
    auto core_log = logging::core::get();
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "MEMCHECK"));
    LogManager::setLog("MEMCHECK");
    MOFEM_LOG_TAG("MEMCHECK", "memcheck");
 
    MOFEM_LOG_FUNCTION();
    auto set_log_attr = [&](auto c) {
      MOFEM_LOG_ATTRIBUTES(c, LogManager::BitLineID | LogManager::BitScope);
    };
    set_log_attr("MEMCHECK");
    MOFEM_LOG_C("MEMCHECK", Sev::inform, "### Input parameter: -my_file %s",
                mesh_file_name);
 
    // Create data structure for crack propagation
    CrackPropagation cp(m_field, 3, 2);
 
    CHKERR MoFEM::Interface::getFileVersion(moab, cp.fileVersion);
    MOFEM_LOG("MEMCHECK", Sev::inform)
        << "File version " << cp.fileVersion.strVersion();
 
    // Configure meshes
    CHKERR cp.getOptions();
    // Register MOFEM discrete data managers
    {
      DMType dm_name = "MOFEM";
      CHKERR DMRegister_MoFEM(dm_name);
    }
 
    auto add_bits = [&]() {
      MoFEMFunctionBegin;
      CHKERR m_field.getInterface<BitRefManager>()
          ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
                                           BitRefLevel().set());
      MoFEMFunctionReturn(0);
    };
    CHKERR add_bits();
 
    std::vector<int> surface_ids;
    surface_ids.push_back(cp.getInterface<CPMeshCut>()->getSkinOfTheBodyId());
    std::vector<std::string> fe_surf_proj_list;
    fe_surf_proj_list.push_back("SURFACE");
 
    // create bit levels, fields and finite elements structures in database
    if (cp.doCutMesh) {
 
      if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
 
        // restart code from file
        CHKERR m_field.build_field("LAMBDA_ARC_LENGTH");
        CHKERR m_field.build_finite_elements("ARC_LENGTH");
 
        Range ents;
        CHKERR m_field.getInterface<BitRefManager>()
            ->getAllEntitiesNotInDatabase(ents);
        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);
 
        auto first_number_string = [](std::string const &str) {
          std::size_t const n = str.find_first_of("0123456789");
          if (n != std::string::npos) {
            std::size_t const m = str.find_first_not_of("0123456789", n);
            return str.substr(n, m != std::string::npos ? m - n : m);
          }
          return std::string();
        };
        std::string str_number =
            first_number_string(std::string(mesh_file_name));
 
        if (!str_number.empty())
          cp.startStep = atoi(str_number.c_str()) + 1;
 
        string cutting_surface_name =
            "cutting_surface_" +
            boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
 
        if (m_field.get_comm_rank() == 0)
          CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
              cp.crackAccelerationFactor, cutting_surface_name, VERBOSE,
              CrackPropagation::debug);
        else
          CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
              cp.crackAccelerationFactor, "", NOISY, CrackPropagation::debug);
 
        CHKERR cp.getInterface<CPMeshCut>()->refineMesh(
            nullptr, false, QUIET, CrackPropagation::debug);
 
      } else {
 
        string cutting_surface_name =
            "cutting_surface_" +
            boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
 
        if (!cp.getInterface<CPMeshCut>()->getCutSurfMeshName().empty())
          CHKERR cp.getInterface<CPMeshCut>()->setCutSurfaceFromFile();
        CHKERR cp.getInterface<CPMeshCut>()->copySurface(cutting_surface_name);
        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);
 
    } else {
 
      Range vol;
      CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
      BitRefLevel bit0 = cp.mapBitLevel["mesh_cut"];
      CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
          vol, bit0.set(BITREFLEVEL_SIZE - 1), true, CrackPropagation::debug);
      CHKERR cp.getInterface<CPMeshCut>()->refineAndSplit(
          QUIET, CrackPropagation::debug);
    }
 
    CHKERR cp.createProblemDataStructures(surface_ids, QUIET,
                                          CrackPropagation::debug);
 
    if (cp.propagateCrack == PETSC_TRUE) {
 
      // set inital coordinates
      CHKERR cp.setMaterialPositionFromCoords();
      CHKERR cp.setSpatialPositionFromCoords();
      if (cp.addSingularity == PETSC_TRUE)
        CHKERR cp.setSingularDofs("SPATIAL_POSITION");
 
      SmartPetscObj<DM> dm_elastic;
      SmartPetscObj<DM> dm_eigen_elastic;
      SmartPetscObj<DM> dm_material;
      SmartPetscObj<DM> dm_crack_propagation;
      SmartPetscObj<DM> dm_material_forces;
      SmartPetscObj<DM> dm_surface_projection;
      SmartPetscObj<DM> dm_crack_srf_area;
      CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
                          dm_crack_propagation, dm_material_forces,
                          dm_surface_projection, dm_crack_srf_area, surface_ids,
                          fe_surf_proj_list);
      MOFEM_LOG_ATTRIBUTES("DMSELF", LogManager::BitScope);
      MOFEM_LOG_ATTRIBUTES("DMWORLD", LogManager::BitScope);
 
      CHKERR cp.setSpatialPositionFromCoords();
      if (cp.addSingularity == PETSC_TRUE) {
        CHKERR cp.setSingularDofs("SPATIAL_POSITION");
      }
      if (string(mesh_file_name).compare(0, 8, "restart_") == 0)
        cp.getLoadScale() = cp.getArcCtx()->getFieldData();
 
      if (cp.solveEigenStressProblem) {
        MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve Eigen ELASTIC problem";
        CHKERR cp.getInterface<CPSolvers>()->solveElastic(
            dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
        CHKERR cp.postProcessDM(dm_eigen_elastic, -2, "ELASTIC", true);
      }
 
      MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve ELASTIC problem";
      CHKERR cp.getInterface<CPSolvers>()->solveElastic(
          dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
      CHKERR cp.postProcessDM(dm_elastic, -1, "ELASTIC", true);
 
      // set finite element instances and user data operators on instances
      CHKERR cp.assembleElasticDM(VERBOSE, false);
      CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
      CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
                                         false);
      CHKERR cp.assembleCouplingForcesDM(PETSC_NULL, VERBOSE, false);
 
      CHKERR cp.calculateElasticEnergy(dm_elastic);
      CHKERR cp.getInterface<CPSolvers>()->calculateGc(
          dm_material_forces, dm_surface_projection, dm_crack_srf_area,
          surface_ids);
 
      if (cp.doCutMesh == PETSC_TRUE) {
        CHKERR cp.getInterface<CPSolvers>()->solveCutMesh(
            dm_crack_propagation, dm_elastic, dm_eigen_elastic, dm_material,
            dm_material_forces, dm_surface_projection, dm_crack_srf_area,
            surface_ids, CrackPropagation::debug);
      } else {
        CHKERR cp.getInterface<CPSolvers>()->solvePropagation(
            dm_crack_propagation, dm_elastic, dm_material, dm_material_forces,
            dm_surface_projection, dm_crack_srf_area, surface_ids);
      }
 
    } else {
 
      // set inital coordinates
      CHKERR cp.setMaterialPositionFromCoords();
      CHKERR cp.setSpatialPositionFromCoords();
      if (cp.addSingularity == PETSC_TRUE)
        CHKERR cp.setSingularDofs("SPATIAL_POSITION");
 
      SmartPetscObj<DM> dm_elastic;
      SmartPetscObj<DM> dm_eigen_elastic;
      SmartPetscObj<DM> dm_material;
      SmartPetscObj<DM> dm_crack_propagation;
      SmartPetscObj<DM> dm_material_forces;
      SmartPetscObj<DM> dm_surface_projection;
      SmartPetscObj<DM> dm_crack_srf_area;
      CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
                          dm_crack_propagation, dm_material_forces,
                          dm_surface_projection, dm_crack_srf_area, surface_ids,
                          fe_surf_proj_list);
 
      SmartPetscObj<DM> dm_bc;
      if (m_field.getInterface<MeshsetsManager>()->checkMeshset(
              "ANALITICAL_DISP"))
        CHKERR cp.createBcDM(dm_bc, "BC_PROBLEM",
                             cp.mapBitLevel["spatial_domain"],
                             BitRefLevel().set());
 
      if (cp.solveEigenStressProblem) {
        MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve Eigen ELASTIC problem";
        CHKERR cp.getInterface<CPSolvers>()->solveElastic(
            dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
        CHKERR cp.postProcessDM(dm_eigen_elastic, -2, "ELASTIC", true);
      }
 
      // set finite element instances and user data operators on instances
      MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve ELASTIC problem";
      CHKERR cp.getInterface<CPSolvers>()->solveElastic(
          dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
 
      // set finite element instances and user data operators on instances
      CHKERR cp.assembleElasticDM(VERBOSE, false);
      CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
      if (cp.doCutMesh)
        CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
                                           false);
 
      CHKERR cp.getInterface<CPSolvers>()->calculateGc(
          dm_material_forces, dm_surface_projection, dm_crack_srf_area,
          surface_ids);
      CHKERR cp.calculateElasticEnergy(dm_elastic);
 
      // post-processing
      CHKERR cp.postProcessDM(dm_elastic, 0, "ELASTIC", true);
      // set vector values from field data
      CHKERR cp.savePositionsOnCrackFrontDM(dm_elastic, PETSC_NULL, 2, false);
      CHKERR cp.writeCrackFont(cp.mapBitLevel["spatial_domain"], 0);
 
      // Run tests
      CHKERR cp.tetsingReleaseEnergyCalculation();
 
      auto refined_ents_ptr = m_field.get_ref_ents();
      v_ref_entity.reserve(refined_ents_ptr->size());
      for (auto e : *refined_ents_ptr)
        v_ref_entity.emplace_back(e);
 
      auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
      v_ref_element.reserve(refined_finite_elements_ptr->size());
      for (auto e : *refined_finite_elements_ptr)
        v_ref_element.emplace_back(e);
 
      auto fields_ptr = m_field.get_fields();
      v_field.reserve(fields_ptr->size());
      for (auto e : *fields_ptr)
        v_field.emplace_back(e);
 
      arc_dof = cp.arcLengthDof;
      mwls_approx = cp.mwlsApprox;
      v_elastic_ele_str.emplace_back(cp.elasticFe);
      v_elastic_ele_str.emplace_back(cp.materialFe);
      vec_fe.emplace_back(cp.feLhs);
      vec_fe.emplace_back(cp.feRhs);
      vec_fe.emplace_back(cp.feMaterialRhs);
      vec_fe.emplace_back(cp.feMaterialLhs);
      vec_fe.emplace_back(cp.feEnergy);
      vec_fe.emplace_back(cp.feCouplingElasticLhs);
      vec_fe.emplace_back(cp.feCouplingMaterialLhs);
      vec_fe.emplace_back(cp.feSmootherRhs);
      vec_fe.emplace_back(cp.feSmootherLhs);
      vec_fe.emplace_back(cp.feGriffithForceRhs);
      vec_fe.emplace_back(cp.feGriffithConstrainsDelta);
      vec_fe.emplace_back(cp.feGriffithConstrainsRhs);
      vec_fe.emplace_back(cp.feGriffithConstrainsLhs);
      vec_fe.emplace_back(cp.feSpringLhsPtr);
      vec_fe.emplace_back(cp.feSpringRhsPtr);
      vec_fe.emplace_back(cp.feRhsSimpleContact);
      vec_fe.emplace_back(cp.feLhsSimpleContact);
      vec_fe.emplace_back(cp.feMaterialAnaliticalTraction);
    }
  }
  CATCH_ERRORS;
 
  if (auto a = arc_dof.lock()) {
    MOFEM_LOG("MEMCHECK", Sev::error)
        << "cp.arcLengthDof.use_count() " << a.use_count();
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
  }
 
  for (auto v : vec_fe)
    if (auto a = v.lock()) {
      MOFEM_LOG("MEMCHECK", Sev::error) << "fe " << a.use_count();
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
    }
 
  for (auto v : v_elastic_ele_str)
    if (auto a = v.lock()) {
      MOFEM_LOG("MEMCHECK", Sev::error)
          << "elastic structure " << a.use_count();
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
    }
 
  if (auto a = mwls_approx.lock()) {
    MOFEM_LOG("MEMCHECK", Sev::error)
        << "cp.mwlsApprox.use_count() " << a.use_count();
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
  }
 
  for (auto e : v_ref_entity) {
    if (auto a = e.lock()) {
      MOFEM_LOG("MEMCHECK", Sev::error) << "v_ref_entity " << a.use_count();
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
    }
  }
 
  for (auto e : v_ref_element) {
    if (auto a = e.lock()) {
      MOFEM_LOG("MEMCHECK", Sev::error) << "v_ref_element " << a.use_count();
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
    }
  }
 
  for (auto e : v_field) {
    if (auto a = e.lock()) {
      MOFEM_LOG("MEMCHECK", Sev::error) << "v_field " << a.use_count();
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
    }
  }
 
  MoFEM::Core::Finalize();
  return 0;
}