SOniCS: Develop intuition on biomechanical systems through interactive error controlled simulations

,


Introduction
SOniCS is an ideal tool to develop more intuitive understanding of biomechanics systems.It enables the development of accelerated surrogate mod-• The material model complexity.Conversely to engineering materials such as steel or copper, the mechanical properties of living organs were only recently quantified [Gibbons, 1934, Payan andOhayon, 2017] (early 90s against 1700 for copper) and show immense variability [Mihai et al., 2017].
For instance, hyperelasticity equations can be written as a minimization of a tensorial function.In most cases, this minimization is solved using gradient-descent algorithms requiring the first and second derivatives of the functional.Therefore, obtaining such high-order nonlinear derivatives is not straightforward and can be prone to manual errors.
• The complexity of the simulation.In addition to the complexity of the material model, the simulation setup itself can be problematic.For instance, the material parameters are patient-specific and require data-driven or inference methods [Han et al., 2011, Urcun et al., 2021a].The simulations can also (partially) involve unknown boundary conditions [Tagliabue et al., 2021], contact with other organs [Courtecuisse et al., 2014, Miguez Pacheco et al., 2014], or surgical tools [Cotin et al., 1999, Lim et al., 2006].The problem can also require multi-physics models (such as FSI [Borazjani, 2013, Bianchi et al., 2017] (Fluid-Structure Interactions)) or incompressibility [Weiss et al., 1996, Mazier et al., 2022], where classical displacement-based finite elements are prone to locking.
• The error control and uncertainty of the solution.When dealing with biomechanical simulations, several uncertainties always arise from the material parameters, loads, geometry, or boundary conditions.This is mainly due to the difficulty in estimating the mechanical properties through exvivo methods or the topology of biological tissues using medical imaging.
Similarly, error control and mesh adaptivity are necessary to ensure homogeneous convergence of the solution over the domain and that the mesh is optimal given a quantity of interest [Allard et al., 2012].Therefore, quantifying the uncertainty or controlling the error on quantities of interests often requires making a very large number of simulations [Rappel et al., 2019a, Rappel et al., 2019b] which is incompatible with surgical timing [Hauseux et al., 2017, Hauseux et al., 2018].Meanwhile, those approaches could be functional in clinical settings by using accelerated simulations [Bui et al., 2018] to build surrogate models or/and machine learning models for faster solutions of those highly non-linear parametric problems.
• The real-time aspect.In addition to the previous point, for clinical environments, the running time of the numerical simulation is crucial.When performing an operation, the surgeon cannot, in general, spend minutes waiting for the model predictions.Eventually, this aspect is also applicable to artificial intelligence.Indeed, in order to build an efficient machine learning model, a significant amount of data is necessary.Using numerical simulations to create synthetic data is now standard and directly dependent on the simulation time [Odot et al., 2022, Deshpande et al., 2022].Consequently, the run-time of the simulation can be considered a principal feature of biomechanics simulations, and indeed, of any non-intuitive nonlinear problems subject to significant uncertainties in loading, boundary and initial conditions, and parameters.Nowadays, gaming engines such as Unity3D [Menard, 2011] or Unreal Engine [Sanders, 2016] offer real-time animation where physics-based algorithms can be included [Comas et al., 2008, Verschoor et al., 2018, Turini et al., 2019].
• The interaction with the user.In a surgical simulation setting, external variables can impact the simulation during the execution.For example, the exact movement of the surgeon's tool influences the simulation during run-time.Thus, the parameters of the simulation must be tuned, "live", to integrate interactions between the user and the simulation [Niroomandi et al., 2013, Wu et al., 2014].
• The visual rendering.Depending on the research field, visualization can play a critical role in the understanding of the results.For instance, in the computer graphics community, photo-realistic visualization is one of the main objectives [Gilles et al., 2011, Malgat et al., 2015].Contrastingly, in the mechanical engineering culture, visualization is a manner of extracting and understanding quantitatively a solution or data set (i.e., stress or displacement fields).For medical simulations, an optimal solution must combine both the accuracy of the results and photorealistic rendering reflecting the clinical ground truth all within clinical time frames [Guo et al., 2021].
According to the state-of-the-art, SOFA [Faure et al., 2012] (Simulation Open Framework Architecture) appears to be a suitable compromise.Indeed, SOFA employs efficient rendering while providing the possibility to interact in realtime with the running simulation.One can note that real-time computing is only possible depending on the complexity of the problem.Indeed, using excessive numbers of degrees of freedom (DOFs) or solving a highly nonlinear problem cannot result in a real-time simulation (without using model order reduction [Chinesta et al., 2014, Goury et al., 2016, Goury and Duriez, 2018] or machine learning [Deshpande et al., 2022]).SOFA can also manage complex simulations through an efficient implementation of contact, for example, or enabling multi-physics coupling.Therefore, only a few material models and elements are available.A similar issue is observed for material models where numerous implementations only focused on a 3-dimensional isotropic behavior.Therefore, coding a new material model or element in SOFA requires advanced C++ skills that may discourage individuals from using the software.
To alleviate the problem of complex material models, FEniCS [Alnaes et al., 2015] seems like an appropriate solution.Indeed, FEniCS may not possess all of SOFA's features, but definitely overcomes SOFA's capabilities for material model complexity.With FEniCS, the user can generate a large number of material models, regardless of the element's geometry or interpolation.Plus, it authorizes an export of the pertinent finite element tensors in C code to be efficiently plugged into SOFA.The benefits of the synergy are considerable.By using SOFA's interactivity and real-time features, the user can easily prototype a real-time simulation.Indeed, by modifying, in "live", various boundary conditions, geometries, or topologies, the user can effortlessly and rapidly verify modeling hypotheses of a specific problem.Combined with the specificities of FEniCS, the user can additionally smoothly prototype complex material models for modeling elaborate scenarios.Such feature has already been used by coupling FEniCS and Acegen but with different objectives [Lengiewicz et al., 2021].
To the authors' knowledge, this paper is the first to use FEniCS code generation capabilities for such an endeavor and is the first coupling between FEniCS and SOFA.
This paper has the following outline.We will first briefly introduce SOFA and FEniCS, highlighting the relative advantages and design choices in each.Section 2 will detail the plugin functionalities and a short tutorial for importing a Saint Venant-Kirchhoff model from FEniCS in SOFA.Then, section 3 will focus on confirming our implementation for various numerical tests.We will use a manufactured solution in 3.1 as validation and compare our solutions for a cantilever beam problem with SOFA in 3.2.The last test consists in implementing a new material model (Mooney-Rivlin) in SOFA and benchmarking it with FEBio [Maas et al., 2012] in 3.3.Finally, in the last section 4, we will use our plugin in a complex haptic simulation that cannot be implemented in FEniCS, using a custom material model inexistent in SOFA.

