LiToTac: An Interactive-Interface Software f or Finite Element Analysis of Multiple Contact Dynamics

: In order to investigate the mechanical behavior of systems with complex architecture and a large number of contacting bodies, a ﬁnite element software, named LiToTac, has been developed by using the object-oriented programming technique. This software, with an interactive graphical user interface, is able to handle highly non-linear problems including multiple contacts and large deformation. More importantly, the contact detection based on a hybrid three-stages methodology can be performed automatically, which is more efﬁcient than the common strategies of pre-deﬁning contact zones in commercial FEM software like ANSYS, ABAQUS, etc. In addition, the contact solver in LiToTac is portable between dynamic and quasi-static codes and can accurately solve contact coupled with friction in a reduced system. Several numerical examples are carried out to illustrate the functionality and capacity of the software package.


Introduction
Contact/impact problem is a very difficult issue in computational mechanics and many engineering fields. The main difficulty comes from the high nonlinearities on the contact surface and the discontinuous nature of motion. Over the years, numerical solutions of such problems accounting for large deformation effects have been intensively studied [Mahmoud, El-Shafei, Abdelrahman et al. (2013); Mlika, Renard and Chouly (2017); Hartmann, Oliver, Weyler et al. (2009) ;Chen, Joli and Feng (2015)]. Consequently, systems composed of only a few interacting bodies may require significant computational efforts. More complex scenarios such as soft textile composite reinforcements [Misra, Dixit and Mali (2014)] and biomechanical systems [Bei and Fregly (2004) ;Lin, Walter, Banks et al. (2010)] are particularly challenging and prone to very long simulating time.
In general, the nonlinear contact system of equations in commercial software can either be solved by the penalty methods or the multiplier Lagrangian methods. As is known, the penalty methods are easy to be implemented, but the contact boundary conditions and friction laws can not be accurately satisfied. Besides, the value of penalty factor is problem-dependent, and depends particularly strongly on the stiffness ratio between the contact objects, which may result in spurious oscillations in contact forces. The multiplier Lagrangian methods are able to enforce the zero-penetration condition exactly by introducing a set of multipliers representing contact forces. Nevertheless, both the generalized coordinates and multipliers are unknown values, thus leading to an increase in the size of the equation systems. In LiToTac, the techniques of bi-potential [De Saxcé and Feng (1998)] are adopted, which shows several favorable properties as follows: (1) The reaction force appearing as the projection of a linear combination of the reaction with the relative velocity onto Coulomb's cones, can lead to a unique variational inequality rather than two minimum principles, as compared to classic augmented Lagrangian methods [Jean and Touzot (1988)].
(2) The unique mathematical operator of projection can can make a certain mathematic reduction, because it does not consider the Lagrangian multipliers as additional unknowns.
(3) The nonlinearities of equations are separated to be solved, which can overcome the computational complexity resulting from the non-differentiable contact bipotential and the inequalities of frictional contact rules. (4) The whole solution process requires only one user-defined parameter: the friction coefficient. This enables users to obtain very good results even though they are not familiar with the contact mechanics theory. The article is organized as follows: in Section 2, the object-oriented approach to design the software is presented. The structure of the developed software and the simulation methods are given in Section 3. Section 4 describes the main functionalities of the software. Section 5 presents several examples to highlight the graphical representation of the numerical solutions. In Section 6, a few concluding remarks are drawn.
2 Object-oriented programming (OOP) LiToTac, as a prototype finite element software of the CAE family, employs the OOP approach for the designing process. Several crucial components have been incorporated into the software using an OOP architecture, these include an interactive graphical user interface, standard mechanical elements and material libraries to assist in the modeling, postprocessing for the manipulation and display of results. The code is implemented with careful software-engineering procedures, management of data, re-use of the programming strategy. The basic features of object-oriented programming include data abstraction, encapsulation, data-hiding, modularity, hierarchy, inheritance, and polymorphism, etc, allowing the finite element code to be easily extended for implementing new ideas and new algorithms. The detailed explanation of these key concepts can be referred to the articles [Fenves (1990); Dubois-Pèlerin and Pegon (1998)]. In this section, we will only briefly present them as  CMES, vol.118, no.1, pp.111-137, 2019 they relate to our own program.

