Bone Model

Bone Model

Background
Benchmark ID	boneModelB010 boneModelB025 boneModelB050 boneModelB120 boneModelBS01 boneModelBS10
Category	oberwolfach
System-Class	LTI-SOS
Parameters
nstates	986703 2159661 4934244 11969688 127224 914898
ninputs	1
noutputs	3
nparameters	0
components	B, C, K, M
Copyright
License	GPLv2
Creator	Christian Himpe
Editor	Christian Himpe Kathryn Lund
Location
NA

Description: Trabecular Bone Micro-Finite Element Models

Figure 1

Figure 2

Three-dimensional serial reconstruction techniques allow us to develop a 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 such an analysis needs to be a dynamic one.

The present benchmark presents six bone models varying in dimension from about two hundred thousand to twelve million equations with the goal to research on scalability of model reduction software. Each model represents a second-order system of the form

$\begin{align} M \ddot{x}(t) + K x(t) &= B u(t), \\ y(t) &= C x(t), \end{align}$

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.

Table 1: Bone micro-finite element models.

	BS01	BS10	B010	B025	B050	B120
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	BoneModel-dim1e6-BS01.dat.gz (305.1 kB)	BoneModel-dim1e7-BS10.dat.gz (2.8 MB)	BoneModel-dim1e7-B010.dat.gz (3.9 MB)	BoneModel-dim1e7-B025.dat.gz (8.6 MB)	BoneModel-dim1e7-B050.dat.gz (19.6 MB)	BoneModel-dim1e8-B120.dat.gz (48.5 MB)

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.

Origin

This benchmark is part of the Oberwolfach Benchmark Collection.

Data

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 number of 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, 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 the first three components of the state vector, that is, three displacements $UX$, $UY$, and $UZ$ for the first node.

The archive BoneModel-assemble.tar.gz 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 compressed 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 .

Files (in .mat format) can be found for the three smaller systems (BS01, BS10, and B010) on SuiteSparse.

Dimensions

System structure:

$\begin{align} M \ddot{x}(t) + K x(t) &= B u(t) \\ y(t) &= C x(t) \end{align}$

System dimensions:

$M, K \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^{n \times 1}$, $C \in \mathbb{R}^{3 \times n}$.

System variants:

BS01: $n = 127224$, BS10: $n = 914898$, B010: $n = 986703$, B025: $n = 2159661$, B050: $n = 4934244$, B120: $n = 11969688$.

Citation

To cite this benchmark, use the following references:

• For the benchmark itself and its data:

The MORwiki Community, Bone Model. MORwiki - Model Order Reduction Wiki, 2018. http://modelreduction.org/index.php/Bone_Model

@MISC{morwiki_bone,
  author =       {{The MORwiki Community}},
  title =        {Bone Model},
  howpublished = {{MORwiki} -- Model Order Reduction Wiki},
  url =          {http://modelreduction.org/index.php/Bone_Model},
  year =         {20XX}
}

• For the background on the benchmark:

@ARTICLE{RieWHetal95,
  author =       {B. van Rietbergen and H. Weinans and R. Huiskes and A. Odgaard},
  title =        {A new method to determine trabecular bone elastic properties and loading using micromechanical finite-elements models},
  journal =      {Journal of Biomechanics},
  volume =       {28},
  number =       {1},
  pages =        {69--81},
  year =         {1995},
  doi =          {10.1016/0021-9290(95)80008-5}
}