SOFA
SOFA was created in 2007 by a joint effort from Inria, CNRS, USTL, UJF, and MGH.Such piece of software aims to provide an efficient framework dedicated to research, prototyping, and the development of physics-based simulations.It is an open-source library distributed under the LGPL license, hosted on GitHub at https://github.com/sofa-framework/sofa,and developed by an international community.SOFA is modular.Users can create public or private plugins to include additional features.
SOFA is a C++ library, including Python wrappers for a user-friendly prototyping interface.It was originally designed for deformable solid mechanics but has been extended to various domains such as robotics, registration, fluid simulations, model-order reduction, and haptic simulations [Duriez et al., 2006, Duriez, 2013].SOFA exhibits many attractive features, but among them, the combination of multi-model representations and mappings differentiate it from other software.Multi-model representation: Most classical FE software uses an identical discretization for the whole model.Consequently, if one user wants to refine the mesh along a contact surface, the FE mesh will undergo the same refinement in the contact region.It can induce slow simulations for solving the FE system while the user was initially only interested in the contact part.Conversely, utilizing a multi-model representation approach, users can split the principal model into three distinct sub-models: deformation, collision, and visual.Thus, the user can decide to have high fidelity deformations with flawed contact detection while maintaining a fine rendering, or vice-versa.Similarly, an object can be made of several deformation models.For example, one can model a muscle by the interaction of FE 3D tetrahedra for the volume and 1D beams for modeling ligaments, using 2 different solvers.Mappings: In SOFA, the "mappings" are responsible for the communication between the different models.The models have parent-child relationships constructing a hierarchy (and a DOF hierarchy by extension).It enables propagating the positions, velocities, accelerations, and forces across the different models.For example, if the contact model calculates a force, it is mapped on the deformation model that will communicate back the computed displacement.
Finally, by combining the multi-model representation with the mappings, SOFA can build complex real-time simulations with high fidelity rendering.In addition to a scenegraph structure and visitors (responsible for going through the model hierarchy) implementation, it can account for interactivity with the users.
Despite the advantages provided by SOFA, some drawbacks have to be acknowledged.First, in terms of solid mechanics simulations, only a few elements and material models are available.Indeed, SOFA only proposes Lagrange linear elements, and the geometries are limited to segments, triangles, quadrangles, tetrahedra, and hexahedra.The following material models are coded: Boyce and Aruda [Arruda and Boyce, 1993], Costa [Costa et al., 2001], isotropic and anisotropic Hookean, Mooney Rivlin (2 invariants) [Mooney, 1940], classical and stabilized Neo-Hookean, Ogden [Ogden, 1972], Saint Venant-Kirchhoff, and Veronda Westman [Veronda and Westmann, 1970].Despite a reasonable number of mechanical models, a few of them are actually implemented for each element type.Secondly, the benefits provided by the mappings can also turn out to be a disadvantage when it comes to implementation.Indeed, the structure of the mappings is usually complex for unexperimented C++ users, and the mechanical tensors such as the Cauchy-Green or Piola-Kirchhoff are rarely computed.The two previous drawbacks are associated to the same flaw: the strong coupling between the material models and the topology of the element assumed within SOFA's architecture.This coupling implies that changing an element's topology or interpolation will involve a new mapping or the rewriting of the material model, even in the case of a similar material model.

