Real-time elastodynamic simulations in Extended Reality (XR) environments show promise for supporting remote emergency medical procedures during disasters. However, determining the optimal combination of energy models and computational frameworks for real-time 3D elastodynamic simulations is complex due to the multitude of existing methods and their potential combinations. This paper introduces a 2D energy-based Finite Element Method (FEM) prototyping framework designed to facilitate comparative analysis of various energy models across diferent computational frameworks. To enable the reproduction of disaster scenarios in 2D space, we propose an algorithm for semi-automatic conversion of arbitrary geometries into 2D FEM-compatible meshes. Using this framework, we simulated a leg trapped under collapsed building debris, computing and visualizing stress distributions to aid medical decision-making. Our framework facilitates the implementation of common energy models and enables eficient comparison of their combinations with various real-time computation frameworks, aiding the selection of optimal approaches for XR-oriented medical elastodynamics simulations.
To demonstrate the practical applicability of our framework, we modeled a case study of a leg trapped under In recent years, multiple applications of XR technologies collapsed building debris. This simulation computed and in medical practices have led to notable advancements visualized Von Mises stress distributions, providing critin remote disaster healthcare [1, 2]. During major dis- ical data that could inform medical decision-making in asters, physical barriers may substantially impact the real-world disaster response situations. The contribuprovision of professional medical treatment for victims. tions of this paper are: One solution is combining elastodynamic simulations with XR technology, which has multiple benefits for med- • A versatile 2D energy-based FEM prototyping ical treatment and medical education training [1]. One framework challenge is that there are several existing energy models • An algorithm for semi-automatic conversion of and solvers, each of which has its own strengths and lim- arbitrary geometries into 2D FEM-compatible itations. Determining the optimal combination of these meshes elements for a specific medical scenario is a non-trivial • A comparative analysis of common energy modtask, requiring thorough comparative analysis and perfor- els and their combinations with various real-time mance evaluation. To address this challenge, we propose computation frameworks a 2D energy-based FEM prototyping framework. This • A case study demonstrating the application of the framework is designed to compare various energy mod- framework in simulating a disaster scenario with els (e.g., Saint Venant-Kirchhof (StVK), Co-rotational, potential medical implications and Neo-Hookean) across diferent computational frameworks (e.g., Explicit FEM, Position-Based Dynamics (PBD) This paper is organized as follows: Section 2 reviews [3], and Extended Position-Based Dynamics (XPBD) [4]). related work in the fields of medical elastodynamics. SecIt serves as a stepping stone towards more complex 3D tion 3 describes the methodology behind our 2D energysimulations. By focusing on 2D representations, we can based FEM framework and the semi-automatic mesh conreduce computational complexity while still capturing es- version algorithm. Section 4 presents the results of our sential physical behaviors, allowing for rapid prototyping case study and comparative analysis. Finally, Section 5 and evaluation of diferent approaches. discusses the limitations of our framework and summarizes future works.
The finite element method has a broad range of applications across various fields, including but not limited to engineering, physics, and biomedical sciences [5, 6]. Due to its ability to accurately simulate the characteristics of deformable objects and represent multiple complex materials and structures, it has been widely used in 3D elastodynamics to simulate brittle and ductile fracture of materials [7, 8, 9].
In the field of medical elastodynamics, FEM has found numerous applications. For example, Sase et al. [10] proposed a GPU-accelerated FEM-based approach for simulating brain fissure opening in surgical procedures. They also introduced a volume embedding method [11] to preserve complex topological structures during surgical simulations, and a penalty based method for rigiddeformable objects coupling [12]. In these cases, all of their methods employed the Corotational FEM [13] for rapid computation of organ deformation simulations. Our study extends this line of research by comparing the performance of the Corotational energy model with St. Venant-Kirchhof (StVK) and Neo-Hookean energy models. Unlike previous studies, we utilize explicit FEM rather than implicit solvers. This approach allows us to simply evaluate the eficacy of diferent energy models in the context of real-time medical elastodynamics simulations.
Position-Based Dynamics (PBD) is a real-time physicsoriented simulation framework initially proposed by Müller et al. in the computer graphics community [3]. Unlike other physics-based simulation methods, PBD primarily focuses on converting physics models into constraint forms and solving them. Although this framework has proven capable of achieving accuracy comparable to physics-based models through multiple iterations of computation, developers primarily use in real-time physical simulations with fewer iterations to achieve satisfactory visual efects. Macklin et al. [ 4] identified limitations in the PBD framework when simulating elastic objects, as the stifness of materials depends on the timestep. In response, they introduced an approximate implicit Euler method called extended PBD (XPBD).
Given the stability and ease of parallelization of the PBD framework, numerous applications have emerged for medical elastodynamics simulations. For instance, Tai et al. [14] developed a virtual surgical training system in Augmented Reality (AR), using XPBD for soft tissue elastodynamic simulations. Camara et al. [15] utilized the PBD-based library NVIDIA FleX to develop a simulation platform for exploring optimal material properties and other parameters. Moreover, Yu et al. [16] employed PBD methods in a Virtual Reality (VR) environment to develop a real-time medical education training system. While these studies demonstrate the direct application of PBD/XPBD in medical simulations, our framework develops energy model-based PBD and XPBD approaches, building upon the FEM-PBD method proposed by Bender et al. [17].
This section will describe the implementation of our developed 2D energy-based FEM prototyping framework. First, we will describe the computation of various energy models when employing explicit FEM methods. Subsequently, we will detail the calculation of diferent energy models within the position-based approach. Finally, we will introduce the semi-automatic mesh conversion algorithm for 2D scenarios.
In continuum mechanics, the deformation of a material can be conceptualized as an afine mapping (· ) from the material space to the deformed world space [5]. This mapping allows us to describe the position of any point in the deformed world space. Specifically, for a material point ¯ in the material space, its position in the world space is given by: =(¯) = ¯ + = =
( ¯ + ) (¯) ¯ ¯ (1) where is the deformation gradient and is a translation matrix.
This formulation provides a fundamental basis for analyzing material deformation in our 2D energy-based FEM framework. It enables us to track the movement and deformation of each point in the material over time, which is crucial for computing various energy models and simulating elastodynamic behavior.
The energy models Ψ used in elastodynamics calculations can be viewed as evaluation functions that measure the degree of material deformation. These can be expressed as [13, 18, 19, 20, 21]: Ψ StVK = ‖‖2F +
(tr())2 2 Ψ Co-rotational = ‖ − ‖2F + (tr( − ))2 Ψ Neo-Hookean = (tr( )− 3)− log(det( )) 2 log2(det( )) (2) 2 1 2 + 2 = ( − ) = = 2(1 + )
(1 + )(1 − 2 ) Frobenius norm. where is a matrix representing the pure rotational part of , which must satisfy = . denotes Green’s strain tensor. tr(· ) represents the trace of a matrix. The Lamé parameters, denoted as and , characterize the elastic properties of the material. These coeficients are directly related to two fundamental material constants: Young’s modulus and Poisson’s ratio . det(· ) denotes the determinant of a matrix. ‖ · ‖ 2F represents the squared
In the explicit FEM framework, we can derive the nodal forces for each independent element based on various energy models Ψ (i.e., elastic potential function ( + ∆ ) is minimized to zero. For an individual constraint, the position adjustment ∆ is derived by resolving the linearized equation. This process can be (3) represented as follows: energy) as follows: = − Ψ Ψ = − = − ( ) Where represents the area of the element in its undeformed configuration. is an energy-independent matrix that can be derived directly from the defined element. ( ) denotes the first Piola-Kirchhof stress tensor, which can be derived from Equation 2.
Based on the derivation of the aforementioned energy models, we can eficiently decompose the calculation of nodal forces for each element in the explicit FEM method into energy-independent and energy-dependent components. This decomposition facilitates the seamless substitution of various energy models within the explicit FEM framework, enabling straightforward validation of their performance and behavior. The algorithm for the 2D explicit FEM is presented in Algorithm I.
To integrate the FEM method into the PBD framework, we primarily developed our framework using the position-based energy reduction method proposed by Bender et al. [17] Additionally, we attempted to implement the energy constraints within the XPBD method [4] and conducted comparative experiments.
The position-based energy reduction method aims to determine position adjustments ∆ such that the energy ( + ∆ ) ≈ () + ∇ ()∆ = 0 ∇ = ∫︁ Ψ
= ∫︁ ( ) d¯ = d¯ ∫︁ Ψ d¯
To constrain ∆ to align with the direction of ∇, a Lagrange multiplier is introduced, which is defined such that: () = − ∑︀ |∇ ()|2 Algorithm I 2D Explicit FEM Algorithm 1: for each finite element do 2: compute 3: compute ( ) 4: compute (Equation (3)) 5: accumulate forces for each vertex in 6: end for 7: ← − 1( + ext) 8: +1 + ∆ 9: +1 ←← + ∆ 10: damp velocities +1 (4) (5)
1 represents the inverse of the particle where = mass, and the position correction for each particle is determined by: ∆ = ∇()
geometries onto a 2D plane, we first establish a camera coordinate system using three vectors: camera position , focal point , and up vector . From these, we derive the (6) orthonormal basis: ˆ = |−− | , ˆ = |×× | , ˆ = ˆ × ˆ
We then construct the camera-to-world transformation matrix :
In contrast to the energy-based FEM framework, the energy-based PBD framework requires a priori conversion of all vertices of each element into a particle-based ⎡x x − x x⎤ ddeaptaenstdreuncttumraes,se,nvseuloricnitgy,thaantdepaocshitvioenrtienxfopromssaetsiosens. Ains- = ⎢⎢⎣yz yz −− yz yz ⎥⎦⎥ (8) for XPBD method can be conceptualized as a solver that 0 0 0 1 approximates an implicit Euler solution. It difers from the PBD framework in its computation of position cor- For a point = [ ] in world coordinates, we convert rections, as illustrated below: it to camera coordinates ′ using: ∆ = − 1∇() ∆ ∆ =
− () − ∇ − 1∇ + ′ = [x 1] · − 1 (9) (7)
Finally, we project ′ = (′, ′, ′, ′) onto the 2D plane as proj = (′, ′). where is a compliance coeficient. The algorithm for the energy-based XPBD is presented in Algorithm II.
ter projecting the 3D geometry to 2D, we obtain a large set of 2D point data that needs to be processed. Our next 3.3. Mesh Conversion task is to reconstruct the polygon shape information To facilitate the simulation of complex and irregular geo- based on these 2D point data. While the most straightmetric shapes in elastodynamics within our developed forward approach typically involves using a convex hull framework, we propose a computational pipeline capa- algorithm to obtain the contour information of these ble of transforming arbitrary geometries into 2D FEM- points, this method struggles to handle concave polycompatible meshes. The following sections elucidate the gon shapes. Therefore, we opted to employ the alpha step-by-step implementation process of this computa- shapes algorithm [22] for 2D polygon shape reconstructional pipeline. tion. An experimental result of applying the alpha shape algorithm can be found in Figure 1. Generation of FEM-Compatible Meshes from 2D evaluate computational performance. The following exPolygon Directly converting the polygons into FEM- perimental parameters were used: a Young’s modulus of compatible meshes is challenging due to the massive 1000 Pa, a Poisson’s ratio of 0.3, a mass of 1 kg, and a time number of vertices in polygons. This is because numer- step of 0.001 s. As illustrated in Figure 2, the simulation ous vertices allow for accurate representation of the orig- results show that there are no significant diferences in viinal geometric object shape, while they complicate the sual efects across the diferent models for the cantilever Delaunay triangulation process, which is essential for beam deformation. Furthermore, to precisely capture the producing FEM-computable meshes. As an example of similarities and diferences among these energy models, complex polygon shapes, Delaunay triangulation often we recorded the variations in the vertical displacement generates numerous small-area triangles. These small of the cantilever beam tip node, as presented in Table triangles can impact both the stability and computational 1. The statistical results indicate that the displacement eficiency of FEM solvers. To overcome this challenge, amounts generated through deformation are relatively we implement a polygon simplification step before ap- similar for the Co-rotational and Neo-Hookean models. plying Delaunay triangulation. The Douglas-Peucker As for the computational performance, the StVK model algorithm primarily forms the basis of the polygon sim- and the Neo-Hookean model require approximately 14 plification process; we choose it for its eficiency in re- seconds to compute 2000 steps, whereas the corotational ducing the number of vertices while preserving essential method necessitates nearly 20 seconds. Additionally, we shape characteristics. For Delaunay triangulation, we conducted an experiment to test the stability of diferent additionally incorporate maximum area and minimum energy models by increasing the Young’s modulus and angle constraints to produce FEM-compatible meshes. reducing the time step. The results of this experiment The complete mesh conversion process is illustrated in demonstrated that all of the energy models can achieve Figure 1. stable simulation for the cantilever beam.
Given these results, the StVK model emerges as the optimal choice for this cantilever beam case study. It not 4. Experimental Results only ensures simulation stability across diferent parameter settings but also demonstrates the highest computational eficiency.
Case A In this case study, we simulate a cantilever beam (5 m long, 1 m high), which is fixed at its left end, and utilize explicit FEM frameworks with diferent energy models to test its deformation under gravity and Algorithm II 2D Energy Based XPBD Algorithm 1: +1 + ∆ − 1 ext 2: +1 ←← + ∆ +1 3: ← 0 4: for iter := 1 to maxIterations do 5: for each energy constraint () do 6: compute ∆ 7: compute ∆ 8: ← + ∆ 9: ← + ∆ 10: end for 11: end for 12: +1 ← Δ1 (+1 − ) 13: damp velocities +1 Case B In this case study, we initially fix both ends of the bar-shaped elastic body. Subsequently, we apply a stretching operation to the nodes on the right side of the elastic body. The fundamental experimental parameters remain consistent with those in Case A. The final experimental results are presented in Figure 3. Examining the longitudinal contraction of the elastic body post-stretching, we observe that the StVK model exhibits more pronounced contraction compared to the corotational and Neo-Hookean models. However, when we attempted to test the stability of various energy models by increasing the Young’s modulus and reducing the time step, the StVK model required a smaller time step than the corotational and Neo-Hookean models to achieve stable simulations. From these observations, we can conclude that the StVK energy model is the optimal choice when simulating elastic bodies with lower Young’s module in 2D scenarios. However, when simulating objects
with higher Young’s module and prioritizing computational performance and stability, the corotational and Neo-Hookean models prove more suitable.
In this paper, we develop an energy model-based 2D FEM framework for XR-oriented medical elastodynamics. The Case C In this case study, we adopt the same exper- main goal of developing this framework is to assist in imental setup as in Case A. However, we will evaluate determining the optimal combination of energy models and compare the performance of the StVK model imple- and computational frameworks for real-time 3D elastomented within three diferent computational frameworks: dynamic simulations. Through comparative experiments explicit FEM, PBD, and XPBD. The simulation results in conducted in 2D space across 4 case studies, we analyze Figure 4 show that even when using the same energy the performance of various energy models and elastomodel, elastic bodies exhibit diferent behaviors under dynamic frameworks, subsequently identifying their apPBD and XPBD frameworks. It should be noted here plicable simulation scenarios. Furthermore, to enable that the deformation under the PBD framework is not our proposed framework to simulate elastodynamics of significant. This is because the material stifness in the more complex geometries, we introduce an algorithm for PBD framework depends on the time step. From this ob- semi-automatic conversion of arbitrary geometries into servation, we should be cautious when using PBD-based 2D FEM-compatible meshes. We successfully employ the frameworks in medical dynamics simulations, especially resulting meshes in dynamic simulations. However, it when simulation accuracy is required. should be noted that our framework does not include parallel versions of the PBD and XPBD approaches, which Case D In this case study, we utilized the mesh con- laidmviatsntoaugreasbiinlitoyutrocdoemmpoanrasttirvaeteepxopteerni mtiaelnptse.rfIonrmthaenfcueversion method described in Section 3.3 to project a 3D ture, we would like to implement parallel algorithms for human body model onto a 2D plane. This projection en- each approach and incorporate an implicit FEM solver abled us to conduct an elastodynamic experiment simu- for comparing the stability of existing methods. lating a simplified scenario where a human leg is trapped Recently, deep learning-based FEM approaches [23, 24, under building debris. For this study, we used the follow- 25, 26] have emerged and demonstrated superior accuing parameters: Young’s modulus of 1.85 MPa, Poisson’s racy and performance compared to classical FEM methratio of 0.3, total mass of 70 kg, and time step of 0.0001 ods, with some implementations achieving real-time pers. We employed the Neo-Hookean energy model and formance. However, our literature review reveals that Explicit FEM framework as our solver. Through visual- these approaches have primarily been validated in lowerization techniques, we were able to illustrate the stress dimensional spaces, with most demonstrations limited distribution on the leg when compressed by the debris. to 1D and a few in 2D scenarios. While authors suggest This approach demonstrates the practical application of that their methods can be extended to 3D applications, our framework in simulating medical scenarios relevant the performance and interactive capabilities in such sceto disaster response. narios remain unverified. In contrast, all energy-based
This work was supported by Innovative Science and Technology Initiative for Security Grant Number JPJ004596, ATLA, Japan. methods examined in our comparative study have been proven efective and stable in 3D implementations, which is not yet the case for deep learning-based approaches. This observation partly motivates our choice of methods for comparison. As future work, we plan to investigate the stability of deep learning-based FEM methods and evaluate their efectiveness in 3D scenarios, particularly for XR-oriented medical dynamics applications. mentation and production practicalities, in: ACM SIGGRAPH 2020 Courses, 2020, pp. 1–182. [22] N. Akkiraju, H. Edelsbrunner, M. Facello, P. Fu, Alpha shapes: definition and software, volume 63, 1995. [23] J. Jung, K. Yoon, P.-S. Lee, Deep learned ifnite elements, Computer Methods in Applied Mechanics and Engineering 372 (2020) 113401. URL: https://www.sciencedirect.com/ science/article/pii/S0045782520305867. doi:https: //doi.org/10.1016/j.cma.2020.113401. [24] R. Phellan, B. Hachem, J. Clin, J.-M. Mac-Thiong, L. Duong, Real-time biomechanics using the finite element method and machine learning: Review and perspective, Medical Physics 48 (2021) 7–18. [25] R. E. Meethal, A. Kodakkal, M. Khalil, A. Ghantasala, B. Obst, K.-U. Bletzinger, R. Wüchner, Finite element method-enhanced neural network for forward and inverse problems, Advanced Modeling and Simulation in Engineering Sciences 10 (2023) 6. [26] M. Lienen, S. Günnemann, Learning the dynamics of physical systems from sparse observations with finite element networks, arXiv preprint arXiv:2203.08852 (2022).