Calculation of the muscle-induced forces on the finite element mesh

In this section, muscle forces are converted to finite-element forces acting on the nodes of the surface mesh of each muscle group. The 3 loading models proposed by Grosse et al. (2007) have been implemented in python with a special care about efficiency.

In the simplest model – the Uniform-Traction model – the total traction force \bar{F} applied by a muscle on the bone is uniformly distributed on the surface mesh of the muscle attachment. For each triangle of this mesh, the traction is assumed to be directed towards a focal node modelling the opposite muscle attachment (see Figure 1).

*Figure 1: Computation of the nodal forces acting on a triangle of the surface mesh of a muscle group. In the simplest case (*Uniform-Traction model), the forces are aligned towards the force focal node \boldsymbol{o} representing the opposite muscle attachment.

If \boldsymbol{o} is the focal node of a given muscle group on the opposite side of the mesh, the force direction \boldsymbol{d}^i related to the i^{\text{th}} triangle corresponds to the unit vector

\boldsymbol{d}^i = \frac{\boldsymbol{o}-\boldsymbol{c}^i}{||\boldsymbol{o}-\boldsymbol{c}^i||} \qquad (i=1,2,...,n_{\text{tri}})

where n_{\text{tri}} is the total number of triangles and \boldsymbol{c}^i is the barycentre of triangle i.

The nodal force corresponding to a unit traction \boldsymbol{d}^i acting on the j^{\text{th}} node of the i^{\text{th}} triangle is then computed as

\boldsymbol{f}_j^i=\frac{A_i}{3}\boldsymbol{d}^i\qquad(j=1,2,3)

where A_i is the area of triangle i.

The forces of all the individual triangles are then assembled into a structural vector \boldsymbol{F}^* whose components are the sum of the contributions of the neighbouring triangles of each node:

\boldsymbol{F}^* = \sum_{i,j}\boldsymbol{f}_j^i

Eventually, this load vector is scaled so that its magnitude corresponds to the prescribed value \bar{F}:

\boldsymbol{F}=\bar{F}\frac{\boldsymbol{F}^*}{||\boldsymbol{F}^*||}

The 2 other methods follow the same procedure except that the direction \boldsymbol{d}^i is computed in a more realistic way. Indeed, if a triangle does not have a direct line of sight to the focal point (Figure 2, right), the direction calculated by the latter method gives a compressive force towards the interior of the bone, although the muscle fibres wrap the surface of the bone. The forces should thus be mainly oriented along a tangent to the surface. This remark leads to a second model of force distribution – the Tangential-Traction model – which is implemented as follows. Let us define \boldsymbol{n}, the outward unit normal to a given triangle of the surface mesh. If \boldsymbol{d}.\boldsymbol{n}>0 (as in Figure 2, left), the previous direction \boldsymbol{d}^i is kept and the forces acting on the triangle are the same as before. However, if \boldsymbol{d}.\boldsymbol{n}<0 (Figure 2, right), the direction is projected onto the plane of the triangle and normalised again. The corrected direction \boldsymbol{d}_T becomes (note that the index i will be omitted in the following expressions to simplify the notations):

\boldsymbol{d}_T = \frac{\boldsymbol{d}_U-(\boldsymbol{d}_U.\boldsymbol{n})\boldsymbol{n}}{||\boldsymbol{d}_U-(\boldsymbol{d}_U.\boldsymbol{n})\boldsymbol{n}||}

where \boldsymbol{d}_U is the direction computed as in the Uniform-Traction model. The new direction \boldsymbol{d}_T is much more realistic and is not expensive in terms of additional calculations.

Figure 2: Projection of the force direction \boldsymbol{d}_U when the triangle does not have a direct line of sight to the focal point. The new direction \boldsymbol{d}_T is tangent to the surface mesh in the plane defined by the focal point \boldsymbol{o} and the normal vector \boldsymbol{n}.

Figure 3: Calculation of the approximation \tilde{\boldsymbol{s} of fibre length \boldsymbol{s} from the origin of the muscle fibre bundle to the centre of the current triangle and the local radius R of curvature at the same point. These 2 lengths must be measured in the plane defined by the focal point \boldsymbol{o} and the normal vector \boldsymbol{n}.

The third model – the Tangent-Plus-Normal model – starts from the previous method and adds a normal component to the tangential contribution which models the compressive action of the muscle fibres that wrap around the bone on top of each other. This additional pressure acts only if the surface is locally convex, and its amplitude is proportional to the length s of the muscle fibre bundle in the plane defined by the focal point \boldsymbol{o} and the normal vector \boldsymbol{n} as schematised in Figure 3. The new direction \boldsymbol{d}_{T+N} including both tangential and normal components is evaluated as

\boldsymbol{d}_{T+N} = \boldsymbol{d}_{T}+\frac{s}{R}\boldsymbol{n}\qquad\text{when } R<0

where R is the local radius of curvature in the plane defined by \boldsymbol{o} and \boldsymbol{n}. Note that \boldsymbol{d}_{T+N} is not a unit vector anymore.

Contrary to the Tangential-Traction model, the evaluation of this expression is extremely costly since it requires to cut the surface mesh by a plane. The cutting operation has been implemented with the help of the VTK library[1] which is commonly used for displaying meshes on computer screens. VTK also comes with a series of very robust algorithms for modifying and manipulating meshes. These algorithms are very efficient since they target real-time visualisation.

For each triangle, the surface mesh of the muscle attachment is first cut by a plane normal to \boldsymbol{d}_{T} and passing through the barycentre \boldsymbol{c} of the triangle. This operation divides the mesh into a front and a back part (see Figure 3) which are needed to identify the part containing the origin of the fibre bundle. Next, the intersection of both parts with the plane defined by \boldsymbol{o} and \boldsymbol{n} is calculated. The result of this operation is 2 sets of piecewise-linear segments which correspond to the mesh displayed in Figure 3. The vertices are sorted and the fibre length \boldsymbol{s} is approximated by \tilde{\boldsymbol{s}}, the sum of the lengths of the line segments coming from the back part of the surface mesh (the blue segments in Figure 3). Then, a smooth high-order polynomial approximation p(t) of the bone surface shape is built in the local axes (\boldsymbol{d}_T, \boldsymbol{n}) from the position of the first 2 vertices found on both sides of the current triangle. The curvature of this polynomial at the barycentre \boldsymbol{c} of the triangle can be easily calculated by

\frac{1}{R} = \frac{p''(0)}{(1+p'(0)^2)^{3/2}}

Where p' and p'' denote the first and second derivatives of p(t).

Once the \boldsymbol{d}_{T+N} vectors have been calculated for every triangle of the muscle attachment surface, they are assembled and scaled as in the Uniform-Traction model. Thanks to VTK, the total CPU time required to create this much more realistic distribution of forces remains reasonable even for very fine meshes (a few seconds).

[1] The Visualization Toolkit : https://vtk.org/