v0.14.0
Loading...
Searching...
No Matches
bone_adaptation.cpp

Main program for running bone adaptation.

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();
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;
}
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);
}
}
return 0;
}
Elastic materials.
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ NODESET
Definition: definitions.h:146
@ DISPLACEMENTSET
Definition: definitions.h:150
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
char mesh_file_name[255]
Definition: enable_if.hpp:6
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
Definition: Remodeling.hpp:132
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
Definition: Remodeling.hpp:133
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
Volume finite element.
Definition: Remodeling.hpp:41
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:79
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:53
Implementation of bone remodeling finite element.
Definition: Remodeling.hpp:36
MoFEMErrorCode addElements()
Set and add finite elements.
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode buildDM()
Set problem and DM.
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:76
Add boundary conditions form block set having 6 attributes.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.