Figure 1: Base abstraction classes of OMT platform
One main characteristic of OOP is to define various classes. A class incorporates the set of objects that shares the same structure and behavior. An object is called an instance of the class, which can be treated as the realization of the data type defined by the class. In our code, the object-oriented architecture is organized around several important base classes of OMT (Objected and Methodological Technology), as illustrated in Fig. 1. The main functions of these classes are summarized as follows: • OmtProductBase: this base class is in charge of managing the version, name, registered modules of a software developed on OMT platform.
• ModuleOmtBase: the user project is organized by splitting the problem into different modules, while each module derives from the base abstract class ModuleOmt-Base.
• OmtDataDef: this base class is used to define specific data class in the user modules, such as the boundary condition, velocity, load, etc.
• OmtSolveSetup: this base class is generated to create different setups for different solver engines.
By using the inheritance concepts, new classes for contact system codes can be directly created from the existing base class. In other words, the derived classes inherit the data and member functions of the base class. And at the same time, some additional methods and attributes relating to it can be implemented. As is shown in Fig. 2, we register the new developed software as "ProductLiToTac". Then some version information about the current developed software would be provided to the OMT platform according to the virtual methods in OmtProductBase class. ModuleContactDetection, ModuleInitVelocity, ModuleLimits, ModuleLoad, ModuleMesh, ModuleMaterial, ModuleSolution are the professional concrete modules derived from the base abstract class ModuleOmtBase. Each module, as its name explains, corresponds to a necessary function for contact dynamic simulation. These modules are created individually, and inherit the basic function members like InitalModule(), ReadModule(), GetModuleActions(), etc from the OmtProductBase class. As a result, they can independently manage their own data and behaviors. It is clear that the modularization by splitting one project into several small modules can be a big advantage of OOP in large scale software design, because each programmer just needs to take charge of one part of them, then the total design time can be significantly reduced. Also, one can easily extend the software to multiple suites by plugging in other functional modules. For instance, ModuleInitialVelocity is in charge of adding initial velocities to nodes before performing time loops. It manages the related date from the geometry and from the interactive interface. If someone wants to use acceleration as the initial condition, ModuleInitalVelocity can be deleted, and then be replaced with a ModuleInitialAcceleration. These kinds of operations do not affect the running of other modules and also have no influence on the final simulation results. The encapsulation of OOP is a way of defining private data members and function members that are hidden from view outside of the object's definition. That means only the object's own method can directly have access to manipulate its field. In our code, the ModuleMesh class has two private data: NumNodes and NumElements. The value of these two items can be obtained by inputting data operation defined by the ModuleMesh, but they are not allowed to be recalled or modified by the methods of other modules, even though these modules have to gather information from mesh geometry. By allowing the developers to limit the inter-dependencies between software components, the encapsulation mechanism can enormously help reduce the system complexity and increase robustness. Polymorphism, which means a variety of forms", is the fundamental technique of OOP. Polymorphism is generally implemented by virtual member functions, which allow redefinition in the subclasses. This characteristic enables objects to respond to messages differently for the same function. Here, we take the operation of saving data function as an example to illustrate this concept in our code. It is very common to click the "save project" button in the menu items when one wants to reuse the current data or simulation operations for the next time. By this operation, the OMT manager will firstly create a project file, and then send messages to each module to call the function "SaveProductData". SaveProduct-Data is a virtual function in the base class of ModuleOmtBase, while its subclasses override this method too. As a result, although all the modules are told to save ProductData, the concrete module only saves its own data inside the project file. Also, through polymorphism, the ModuleMesh can draw any geometry components (selected nodes, elements, surface sets, etc by the operation on the interface) in run-time without pretreatment of geometrical types before compiling. This is so-called dynamic polymorphism, which ensures multiple methods can be implemented by one connector.
3 General structure of the software and computational methods  Once the preprocessing and the contact detection are finished, an input file will be produced for the contact detection engine and solver engine. Within a time step, the contact detection engine based on a hybrid methodology consisting of Octree and bounding volume hierarchy (BVH), can calculate the contact information (actual contact area, nodes, penetration, etc) and transform these data to the finite element solver engine. The Octree and BVH need to be updated based on the current coordinates, thus they have to consider the solved displacements for each time step. The solver engine includes three analysis types: static solver, implicit dynamic solver, and semi-explicit dynamic solver. The final computed results are written in an output file and displayed by the post-processor.

