This module allows for computational simulation of bone functional adaptation by means of a higher order finite element approach based on open system thermodynamics (mainly this paper). Calculations can be enriched with initial bone density aquired from CT scan images. For more instructions, tutorials visit MoFEM page. 
Install bone remodeling module on your MoFEM build
Go to your build (users modules)
cd $MOFEM_INSTALL_DIR/um/users_modules/users_modules/
git clone https://karol41@bitbucket.org/likask/mofem_um_bone_remodelling.git bone_remodelling
cd ..
cmake -DCMAKE_BUILD_TYPE=Release -DWITH_METAIO=1 -DCMAKE_C_FLAGS="-Wall" \
-DCMAKE_CXX_FLAGS="-Wall -Wno-bind-to-temporary-copy -Wno-overloaded-virtual" \
-DCMAKE_EXE_LINKER_FLAGS="-L$MOFEM_INSTALL_DIR/local/lib" users_modules
cd bone_remodelling
make
On Docker cmake command is slightly different:
cmake -DWITH_METAIO=1 -DBUILD_SHARED_LIBS=yes -DCMAKE_BUILD_TYPE=Release users_modules
Bone remodelling, unlike the other modules, requires compilation with MetaIO library for handling CT scan images. Entire installation process can be also found on Youtube video.
Check before run
- Check where is your mesh file
- Check where is your mhd (CT) file
- Check where is your loading history data file
- Check where is your config file (if required)
Density mapping
mpirun -np 4 ./mass_calculation \
-meta_image "file.mhd" -my_file "file.h5m" \
-param_density 1,0 \
-cube_size 3 \
-dist_tol 1 \
-my_order 2 \
-lambda 0 \
-scale 1 \
-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
&& mbconvert out.h5m out.vtk
Bone adaptation, 2D femur example
First, we partition the mesh: (we can skip this step if only one processor is used)
../tools/mofem_part -my_file ./examples/femur_2d.h5m -my_nparts 4
cp out.h5m femur_2d_on_4parts.h5m
Next, we start analysis on a partitioned output mesh file with the following:
mpirun -np 4 ./bone_adaptation \
-my_file femur_2d_on_4parts.h5m -my_order 2 -my_output_prt 1 \
-ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package mumps \
-ksp_atol 1e-12 -ksp_rtol 1e-12 -snes_monitor -snes_type newtonls \
-snes_linesearch_type basic \
-snes_max_it 1000 -snes_atol 1e-6 -snes_rtol 1e-8 \
-ts_type beuler -ts_max_snes_failures 1 -ts_monitor \
-ts_dt 0.1 -ts_final_time 20 -my_load_history ./examples/load_history2.in \
-mass_postproc -ts_monitor -young_modulus 500 -poisson_ratio 0.2 -c 1 \
-rho_ref 1.2 -psi_ref 0.01 -m 3 -n 2 | tee log
Optional parameters:
-analysis_mesh #saves analysis meshes at each time step
-less_post_proc #reduces output files size
-log_view #provides summary and diagnostic information
-my_config #loads configuration file
-with_adolc #calculate the material tangent with automatic differentiation (for testing)
-b 0 -rho_max 2.0 -rho_min 0.3 #parameters for 'bell function'
-ts_adapt_type basic -ts_adapt_basic_safety 0.8 #adaptive time stepping
-ts_adapt_dt_max 5.0 -ts_adapt_dt_min 0.005
-ts_rtol 0.01 -ts_atol 0.01 #adaptive param
For bigger problems (up to 18 cores) we use field split preconditioner with mumps on each sub-problem:
mpirun -np 4 ./bone_adaptation \
-my_file femur_2d_on_4parts.h5m -my_order 6 -my_output_prt 1 \
-snes_monitor -snes_type newtonls \
-snes_linesearch_type basic \
-snes_max_it 100 -snes_atol 1e-6 -snes_rtol -snes_max_linear_solve_fail 1000 1e-8 -ts_type beuler -ts_max_snes_failures 1 \
-ts_dt 0.1 -ts_final_time 20 -my_load_history ./examples/load_history2.in \
-ts_monitor -young_modulus 500 -poisson_ratio 0.2 -c 1 \
-rho_ref 1.2 -psi_ref 0.01 -m 3 -n 2 -mass_postproc\
-ksp_type gmres -pc_type fieldsplit \
-ksp_atol 1e-8 -ksp_rtol 1e-8 \
-ksp_max_it 1000 \
-fieldsplit_0_ksp_type preonly \
-fieldsplit_0_pc_type lu \
-fieldsplit_0_pc_factor_mat_solver_package mumps \
-fieldsplit_1_ksp_type preonly \
-fieldsplit_1_pc_type lu \
-fieldsplit_1_pc_factor_mat_solver_package mumps \
-pc_fieldsplit_type multiplicative -ksp_monitor | tee log
If we want to run analysis for even larger problem (on cluster) block jacobi precoditioners can be utilised:
mpirun -np 64 ./bone_adaptation -my_file out.h5m -my_order 3 -my_output_prt 100 \
-ksp_atol 1e-12 -ksp_rtol 1e-12 -snes_monitor -snes_type newtonls -snes_linesearch_type basic -snes_max_it 30 -snes_atol 1e-6 -snes_rtol 1e-8 \
-ts_type beuler -ts_max_snes_failures 1 -snes_max_linear_solve_fail 1000 -ts_dt 0.5 -ts_final_time 5 \
-my_load_history load_history3.in -mass_postproc -ts_monitor -less_post_proc \
-young_modulus 900 -poisson_ratio 0.3 -c 1 -rho_ref 1 -psi_ref 0.0025 -m 3 -n 2 -less_post_proc -equilibrium_stop_rate 0.001 -log_view \
-ksp_type fgmres -pc_type fieldsplit \
-ksp_atol 1e-8 -ksp_rtol 1e-8 \
-ksp_max_it 1000 \
-ksp_monitor \
-fieldsplit_0_ksp_type gmres \
-fieldsplit_0_ksp_max_it 100 \
-fieldsplit_0_pc_type bjacobi \
-fieldsplit_0_sub_pc_type ilu \
-fieldsplit_0_sub_pc_factor_levels 2 \
-fieldsplit_1_ksp_type gmres \
-fieldsplit_1_ksp_max_it 100 \
-fieldsplit_1_pc_type bjacobi \
-fieldsplit_1_sub_pc_type ilu \
-fieldsplit_1_sub_pc_factor_levels 2 \
-pc_fieldsplit_type multiplicative
Do output VTK script
Produce mass, energy graphs in gluplot
grep "Elastic" log | awk '{print $2,$4,$7}' | tee log_upd
gnuplot
plot "log_upd" using 1:2 title 'mass' with linespoints
plot "log_upd" using 1:3 title 'psi' with linespoints
Acknowledgements
This user module has been co-developed by the School of Engineering and the Weipers Centre Equine Hospital, School of Veterinary Medicine at the University of Glasgow. The authors are particularly grateful to the staff of the Small Animal Hospital, where the CT scans all required bones were performed. Kelvin Smith PhD Scholarship is acknowledged for providing financial support for this research.
