bone_adaptation.cpp File Reference

#include <boost/program_options.hpp>
#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <Remodeling.hpp>

Go to the source code of this file.

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

Variables
static char	help [] = "\n"

Function Documentation

main()

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

Examples

bone_adaptation.cpp.

Definition at line 41 of file bone_adaptation.cpp.

                                 {
 
  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();
  return 0;
}

Variable Documentation

help

char help[] = "\n"

static

Examples

bone_adaptation.cpp.

Definition at line 35 of file bone_adaptation.cpp.