Main program for running bone adaptation.
#include <boost/program_options.hpp>
namespace po = boost::program_options;
static char help[] =
"\n";
int main(
int argc,
char *argv[]) {
try {
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",
"-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;
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
MeshsetsManager *meshsets_manager_ptr;
CHKERR meshsets_manager_ptr->setMeshsetFromFile();
PetscPrintf(PETSC_COMM_WORLD, "Read meshset. Added meshsets for bc.cfg\n");
PetscPrintf(PETSC_COMM_WORLD, "%s",
static_cast<std::ostringstream &>(
std::ostringstream().seekp(0) << *it << endl)
.str()
.c_str());
cerr << *it << endl;
}
common_data.
preProcRhs = boost::shared_ptr<Remodeling::FePrePostProcessRhs>(
common_data.
preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
if (!flg_elastic_only) {
bool flag_cubit_disp = false;
flag_cubit_disp = true;
}
boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
if (!flag_cubit_disp) {
dirihlet_bc_ptr =
m_field, "DISPLACEMENTS", "DISPLACEMENT"));
} else {
dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
}
CHKERR DMMoFEMTSSetIFunction(common_data.
dm,
"FE_REMODELLING", 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);
CHKERR DMMoFEMTSSetIFunction(common_data.
dm,
"ELASTIC",
NULL, NULL);
}
for (boost::ptr_map<string, NeummanForcesSurface>::iterator fit =
CHKERR DMMoFEMTSSetIFunction(common_data.
dm, fit->first.c_str(),
&fit->second->getLoopFe(), NULL, NULL);
}
for (boost::ptr_map<std::string, EdgeForce>::iterator fit =
CHKERR DMMoFEMTSSetIFunction(common_data.
dm, fit->first.c_str(),
&fit->second->getLoopFe(), NULL, NULL);
}
for (boost::ptr_map<std::string, NodalForce>::iterator 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,
CHKERR DMMoFEMTSSetIJacobian(common_data.
dm,
"FE_REMODELLING", 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);
CHKERR DMMoFEMTSSetIJacobian(common_data.
dm,
"ELASTIC",
NULL, NULL);
}
CHKERR DMMoFEMTSSetIJacobian(common_data.
dm,
"FE_REMODELLING", NULL, NULL,
dirihlet_bc_ptr.get());
CHKERR DMMoFEMTSSetIJacobian(common_data.
dm,
"FE_REMODELLING", NULL, NULL,
CHKERR DMCreateMatrix(common_data.
dm, &common_data.
A);
CHKERR DMCreateGlobalVector(common_data.
dm, &common_data.
D);
CHKERR DMCreateGlobalVector(common_data.
dm, &common_data.
F);
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);
}
}
return 0;
}
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERR
Inline error check.
#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.
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
MoFEMErrorCode getParameters()
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
DM dm
Discretization manager.
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Implementation of bone remodeling finite element.
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.
boost::ptr_vector< MethodForForceScaling > methodsOp
Add boundary conditions form block set having 6 attributes.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.