The modular mechanics plugin (Caribou)
The initial goal of the plugin (called Caribou at the time of writing, https: //github.com/mimesis-inria/caribou)was to quickly implement new shape functions and their derivatives for different Immersed-Boundary and meshless domain discretization while keeping the compatibility with the existing SOFA surgical simulations [Brunet, 2020].Besides, the plugin enabled to effortlessly implement different volumetric quadrature schemes and several hyperelastic material models.Hence, the software design had to be generic enough to combine all the previous requirements.It also had to be efficient enough to avoid the creation of a bottleneck that would prevent the biomechanical model from meeting its computational speed requirement.Hence, the plugin was made as an extension to SOFA, bringing a redesigned software architecture.
In the plugin, the authors implemented a compile-time polymorphism design using generic C++ template programming.The idea is to write the code as close as possible to equations found in traditional FE books.Then, the C++ compiler optimizes the set of operations executed during the simulation while keeping an object-oriented code.In this design, the "Element" concept was created as a generic computational class that would be inherited by all element types.Similarly to OpenXFEM++ [Bordas et al., 2006], it provides a flexible implementation to add interpolation and quadrature numerical procedures quickly.Since standard isoparametric elements have a number of nodes, quadrature points, and shape functions already known at compile-time, most modern compilers will be able to aggressively inline the code to optimize the computation.
Finally, the plugin allows the creation of additional material models by simply defining three methods per material: the strain energy density function, the second Piola-Kirchhoff stress tensor function, and its derivative functions.These three functions are evaluated at a given integration point automatically provided by the plugin.This design delivers an undeniable advantage: writing a new material model is now independent of the topology and integration scheme.However, it comes with a non-negligible cost.The author of the new material model has to manually differentiate the strain energy twice and write it in C++.This manual intervention is error prone and can quickly become a substantial drawback for complex materials.

FEniCS
The FEniCS Project (FEniCS) [Alnaes et al., 2015] is a collection of tools for the automated solution of partial differential equations using the finite element method.Like SOFA, the FEniCS components are distributed under open-source licenses (LGPL v3 or later, and MIT) and development is hosted on GitHub at https://github.com/fenics.
A distinguishing feature of FEniCS is the ability to allow the user to write variational of weak formulations of finite element methods in a high-level Pythonbased domain specific language (DSL), the Unified Form Language (UFL) [Alnaes et al., 2014].Subsequently, that high-level description can be compiled/transformed using the FEniCS Form Compiler (FFC) [Kirby and Logg, 2006] into low-level and highperformance kernels.These kernels can calculate the corresponding local finite element tensor for a given cell in the mesh.UFL is also used by other finite element solvers with independently developed automatic code generation capabilities, notably Firedrake [Rathgeber et al., 2016] and Dune [Bastian et al., 2021].Compared with SOFA, FEniCS is limited in scope; its primary focus is the specification and solution of partial differential equations via the finite element method, leaving difficult problems like mesh generation, post-processing and visualisation to leading third-party packages such as Gmsh [Geuzaine and Remacle, 2009] and Paraview [Ahrens et al., 2005].
In the context of implementing finite element models of hyperelastic materials, this automatic approach has a number of advantages over the traditional route used by most finite element codes (including, to some extent, SOFA); differentiating analytical expressions for the residual (first derivative) and Jacobian (second derivative) of the energy functional for the hyperelastic model, picking a suitable finite element basis and then hand-coding the corresponding finite element kernels in a low-level language (C, Fortran, C++) for performance.Specifically: 1.The symbolic residual and Jacobian can be derived automatically using the symbolic differentiation capabilities of UFL.By contrast taking these derivatives by hand can be tedious and error-prone.
2. The compilation of the UFL description of the problem by FFC into the associated low level kernels is entirely automated.Again, this step is often time consuming and difficult to perform manually.
3. Because of the high-level description of the problem it is possible to experiment quickly with different concrete finite element formulations (material models, basis functions, element topology etc.) without manually modifying low-level kernel code.
The potential of this high-level approach for solid mechanics were recognised early on in the development of FEniCS, with two chapters in the FEniCS Book [Ølgaard andWells, 2012, Narayanan, 2012] promoting this direction.Since then, FEniCS has been used in a large number of publications on the topic of hyperelastic large-deformation elasticity e.g.[Baroli et al., 2012, Phunpeng and Baiz, 2015, Weis et al., 2017, Nguyen-Thanh et al., 2019, Patte et al., 2022].
Recently the FEniCS Project has undergone a major redevelopment, resulting in the new FEniCSx components; DOLFINx (the finite element problem solving environment, replacing DOLFIN), FFCx (the FEniCSx Form Compiler, replacing FFC) and Basix [Scroggs et al., 2022] (a finite element basis function tabulator, replacing FIAT [Kirby, 2004]).UFL is largely unchanged from the version used in the old FEniCS components and Firedrake.
In this work we do not use DOLFINx.DOLFINx contains the basic finite element data structures and algorithms (e.g.meshes, function spaces, assembly, interfacing with linear algebra data structures in e.g.PETSc []) and therefore there is a significant overlap with the functionality already available in SOFA.Directly interfacing DOLFINx and SOFA at the Application Programming Interface (API) level would be a significant technical challenge due to the substantial differences in their internal data structures.Dedicated weak coupling libraries such as PreCICE [Rodenberg et al., 2021] could be an interesting alternative to API coupling, but it is not a path that we explore in this work.
Instead, the approach taken by SOniCS is to only use UFL and FFCx (which in turn depends on Basix) to convert the high-level description of the finite element problem into low-level C code, which are then called using SOFA's existing C++ finite element data structures and algorithms.Compared with coupling DOLFINx and SOFA directly, our approach creates a relatively light compile-time coupling between FEniCS (specifically, the generated C code) and SOFA (a large complex C++ code with many dependencies).Consequently, no additional runtime dependencies required for SOFA.This methodology will be familiar to users of the DOLFINx C++ inteface where C finite element kernels are generated in a first step using UFL and FFCx and are then integrated into the DOLFINx solver in a second step through a standard compile/include/link approach.Without going into excessive detail, two changes in the redeveloped FEniCSx components have made SOniCS significantly easier to realise: 1. Basix and FFCx have full support for Serendipity finite elements of arbitrary polynomial order following the construction of Arnold and Awanou [Arnold and Awanou, 2011].Serendipity elements are used in SOFA and there was a desire to continue supporting Serendipity basis function due to their lower number of degrees of freedom per cell and generally lower number of local computations compared with standard tensor-product Lagrange elements.We remark that despite the widespread use of Serendipity elements in many solvers, they can only obtain optimal order convergence on affinely-mapped meshes, see e.g.[Arbogast et al., 2022] for more details.
2. FFCx outputs C99 compliant code according to the UFCx interface, which is specified as a C header file included with FFCx.This is in contrast with FFC, which outputs C++03 compliant code conforming to an interface specified with a C++ header file.This switch makes it significantly easier to call FFCx generated kernels from libraries with a C Foreign Function Interface (FFI) such as Python and Julia, or any language which can easily call functions with a C ABI (e.g.Fortran).Although SOFA is a C++ libraries and could certainly call C++ generated kernels, the C interface is simpler to use, consisting only of structs containing basic native data types and functions.

SOniCS
In this section, we present more in-depth the SOniCS (SOFA + FEniCS) plugin.We first introduce the procedure for defining the material model, the element geometry, and the quadrature rule or degree using the UFL (Python) syntax.
For simplicity, we only focused on a Saint Venant-Kirchhoff model.But the method can be generalized to all element types following the pipeline shown in figure 1.Secondly, we explain the methodology for converting the UFL script into efficient C kernels.Finally, we show the interface between the SOniCS plugin and SOFA, stating the conceptual and coding differences.For simplicity and as the modular mechanics plugin's name (Caribou) might change, we use the name SOFA to denote the combination of SOFA and the modular mechanics plugin.
Figure 1: Description of the SOniCS pipeline (on the right) and differences with SOFA(on the left).In SOFA, each element has to be defined, embedding its geometry, shape functions (including derivatives), and the quadrature scheme and degree.In SOniCS, it has been replaced by two Python lines of code for describing the element and its quadrature.The same benefit goes for the material model description.In SOFA, each material has to be created in a separated file stating its strain energy, derivating by the hand the second Piola Kirchhoff tensor (S) and its Jacobian.It was replaced in SOniCS by only defining the strain energy of the desired material model in UFL.The derivative of the strain energy will then be automatically calculated using the ffcx module.Finally, both plugins share the same Forcefield methods for assembling the global residual vector (R) and stiffness matrix (K).

UFL: from FE model to Python code
The first step is to define the FE model using the UFL syntax in a Python script.As an example, in listing 1, we describe the 3D simplest hyperelastic model and element: Saint Venant-Kirchhoff with linear Lagrange finite elements on tetrahedra.In the context of hyperelastic simulations, we will exclusively describe features that users could be interested in customizing.Listing 1: Python code example (material.py) of a Saint Venant-Kirchhoff material model using Lagrange linear tetrahedron.
Element: After importing the necessary packages, we can define the element geometry.In listing 1, on line 10, we used a linear Lagrange tetrahedron.The user can easily modify different parameters of the element, such as the geometry, the family type, or the interpolation degree.For example, by only changing line 10 to element = VectorElement (" Serendipity " , hexahedron , 2) the element is now a quadratic Serendipity hexahedron.Note that VectorElement creates by default a function space of vector field equal to the spatial dimension.A complete list of the element available in the Basix documentation [Scroggs et al., 2022].
Material model: By definition, boundary value problems for hyperelastic media can be expressed as minimization problems.For a domain Ω ⊂ R 3 , the goal is to find the displacement field u : Ω → R 3 that minimizes the total potential energy Π.The potential energy is given by where ψ is the elastic stored energy density, B is a body force (per unit reference volume), and T is a traction force (per unit reference area) prescribed on a boundary ∂Ω (of measure ds) of the domain Ω (of measure dx).
In listing 1, at line 37, we define the strain energy of a Saint Venant-Kirchhoff material model as: where λ and µ are Lamé material constants while E is the Green Lagrange strain tensor.Therefore, we observe that tensor E is expressed with respect to the displacement u, which is the unknown displacement field.Moreover, λ and µ are a function of the Young's modulus and of Poisson's ratio that are assumed to be known constants.
If the user would like to change the material model for Neo-Hookean with the strain energy density Quadrature rule: In listing 1, at line 37, when calculating the total potential energy it is also possible to choose the quadrature rule and the degree.In our example, we selected a quadrature degree of 1, triggering by default the Zienkiewicz and Taylor scheme [Zienkiewicz et al., 2014] for tetrahedra.Hence, the user could also choose to use a Gauss-Jacobi quadrature of degree 2 [Ralston and Rabinowitz, 2001] by replacing line 37 with: Pi = psi * dx ( degree =2 , scheme =" Gauss -Jacobi ") -inner (B , u ) * dx ( degree =2 , scheme =" Gauss -Jacobi ") -inner (T , u ) * ds ( degree =2 , scheme =" Gauss -Jacobi ") .

FFCx: from Python code to efficient C kernels
In the particular case of static hyperelastic simulations, we solve the following non-linear system of equation The R tensor is called the residual vector and is defined as the Gâteaux derivative of the total potential energy Π with respect to change in the displacement The tensor K is the Jacobian (also called stiffness in the context of mechanics) matrix and corresponds to the derivative of R Solving the non linear system in equation 4 can be achieved using the Newton-Raphson algorithm that will iteratively solve a set of linear systems, assuming an initial guess u n u n+1 = u n − du. (7) The two tensors can be derived symbolically and exported using the UFL syntax with the function derivative, as shown at lines 40 and 43 in listing 1.A simple call to ffcx will create a .cand .hfiles containing the code for generating the local R and K tensors.

Integration in SOniCS
In SOFA, the definitions of the residual and stiffness tensors are carried out within a C++ file, HyperelasticForcefield.cpp.Each material model is in a separate file.So far, only Saint Venant-Kirchhoff and Neo-Hookean models have been implemented.The users can easily access those functionalities through Python wrappers: node .addObject ( " S a i n t V e n a n t K i r c h h o f f M a t e r i a l " , young_modulus =E , poisson_ratio = nu ) node .addObject ( ' H y p e r e l a s t i c F o r c e f i e l d ') Listing 2: Python definition of a hyperealstic forcefield in SOFA using a Saint Venant-Kirchhoff material model.
The HyperelasticForcefield contains several functions, but we are particularly focusing on two of them.The addForce and assemble_stiffness functions are assembling the global residual and stiffness tensors respectively.Algorithm 1 details the addForce function, while algorithm 2 presents our reimplementation of the procedure.We did not detail the assemble_stiffness as it involves the exact same differences between the two implementations.
Algorithm 1 SOFA addForce function.The addForce function is in charge of assembling the global residual vector.for i in range(0, NumberOfNodesPerElement) do i indicates the element node index while global(i) denotes the global node index 20: end for 21: end for The structure is similar but we can still observe a few differences.
Algorithm 2 SOniCS addForce function.The addForce function is in charge of assembling the global residual vector.for i in range(0, NumberOfNodesPerElement) do end for 16: end for • The new implementation of algorithm 2 is more concise and involves less visible tensorial operations because all those operations are efficiently hard coded in the C file provided by FFCx.For example, in algorithm 1, lines 5 to 17 were replaced in the new algorithm 2 by solely line 11.
• Algorithm 2 needs to have access to the initial position of the object and to the displacement vector.This was indeed not needed in the previous implementation since the modular mechanics plugin takes advantage of writing the deformation gradient only based on the current nodal coordinates x: ∇ Ω0 and x 0 respectively denote the gradient and the nodal coordinates in the initial configuration, thus saving one extra vector operation.
• In the SOFA implementation, the boundary conditions and body forces are treated in separate files.In the new implementation of the Forcefield 1, the boundary conditions and body forces are now directly carried out in the Forcefield on lines 11 and 12.It avoids calling another function to loop again through every element of the object, thus speeding up the assembly of the residual vector.
Based on this new implementation, we created a new forcefield HyperelasticForcefield_FEniCS as close as possible to the existing syntax of HyperelasticForcefield.We also needed to tune the existing material definition to replace unnecessary calculations and allow us to read the corresponding .cfile.
node .addObject ( " F En iC S_ Mat er ia l " , material = " S a i n t V e n a n t K i r c h h o f f " , young_modulus =E , poisson_ratio = nu ) node .addObject ( ' H y p e r e l a s t i c F o r c e f i e l d _ F E n i C S ') Listing 3: Python definition of a hyperelastic forcefield in SOniCS using a Saint Venant-Kirchhoff material model.
Finally, the last hurdle was the element definitions.Indeed, SOFA and FEniCS do not use the same vertices, edges, and facets ordering (as shown in figure 2).To avoid any conflict with the existing users of SOFA and solve the ordering issue, we proposed rearranging the topology indices, edges, and vertices and create new elements.It ensured an accurate integration over the elements (especially for quadratic Serendipity integrating over the edges) and preserved an appropriate visualization.Those elements have been interfaced with the existing topology named CaribouTopology.[4 , 5 , 0 , 1 , 7 , 6 , 3 , 2]]) Listing 4: Python definition of hexahedron topology in SOFA and SOniCS.

