Trabecular Bone Micro-Finite Element Models
B. van Rietbergen, E. B. Rudnyi, J. G. Korvink
Three-dimensional serial reconstruction techniques allow us to develop very detailed micro-finite element (micro-FE) model of bones that can very accurately represent the porous bone micro-architecture. Fig. 1 sketches the micro finite element analysis . Micro computed tomography (CT) is employed to make 3D high-resolution images (~50 microns) of a bone. Then the 3D reconstruction is directly transformed into an equally shaped micro finite element model by simply converting all bone voxels to equally sized 8-node brick elements. This results in finite element (FE) models with a very large number of elements. Such models can be used, for example to study differences in bone tissue loading between healthy and osteoporotic human bones during quasi static loading .
There is increasing evidence, however, that bone responds in particular to dynamic loads . It has been shown that the application of high-frequency, very low magnitude strains to a bone can prevent bone loss due to osteoporosis and can even result in increased bone strength in bones that are already osteoporotic. In order to better understand this phenomenon, it is necessary to determine the strain as sensed by the bone cells due to this loading. This would be possible with the micro-FE analysis, but then such an analysis need to be a dynamic one.
The present benchmark presents 6 bone models varying in dimension from about 200 thousand to 12 millions equations with the goal to research on scalability of model reduction software. Each model represent a second order system in the form
M d2x/dt2 + K x = B
y = Cx
where the matrices M and K are symmetric and positive definite. The goal of model reduction is to speed up harmonic response analysis in the frequency range 1-100 Hz.
The matrix properties are given in Table 1 below.
|number of elements||20098||192539||278259||606253||1378782||3387547|
|number of nodes||42508||305066||329001||719987||1644848||3989996|
|number of DoFs||127224||914898||986703||2159661||4934244||11969688|
|nnz in half M||1182804||9702186||12437739||27150810||61866069||151251738|
|nnz in half K||3421188||28191660||36326514||79292769||180663963||441785526|
|file||(File 1), 312402 Bytes||(File 2), 2910348 Bytes||(File 3), 4100941 Bytes||(File 4), 8994863 Bytes||(File 5), 20588663 Bytes||(File 6), 50830771 Bytes|
It should be stressed that the first two models have been obtained differently and they are much simpler to deal with than the last four. The connectivity in the last four models is about four times higher. This can be seen by comparing models BS10 and B010. Although models look similar by number of nonzeros in the system matrices, the model B010 is much harder to solve: the number of nonzero elements in the factor for model B010 is about four times more than for BS10.
The method allows for the compact representation of the models, as the element mass and stiffness matrices are the same for all elements. As a result, a file describing the node indices for each element is enough to assemble the global matrix. Each node has three degrees of freedom (UX, UY, UZ) and it contributes three consecutive entries to the state vector. The node numbering is natural from the first to the last. The assembly procedure as a pseudo-code is presented below (indices start from one). It is assumed that the last 300 degrees of freedom are fixed as zero Dirichlet boundary conditions. For simplicity, the pseudo-code does not take into account that the matrix is symmetric.
The data file for each model contains the number of elements, nel, and the number of nodes, nnod, in the first line and then nel lines with eight numbers for node indices in each line.
1)Read the element stiffness matrix elemK(24,24), 8 nodes * 3
degrees of freedom per node.
2)Read the number of elements, nel, and number of node, nnod, from
the first line of the data file.
3)Number of degrees of freedom, ndof = nnod *3 - 300.
4)Allocate space for the sparse global matrix, matK(ndof,
5)Assembly: Do (k = 1, nel) Read eight node numbers from the k-th line, nodeindex(8); Construct the index for degrees of freedom, dofindex(24): Do (i = 1, 8) Do (j = 1, 3) dofindex((i - 1)*3 + j)= (nodeindex(i) - 1)*3 + j; Use dofindex to assemble the element matrix elemK: Do (i = 1, 24) Do (j = 1, 24) If (dofindex(i) < ndof AND dofindex(j) < ndof) matK(dofindex(i), dofindex(j)) += elemK(i, j).
The input matrix contains a single column with B(1) = 1. The output matrix takes first three components of the state vector, that is, three displacements UX, UY and UZ for the first node.
The archive File 7 contains the element mass and stiffness matrices as well as the sample code in C++ to assemble the dynamic system. The code can write the dynamic system in the Matrix Market format or can be used as a hook to transform the global matrices to an appropriate format. The gzipped data files for element assembly as described above can be downloaded from Table 1.
Model reduction for models BS010 and BS10 was performed in . The benchmarking of the parallel MUMPS direct solver  for the stiffness matrices is described in .
 B. van Rietbergen, H. Weinans, R. Huiskes, A. Odgaard, A new method to determine trabecular bone elastic properties and loading using micromechanical finite-elements models. J. Biomechanics, v. 28, N 1, p. 69-81, 1995.
 B. van Rietbergen, R. Huiskes, F. Eckstein, P. Rueegsegger, Trabecular Bone Tissue Strains in the Healthy and Osteoporotic Human Femur, J. Bone Mineral Research, v. 18, N 10, p. 1781-1787, 2003.
 L. E. Lanyon, C. T. Rubin, Static versus dynamic loads as an influence on bone remodelling, Journal of Biomechanics, v. 17, p. 897-906, 1984.
 E. B. Rudnyi, B. van Rietbergen, J. G. Korvink. Efficient Harmonic Simulation of a Trabecular Bone Finite Element Model by means of Model Reduction. 12th Workshop "The Finite Element Method in Biomedical Engineering, Biomechanics and Related Fields", University of Ulm, 20-21 July 2005. Proceedings of the 12th FEM Workshop, p. 61-68. ISBN: 3-9806183-8-2.
 P. R. Amestoy and A. Guermouche and J.-Y. L'Excellent, S. Pralet, Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, v. 32, N 2, pp. 136-156 (2006). http://graal.ens-lyon.fr/MUMPS/
 E. B. Rudnyi, B. van Rietbergen, J. G. Korvink. Model Reduction for High Dimensional Micro-FE Models. TAM'06, The Third HPC-Europa Transnational Access Meeting, Barcelona, 14 - 16 June 2006.