On the approximation in the smoothed finite element method (SFEM)

This letter aims at resolving the issues raised in the recent short communication (Int. J. Numer. Meth. Engng 2008; 76 (8):1285–1295. DOI: 10.1002/nme.2460 ) and answered by (Int. J. Numer. Meth. Engng 2009; DOI: 10.1002/nme.2587 ) by proposing a systematic approximation scheme based on non‐mapped shape functions, which both allows to fully exploit the unique advantages of the smoothed finite element method (SFEM) (Comput. Mech. 2007; 39 (6):859–877. DOI: 10.1007/s00466‐006‐0075‐4 ; Commun. Numer. Meth. Engng 2009; 25 (1):19–34. DOI: 10.1002/cnm.1098 ; Int. J. Numer. Meth. Engng 2007; 71 (8):902–930; Comput. Meth. Appl. Mech. Engng 2008; 198 (2):165–177. DOI: 10.1016/j.cma.2008.05.029 ; Comput. Meth. Appl. Mech. Engng 2007; submitted; Int. J. Numer. Meth. Engng 2008; 74 (2):175–208. DOI: 10.1002/nme.2146 ; Comput. Meth. Appl. Mech. Engng 2008; 197 (13–16):1184–1203. DOI: 10.1016/j.cma.2007.10.008 ) and resolve the existence, linearity and positivity deficiencies pointed out in (Int. J. Numer. Meth. Engng 2008; 76 (8):1285–1295). We show that Wachspress interpolants (A Rational Basis for Function Approximation. Academic Press, Inc.: New York, 1975) computed in the physical coordinate system are very well suited to the SFEM, especially when elements are heavily distorted (obtuse interior angles). The proposed approximation leads to results that are almost identical to those of the SFEM initially proposed in (Comput. Mech. 2007; 39 (6):859–877. DOI: 10.1007/s00466‐006‐0075‐4 ). These results suggest that the proposed approximation scheme forms a strong and rigorous basis for the construction of SFEMs. Copyright © 2009 John Wiley & Sons, Ltd.     

Assessment of the cost and accuracy of the generalized FEM

In this paper, we address the cost versus accuracy capabilities for the generalized FEM (GFEM) which was developed in (Comput. Methods Appl. Mech. Eng. 2003; 192 :3109–3161, Int. J. Numer. Meth. Engng 2004; 60 :1639–1672, Ph.D. Thesis, Texas A&M University, College Station, TX, August 2003 (Advisor: T. Strouboulis)), and also the construction of two‐sided a posteriori error estimates, which can be used to assess the achieved accuracy at all levels of the method. Copyright © 2006 John Wiley & Sons, Ltd.     

On the essence and the evaluation of the shape functions for the smoothed finite element method (SFEM)

This paper is written in response to the recently published paper (Int. J. Numer. Meth. Engng 2008; 76 :1285–1295) at IJNME entitled ‘On the smoothed finite element method’ (SFEM) by Zhang HH, Liu SJ, Li LX. In this paper we

• (1) repeat briefly the important essence of the original SFEM presented in (Comp. Mech. 2007; 39 : 859–877; Int. J. Numer. Meth. Engng 2007; 71 :902–930; Int. J. Numer. Meth. Engng 2008; 74 :175–208; Finite Elem. Anal. Des. 2007; 43 :847–860; J. Sound Vib. 2007; 301 :803–820), and
• (2) examine further issues in the evaluation of the shape functions used in the SFEM.

It will be shown that the ‘SFEM’ presented in paper (Int. J. Numer. Meth. Engng 2008; 76 :1285–1295) is not at all our original SFEM presented in (Comp. Mech. 2007; 39 :859–877; Int. J. Numer. Meth. Engng 2007; 71 :902–930; Int. J. Numer. Meth. Engng 2008; 74 :175–208; Finite Elem. Anal. Des. 2007; 43 :847–860; J. Sound Vib. 2007; 301 :803–820). Therefore, all these ‘Theorems’, ‘Corollaries’ and ‘Remarks’ presented in paper (Int. J. Numer. Meth. Engng 2008; 76 :1285–1295) have nothing to do with our original SFEM. The properties of the original SFEM stand as they were presented in our original papers (Comp. Mech. 2007; 39 :859–877; Int. J. Numer. Meth. Engng 2007; 71 :902–930; Int. J. Numer. Meth. Engng 2008; 74 :175–208; Finite Elem. Anal. Des. 2007; 43 :847–860; J. Sound Vib. 2007; 301 :803–820). Finally, we brief on our advancements made far beyond our original SFEM and our visions on future numerical methods. Copyright © 2009 John Wiley & Sons, Ltd.     