Numerical examples
In this section we describe three numerical examples used for the validation of our SOniCS implementation.We use the same domain description for each example while varying the boundary conditions and material parameters of each simulation.
Let Ω be a domain represented by a squared-section beam of dimensions 80×15×15m 3 , considered fixed on the right side (u = 0 on Γ D ) while Neumann boundary conditions are applied on the left side (Γ N ), as shown in figure 3a and  3b.
Ω was discretized using two different geometrical elements using linear and quadratic interpolations.P1 and P2 elements stand for linear or quadratic tetrahedra, while Q1 and Q2 denote linear and quadratic hexahedra.
To solve each hyperlastic formulation described in equation 4 we used an identical implementation of the classical Newton Raphson (NR) solver for SOn-iCS, SOFA and FEBio.The solver had the following parameters: a maximum of 25 iterations with a residual and displacement tolerance of 10 −10 .In order to compare the running time of the two implementations, we, therefore, introduced the mean NR iteration time.We defined the NR iteration time as the duration for assembling and factorizing the system matrix, solving and propagating the unknown increment, updating, and computing the force and displacement residual.After checking that the same number of iterations have been achieved, we averaged the total time over the number of iterations needed for the solver convergence.All calculations were performed using an Intel® Core™ i5-6300HQ CPU @ 2.30GHz × 4 processor with a 16GiB memory and a NV117 / Mesa Intel® HD Graphics 530 (SKL GT2) graphics card.
We evaluated the soundness of the SOniCS solution using SOFA, FEBio or a manufactured solution as the reference solution and computed the Euclidean relative L 2 error (e(u, v)) [Odot et al., 2022].
where u and v are the calculated displacements for SOniCS and the reference implementation, respectively.
In this section, we first compare our solution with an analytical one: the manufactured solution.Then, we consider a clamped cantilever beam subject to Neumann boundary conditions and compare its deformation with the SOFA solution.Finally, using the same cantilever beam, we implemented a Mooney Rivlin model (uncoded in SOFA) using SOniCS and compared the solution with FEBio.

