bone_adaptation.cpp

Main program for running bone adaptation

/**
 * \file bone_adaptation.cpp
 * \example bone_adaptation.cpp
 * \ingroup bone_remodelling
 * 
 * \brief Main program for running bone adaptation
 *
 */

/* 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 <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;

#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <Remodeling.hpp>

using namespace BoneRemodeling;

static char help[] = "\n";

// #ifndef WITH_METAIO
//   #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
// #endif

int main(int argc, char *argv[]) {

  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);

  try {

    char mesh_file_name[255];
    PetscBool flg_file = PETSC_FALSE;
    PetscBool flg_elastic_only = PETSC_FALSE;

    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Main set up", "none");
    CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_file);

    CHKERR PetscOptionsBool(
        "-elastic_only", "is used for testing, calculates only elastic problem",
        "", PETSC_FALSE, &flg_elastic_only, PETSC_NULL);

    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);

    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;

    // Read parameters from line command
    if (flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
    }
    const char *option;

    // option = "";
    option = "PARALLEL=READ_PART;"
             "PARALLEL_RESOLVE_SHARED_ENTS;"
             "PARTITION=PARALLEL_PARTITION;";
    CHKERR moab.load_file(mesh_file_name, 0, option);
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);

    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;

    // Add meshsets with material and boundary conditions
    MeshsetsManager *meshsets_manager_ptr;
    CHKERR m_field.getInterface<MeshsetsManager>(meshsets_manager_ptr);
    CHKERR meshsets_manager_ptr->setMeshsetFromFile();

    PetscPrintf(PETSC_COMM_WORLD, "Read meshset. Added meshsets for bc.cfg\n");
    for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
      PetscPrintf(PETSC_COMM_WORLD, "%s",
                  static_cast<std::ostringstream &>(
                      std::ostringstream().seekp(0) << *it << endl)
                      .str()
                      .c_str());
      cerr << *it << endl;
    }

    Remodeling::CommonData common_data;
    common_data.feRhs =
        boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
    common_data.feLhs =
        boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
    common_data.preProcRhs = boost::shared_ptr<Remodeling::FePrePostProcessRhs>(
        new Remodeling::FePrePostProcessRhs());
    common_data.preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
        new Remodeling::FePrePostProcessLhs());

    CHKERR common_data.getParameters();

    Remodeling remodelling(m_field, common_data);

    // Remodeling
    if (!flg_elastic_only) {
      remodelling.getParameters();
      remodelling.addFields();
      remodelling.addMomentumFluxes();
      remodelling.addElements();
      remodelling.buildDM();

      bool flag_cubit_disp = false;
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
               m_field, NODESET | DISPLACEMENTSET, it)) {
        flag_cubit_disp = true;
      }

      boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
      // if normally defined boundary conditions are not found, try to use
      // DISPLACEMENT blockset
      if (!flag_cubit_disp) {
        dirihlet_bc_ptr =
            boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
                m_field, "DISPLACEMENTS", "DISPLACEMENT"));
        dynamic_cast<DirichletSetFieldFromBlockWithFlags &>(*dirihlet_bc_ptr)
            .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
      } else {
        dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
            new DirichletDisplacementBc(m_field, "DISPLACEMENTS"));
        dynamic_cast<DirichletDisplacementBc &>(*dirihlet_bc_ptr)
            .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
      }

      // Add elements Time Solver via DM interface.
      // This part is to calculate residual vector.
      CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
                                   common_data.preProcRhs.get(), NULL);
      CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
                                   dirihlet_bc_ptr.get(), NULL);
      CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING",
                                   common_data.feRhs.get(), NULL, NULL);
      if (common_data.nOremodellingBlock) {
        CHKERR DMMoFEMTSSetIFunction(common_data.dm, "ELASTIC",
                                     &(common_data.elasticPtr->getLoopFeRhs()),
                                     NULL, NULL);
      }

      // Add surface forces
      for (boost::ptr_map<string, NeummanForcesSurface>::iterator fit =
               common_data.neumannForces.begin();
           fit != common_data.neumannForces.end(); fit++) {
        CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
                                     &fit->second->getLoopFe(), NULL, NULL);
      }

      // Add edge forces
      for (boost::ptr_map<std::string, EdgeForce>::iterator fit =
               common_data.edgeForces.begin();
           fit != common_data.edgeForces.end(); fit++) {
        CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
                                     &fit->second->getLoopFe(), NULL, NULL);
      }

      // Add nodal forces
      for (boost::ptr_map<std::string, NodalForce>::iterator fit =
               common_data.nodalForces.begin();
           fit != common_data.nodalForces.end(); fit++) {
        CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
                                     &fit->second->getLoopFe(), NULL, NULL);
      }

      CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
                                   dirihlet_bc_ptr.get());
      CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
                                   common_data.preProcRhs.get());

      // Add elements Time Solver via DM interface.
      // This part is to calculate tangent matrix.
      CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
                                   common_data.preProcLhs.get(), NULL);
      CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
                                   dirihlet_bc_ptr.get(), NULL);
      CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING",
                                   common_data.feLhs.get(), NULL, NULL);
      if (common_data.nOremodellingBlock) {
        CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "ELASTIC",
                                     &common_data.elasticPtr->getLoopFeLhs(),
                                     NULL, NULL);
      }
      CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
                                   dirihlet_bc_ptr.get());
      CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
                                   common_data.preProcLhs.get());

      CHKERR DMCreateMatrix(common_data.dm, &common_data.A);
      CHKERR DMCreateGlobalVector(common_data.dm, &common_data.D);
      CHKERR DMCreateGlobalVector(common_data.dm, &common_data.F);

      // Initialize vectors with initial data from mesh
      CHKERR DMoFEMMeshToLocalVector(common_data.dm, common_data.D,
                                     INSERT_VALUES, SCATTER_FORWARD);

      CHKERR TSCreate(PETSC_COMM_WORLD, &common_data.ts);
      CHKERR TSSetType(common_data.ts, TSBEULER);

      CHKERR remodelling.solveDM();

      CHKERR MatDestroy(&common_data.A);
      CHKERR VecDestroy(&common_data.D);
      CHKERR VecDestroy(&common_data.F);
      CHKERR TSDestroy(&common_data.ts);
    }

  }
  CATCH_ERRORS;

  MoFEM::Core::Finalize();
  PetscFunctionReturn(0);
}