v0.8.23
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 <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();
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;
}
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>(
common_data.preProcLhs = boost::shared_ptr<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;
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);
}
}
PetscFunctionReturn(0);
}