Manufactured solution
Aiming at code verification, the method of the manufactured solution consists in choosing an exact solution to the problem as an analytical expression [Chamberland et al., 2010].The chosen analytical expression is then inserted into the Partial Differential Equation (PDE) under consideration to find the conditions that lead to this solution.In general, the manufactured chosen solution is expressed in simple primitive functions like sin(), exp(), tanh(), etc...In the context of hyperelastic equations, we considered the following manufactured solution for the displacement Starting from the above chosen displacement and using continuum mechanics laws, the relative analytical forces are applied as Neumann boundary conditions and deduced as follows Where F is the deformation gradient, I d is the identity matrix of dimension d, P is the first Piola-Kirchhoff stress tensor and W is the strain energy density depending on the material model constitutive law.We used a Saint Venant-Kirchhoff material (2) with a Young's Modulus of 3 kPa and a Poisson's ratio of 0.3, the computation of this solution was performed using Python Sympy package [Meurer et al., 2017].
In this experiment, we generated 8 and 6 discretizations of P1 and P2 elements, respectively, with a decreasing element size in both scenarios.For each discretization, we applied the relative analytical forces deduced from the manufactured solution (10) and used SOniCS Saint Venant-Kirchhoff material implementation with the same parameters as Sympy's to fill the domain.The displacements obtained were compared to the chosen analytical solution in equation (10) for each discretization.The results are presented in figure 4, the error metric is the relative mean error in equation ( 9).
Figure 4: Plot of the mesh convergence analysis of the manufactured solution.The L 2 error between the analytical and the SOniCS simulation is calculated for different mesh refinement with fixed parameters (E=3 kPa and ν=0.3) for P1 linear tetrahedra (blue) and P2 quadratic tetrahedra (red) elements.