Automatic contact detection
Contact detection is usually considered as a prerequisite for solving contact forces. It can be found in applications of virtual reality, game development, and robotics, etc. In the literature, the range of contact detection algorithms proposed for the finite element modeling is various. These include methods based on the slave-master surface [Hallquist, Goudreau and Benson (1985)], the nested bucket [Benson and Hallquist (1990)], the sweep and prune [Cohen, Lin, Manocha et al. (1995)], the closest features [Gilbert, Johnson and Keerthi (1988); Cameron (1997); Lin and Canny (1991)], the LC-grid [Chen, Lei and Zang (2014)], etc. For survey articles on contact detection, the reader is referred to [Jiménez, Thomas and Torras (2001); Lin and Gottschalk (1998)].
In the multiple contact systems with a large number of bodies, contact detection becomes quite complicated. For one thing, bodies within the scene can move and deform unpredictably from the previous loading steps. For another thing, highly frequent occurrences of contact events can significantly increase the computational cost. The most widely used methods in CAE software such as ANSYS and ABAQUS, are based on slave-master approaches, in which the user should pre-define contact zones. This may result in a heavy burden for users to find potential contact nodes and surfaces in order to establish a list of contact zones. Moreover, for systems with very complex geometry architecture, contact zones may be inside the body, which is generally impossible to be predefined by simple operations like picking nodes or elements from the interface. In our code, these problems can be solved by performing automatic contact detection based on a three-stages methodology. The main idea behind each stage will be briefly introduced in the following subsections.  The top-stage inspections implemented in ModuleContactDetection is based on a pointerbased Octree data structure [Wilhelms and Van Gelder (1992)], which can quickly enumerate pairs of objects that may potentially collide. The Octree (Fig. 4) represents the 3D volume as a hierarchy of discrete octants, in which each parent octant is recursively subdivided into eight child octants. The octant can be added or removed as required, thus is well suited to dynamic problems that update frequently. The definition of a typical Octree node is shown in Fig. 5. Each octant stores pointers to its eight children. But only these non-empty octants are considered to take memory, while others are defined Null. This allows all the search paths to be tested with a lower memory 119 storage. The recursive location of objects terminates when the branch reaches the tree depth, or the number of objects per octant below the predefined value. The Octree structure is searched from the root down to determine collisions. Only objects that share the same octants are taken to have potential contacts. Note: each object in an Octree hierarchy should be firstly culled by an approximated bounding box for a rough test. To locate the real contact nodes or elements in the following stages, each object is also decomposed of a finite element mesh.

