Survey of Finite Element Method-Based Real-Time Simulations

: The ﬁnite element method (FEM) has deservedly gained the reputation of the most powerful, highly e ﬃ cient, and versatile numerical method in the ﬁeld of structural analysis. Though typical application of FE programs implies the so-called “o ﬀ -line” computations, the rapid pace of hardware development over the past couple of decades was the major impetus for numerous researchers to consider the possibility of real-time simulation based on FE models. Limitations of available hardware components in various phases of developments demanded remarkable innovativeness in the quest for suitable solutions to the challenge. Di ﬀ erent approaches have been proposed depending on the demands of the speciﬁc ﬁeld of application. Though it is still a relatively young ﬁeld of work in global terms, an immense amount of work has already been done calling for a representative survey. This paper aims to provide such a survey, which of course cannot be exhaustive.


Introduction
Model building belongs to basic but most important engineering activities in which various options are explored to reach the solution that offers a compelling compromise between the computational effort and reached accuracy of performed simulations.This procedure has progressed from rather simple sketches in the sand, via manual computations and graphical solutions performed on the paper by means of simplified models, to high fidelity solutions performed on modern computers based on models of various levels of complexity and involving significant computational effort.Engineers combine the available knowledge of physical processes with innovation and creativity to shape the world in ways that are sometimes breathtaking.The modern engineering tools, hardware, and software are the key factors that enable highly demanding computations that engineers could only dream of just a couple of decades ago.The power of modern hardware components allows to analyze different physical processes at practically any level of detail.
Focusing on the field of structural analysis, the finite element method (FEM) has imposed itself as currently the most powerful tool that enables highly efficient modeling and simulation of structures characterized by complex geometry and exposed to arbitrary boundary and initial conditions.As a numerical method, it substantially relies on the usage of modern hardware.In typical engineering tasks, the FEA computations are performed as "off-line" computations.This practically implies a sequenced use of all three essential parts of a FEA program-a model is prepared in the preprocessor; the solver provides the solution of the problem with all required quantities; finally, the postprocessor offers tools to observe and analyze the obtained solution.Hence, setting external excitations (done in the preprocessor) and observing the structural response (done in the postprocessor) are two processes separated by a substantial amount of time.
On the other hand, the remarkable pace of hardware development in the past couple of decades was the crucial impetus for a number of researchers to raise the question about the possibility of real-time simulations based on FE models.Such an intriguing option opens up the door to a whole new virtual world, thus offering unforeseen possibilities in various fields of work ranging from engineering tasks, via training simulators, to entertainment applications.It is no wonder that it has attracted a significant number of researchers to provide their ideas on how to make FEM-based real-time simulation feasible.
The authors believe that, although still relatively young, the field has reached a level of maturity that deserves a not necessarily exhaustive but representative overview of the work conducted in the field, with the primary idea to address various directions of development.Prior to that, the term "real-time simulation" will be addressed together with various means to achieve it.