Benchmark with SOFA
The cantilever beam deflection is a classical mechanical test case, as you can smoothly refine the mesh due to the simplicity of the geometry or modify its boundary conditions to fit real-life experiments.In this context, the beam was still clamped on the right side (natural Dirichlet condition on Γ D ) while Neumann boundary conditions were applied on the left side (Γ N ).To compare our SOniCS implementation, we model the deformation of the beam with two hyperelastic material models: Saint Venant-Kirchhoff and Neo-Hookean.We fixed the mechanical parameters and the Neumann boundary conditions equal to −10 Pa in the y direction, until reaching sufficient large deformations with the same parameters as before: E = 3 kPa and ν = 0.3.
This study aims at comparing the finite element solutions provided by SOn-iCS and SOFA under the same constraints in terms of computational and running time performances.To do so, we computed the mean relative L 2 error (e(u SOniCS , u SOFA )) between the SOniCS solution using SOFA as the reference solution.
The results obtained are presented in tables 1 and 2 for Saint Venant-Kirchhoff and Neo-Hookean materials, respectively.

Benchmark with FEBio
FEBio is a open-source finite element package specifically designed for biomechanical applications.It offers modeling scenarios, a wide range of constitutive material models, and boundary conditions relevant to numerous research areas in biomechanics.In this section, FEBio was used to compute the same scenarios as in section 3.2 to evaluate the trustworthiness of SOniCS.A more advanced constitutive material model, Mooney Rivlin, was introduced for this purpose Where C 01 , C 10 , and K are the material constants in addition to the modified invariants Tables 3, 4, 5 show the results obtained for Saint Venant-Kirchhoff, Neo-Hookean and Mooney Rivlin material models and considering the four discretizations implemented so far in SOniCS.The error evaluation is still based on the mean relative error defined in equation 9 using FEBio as the reference while using a Newton Raphson solver with the same characteristics in both cases.4 Hozapfel and Ogden anisotropic material model coupled with haptic simulation In the context of numerical surgical simulations, a robot haptic feedback has been shown to be a consistent tool for drastically improving user-interactions and opening up countless applications.Among them, haptic devices have mainly been used as a training tool for surgeons.Indeed prior to surgery, under the assumption of known geometry and mechanical properties of the patient's organ, a surgeon would be able to plan and better choose between specific surgical paths/approaches.In this paper, we used the 3D between the instrument and the simulations.For this hypothetical simulation, we virtually simulate the contact between a surgical tool and a liver during surgery.The liver was described by an anisotropic Holzapfel Ogden model [Holzapfel andOgden, 2009, Pezzuto et al., 2014], with an existing FEn-iCS implementation [Hauseux et al., 2018].
ψ iso and ψ vol are the isochoric and volumetric part of the strain energy density function respectively.The volumetric part can be evaluated as a function of the bulk modulus κ of the material and J with The transversely isotropic behavior can be obtained by removing the parameters a f s , b f s , a s and b s , while the isotropic behavior is obtained by also suppressing the two parameters a f and b f .This kind of model is frequently used to model orthotropic materials (e.g.muscle with fibers or tendons).Vectors f 0 , s 0 are the unit base vectors normal to the planes of symmetry.For our application, we selected the same material properties as chosen in [Hauseux et al., 2018].
The instrument is assumed to be a rigid body kinematically constrained by the haptic device position at each time step.The contact forces generated by the collision between the instrument and the liver are calculated using frictional contact and an implicit Euler scheme to solve the dynamic system [Courtecuisse et al., 2013].Finally, the contact forces are transmitted back to the user's hand through the haptic device.The result is a simulation running at 100 FPS (Frames Per Seconds) on average displayed in figure 5.A full video of the interaction with the haptic device is available in the Supplementary Materials.The real-time performance has been obtained by using a small number of DOFs (543).A higher number of DOFs could be used once a GPU implementation is available.
Figure 5: Numerical simulation of a liver in contact with a surgical tool connected to a haptic device.On the right, the 3D Systems Touch Haptic Device is controlled by the user.On the left, the liver is modeled using a Holzapfel Ogden anisotropic material in contact with the surgical tool (in red) guided by the user.In case of contact detection, the contact forces are transmitted to the user through the haptic device.A video of the simulation is available as supplementary materials.

