Inserting needles into tissues is an often performed clinical procedure. Tissue deformations influence the outcome of an insertion, and therefore such insertions are hard to do. Training with computer simulations could improve accuracy of such insertions.

DiMaio and Salcudean (2002) show needle insertion in 2D soft material, simulated with the FEM on a semi-regular grid. By using precomputations, they are able to reach update rates of 500 Hz. The simulation drives a force-feedback device, and allows interactive manipulation of the soft tissue. The model is validated by measuring 2D deformations of a block of tissue during an automated needle insertion.

Alterovitz et.al. (2003) show needle insertion in 2D, on a regular grid using FEM. By using a dynamic method, they are able to reach update rates of 50 Hz (1250 elements.) Both approaches are 2D, and DiMaio's is fundamentally linear, due to the precomputations involved. In this work we explore to what extent needle insertion can be generalized to 3 dimensions and more complex material models.

In 2D, the needle is a line segment (or curve). Since line segments are
*n-1* dimensional in 2D, applying forces and setting boundary
conditions on them is well-defined, and yields an approximation to a
mechanical problem.

For mechanically valid problem, forces in 3D should be applied to surfaces. Therefore, the needle must be modeled as an object with a surface. For an object of 0.1 m in diameter, and a needle of 0.001 m, this implies that the elements should have sizes of approximately 0.001 m in the vicinity of the inserted needle. This implies

- A huge disparity between element size and object size: in other words small elements, and a large of small elements would be required. If we assume that H = 0.1 m (size of the object), and h = 0.001 m (size of the elements), then a uniform mesh of a cube will require 100 nodes to each side, and (assuming tets) 6M elements.
- As a consequence, precomputations become infeasible: condensing a 1M system requires 1G storage and 1 petaflop of computation power. A single force response costs 60M elements * 100 flops/element = 6 GFLOP in CPU cycles. During the simulation, 10k nodes are visible, so updates to the condensation matrices require 100 mflop per update
- As a consequence, very small time steps will be required for a dynamic integration: for h = 0.001, and wave speed c = 3 m/s, the critical time step are around 0.3 ms. Real-time performance requires 0.3 ms for computing a force response. This works out to 18 tflop/second computation power (Assuming linear elasticity).

When adaptive meshing is used, it is not possible to use condensation. In this scenario, a force computation costs 30k elements = 3 mflop. For 25 Hz, this requires 75 mflop/second, which is feasible on todays computers. However, for real-time performance, we require 3000 Hz (assuming a uniform time step for the entire mesh). This implies 9 gflop/second, which is not feasible for today's computers, but might be in the future.

The relaxation chapter of my thesis shows that static optimization and dynamic relaxation are on par when it comes to speed. This suggests that both multi resolution methods will probably also be as quick when properly implemented, and that differences amount to tweaking (dynamic methods need parameter choices) and physical realism (in a static optimization, there is no realistic force evolution and only conservative external forces can be handled.)

It uses simplex refinement meshing on a regular grid to achieve adaptivity. It seems to work, but now the results should be validated. We could have used Delaunay meshing, but this does not generalize well to 3 dimensions. For 3D tetrahedral meshing, our MICCAI paper shows that relocating nodes by projecting them is a risky business: it leads to degenerate elements when relocating nodes. Therefore, we try to see if we can do without relocating.

For residual tolerance 0.1, a boundary change requires 5 to 20 iterations of the CG relaxation. This allows haptic rates for the insertion.