Computational method of inverse elastostatics for anisotropic hyperelastic solids

The paper presents a computational method for predicting the initial geometry of a finitely deforming anisotropic elastic body from a given deformed state. The method is imperative for a class of problem in stress analysis, particularly in biomechanical applications. While the basic idea has been established elsewhere Comput. Methods Appl. Mech. Eng. 1996; 136 :47–57; Int. J. Numer. Meth. Engng 1998; 43 : 821–838), the implementation in general anisotropic solids is not a trivial exercise, but comes after a systematic development of Eulerian representations of constitutive equations. In this paper, we discuss the general representation in the context of fibrous hyperelastic solids, and provide explicit stress functions for some commonly used soft tissue models including the Fung model and the Holzapfel model. A three‐field mixed formulation is introduced to enforce quasi‐incompressibility constraints. The practical utility of this method is demonstrated using an example of aneurysm stress analysis. Copyright © 2006 John Wiley & Sons, Ltd.     

A comparative study between two smoothing strategies for the simulation of contact with large sliding

The numerical simulation of contact problems is still a delicate matter especially when large transformations are involved. In that case, relative large slidings can occur between contact surfaces and the discretization error induced by usual finite elements may not be satisfactory. In particular, usual elements lead to a facetization of the contact surface, meaning an unavoidable discontinuity of the normal vector to this surface. Uncertainty over the precision of the results, irregularity of the displacement of the contact nodes and even numerical oscillations of contact reaction force may result of such discontinuity. Among the existing methods for tackling such issue, one may consider mortar elements (Fischer and Wriggers, Comput Methods Appl Mech Eng 195:5020–5036, 2006; McDevitt and Laursen, Int J Numer Methods Eng 48:1525–1547, 2000; Puso and Laursen, Comput Methods Appl Mech Eng 93:601–629, 2004), smoothing of the contact surfaces with additional geometrical entity (B-splines or NURBS) (Belytschko et al., Int J Numer Methods Eng 55:101–125, 2002; Kikuchi, Penalty/finite element approximations of a class of unilateral contact problems. Penalty method and finite element method, ASME, New York, 1982; Legrand, Modèles de prediction de l’interaction rotor/stator dans un moteur d’avion Thèse de doctorat. PhD thesis, École Centrale de Nantes, Nantes, 2005; Muñoz, Comput Methods Appl Mech Eng 197:979–993, 2008; Wriggers and Krstulovic-Opara, J Appl Math Mech (ZAMM) 80:77–80, 2000) and, the use of isogeometric analysis (Temizer et al., Comput Methods Appl Mech Eng 200:1100–1112, 2011; Hughes et al., Comput Methods Appl Mech Eng 194:4135–4195, 2005; de Lorenzis et al., Int J Numer Meth Eng, in press, 2011). In the present paper, we focus on these last two methods which are combined with a finite element code using the bi-potential method for contact management (Feng et al., Comput Mech 36:375–383, 2005). A comparative study focusing on the pros and cons of each method regarding geometrical precision and numerical stability for contact solution is proposed. The scope of this study is limited to 2D contact problems for which we consider several types of finite elements. Test cases are given in order to illustrate this comparative study.     

High‐quality surface remeshing using harmonic maps—Part II: Surfaces with high genus and of large aspect ratio