Real-Time Simulation and Means to Achieve It
In different fields, different interpretations of the term "real-time simulation" may be encountered, and, hence, it is worthwhile to address this aspect at the very beginning.According to the standard DIN ISO/IEC 2382 [1], the term denotes such an operation of a computing system that provides the results within a prescribed period of time.In contrast to the most frequent, intuitive understanding of the term, the processing time can actually differ from real (i.e., wall clock time_.The demanded processing time could be longer, but also in certain cases even shorter than real time, and this primarily depends on the application requirements.The intuitive understanding of the term would imply that the time needed for computation and presentation of the results (usually on a screen) is the same as the wall clock time, or shorter since the process can be put on idle in order to get synchronized with the wall clock time.However, it should be noticed that in many applications the simulation would still serve the purpose even if it performed somewhat slower than wall clock time.
Speaking of FEM-based real-time simulation of deformable structures, there are several means to influence the speed of simulation:

•
Use of powerful computer components, including the possibility of changing the strategy how the computer components are used; • Faster algorithms, code optimization, and particularly parallelization of computation; • Development of various formalisms and formulations for structural analysis.
Although the influence of each of the above mentioned means to speed up the computation is relatively obvious, the authors still find a parallel with the objective of trying to set a faster time in a racing car around a circuit rather neat and illustrative.In this parallel, better hardware obviously means providing the driver with a faster car.This is a kind of "brute force" approach and, despite the Moore's law of hardware development [2], it certainly has its limits with the current technology.Software optimization would be equivalent to improving the driver's skills.Now, the question is where the analogy with the third group of means exactly resides.This one is actually most interesting.The analogy implies that this group of means bends the classical rules so that the circuit start and finish line remain the same, but not necessarily the circuit itself (i.e., one may construct his/her own circuit).If one could come up with a shorter circuit containing less/milder curves that still leads from the same start (input data) to the same finish line (results), then this could easily turn out to be the decisive advantage.This survey will not consider the developments belonging to the software optimization.Regarding hardware development, only one important and novel aspect will be mentioned at this point, without detailed discussion.Namely, despite the fact that the hardware development resulted in striking computational power, a very promising approach actually consists in adjusting the way of how the hardware is used.Instead of the standard approach that implies the use of central processing unit (CPU) for computation, the tendency is to use graphics processing units (GPU) for this purpose.The reason is to be sought in the fact that modern GPUs contain a number of stream processors that offer massive computational power [3,4].This approach demands to resolve problems such as transferring the data necessary for computation to the GPU [5].The focus of this survey will be on the third group.The objective is to give an overview of different directions stemming from various ideas how to approach the problem and the resulting art of performing a suitable trade-off between the simulation accuracy and numerical efficiency.The proposed ways to reduce the computational power, the associated advantages and drawbacks, together with possible fields of application will be discussed.

Mass-Spring Systems
Mass-spring systems represent relatively simple models, in which concentrated masses are interconnected by a set of springs.This heuristic approach actually implies that a 3D (continuum) body is replaced by a truss structure.The idea is relatively old and was proposed by Hrennikoff in 1941 [6], who; however, never aimed at real-time simulation but simply at a computational method applicable in the field of structural mechanics.The obvious advantage of the approach resides in the rather simple formulation of the spring (truss) element that transfers the force only in the length direction, thus significantly reducing the number of degrees of freedom.Additionally, if the geometrically nonlinear effects are to be accounted for, it is done in a straightforward manner as the force is always transferred between the current positions of the spring nodes.Figure 1 shows a few illustrative screen-shots from an interactive simulation involving a simple cloth-like structure formed of a mass-spring system and exposed to large, geometrically nonlinear deformations [7].Hence, the major advantages of the approach are the simplicity of formulation and computational efficiency.Those properties rendered the approach extremely attractive in the 1990s and the beginning of the new millennium as the limited power of available hardware components at that time called for such solutions.In that period, the mass-springs systems were kind of a workhorse for real-time simulations in a number of fields, ranging from entertainment industry up to biomechanics and surgery simulators.The focus of this survey will be on the third group.The objective is to give an overview of different directions stemming from various ideas how to approach the problem and the resulting art of performing a suitable trade-off between the simulation accuracy and numerical efficiency.The proposed ways to reduce the computational power, the associated advantages and drawbacks, together with possible fields of application will be discussed.

Mass-Spring Systems
Mass-spring systems represent relatively simple models, in which concentrated masses are interconnected by a set of springs.This heuristic approach actually implies that a 3D (continuum) body is replaced by a truss structure.The idea is relatively old and was proposed by Hrennikoff in 1941 [6], who; however, never aimed at real-time simulation but simply at a computational method applicable in the field of structural mechanics.The obvious advantage of the approach resides in the rather simple formulation of the spring (truss) element that transfers the force only in the length direction, thus significantly reducing the number of degrees of freedom.Additionally, if the geometrically nonlinear effects are to be accounted for, it is done in a straightforward manner as the force is always transferred between the current positions of the spring nodes.Figure 1 shows a few illustrative screen-shots from an interactive simulation involving a simple cloth-like structure formed of a mass-spring system and exposed to large, geometrically nonlinear deformations [7].Hence, the major advantages of the approach are the simplicity of formulation and computational efficiency.Those properties rendered the approach extremely attractive in the 1990s and the beginning of the new millennium as the limited power of available hardware components at that time called for such solutions.In that period, the mass-springs systems were kind of a workhorse for real-time simulations in a number of fields, ranging from entertainment industry up to biomechanics and surgery simulators.Mass-springs systems found their application in numerous medical applications for simulation of soft tissue deformations.This is basically due to the fact that the demand for accuracy in this type of application is reduced to "plausible physical behavior."Hence, in most of such applications, the simulation simply needs to appear realistic to human eye and that kind of behavior can easily be reproduced by mass-spring models even in geometrically nonlinear regimes, despite the significant simplifications implemented.Models with all their specifics were proposed in the field of orthodontics [8], robotic surgery [9], real-time simulation of muscles [10], and various surgery simulations [11][12][13], some of which were proposed with GPU accelerated computations [14].
In virtual reality simulators and entertainment products, the use of mass-spring systems is particularly versatile.For instance, they come in handy in modeling very thin structures whose physical behavior is dominated by the membrane stiffness, while the bending stiffness is practically negligible.A typical example would be cloth.O'Connor and Stevens [15] dedicated their work to find proper spring stiffness that enables realistic real-time simulation of collision between a cloth and other objects.In order to achieve enhanced behavior, beside classical springs that provide structural basis of the cloth, they introduced shear and bend springs, which in combination with a specific pattern of interconnections produce more realistic cloth behavior.Another interesting example of mass-spring models application is in the simulation of facial expressions.For this purpose, Zhang et Mass-springs systems found their application in numerous medical applications for simulation of soft tissue deformations.This is basically due to the fact that the demand for accuracy in this type of application is reduced to "plausible physical behavior."Hence, in most of such applications, the simulation simply needs to appear realistic to human eye and that kind of behavior can easily be reproduced by mass-spring models even in geometrically nonlinear regimes, despite the significant simplifications implemented.Models with all their specifics were proposed in the field of orthodontics [8], robotic surgery [9], real-time simulation of muscles [10], and various surgery simulations [11][12][13], some of which were proposed with GPU accelerated computations [14].
In virtual reality simulators and entertainment products, the use of mass-spring systems is particularly versatile.For instance, they come in handy in modeling very thin structures whose physical behavior is dominated by the membrane stiffness, while the bending stiffness is practically negligible.A typical example would be cloth.O'Connor and Stevens [15] dedicated their work to find proper spring stiffness that enables realistic real-time simulation of collision between a cloth and other objects.In order to achieve enhanced behavior, beside classical springs that provide structural basis of the cloth, they introduced shear and bend springs, which in combination with a specific pattern of interconnections produce more realistic cloth behavior.Another interesting example of mass-spring models application is in the simulation of facial expressions.For this purpose, Zhang et al. [16] proposed to model facial skin and a set of facial muscles by means of mass-spring systems.Their models incorporate materially nonlinear springs.
Distinctive works in this field addressed the problem of ambiguity of model building.Namely, the major problem with the mass-spring models is how to determine the distribution of point masses and the topology of their connections so as to have acceptable description of the actual physical behavior.There are obviously a large number of possibilities for any considered physical object.Bianchi et al. [17] proposed the idea to optimize the topology of interconnections by means of genetic algorithms based on the results obtained using a known reference model.They assumed constant mass and stiffness distribution in order to reduce the problem to topology identification only.While in this work the reference model was just a rather fine mass-spring system, in a further work by Bianchi et al. [18] the approach was extended so as to use the results by the continuum-based finite element (FE) model as a reference in order to identify the stiffness distribution besides the topology.The approach can be successful for deformations in the linear domain, but the complexity of deformations in the geometrically nonlinear regime significantly deteriorates success of the approach.A very important conclusion by the authors is that homogenous stiffness parameters (i.e., the same stiffness of all springs) fail to properly simulate linear elastic material and one has to distinguish at least between diagonal and straight springs in the mesh when setting their stiffness.Bourguignon and Cani [19] proposed a method for controlling anisotropy in mass-spring models.To do that, they actually used additional elements only to determine forces that act onto the nodes of the mass-spring model.Additionally, they suggested the use of "extra forces" in the simulation that would enable constant volume deformations.

Model Order Reduction Techniques
FE-models are characterized by sufficiently fine discretization of geometry that leads to a high-dimensional system of equations.Solving the resulting system of equations is actually the most demanding and time-consuming part of the overall FEM solution.This fact was the motivation for researchers to develop ideas for performing model order reduction by finding a low-dimensional but still adequate approximation of the initial system of equations.
The overall motion in multi-body systems (MBS) is represented as a superposition of rigid-body motion and small deformational motion that is conveniently described in a local (body-fixed) coordinate system.The consideration of deformation in MBS dynamics together with all types of nonlinearities, multi-physics (coupled-field) effects, etc., can be successfully done by means of MBS-FEM co-simulation [20][21][22][23].The great fidelity of the approach is certainly its undisputable advantage, but it comes at the price of an extreme computational effort combined with the additionally time-consuming data transfer between the two systems.To reduce the effort and therewith significantly improve the simulation speed, model order reduction is extensively used in MBS dynamics involving elastic structures.It should be emphasized that real-time simulation is usually not a primary requirement in this field.Nevertheless, computational efficiency is of great importance in MBS systems, and since they are mostly used for engineering tasks, the accuracy is a demand of equal or even greater importance.
Among different approaches to model order reduction, the one typically used in the MBS systems is based on modal vectors and coordinates.The essence of the idea is to perform a great deal of computation prior to simulation.In a step performed prior to simulation, orthogonal mode shapes are determined and they become the degrees of freedom further used to determine the deformational part of the motion.In this manner, FE models with immensely large numbers of degrees of freedom (~10 6 ) could be reduced down to a few dozens of degrees of freedom (DOFs) in modal space.The fact that orthogonal mode shapes are used leads to further reduction of computational effort as the resulting equations describing deformation are decoupled, provided modal damping coefficients are used.In fact, this general idea is already known from FEM, where the modal superposition technique is often used to determine structural transient response [24].In FEM the so-called normal modes are used, which are determined by performing modal analysis for the known kinematic boundary conditions.They represent the patterns of vibration of a system at specific frequencies (eigenfrequency).Figure 2 shows screen-shots from an interactive simulation of a car rear-axle.The original (full extent) FE model contains 300,000 DOFs, while the model reduced in modal space only 10 DOFs and easily enables real-time simulation with suitable accuracy in the lower frequency range (up to the frequency of the highest mode used in the simulation).For improved versatility, MBS systems typically use the component mode synthesis (CMS) technique, particularly the Craig-Bampton method [25].The idea is to divide the DOFs of the elastic structure into boundary and interior DOFs.The boundary DOFs belong to those nodes of the FE model that are kept in the simulation model for the purpose of interaction with the surrounding (i.e., for the purpose of kinematic and dynamic boundary conditions).In the next step two sets of modes are determined: constraint modes and fixed boundary normal modes.The constraint modes are static modes and their number corresponds to the number of boundary DOFs.A constraint mode is obtained by assigning a unit displacement to one of the boundary DOFs while all remaining boundary DOFs are fixed.The fixed-boundary normal modes are obtained by classical modal analysis whereby all boundary DOFs are fixed.Such a combination of modes does not make an orthogonal set of modes.In order to keep the numerical advantages of a decoupled system, the soobtained set of modes is orthogonalized prior to simulation.However, the disadvantage of the orthogonalization is in the fact that the obtained orthogonal modes lack physical interpretation, which, for instance, prevents assignment of (experimentally determined) modal damping coefficients.
Regardless of the approach used to obtain the mode shapes, they always characterize the structure in the undeformed (i.e., initial configuration).Hence, the method is intrinsically meant for linear analysis and it is; therefore, suitable for small deformations.However, different approaches were proposed recently in order to extent the approach at least into the realm of moderate geometric nonlinearities.The approach based on the introduction of the geometric stiffness matrix was first introduced in the commercial software package SIMPACK [26].In that approach it is essential to differentiate between the forces which easily produce relatively large displacements but small strains (i.e., the structure is rather flexible with respect to them) and the forces that give rise to relatively small deformations but large stresses in the structure [27].Furthermore, the geometric stiffness matrix is computed for each force belonging to the latter group and assuming its magnitude equal to one.Over the course of simulation, the geometric stiffness matrix for each of the forces from the second group is scaled by using the current magnitude of the force.The overall stiffness matrix is obtained by superposition of the linear stiffness matrix and the current values of all geometric stiffness matrices.Thus, the introduced geometric stiffness matrix is force-scaled.In this manner, the exquisite numerical efficiency is retained, while the stress stiffening effects are included for a carefully selected group of forces.Marinkovic and Zehn [28,29] proposed to introduce a geometric stiffness for each mode shape.It is computed prior to simulation using geometrically nonlinear analysis in which the structure is deformed according to the mode shape.In doing that, the intensity of deformation corresponds to the modal coefficient equal to one.Hence, during the simulation, the modal geometric stiffness coefficients are updated by using the current values of the modal coefficients.In this case, the geometric stiffness matrix is deformation-scaled.
Frequently the geometric nonlinearities are caused by finite local rotations of structural subdomains.Techniques denoted as modal warping were proposed as a remedy that would again extend the applicability of modal superposition to moderate geometrically nonlinear deformations.In order to keep the numerical advantages of a decoupled system, the so-obtained set of modes is orthogonalized prior to simulation.However, the disadvantage of the orthogonalization is in the fact that the obtained orthogonal modes lack physical interpretation, which, for instance, prevents assignment of (experimentally determined) modal damping coefficients.
Regardless of the approach used to obtain the mode shapes, they always characterize the structure in the undeformed (i.e., initial configuration).Hence, the method is intrinsically meant for linear analysis and it is; therefore, suitable for small deformations.However, different approaches were proposed recently in order to extent the approach at least into the realm of moderate geometric nonlinearities.The approach based on the introduction of the geometric stiffness matrix was first introduced in the commercial software package SIMPACK [26].In that approach it is essential to differentiate between the forces which easily produce relatively large displacements but small strains (i.e., the structure is rather flexible with respect to them) and the forces that give rise to relatively small deformations but large stresses in the structure [27].Furthermore, the geometric stiffness matrix is computed for each force belonging to the latter group and assuming its magnitude equal to one.Over the course of simulation, the geometric stiffness matrix for each of the forces from the second group is scaled by using the current magnitude of the force.The overall stiffness matrix is obtained by superposition of the linear stiffness matrix and the current values of all geometric stiffness matrices.Thus, the introduced geometric stiffness matrix is force-scaled.In this manner, the exquisite numerical efficiency is retained, while the stress stiffening effects are included for a carefully selected group of forces.Marinkovic and Zehn [28,29] proposed to introduce a geometric stiffness for each mode shape.It is computed prior to simulation using geometrically nonlinear analysis in which the structure is deformed according to the mode shape.In doing that, the intensity of deformation corresponds to the modal coefficient equal to one.Hence, during the simulation, the modal geometric stiffness coefficients are updated by using the current values of the modal coefficients.In this case, the geometric stiffness matrix is deformation-scaled.
Frequently the geometric nonlinearities are caused by finite local rotations of structural subdomains.Techniques denoted as modal warping were proposed as a remedy that would again extend the applicability of modal superposition to moderate geometrically nonlinear deformations.Choi and Ko [30] keep the track of local rotations that occur during the deformation computed by the linear model.In the next step, they warp the computed modal basis according to the local rotations of the FE nodes.Similarly, Marinkovic and Zehn [29,31] first identify substructures that can perform relatively large local rotations compared to the rest of the structure.They compute the overall structural response by means of linear mode shapes.Based on the deformed state, the average local rotations of substructures are determined and used to rotate ("warp") the nodal displacements previously computed in the classical manner using linear mode shapes.The improvement in accuracy of determining the deformed shape was demonstrated in term of numbers for a simple cantilever beam and a car rear-axle.It should be emphasized that the procedure of "warping" the displacements introduces a rather small additional numerical effort so that the numerical efficiency of the reduced model based on modal space remains extremely high.Another approach can be applied in cases where the nonlinear behavior is restricted to a relatively small part of a structure.In such a case, as suggested by Zehn and Marinkovic [32], model order reduction can be performed for the large portion of the structure that can be adequately described by a linear model, whereas the small part characterized by nonlinearities remains a full-extent (not reduced) FE model with all necessary nonlinearities (geometrical and/or material) accounted for.
The idea of using modal space based solutions in order to speed up the simulation and reach interactive frame rates was also successfully applied in computer graphics.Upon its introduction to the field by Pentland and Williams [33], it was versatilely used for a number of tasks, including interactive simulation of deformable models with constraints such as collision [34] simulation of movements of trees exposed to turbulent wind [35], bone-based character animation for interactive simulations [36], etc.

FEM-Based Training of Neural-Networks
A rather specific direction of the development that could be interpreted as an "indirect" model order reduction technique relies on the use of FEM results to train the neural networks for fast computations.Additionally, in this case, the idea is to perform rather cumbersome and time-consuming computations prior to simulation.The fact that the approach allows to tackle the problem of nonlinear structural deformations, whereby both geometrical and material nonlinearities can be addressed, makes it particularly interesting.Hence, the (linear or nonlinear) FEM computations provide a set of results used in the training phase of neural networks.
The approach was applied by Hambli et al. [37] for real-time simulation of tennis racket and ball deformation.As input data for the artificial neural network, they used a set of randomly generated parameters, such as the impact angle, ball impact velocity and impact zone, while the deformation of the ball and strings together with the impact force made the output data set.This development aimed at a virtual reality simulator including a haptic device with force-feedback.Ordaz-Hernandez et al. [38] used a multi-layer feed-forward neural network to produce a reduced model for geometrically nonlinear deformations of a cantilever beam.Ćojbašić et al. [39] reported on a similar work for shell structures.Torano et al. [40] used fuzzy logic, neural networks and three-dimensional (3D) finite element calculations to develop a numerical model that allows fast predictions of the response of longwall coal mining installation to changing operation conditions.The results were further represented to the user by means of virtual reality modelling language (VRML) tools.Morooka et al. [41] presented a method for simulating the deformation of organ models that relies on usage of neural networks to find a relation between external forced and model deformations.They confirmed the possibility of real-time simulations whereby the results were comparable with those yielded by the nonlinear FEM.Tzong-Ming et al. [42] used Ansys for FE computations and Matlab to perform neural training.Neural networks memorized only surface displacements of the models and the parametric deformation matrices were converted into "behavioral modules" for the virtual reality engine.
Systems for real-time control can particularly benefit from strongly reduced models that offer adequate accuracy in simulating the system to be controlled.Such models can be used to improve the predictor part of the controller as Kalkkuhl et al. [43] demonstrated by implementing the FEM-based neural-network approach for nonlinear longitudinal dynamics of an experimental lorry.Runge et al. [44] used neural networks to capture the finest nonlinearities based on FE and needed for simulation and control of modular soft robots.

Full Extent FE Models in Real-Time Applications
Though keeping FE models in their full extent (all nodal DOFs kept, which means no reduction performed) offers great fidelity in modeling, it is, on the other hand, a rather challenging task to perform real-time simulations based on such models.Even with a relatively moderate number of DOFs, the resulting systems of equations are usually prohibitively large for a successful real-time simulation based on rigorous FE formulations and using conventional hardware components.This is particularly valid if nonlinearities are involved in the simulation.Hence, various ideas on how to simplify the way of formulating the problem and finding the solution and thus perform a reasonable trade-off between accuracy and numerical efficiency were proposed in order to meet the objective of real-time simulation.The solutions depend on the specific field of application.
In the early days of the FEM-based virtual reality and augmented reality (AR) concept, the hardware components could offer significantly less computational power compared to what is available at present.Therefore, at that time, if the full extent FE models were to be used, it was necessary to restrict their use to linear regimes as the assumption of linearity allows for significant savings in numerical effort.However, this approach remained until the present days in applications in which linear models do meet the simulation requirements.Huang et al. [45] aimed at enhancements in structural analysis by means of AR technologies.For this purpose, they developed a system that integrates sensor measurement and real-time FEA simulation into an AR-based environment.A network of wireless sensors acquired spatially distributed loads, while FEA generated fast solutions for the sensed loads.The speed of computation was enabled by introducing the assumption of linear and quasi-static behavior, which allowed to pre-compute the inverse of the stiffness matrix.Finally, the obtained FEA results were superimposed directly onto real-world objects for enhanced data exploration.A rather similar approach was seen in the work by Fiorentino et al. [46].Additionally, in this work the general idea was to visualize the FE results by overlaying them over the real model.The user may interactively define the boundary conditions on the observed model and based on the FE results assess the stress/strain distribution.A cantilever test was used for demonstration.Cerracchio et al. [47] used the linear inverse FEM to perform real time full-field reconstruction of displacements based on limited set of discrete strain data.To obtain the strain data, in situ strain measurements of the deformed structural shape were done.They demonstrated the approach on a composite stiffened structures for static, dynamic and thermal loadings.
Linear FE models were also seen as a trade-off between numerical efficiency and biomechanical realism in the field of surgery simulators that demand modeling soft tissues.In one of the pioneer works in 90's, Bro-Nielsen and Cotin [48] used linear elasticity combined with simple tetrahedral elements.They further reduced the numerical effort by performing condensation so that only the surface nodes remained during the simulation available for interaction.In addition, the condensed stiffness matrix was inverted before the simulation.Cotin et al. [49] proposed physical models based on FEM and linear elasticity, in which precomputed deformations were combined to enable real-time simulations.Pincinbono et al. [50] developed several biomechanical models based on linear elasticity and FEM with anisotropic mechanical behavior included.They implemented transversely isotropic properties and added external elastic membrane to the models to represent relatively stiff skin of internal organs ("surface anisotropy").Furthermore, their development addressed the problem of contact with an instrument and determination of reaction forces for the need of force-feedback devices and haptic rendering.The attractiveness of this field of application in this period is reflected in the fact that survey papers presenting developments in real-time simulation for surgery simulation were already published at the end of 1990s [51] and the beginning of 2000s [52,53].
Inclusion of nonlinear effects in real-time simulations is actually of crucial importance in many applications, but it also implies considerably larger numerical effort, in particular when full extent FE models are used.Commercially available FEA programs mostly implement the total or updated Lagrangian formulations [24] as standard solutions with the basic difference between them in the choice of the reference configuration.Those formulations are characterized by high accuracy, but also by significant numerical effort and convergence issues.Székely et al. [54] applied total Lagrangian formulation in combination with material nonlinearities in modeling soft tissue deformations for laparoscopic surgery simulation.They addressed the high computational demand of the approach by using a "brute force" method (i.e., by parallelizing the computation on a large, three-dimensional network of high-performance processor units).Despite of the fact that they report on a not highly realistic simulation of uterus tissue deformation.On the other hand, various ideas were implemented to keep the numerical effort in acceptable limits despite the use of rigorous FE formulations.For instance, Heng et al. [55] proposed an approach denoted as a hybrid condensed FE model.The idea is to partition the model into two regions-a part of the model that is being interacted with and the rest of the model.Hence, a complex nonlinear FE model that can also deal with topological change was used to model a small-scale FE model that represents the operational region, while a linear and topology-fixed FE model was used to model the large-scale nonoperational region.The development was done for a virtual training system for knee arthroscopic surgery.In another field of application, Dulong et al. [56] focused their work on development of real-time interaction between a designer and a virtual prototype as a promising way to perform optimization of parts design.In the first step, which they refer to as the "training phase," they performed a set of nonlinear FE computations.Over the course of simulation, deformations for a current load case is determined by interpolating between the results obtained in the "training phase."One may actually notice here a certain parallel with the approach based on neural networks.
The complexity of simulations based on rigorous nonlinear FE computations drove the researchers to propose innovative formulations for nonlinear FE computations but also to consider different options for solvers.The consideration of solvers is closely related to time integration schemes.A number of authors emphasize the low computational effort needed to resolve a time-step (decoupled equations, provided a lumped mass matrix is used) and no need for iterations due to nonlinearities as major advantages of an explicit time-integration scheme [57,58].However, those advantages come together with the disadvantage of a limited and typically rather short time-step due to the conditional stability of the method.Soft materials, such as tissues often involved in real-time simulations, would allow for somewhat larger time-steps as a consequence of their low stiffness.On the other hand, high internal damping is also intrinsic for such materials and this is a parameter that further decreases the time-step.Since real-time simulations are of primary interest for processes that are observable to human perception (i.e., relatively slow dynamics is observed), an adequate time-step of numerical integration could easily be in the order of magnitude of 10 −2 s.Such a large time-step calls for implicit time integration schemes that are unconditionally stable [59,60].This implies that the system of equations to be resolved remains coupled, which invokes significant numerical effort.The standard solution in commercially available FEA programs is a direct solver.It offers a great stability of solution procedure for various types of structures and engineering accuracy, but lacks options for a trade-off between accuracy and computational speed.Iterative solvers, on the other hand, are rather advantageous with respect to the last aspect, particularly if used with adequately selected preconditioners to improve the convergence rate [61].One may simply limit the number of iterations and thus obtain a solution that is not necessarily very accurate but still sufficiently good for the purpose of a real-time simulation [62].In this manner, the several orders of magnitude larger time-step used in an implicit time integration scheme can easily make up for a larger numerical effort needed to resolve the coupled system of equations compared to the numerical effort needed for a time-step within an explicit time integration scheme.
More interesting is the innovativeness in proposing nonlinear FE formulations for real-time simulations.In doing that, the aim is to reduce the numerical effort by neglecting certain nonlinear effects and thus simplify the rigorous FE formulations.The developed formulations are supposed to combine the advantages of linear and geometrically nonlinear FE formulations.In other words, they should be stable and as efficient as possible, while covering the geometrically nonlinear effects to a large extent.Probably the most interesting approach is based on the co-rotational FEM formulation [60,63,64].The basic idea is to account for large local rigid-body rotations by means of local coordinate frames that are attached to the structure and perform the same rigid-body rotation as the local area of the structure.The usual simplification implies that the local deformations are considered to be linear with respect to the local coordinate frames.This means that the influence of local deformations and stress state ("stress stiffening") onto the tangent stiffness matrix are neglected.The stiffness matrices of the local regions are precomputed and only rotated during the simulation based on the current and initial positions of the local coordinate frames.Capell et al. [65] presented a framework for the skeleton-driven animation of elastically deformable characters.They associated regions of FE meshes with particular bones.Those regions were structural sub-domains whose rigid-body rotations were accounted for by using local coordinate frames, while the FE computations were linearized with respect to them.A similar approach, but extended in order to provide finer resolution of accounting for local rotations was proposed by Müller et al. [66].They implemented local coordinate frames at the nodes and compute a tensor field at every time-step to describe the local rotations of all the vertices in the mesh.Etzmuss et al. [67] presented development for viscoelastic, highly flexible surfaces.They attached local coordinate frames to triangular elements to cover local rigid-body rotations on the element level.Some researchers extended the approach to 3D elements (tetrahedrons) and also referred to this formulation as the "warped stiffness" approach [68].Marinkovic and Zehn analyzed the accuracy in terms of numbers provided by this simplified co-rotational formulation for both solid [69] and shell [70] elements, and showed that if nonlinearities in deformational behavior are dominated by the local rotations, the accuracy for predicted displacements lies within a few percent of the results by rigorous geometrically nonlinear FE formulations.They also demonstrated that FE models of up to several thousand elements are applicable in real-time simulations on conventional hardware components without parallelization and presented this development for the field of surgery simulation [71].
Since several thousand elements cannot provide an appealing representation of complex geometries, a coupled mesh technique has been proposed [66,69,71,72].The idea is to use a relatively coarse FE mesh for computations, while a rather fine surface mesh is used to represent the actual geometry.Coupling of the two meshes is done by assigning the surface vertices of the geometry mesh to the closest finite element in the initial configuration.In most cases, a vertex is actually positioned within an element.If this is not the case, the minimum absolute value of the sum of all three local vertex coordinates with respect to each finite element is decisive.Figures 3-5 show a human head, heart, and a dog model, respectively, that implement the coupled mesh technique and the deformations during an interactive real-time simulation based on the above addressed simplified co-rotational FE formulation.
to the closest finite element in the initial configuration.In most cases, a vertex is actually positioned within an element.If this is not the case, the minimum absolute value of the sum of all three local vertex coordinates with respect to each finite element is decisive.Figures 3-5 show a human head, heart, and a dog model, respectively, that implement the coupled mesh technique and the deformations during an interactive real-time simulation based on the above addressed simplified corotational FE formulation.(a)

Conclusions
FEM-based real-time simulation has entered numerous fields ranging from entertainment to engineering.Although we are witnessing the early days of its development, an immense amount of work has already been done to respond to the challenges of performing real-time simulations using currently available computer hardware technology.This paper aimed at a not exhaustive but representative survey of the development in the field.It would certainly take a prohibitively large space to mention all or most of the work done and, hence, the objective was to represent different

Conclusions
FEM-based real-time simulation has entered numerous fields ranging from entertainment to engineering.Although we are witnessing the early days of its development, an immense amount of work has already been done to respond to the challenges of performing real-time simulations using currently available computer hardware technology.This paper aimed at a not exhaustive but representative survey of the development in the field.It would certainly take a prohibitively large space to mention all or most of the work done and, hence, the objective was to represent different directions and ideas on how to balance between conflicting objectives of high numerical efficiency and acceptable accuracy required for real-time simulations.
The solution in the form of mass-spring systems obviously offers attractive simplicity and efficiency, which come at the price of rather problematic accuracy and ambiguity of model building.Though considered obsolete today, the approach still represents an option in the fields where the accuracy is not a critical demand and can also be used for a rough approximation of deformable environment of the objects the focus of interaction is put on.
Model order reduction techniques have already proven to be extremely useful in the field of MBS dynamics for engineering tasks, but also found their applications in various interactive simulations in the field of computer graphics.Different attempts have been made to tackle the most serious shortage of the approach, namely its intrinsic limitation to linear deformations.Implementation of the geometric stiffness and modal warping were mentioned as moderate remedies, but interestingly, so far, no report could be found that combines both approaches, which should offer a new quality due to their complementariness.FEM-based training of neural networks was addressed as a sort of indirect approach to model reduction.It is a promising approach that may address nonlinearities in structural behavior as well, but its success strongly depends on the "density" of the deformation set used in training the neural networks.
Finally, the most daring solutions keep the FE models in their full extent and, additionally, involve nonlinear deformational behavior.The hardware developments in the past two decades encouraged the researchers to explore new ideas in order to address this particular challenge.Currently, the approach based on the simplified co-rotational FE formulation, or "warped stiffness" as denoted by some authors, in combination with the coupled mesh technique appears to be rather promising for applications on conventional hardware components.The formulation implements simplifications that allow to use the precomputed linear stiffness matrices of each element, while the coupled mesh technique enables the use of relatively coarse FE meshes for complex geometries.
With the hindsight on development done, it is reasonable to expect that new brilliant ideas are probably just around the corner.In addition, having in mind the Moore's law of hardware development and the already existing new strategies of hardware usage that bring significant increase in computational power, it is certain that the field of FEM-based real-time simulations has a very promising future in the field of structural analysis.In addition, it should be emphasized that related fields, such as computational fluid dynamic (CFD) and particle physics, encounter similar challenges and the interest in real-time simulations by means of conventional and state-of-the-art hardware components in those fields is equally large.The developments in all those fields already began to merge in order to enable real-time simulation of coupled-field effects, which represents a quantum leap toward high-fidelity simulation.

Figure 1 .
Figure1.Cloth-like structure as a mass-spring system in an interactive simulation[7].

Figure 1 .
Figure1.Cloth-like structure as a mass-spring system in an interactive simulation[7].
Appl.Sci.2019, 9, 2775 5 of 15 Appl.Sci.2019, 9, x 5 of 15 (eigenfrequency).Figure2shows screen-shots from an interactive simulation of a car rear-axle.The original (full extent) FE model contains 300,000 DOFs, while the model reduced in modal space only 10 DOFs and easily enables real-time simulation with suitable accuracy in the lower frequency range (up to the frequency of the highest mode used in the simulation).

Figure 2 .
Figure 2. Modal space reduced finite element (FE) model: Screen-shots from an interactive simulation of a car rear-axle.

Figure 2 .
Figure 2. Modal space reduced finite element (FE) model: Screen-shots from an interactive simulation of a car rear-axle.For improved versatility, MBS systems typically use the component mode synthesis (CMS) technique, particularly the Craig-Bampton method[25].The idea is to divide the DOFs of the elastic structure into boundary and interior DOFs.The boundary DOFs belong to those nodes of the FE model that are kept in the simulation model for the purpose of interaction with the surrounding (i.e., for the purpose of kinematic and dynamic boundary conditions).In the next step two sets of modes are determined: constraint modes and fixed boundary normal modes.The constraint modes are static modes and their number corresponds to the number of boundary DOFs.A constraint mode is obtained by assigning a unit displacement to one of the boundary DOFs while all remaining boundary DOFs are fixed.The fixed-boundary normal modes are obtained by classical modal analysis whereby all boundary DOFs are fixed.Such a combination of modes does not make an orthogonal set of modes.In order to keep the numerical advantages of a decoupled system, the so-obtained set of modes is orthogonalized prior to simulation.However, the disadvantage of the orthogonalization is in the fact that the obtained orthogonal modes lack physical interpretation, which, for instance, prevents assignment of (experimentally determined) modal damping coefficients.Regardless of the approach used to obtain the mode shapes, they always characterize the structure in the undeformed (i.e., initial configuration).Hence, the method is intrinsically meant for linear analysis and it is; therefore, suitable for small deformations.However, different approaches were proposed recently in order to extent the approach at least into the realm of moderate geometric nonlinearities.The approach based on the introduction of the geometric stiffness matrix was first introduced in the commercial software package SIMPACK[26].In that approach it is essential to differentiate between the forces which easily produce relatively large displacements but small strains (i.e., the structure is rather flexible with respect to them) and the forces that give rise to relatively small deformations but large stresses in the structure[27].Furthermore, the geometric stiffness matrix is computed for each force belonging to the latter group and assuming its magnitude equal to one.Over the course of simulation, the geometric stiffness matrix for each of the forces from the second group is scaled by using the current magnitude of the force.The overall stiffness matrix is obtained by superposition of the linear stiffness matrix and the current values of all geometric stiffness matrices.Thus, the introduced geometric stiffness matrix is force-scaled.In this manner, the exquisite numerical efficiency is retained, while the stress stiffening effects are included for a carefully selected group of forces.Marinkovic and Zehn[28,29] proposed to introduce a geometric stiffness for each mode shape.It is computed prior to simulation using geometrically nonlinear analysis in which the structure is deformed according to the mode shape.In doing that, the intensity of deformation corresponds to the modal coefficient equal to one.Hence, during the simulation, the modal geometric stiffness coefficients are updated by using the current values of the modal coefficients.In this case, the geometric stiffness matrix is deformation-scaled.

Figure 3 . 15 Figure 3 .Figure 4 .
Figure 3. Human head model: (a) Geometry and the surface representation of the FE mesh; (b) screenshots from an interactive simulation [72].