Discussion
We shows that the SOniCS plugin is an efficient implementation of material models for hyperelastic simulations and enables the user to develop an intuitive understanding of the impact of modelling choices on accuracy reliability of the predictions.We first demonstrated a convergence study for Saint Venant-Kirchhoff material using P1 and P2 elements with the manufactured solution in section 3.1.Indeed, as expected, by refining the mesh, the L 2 relative error almost linearly decreases when increasing the number of DOFs (on a log-log plot), showing the stability of our method.
Then, from section 3.2, two main results are noteworthy from tables 1 and 2. First, the relative error of the displacement between SOniCS and SOFA, for both material models, is close to machine precision for P1, P2, Q1, and Q2 discretizations.Finally, the last comment concerns the mean NR iteration time.We observe that the SOniCS implementation is slightly faster than SOFA.The difference increases when using more DOFs, e.g., 2.40 ms difference for P1 against 5.17 ms difference for P2 elements for the Saint Venant-Kirchhoff model.The reason for this difference is the need for SOFA to compute the shape functions and derivatives, then calculate the local residual and Jacobian using multiple tensor operations.Conversely, SOniCS has all those operations efficiently hard-coded in C kernels, thus performing faster than SOFA.
We compared the SOniCS and FEBio simulations for several material models: Saint Venant-Kirchhoff, Neo-Hookean, and Mooney-Rivlin in section 3.3.For P1 elements, the errors are close to machine precision for all 3 models.For P2 and Q2 elements, the 3 models display similar errors that still represent a minor error (less than 0.08 and 0.1, respectively) which is of same magnitude between SOFA and FEBio, as well as FEniCS and FEBio.The reasons for those minor errors could be the difference in the implementation of the elements or in the choice of solver parameters.For Q1 elements, we reached an accuracy near machine precision for Saint Venant-Kirchhoff and Neo-Hookean.However, the relative error rose close to 11 for the Mooney-Rivlin model.To further understand this divergence, we used a third open-source software AceGen.AceGen also uses an automatic code generation package for the symbolic generation of new finite elements.The cantilever beam scenario presented in section 3.3 was reproduced using AceGen under the same conditions and with a Q1 discretization of the domain.Using the same metric as defined in equation 9, the results are the following: e(u SOniCS , u Acegen ) = 3.33 and e(u FEBio , u Acegen ) = 14.19.Even if SOniCS and Acegen showed similar results, a more in-depth study would be needed to confirm the soundness of our solution.Despite using a finer mesh, the error was getting slightly lower but was still noticeable.Eventually, FEBio has shown difficulty in converging with trivial parameter sets or when increasing the number of DOFs.Thus, as shown in tables 3,4, and 5, on average, SOn-iCS is solving the equation system faster than FEBio.Several reasons could explain those differences, such as the number of quadrature points used, the implementation of the Newton-Raphson scheme, or the solver parameters.
Finally, we showed the capabilities of the SOniCS plugin in simulating com-plex material models such as the Holzapfel Ogden anisotropic model coupled with a haptic device in section 4. The material model was effortlessly implemented for several elements (P1, P2, Q1, and Q2) without needing any manual derivation or coding.The result is a real-time simulator functional for surgeons' training or any other biomechanics simulation replicating the behavior of a liver in contact with a surgical tool.It was not possible to verify our simulation as we did not find any similar implementation of such material in any open-source software.