This paper follows a previous one that was dealing with high‐quality surface remeshing using harmonic maps (Int. J. Numer. Meth. Engng 2010; 83 :403–425). In (Int. J. Numer. Meth. Engng 2010; 83 :403–425), it has been demonstrated that harmonic parametrizations can be used as input for surface meshers to produce high‐quality triangulations. However, two important limitations were pointed out, namely surfaces with high genus and/or of large aspect ratio. This paper addresses those two issues. We first develop a multiscale version of the harmonic parametrization of (Int. J. Numer. Meth. Engng 2010; 83 :403–425) and then combine it with a multilevel partitioning algorithm to come up with an automatic remeshing algorithm that overcomes the above‐mentioned limitations of harmonic maps. The overall procedure is implemented in the open‐source mesh generator Gmsh (Int. J. Numer. Meth. Engng 2009; 79 (11): 1309–1331). Copyright © 2011 John Wiley & Sons, Ltd.     

On the numerical integration of an advanced Gurson model

This article is focused on a new extended version of Gurson's model (J. Eng. Mater. Technol. 1977; 99 :2–15), its numerical integration scheme and its consistent tangent matrix being within an FE code. First, this new advanced Gurson model is proposed, which is an extension of the original to take into account plastic anisotropy and mixed (isotropic+kinematic) hardening. In this paper, only the growth phase of cavities is considered (the nucleation of new voids is ignored). Second, a new numerical algorithm for the integration of this new Gurson model is presented. The algorithm is implicit in all variables and is unconditionally stable. This algorithm is generic and could be used for other anisotropic yield functions and other hardening laws. Third, the consistent tangent matrix is computed in an explicit way by exact linearization of the constitutive equations. To check its efficiency and robustness, the proposed integration algorithm is compared, under some simplified assumptions and choices, with the algorithms of Aravas (Int. J. Numer. Meth. Engng 1987; 24 :1395–1416) and Kojic (Int. J. Numer. Meth. Engng 2002; 53 (12):2701–2720). The performance of the developed consistent modulus, compared to other techniques for the computation of the tangent matrix is assessed. The paper ends with numerical simulations of tensile tests on homogeneous and notched specimens. Copyright © 2010 John Wiley & Sons, Ltd.     

Deflated preconditioned conjugate gradient solvers for linear elasticity

Extensions of deflation techniques previously developed for the Poisson equation to static elasticity are presented. Compared to the (scalar) Poisson equation (J. Comput. Phys. 2008; 227 (24):10196–10208; Int. J. Numer. Meth. Engng 2010; DOI: 10.1002/nme.2932 ; Int. J. Numer. Meth. Biomed. Engng 2010; 26 (1):73–85), the elasticity equations represent a system of equations, giving rise to more complex low‐frequency modes (Multigrid. Elsevier: Amsterdam, 2000). In particular, the straightforward extension from the scalar case does not provide generally satisfactory convergence. However, a simple modification allows to recover the remarkable acceleration in convergence and CPU time reached in the scalar case. Numerous examples and timings are provided in a serial and a parallel context and show the dramatic improvements of up to two orders of magnitude in CPU time for grids with moderate graph depths compared to the non‐deflated version. Furthermore, a monotonic decrease of iterations with increasing subdomains, as well as a remarkable acceleration for very few subdomains are also observed if all the rigid body modes are included. Copyright © 2011 John Wiley & Sons, Ltd.     

Implementation of an algorithm for general material behavior of steel taking interaction between plasticity and transformation‐induced plasticity into account

In this article a semi‐implicit algorithm (predictor–corrector approach) for incorporating the interaction between plasticity and transformation‐induced plasticity (TRIP) in steel is developed. Contrary to the usual elasto‐plasticity, the underlying model of material behavior of steel is far more complex. The interaction between plasticity and TRIP requires extensions of algorithms developed in Doghri (Int. J. Numer. Meth. Engng 1993; 36 :3915–3932) and in Mahnken (Commun. Numer. Meth. Engng 1999; 15 :745–754). A particular feature of the algorithm is that the inner iteration can be reduced to a single scalar equation. Numerical examples illustrate the algorithm's capabilities. Copyright © 2011 John Wiley & Sons, Ltd.     

Tensor‐based Gauss–Jacobi numerical integration for high‐order mass and stiffness matrices