Middle stage for interactions between two FEM meshes
After the filtering of the first stage, the middle stage lies in how to precisely detect the contact elementary pairs (e.g., node-segment(2D), node-triangle(3D), node-quad(3D)) be-tween two meshes. Accordingly, the technique of bounding volume hierarchies (BVH) is adopted for primitive queries. Types of BVs include spheres, axis aligned bounding boxes (AABBs), discretely oriented polytopes (K-DOPs), or a hybrid of them.
In this work, we employ the popular hierarchy of AABBs as the fundamental BVHs, because they provide a good trade-off between tightness of fit and computation cost. We create bounding volumes based on specific geometric components, i.e., vertices, edges and faces, see Fig. 6(b). This strategy can largely reduce duplications of elementary queries compared to traditional methods by merely using triangle elements ( Fig. 6(a)). For example, the vertex A ( Fig. 6(a)) is incident to eight triangles (t 1 ",t 8 ), and A comes into contact with triangle t 9 , which will produce eight times triangle-triangle tests, and 24 times vertex-triangle tests. But in Fig. 6(b), each bounding volume enclosing a vertex is represented only once in the hierarchy, thus the overlapping between vertex A and triangle t 9 would be unique. Furthermore, the separate hierarchies can improve the culling efficiency of elementary testing, because the volume of a vertex can be much smaller than any other type of geometric component.  Figure 7: (1) Triangle-based BVH and Vertex-based BVH, (2) Hierarchy traversal between two meshes A triangle-based hierarchy or a vertex-based hierarchy can be built in a typical top-down manner [Bergen (1997)], see Fig. 7(1). It starts out by enclosing the hull of all finite element primitives with a tight Axis-aligned bounding box (AABB). Then these primitives are successfully sorted into two smaller subsets from one level to a deeper lever, thus producing a binary tree hierarchy. The recursion stops when the leaf nodes of the tree consist of only a single primitive. For the constructing process, the most important part is to find a suitable splitting plane to partition the primitives. To ensure a balanced tree, we position the splitting plane orthogonal to the longest AABB axis, which should pass through the median of the centroid coordinates. Any primitive is classified depending on its location with respect to the splitting plane. In Fig. 7(2), all the contacting elementary pairs between object 0 and a can be searched by a recursive traversal algorithm. The two trees are traversed down from the root, and recursively descend to deeper levels. The primitive contact tests are finally performed between the leaf BV 4 and BV d. This tree search complexity is in logarithmic order, which is particularly efficient f or c ases w ith very f ew c ontact e lements, but a very large data system. Note that, for finite element contact calculation, the simplest linear triangle element can be replaced by other more complicated types in the context of finite element analysis, e.g., 6-node triangle, 4-node quad, 8-node quad, etc.

Bottom stage for local calculation
Finally, the bottom stage functions to discern whether these pairs detected by middle stage are really in penetrated condition. Consider a surface ϕ, any point on the surface can be presented by a parametric form: S(ε, η) = ∑ N i (ε, η)x i , where N i denotes the Lagrangian shape function associated with node position x i . Let point y be the potential counterpart of the contact surface ϕ as illustrated in Fig. 8. Therefore, the closest point from surface ϕ to an arbitrary point y can be obtained after solving the following system of equations: where J(l (k) ) = The iterations usually converge in four or less iterations. The final solution (ε * , η * ) should satisfy the conditions: −1 ≤ ε * ≤ 1, −1 ≤ η * ≤ 1. Herein, the outward normal vector of point (ε * , η * ) is determined by the cross product of r ε and r η The gap vector is given by

LiToTac: An Interactive-Interface Software for Finite Element Analysis
where ,η ) . The problem stated by Eq.
(1) may be nonlinear when the surface ϕ is curved, for example, the quadrilateral face. Therefore, the solution can be obtained by the Newton-Raphson iterative process:

Finite element computation
The solver in ModuleSolution is based on the finite element method for multiple contact dynamic analysis. Typically, the finite element formulation of this problem in the discrete form can be written as: where F ext the vector of external loads, F int the internal forces, R the contact reaction vector, M the mass matrix, A the damping matrix,u the velocity vector andü the acceleration vector. It is worth noting that the internal forces F int are derived from the physics-based deformable model related to the stiffness effect. The variables such as velocity and stress inside the continuum field are not known. The shape interpolating function N is needed to approximate them by the nodal values of the field variable. The shape functions N are generally written in the form of polynomials depending on the element's shape. If the displacement vector u in an element Ω e is approximated by u = N · u e , where u e represents the displacements at nodes. Then the mass matrix M e , the damping matrix A e can be defined as: where ρ e is the mass density, κ e the damping parameter, V e the volume of Ω e . The solution to Eq. (5) differs between explicit and implicit approaches. The explicit method appears to be efficient, but unstable, and without checking convergence. The implicit method is supposed to be accurate. Nevertheless, the computation is much more expen-sive for large deformations. Currently, our code offers three solvers for different types of analysis: (1) the implicit solver [Feng, Joli, Cros et al. (2005); Feng, Magnain and Cros (2006)] to address low-velocity or quasi static cases; (2) the semiexplicit solver [Peng, Feng and Joli (2018)] to handle high-velocity impact problems; (3) the static solver [Feng, Hjiaj, De Saxcè et al. (2006)] to deal with static frictional cases.

