Cartesian grid FEM for direct creation of patient specific models and implant simulation

Autores UPV
CONGRESO Cartesian grid FEM for direct creation of patient specific models and implant simulation


Patient speci c simulation usually requires numerical models obtained from medical image data. Two common approaches can be distinguished. The rst one is based on the creation of volumes by interpolating the bitmaps with closed geometrical surfaces before a Finite Element (FE) mesh can be created. A great segmentation effort is required to obtain realistic geometries. This is typically hard to automatize and usually consumes a great amount of specialized human work. This method is the one commonly used in pre-operative prosthesis simulations. Mixing the geometrical model from the image and the CAD of the prosthetic device causes additional cumbersomeness and hence additional human work, especially in the process of meshing the resulting complex geometries. The second approach directly assigns one hexahedral element to each pixel. The resulting mesh is structured, all its elements have the same good aspect ratio and dimension, allowing efficient information treatment, but it usually involves a high number of degrees of freedom which results in high computational cost of the numerical analysis. In this paper a procedure based on the use of Cartesian grid Finite Element Method (cg-FEM) is proposed for the direct creation of patient speci c FE models from medical images. This procedure simultaneously lightens the segmentation stage and reduces the computational cost, providing a tool to easily include geometries of prosthetic devices in the image-based FE model for efficient analysis of interaction between bone and prosthesis in a preoperative framework. The medical image is directly immersed in a structure of nested hierarchical Cartesian grids so that an integer number of pixels is contained in the elements of each level. To take into account the heterogeneity of biological tissue, the scale the pixel values are expressed in is related with corresponding material properties. All the pixels contained in each element are used for the calculation of its stiffness matrix. For this purpose a special quadrature rule based on the Reimann sum is used in which integration points are located at the centre of each pixel, as a result the elastic behaviour is homogenized in each element. To enhance accuracy in areas where the pixel value eld varies steeply, its distribution is evaluated in each element and used to guide an h-adaptive process. This allows the mesh to be ner near the body contour implicitly de ned in the bitmap. Multi-point constraints are used to enforce C0 continuity in the non conforming mesh obtained from the h-adaptive re nement. The introduction of geometrically de ned prostheses in the model is carried out by directly introducing their CAD into the Cartesian grid hierarchical structure and intersecting the mesh with the geometry. The mesh elements can be divided in three classes, characterized by different procedures for their stiffness matrix calculation, Ke. In the elements completely inside the geometry all the pixels are deactivated and a simple Gauss quadrature is used. In a 2D h-adapted cg-FEM all the inner elements share the same Ke and, as a result, it can be calculated only once. In the elements completely on the image the Riemann sum-based quadrature rule is used. Part of the computational efficiency is kept because the elements belonging to the same level share the same B matrix of the shape function derivatives. The elements on the contour are partially on the geometry and on the bitmap. Only the pixels whose center lays outside the geometry are active and integrated according to the quadrature adopted for the image while the part on the geometry is divided in triangular integration sub-domains for which the proper Gauss integration rule is used. The Dirichlet boundary conditions are enforced by using Lagrange multipliers on curves belonging to the prosthesis geometrical model or created ad hoc for this purpose. Convergence studies and numerical examples are presented and discussed for the 2D case.