In this work, we choose the points and weights of the Gauss–Jacobi, Gauss–Radau–Jacobi and Gauss–Lobatto–Jacobi quadrature rules that optimize the number of operations for the mass and stiffness matrices of the high‐order finite element method. The procedure is particularly applied to the mass and stiffness matrices using the tensor‐based nodal and modal shape functions given in (Int. J. Numer. Meth. Engng 2007; 71 (5):529–563). For square and hexahedron elements, we show that it is possible to use tensor product of the 1D mass and stiffness matrices for the Poisson and elasticity problem. For the triangular and tetrahedron elements, an analogous analysis given in (Int. J. Numer. Meth. Engng 2005; 63 (2):1530–1558) was considered for the selection of the optimal points and weights for the stiffness matrix coefficients for triangles and mass and stiffness matrices for tetrahedra. Copyright © 2009 John Wiley & Sons, Ltd.     

A fourth‐order compact scheme for the Helmholtz equation: Alpha‐interpolation of FEM and FDM stencils

We propose a fourth‐order compact scheme on structured meshes for the Helmholtz equation given by R(φ):=f( x )+Δφ+ξ²φ=0. The scheme consists of taking the alpha‐interpolation of the Galerkin finite element method and the classical central finite difference method. In 1D, this scheme is identical to the alpha‐interpolation method (J. Comput. Appl. Math. 1982; 8 (1):15–19) and in 2D making the choice α=0.5 we recover the generalized fourth‐order compact Padé approximation (J. Comput. Phys. 1995; 119 :252–270; Comput. Meth. Appl. Mech. Engrg 1998; 163 :343–358) (therein using the parameter γ=2). We follow (SIAM Rev. 2000; 42 (3):451–484; Comput. Meth. Appl. Mech. Engrg 1995; 128 :325–359) for the analysis of this scheme and its performance on square meshes is compared with that of the quasi‐stabilized FEM (Comput. Meth. Appl. Mech. Engrg 1995; 128 :325–359). In particular, we show that the relative phase error of the numerical solution and the local truncation error of this scheme for plane wave solutions diminish at the rate O((ξ?)⁴), where ξ, ? represent the wavenumber and the mesh size, respectively. An expression for the parameter α is given that minimizes the maximum relative phase error in a sense that will be explained in Section 4.5. Convergence studies of the error in the L² norm, the H¹ semi‐norm and the l^∞ Euclidean norm are done and the pollution effect is found to be small. Copyright © 2010 John Wiley & Sons, Ltd.     

Algorithms by design with illustrations to solid and structural mechanics/dynamics

A novel procedure, concepts, and new ideas to tailor and design time operators under the notion of algorithms by design is formulated in this exposition with emphasis on applications to the broad area of computational mechanics, but with focus on solid and structural mechanics/dynamics as an illustration. The algorithms by design concepts capitalize upon: (i) the recently developed unified theory underlying computational algorithms (Int. J. Numer. Meth. Engng 2004; 59 :597–668), and (ii) newly established design spaces and algorithmic measures for evaluating the quality of computational algorithms (Int. J. Numer. Meth. Engng 2005; 64 :1841–1870). As a step in the forward direction, in this exposition we embark upon some challenging tasks with the objective to advance, tailor, and foster the design of computational algorithms for time‐dependent problems with desired and/or improved algorithmic attributes in the sense of accuracy, stability and other characteristics including algorithmic complexity in a well educated manner. The design process for computational algorithms is explained in the sense of the algorithms by design concepts via selected numerical illustrations of practical scenarios encountered in solid and structural mechanics/dynamics applications. Copyright © 2006 John Wiley & Sons, Ltd.     

Boundary element solution of the two-dimensional sine-Gordon equation using continuous linear elements