Computation of contact forces
Our contact solver is based on the bi-potential method [De Saxcé and Feng (1998)], in which a formulation extended by the augmented Lagrangian method is provided. Considering a local reference frame (normal vector n and tangential vector τ τ τ), there exists an implicit relationship between the constraint displacement x and the contact reaction forces r: x = Wr +x (7) CMES, vol.118, no.1, pp.111-137, 2019 where W andx represent respectively the global compliance matrix and the free displacement. For each contact point α, among N c instantaneous contact points, the relationship can be written as: where W αβ is an influence matrix that takes account of the coupling between contact points α and β . With regard to point α, Eq. (8) can be treated by considering other contact points (α = β ) as "frozen".
In the work of [De Saxcé and Feng (1998)], the contact bi-potential is formulated as follows: where ℜ − is the set of negative and null real numbers, K µ the so-called Coulomb cone. ∪ K µ (r) denotes the indicator function of the closed convex set K µ .
In order to avoid non-differentiable potentials, it is convenient to use the Augmented La-grangian method [De Saxcé and Feng (1991)], Eq. (9) can be equal to and where r * is the so-called augmented contact forces, and ρ is a positive real parameter, µ the friction coefficient. Proj K µ (r * ) means that r is the projection of r * onto the closed convex Coulomb cone. The Uzawa technique can be used to solve the implicit Eq. ( 8), which leads to a predictioncorrection procedure as: where the correction is explicitly carried out with respect to three possible contact statues as: if

124
Copyright © 2019 Tech Science Press Figure 9: GUI of LiToTac After the prediction-correction procedure converges, the global contact forces R can be obtained by where H(u) is the mapping matrix from the local frame to the global frame.

Functionalities of the software
Qt is an efficient application framework for developing multi-platform applications, by which the code can be run on various systems and hardware platforms with little or even no change. In the current project, Qt has been largely used to design the graphical user interface (GUI) objects such as Dialog, Menu, Toolbar, Label, Icon, etc. Also, to render the scene with high fidelity, the OpenGL library is applied to draw graphic primitives like points, lines, surfaces, volumes, etc. As mentioned, the software LiToTac is an integrated environment which consists of four main parts: the preprocessor, the contact detector, the solver, and the post-processor. The GUI of the program is shown in Fig. 9. Some main characteristics of LiToTac are presented as follows.

Preprocessor
• Mesh: to allow the user to make certain operations on the model geometry such as creating, deleting, modifying the sets of primitives (nodes, elements, surfaces, etc). , vol.118, no.1, pp.111-137, 2019 In addition, it helps to render the view environment based on user's habits. The user can achieve some special effects by defining light, texture, transparency, color, and so on.

CMES
• Material: to set the parameters for different material models. Several material models including linear elasticity, and non-linear hyperelasticity have been implemented in our code. Particularly, users are allowed to quickly add their own material models into the library in very few steps, which provides better convenience over general commercial softwares. • Create input files: after finishing the procedures as described above, an input file with only necessary information is produced. This file can also be generated for supporting ANSYS or other commercial software.

Contact detection
• The contact detection is performed automatically to gather the basic data for calculating contact forces. The construction of Octrees and bounding volume hierarchies should rely on the current coordinates provided by the initial mesh as well as the motion displacements. • For some particular materials, different areas within a body may produce several different frictional coefficients when interacting with another kind of material. This may add complexity to the computation of frictional force and the management of data. In our software, the user can well handle this issue by defining frictional coefficients between two kinds of groups, wherein, one is the node groups, the other is the element edge groups. The value of coefficients will be recorded in a 1D array. Consequently, once a contact detection query is solved, each of the resulting contact pairs can have a individual frictional coefficient based on the node ID and its counterpart element ID. This methodology allows users to make further researches on inhomogeneous friction and anisotropic friction problems.