Conclusion and outlook
We performed several numerical experiments to develop intuitive understanding of new material models in SOFA using the SOniCS plugin.First, we validated the most common hyperelastic material models: Saint Venant-Kirchhoff and Neo-Hookean using a manufactured solution.Then, utilizing FEBio as a reference, we validated our implementation of a Mooney-Rivlin material model using P1, P2, Q1, and Q2 elements.The final application employed a haptic device to interact with an anisotropic Holzapfel Ogden liver model in real-time.
The study used our SOniCS plugin to generate optimized C code for complex material models compatible with SOFA.On one side, we benefited from FEniCS automatic differentiation and code generation capabilities to bypass the difficulties of deriving and implementing the consistent Jacobian in SOFA.On the user side, we implemented compatible and user-friendly SOFA forcefields to use the FEniCS C kernels.A SOFA user can now easily define a new material model by specifying its strain energy function, element geometry or family and the quadrature scheme and degree only in Python.We made the open-source code and all data and test cases available as the supplementary material.
In future work, we intend to apply the SOniCS plugin for solving more complex mechanical phenomena.For example, mixed formulations for solving incompressible materials, viscous or plastic effects, and multi-material systems.In this paper, we only utilized the plugin for solving hyperelasticity equations, but it would be interesting to tackle multi-physics problems such as thermomechanics or magnetomechanics.We only used Lagrangian P1, P2, Q1, and Q2 elements while more geometries such as prisms could be relevant for the SOFA and FEniCS community.Finally, our validation only focused on the deformation field of the structures, while a more in-depth study would also enable us to validate an analytical stress field.
Some work remains to improve user-friendliness of our package.Indeed, even if writing the Python file describing the material model and generating the associated C file is mostly effortless and automated, some steps are still manual.Indeed, the FEniCS_Material C++ class must have knowledge of every new C file created at compile time.Hence, for the moment, the user still has to manually specify two C++ functions in the code and recompile the whole plugin.Even if this step is manageable, it still requires diving into the C++ code, which can discourage a few users.Meanwhile, to mitigate this effect, actions have been taken by providing detailed documentation and tutorials for this crucial step.In addition, future works aim at improving this stage by directly providing the path of the C file in the Python code to trigger just-in-time compilation for new materials.

S
1: for element in elements do 2: X ← element.positionsreturn the current positions of the element 3: R global ← 0 zero the global residual vector of dimension (DOFs ×3) 4: R local ← 0 zero the local residual vector of dimension (element DOFs ×3) ← quadrature.nodes.shapefunctions derivatives return the derivatives of the shape functions of the quadrature nodes ← f(C, MaterialParameters) return the second Piola-Kirchhoff depending on the material parameters and kinematics tensors 13: for i in range(0, NumberOfNodesPerElement) do 14: dx ← dN [i] T 15: R local [i] ← (detJ • w) • F • S • dx allocate the result in the local residual vector

Figure 2 :
Figure 2: Local numbering of element vertices and edges in both FEniCS and SOFA

Figure 3 :
Figure 3: Cantilever beam domain discretization and displacement field.Ω is a domain represented by a squared-section beam of dimensions 80 × 15 × 15m 3 , considered fixed on the right side (u = 0 on Γ D ) while Neumann boundary conditions are applied on the left side (Γ N ).
based on the classic invariants I C = tr(C), II C = 1 2 (tr(C)) 2 − tr C 2 .In order to obtain sufficiently large deformations, we chose the following material parameters: C 01 = 2000 Pa, C 10 = 100 Pa, and K = 1000 Pa.
the element node index while global(i) denotes the global node index

Table 1 :
Mean relative error (e(u SOniCS , u SOFA ) defined in equation 9) and mean NR (Newton-Raphson) iteration time between SOniCS and SOFA for different element geometries and interpolation schemes using Saint Venant-Kirchhoff material model.P1 and P2 elements stand for linear or quadratic tetrahedra, while Q1 and Q2 denote linear and quadratic hexahedra.

Table 2 :
Mean relative error (e(u SOniCS , u SOFA ) defined in equation 9) and mean NR (Newton-Raphson) iteration time between SOniCS and SOFA for different element geometries and interpolation schemes using Neo-Hookean material model.P1 and P2 elements stand for linear or quadratic tetrahedra, while Q1 and Q2 denote linear and quadratic hexahedra.

Table 4 :
Mean relative error (e(u SOniCS , u FEBio ) defined in equation 9) and mean NR (Newton-Raphson) iteration time between SOniCS and FEBio for different element geometries and interpolation schemes using Neo-Hookean material model.P1 and P2 elements stand for linear or quadratic tetrahedra, while Q1 and Q2 denote linear and quadratic hexahedra.

Table 5 :
Mean relative error (e(u SOniCS , u FEBio ) defined in equation 9) and mean NR (Newton-Raphson) iteration time between SOniCS and FEBio for different element geometries and interpolation schemes using Mooney Rivlin material model.P1 and P2 elements stand for linear or quadratic tetrahedra, while Q1 and Q2 denote linear and quadratic hexahedra.