This article studies the boundary element solution of two-dimensional sine-Gordon (SG) equation using continuous linear elements approximation. Non-linear and in-homogenous terms are converted to the boundary by the dual reciprocity method and a predictor–corrector scheme is employed to eliminate the non-linearity. The procedure developed in this paper, is applied to various problems involving line and ring solitons where considered in references [Argyris J, Haase M, Heinrich JC. Finite element approximation to two-dimensional sine-Gordon solitons. Comput Methods Appl Mech Eng 1991;86:1–26; Bratsos AG. An explicit numerical scheme for the sine-Gordon equation in 2+1 dimensions. Appl Numer Anal Comput Math 2005;2(2):189–211, Bratsos AG. A modified predictor–corrector scheme for the two-dimensional sine-Gordon equation. Numer Algorithms 2006;43:295–308; Bratsos AG. The solution of the two-dimensional sine-Gordon equation using the method of lines. J Comput Appl Math 2007;206:251–77; Bratsos AG. A third order numerical scheme for the two-dimensional sine-Gordon equation. Math Comput Simul 2007;76:271–8; Christiansen PL, Lomdahl PS. Numerical solutions of 2+1 dimensional sine-Gordon solitons. Physica D: Nonlinear Phenom 1981;2(3):482–94; Djidjeli K, Price WG, Twizell EH. Numerical solutions of a damped sine-Gordon equation in two space variables. J Eng Math 1995;29:347–69; Dehghan M, Mirzaei D. The dual reciprocity boundary element method (DRBEM) for two-dimensional sine-Gordon equation. Comput Methods Appl Mech Eng 2008;197:476–86]. Using continuous linear elements approximation produces more accurate results than constant ones. By using this approach all cases associated to SG equation, which exist in literature, are investigated.     

An accelerated FFT algorithm for thermoelastic and non‐linear composites

A fast numerical algorithm to compute the local and overall responses of non‐linear composite materials is developed. This alternative formulation allows us to improve the convergence of the existing method of Moulinec and Suquet (e.g. Comput. Meth. Appl. Mech. Eng. 1998; 157 (1–2):69–94). In the present method, a non‐linear elastic (or conducting) material is replaced by infinitely many locally linear thermoelastic materials with moduli that depend on the values of the local fields. This makes it possible to use the advantages of an algorithm developed by Eyre and Milton (Eur. Phys. J. Appl. Phys. 1999; 6 (1):41–47), which has faster convergence. The method is applied to compute the local fields as well as the effective response of non‐linear conducting and elastic periodic composites. Copyright © 2008 John Wiley & Sons, Ltd.     

Minimum theorems in 3D incremental linear elastic fracture mechanics

The crack propagation problem for linear elastic fracture mechanics has been studied by several authors exploiting its analogy with standard dissipative systems theory (see e.g. Nguyen in Appl Mech Rev 47, 1994, Stability and nonlinear solid mechanics. Wiley, New York, 2000; Mielke in Handbook of differential equations, evolutionary equations. Elsevier, Amsterdam, 2005; Bourdin et al. in The variational approach to fracture. Springer, Berlin, 2008). In a recent publication (Salvadori and Carini in Int J Solids Struct 48:1362–1369, 2011) minimum theorems were derived in terms of crack tip “quasi static velocity” for two-dimensional fracture mechanics. They were reminiscent of Ceradini’s theorem (Ceradini in Rendiconti Istituto Lombardo di Scienze e Lettere A99, 1965, Meccanica 1:77–82, 1966) in plasticity. Following the cornerstone work of Rice (1989) on weight function theories, Leblond et al. (Leblond in Int J Solids Struct 36:79–103, 1999; Leblond et al. in Int J Solids Struct 36:105–142, 1999) proposed asymptotic expansions for stress intensity factors in three dimensions—see also Lazarus (J Mech Phys Solids 59:121–144, 2011). As formerly in 2D, expansions can be given a Colonnetti’s decomposition (Colonnetti in Rend Accad Lincei 5, 1918, Quart Appl Math 7:353–362, 1950) interpretation. In view of the expression of the expansions proposed in Leblond (Int J Solids Struct 36:79–103, 1999), Leblond et al. (Int J Solids Struct 36:105–142, 1999) however, symmetry of Ceradini’s theorem operators was not evident and the extension of outcomes proposed in Salvadori and Carini (Int J Solids Struct 48:1362–1369, 2011) not straightforward. Following a different path of reasoning, minimum theorems have been finally derived.     

Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation

The modeling of failure mechanisms in solids due to fracture based on sharp crack discontinuities suffers in situations of complex crack topologies including branching. This drawback can be overcome by a diffusive crack modeling based on the introduction of a crack phase field as proposed in Miehe et?al. (Comput Methods Appl Mech Eng 19:2765?C2778, 2010a; Int J Numer Meth Eng 83:1273?C1311, 2010b), Hofacker and Miehe (Int J Numer Meth Eng, 2012). In this work, we summarize basic ingredients of a thermodynamically consistent, variational-based model of diffusive crack propagation under quasi-static and dynamic conditions. It is shown that all coupled field equations, in particular the balance of momentum and the gradient-type evolution equation for the crack phase field, follow as the Euler equations of a mixed rate-type variational principle that includes the fracture driving force as the mixed field variable. This principle makes the proposed formulation extremely compact and provides a perfect basis for the finite element implementation. We then introduce a local history field that contains a maximum energetic crack source obtained in the deformation history. It drives the evolution of the crack phase field. This allows for the construction of an extremely robust operator split scheme that updates in a typical time step the history field, the crack phase field and finally the displacement field. We demonstrate the performance of the phase field formulation of fracture by means of representative numerical examples, which show the evolution of complex crack patterns under dynamic loading.     

Level‐cut inhomogeneous filtered Poisson field for two‐phase microstructures

Level‐cut homogeneous filtered Poisson fields developed in (J. Appl. Phys. 2003; 94 (6):3762–3770) to model two‐phase microstructures are defined, and their properties are briefly reviewed. Filtered Poisson fields are sums of randomly scaled and oriented kernels that are centered at the points of homogeneous Poisson fields. The cuts of these fields above specified thresholds are called level‐cut homogeneous filtered Poisson fields. It is shown that an arbitrary inhomogeneous Poisson field becomes homogeneous if observed in new coordinates, and that the mapping relating inhomogeneous and homogeneous Poisson fields can be constructed in a simple manner. This mapping and the model in (J. Appl. Phys. 2003; 94 (6): 3762–3770) provide an efficient algorithm for generating arbitrary inhomogeneous two‐phase microstructures. Developments in (Int. J. Numer. Meth. Engng 2008; DOI: 10.1002/nme.2340 ), using arguments essentially identical to those in (J. Appl. Phys. 2003; 94 (6):3762–3770) to define and generate inhomogeneous Poisson fields, overlook the natural extension of results in (J. Appl. Phys. 2003; 94 (6): 3762–3770) to these fields provided by the mapping constructed in this paper. Copyright © 2008 John Wiley & Sons, Ltd.     

Comments on ‘The computation of wavelet‐Galerkin approximation on a bounded interval’

In this letter, we identify and correct the errors in (Int. J. Numer. Meth. Eng. 1996; 39 :2921–2944). And we also develop clearer procedures for the computation of the connection coefficients from the wavelet‐Galerkin scheme. Copyright © 2007 John Wiley & Sons, Ltd.     

A time‐integration method for the viscoelastic–viscoplastic analyses of polymers and finite element implementation

The present study introduces a time‐integration algorithm for solving a non‐linear viscoelastic–viscoplastic (VE–VP) constitutive equation of isotropic polymers. The material parameters in the constitutive models are stress dependent. The algorithm is derived based on an implicit time‐integration method (Computational Inelasticity. Springer: New York, 1998) within a general displacement‐based finite element (FE) analysis and suitable for small deformation gradient problems. Schapery's integral model is used for the VE responses, while the VP component follows the Perzyna model having an overstress function. A recursive‐iterative method (Int. J. Numer. Meth. Engng 2004; 59 :25–45) is employed and modified to solve the VE–VP constitutive equation. An iterative procedure with predictor–corrector steps is added to the recursive integration method. A residual vector is defined for the incremental total strain and the magnitude of the incremental VP strain. A consistent tangent stiffness matrix, as previously discussed in Ju (J. Eng. Mech. 1990; 116 :1764–1779) and Simo and Hughes (Computational Inelasticity. Springer: New York, 1998), is also formulated to improve convergence and avoid divergence. Available experimental data on time‐dependent and inelastic responses of high‐density polyethylene are used to verify the current numerical algorithm. The time‐integration scheme is examined in terms of its computational efficiency and accuracy. Numerical FE analyses of microstructural responses of polyethylene reinforced with elastic particle are also presented. Copyright © 2009 John Wiley & Sons, Ltd.