Solver
The solver consists of the following tasks: • to read the input file from the preprocessor.
• to exchange data with the contact detection engine.
• to solve contact problems via implicit, semi-explicit, and static methods.
• to perform hourglass control for semi-explicit simulation, when one point Gauss integration is applied to accelerate the computing speed. • to create output files of the solutions for the purpose of post-processing.

Post-processor
The roles of the post-processor can be summarized as follows: • to read the result file produced by the solver engine.
• to import the stress analysis solutions. • to display the structure of a solid body by scanning.
In order to test and validate the functions and the capability of LiToTac, several specific examples have been carried out. It is noted that the following analysis are performed on a PC (i7-8550U, 1.80 GHz). , vol.118, no.1, pp.111-137, 2019 5.1 Example 1: 2D quasi-static contact simulation A B C The first model consists of 16 hyperelastic bodies (spheres) and 4 rigid bodies (rectangles), as is shown in the left side of Fig. 10. The highlighted yellow points on the rigid bodies are fixed, while the left rectangle body is imposed by constant force loading. The plots of Von-Mises stress at the final step are presented in the right side of Fig. 10. This example is a very special quasi-static case, which can not be modeled by a traditional static solution. Because several spheres inside the system are enforced without displacement limits, but only with contact constraints, which leads to a singularity of the global stiffness matrix. It is also difficult to apply a general dynamic solver to perform this simulation. Since quasistatic simulations requires extremely accurate and robust computation of dynamic contact forces. Whereas, most dynamic solvers in commercial softwares may produce numerical oscillations when taking account of impact effects. By using the implicit method [Feng, Joli, Cros et al. (2005)], LiToTac can well address this case. As is shown in the left side of Fig. 11, the hyperelastic bodies can deform stably in very low velocity, since it generally exists no obvious oscillations in the evolution of Von-Mises versus time. • to create reports about the energy including total energy, kinematic energy, and elastic energy with change of time. By means of the powerful post-processing capabilities of LiToTac, users can import the results of different situations and plot the reports together to view the influences by applying different computing algorithms. For example, an analysis report of contact searching is given in the right side of Fig. 11. The result shows the computational time of the automatic contact detection (Curve1) is significantly lower than the brute force (Curve2). Therefore, the efficiency of automatic contact detection is confirmed. Additionally, it is worthy noting a multi-window function of LiToTac. In Fig. 9 or Fig. 11, the mesh window and plots window are able to be viewed at the same time. This function can provide great convenience for users when they need to gather information from different views.

LiToTac: An Interactive-Interface Software for Finite Element Analysis
This 3D example is a typical "1+6" structure [Ghoreishi, Davies, Cartraud et al. (2007)] composed of seven hyperelastic fibers. Large deformation and sliding are analyzed by using the static solver. The total torsion angle is θ = 95 • , the bottom surface is fixed, and the displacements of Z direction on the top surface are set to zero. 25 load steps are performed for this problem, so a pure torsion of 3 • 48 is applied to each step. The initial mesh and the final state of distribution of Von-Mises stress are shown in Fig. 12.
From the obtained results, it can be seen that the stress distribution is symmetrical to the center point. The vector quantity of displacements and reaction forces at Step 20 are given in Fig. 13. Fig. 14 depicts the contour plots of contact forces at real contacting nodes at different loading steps, which allows us to observe the inter fiber-fiber interactions. Numerical modeling of fiber systems still remains a very difficult issue in FEM field. The inter contacts between fibers may become unpredicted because of its complicated architecture and large number of bodies. Example 2 presented here aims at illustrating the capacity of LiToTac, thus no detailed analysis of mechanical behavior is performed. The functionalities mentioned such as multi-windows, cut plane of inter structure, animation of real contact areas, view plotting of contact force contour, etc are quite convenient, which permits LiToTac to handle complicated fiber-fiber interacting simulations with much less operations than commercial software.

Example 3: 2D comparison between LiToTac and commercial software
To illustrate the computational efficiency of LiToTac, we consider here one example to make a comparison with the commercial software ANSYS and ABAQUS. This example is inspired by the work of Wriggers et al. [Wriggers, Van and Stein (1990)], which involves a hyperelastic cylinder impacting upon two oblique rigid symmetric surfaces, see Fig. 15. During the pro-cess, there exists no damping except for Coulomb friction between contact surfaces. The Yeoh material model is considered to describe the hyperelasticity behavior. Firstly, the distribution of Von-Mises stress contour at time t = 0.002 s is displayed in Fig. 16. From the results, we can observe that the maximum value of Von-Mises stress in LiToTac is slightly lower than ANSYS and ABAQUS. In Tab. 1, N1, N2, N3 are three selected nodes that interpenetrate with surface AD at time t = 0.002 s. As is shown, the contact penetrations δ of these nodes in LiToTac can be controlled in the order of 10 −16 , while the value just approaches to 10 -5 in ANSYS, and 10 -6 in ABAQUS. This difference indicates that LiToTac can be more accurate than ANSYS and ABAQUS to satisfy the contact impenetrability condition. The performance of computational time is reported in Tab. 2. Through the result, we can find that LiToTac possesses a small performance difference with ABAQUS, but is much faster than ANSYS.  It is also interesting to compare the energy evolution of the two softwares. Fig. 17 plots the total energy as a function of time in the case of frictionless contact. We can see that the total energy for these three softwares is conserved within an acceptable range of error (less than 0.5%). However, the change of energy curves confirms a better robustness of LiToTac over ANSYS and ABAQUS.

Example 4: 3D comparison between LiToTac and commercial software
As a comparison, a 3D example consisting of two deformable blocks is also analyzed here, see  Fig. 19 shows the plots of the Von-Mises stress of Node A versus time. We can observe that the curves produced by LiToTac and ABAQUS are close to each other, and are more robust than the case of ANSYS. Fig. 20 indicates the evolution of total energy as a function of time. The results confirms that LiToTac allows a better conservation of energy than ABAQUS for 3D frictionless contact problem. Whereas, for ANSYS, an increase of total energy can be produced because of numerical instability.

LiToTac: An Interactive-Interface Software for Finite Element Analysis
To investigate the frictional effects on the algorithm performance, we change the frictional coefficient to µ = 0 .5. Fig. 21 gives the plots of whole kinetic energy versus t ime. In L-iToTac and ABAQUS, the kinetic energy decreases as expected by frictional contact effects after time t = 0.05 when the impact occurs. However, in ANSYS, a jump takes places at time t = 0.065. The performance of the applied algorithms in each software is summarized in Tab. 3. The Lagrange-multiplier method in ABAQUS is replaced by penalty method because of its non convergence in dealing with this frictional case. Though the results obtained by computing a 3D example, we can see that LiToTac is faster than ANSYS, but a little bit slower than ABAQUS. First-order implicit+Bi-potential 9.3 9.6

Conclusion
In this paper, we have presented a new finite element software for the modeling and analysis of multiple contact dynamics. This software is able to handle strong non-linearities including multiple contacts and large deformations, and it can also perform automatic contact detections without predefining contact zones. The program has been designed by using the object-oriented principles and the OMT technology. These techniques allow us to simplify the architecture of the program and to incorporate quite easily new features to other specific suites. Through several complicated numerical examples, the new function of automatic contact detection confirms a good convenience over slave-master approaches. Also, by comparing with the commercial software ANSYS and ABAQUS, the bi-potential techniques in LiToTac are demonstrated to possess higher accuracy and better robustness than penalty and Lagrange multiplier algorithms. The software presented can be further extended to study local phenomenon occurring in polymer chains, and DNA, or in the domain of computer graphics to improve the rendering of hairs. It does not need to completely redefine the software architecture. Some optimization methods can also be added to perform parallel computation of large-scale problems.