The year 2021 marks the eightieth anniversary of the invention of the finite element method (FEM), which has become the computational workhorse for engineering design analysis and scientific modeling of a wide range of physical processes, including material and structural mechanics, fluid flow and heat conduction, various biological processes for medical diagnosis and surgery planning, electromagnetics and semi-conductor circuit and chip design and analysis, additive manufacturing, and in general every conceivable problem that can be described by partial differential equations (PDEs). The FEM has fundamentally revolutionized the way we do scientific modeling and engineering design, ranging from automobiles, aircraft, marine structures, bridges, highways, and high-rise buildings. Associated with the development of FEMs has been the concurrent development of an engineering science discipline called computational mechanics, or computational science and engineering.

In this paper, we present a historical perspective on the developments of finite element methods mainly focusing on its applications and related developments in solid and structural mechanics, with limited discussions to other fields in which it has made significant impact, such as fluid mechanics, heat transfer, and fluid–structure interaction. To have a complete storyline, we divide the development of the finite element method into four time periods: I. (1941–1965) Early years of FEM; II. (1966–1991) Golden age of FEM; III. (1992–2017) Large scale, industrial applications of FEM and development of material modeling, and IV (2018–) the state-of-the-art FEM technology for the current and future eras of FEM research. Note that this paper may not strictly follow the chronological order of FEM developments, because often time these developments were interwoven across different time periods.

1 (1941–1965) The Birth of the Finite Element Method

The origin of the FEM as a numerical modeling paradigm may be traced back to the early 1940s. In 1941, A. Hrennikof, a Russian-Canadian structural engineer at the University of British Columbia, published a paper in the ASME Journal of Applied Mechanics on his membrane and plate model as a lattice framework [1]. This paper is now generally regarded as a turning point that led to the birth of the FEM. In the paper, he discretized the solution domain into a mesh of lattice structure, which was the earliest form of a mesh discretization.

On May 3rd, 1941, the same year that Hrennikoff published his paper, R. Courant of New York University delivered an invited lecture at a meeting of the American Mathematical Society held in Washington D.C. on his numerical treatment using a variational method to solve a second order PDE, which arises from Saint–Venant’s torsion problem of a cylinder. In this work, Courant systematically used the Rayleigh–Ritz method with a trial function defined on finite triangle subdomains, which is a primitive form of of the finite element method [2].

Courant’s presentation was later published as a paper in 1943 [3]. Similar works of discretization and variational formulations were also reported in the literature, including McHenry [4], Prager and Synge [5], and Synge and Rheinboldt [6]. As Ray Clough commented in his 1980 paper, “One aspect of the FEM, mathematical modeling of continua by discrete elements, can be related to work done independently in the 1940s by McHenry and Hrennikoff-formulating bar element assemblages to simulate plane stress systems. Indeed, I spent the summer of 1952 at the Boeing Airplane Company trying to adapt this procedure to the analysis of a delta airplane wing, the problem which eventually led to the FEM”.

By the early 1950s, several engineers and academics had further developed and extended these early approaches to solve real engineering problems in aeronautical and civil engineering. In parallel but with different emphases, Argyris at the Imperial College London, and M. J. Turner (1950–1956) at Boeing Company, who was later joined by R.W. Clough of UC Berkeley and H.C. Martin of Washington University, developed what we know today as the earliest form of the finite element method (1954), which was called the Matrix Stiffness Method at the time. In a paper published in 1960 [14], R.W. Clough coined the phrase Finite Element Method, and this unassuming and right-to-the-point phrase was an instant hit, bringing out the essence of the method.

With his deep insights and profound vision into Courant’s variational approach, J.H. Argyris developed the energy method for engineering structures [6,7,9], a foundational development enabling FEM for solids. While importantly, Turner, Clough, Martin, and Topp successfully developed FEM interpolants for triangular elements [10], which is suitable for structural parts with arbitrary shape. In some senses, the invention of the triangle element was a “quantum leap”, and hence for a large spectrum of the engineering community, the inception of the FEM is the publication of the landmark paper by Turner et al. [10]. The following is an excerpt from a 2014 document that celebrates the 50th anniversary of the formation of National Academy of Engineering, which is an official account of that part of FEM history: To ensure safety and avoid costly modifications after planes entered production, engineers needed a reliable method for determine in advance whether

their designs could withstand the stresses of fight. M. Jon Turner, head of Boeing’s Structural Dynamics Unit, addressed that problem in the early 1950s by bringing civil engineering professor Ray Clough of the University of California, Berkeley, and Harold Martin of the University of Washington to Boeing for summer ``faculty internships,’’ Collectively, they created a method of structural analysis that Turner applied at Boeing using computers to perform the myriad calculations need to predict real-world performance. That fruitful collaboration led to Clough’s development a few years later of what he named the finite element method (FEM). Clough formed a research group at UC Berkeley that used FEM in a host of analytical and experimental activities, from designing buildings and structures to withstand nuclear blasts or earthquakes to analyzing structural requirements for spacecraft and deep-water offshore drilling. By revolutionizing the application of computer technologies in engineering, FEM continues to help engineers design to this day all sorts of durable, cost-effective structures. Meanwhile, Turner’s efforts at Boeing contributed to the success of its renowned line of commercial jets, beginning in 1958 with the 707 and continuing in 1964 with the 727, which could land on shorter runways and serve more airports. Equipped with three fuel-efficient turbofan engines, the 727 became the workhorse of commercial aviation and helped achieve a threefold increase in U.S. passenger air traffic in the’60s.

Independently and separately, in the early 1960s, Kang Feng of the Chinese Academy of Science also proposed a discretization-based numerical method for variational principles for solving elliptic partial differential equations [11]. As Peter Lax [11, 12] commented, “Independently of parallel developments in the West, he (Feng) created a theory of the finite element method. He was instrumental in both the implementation of the method and the creation of its theoretical foundation using estimates in Sobolev spaces….”, which was one of the first convergence studies of FEMs.

During this period, several great engineering minds were focusing on developing FEMs. J.H. Argyris with his co-workers at the University of Stuttgart; R. Clough and colleagues such as E. L. Wilson and R.L. Taylor at the University of California, Berkeley; O.C. Zienkiewicz with his colleagues such as E. Hinton and B. Irons at Swansea University; P. G. Ciarlet at the University of Paris XI; R. Gallager and his group at Cornell University, R. Melosh at Philco Corporation, B. Fraeijs de Veubeke at the Université de Liège, and J. S. Przemieniecki at the Air Force Institute of Technology had made some important and significant contributions to early developments of finite element methods.

To understand what happened sixty years ago, we quote an excerpt from a FEM history paper by Clough and Wilson [13], in which they recalled:

When Clough presented the first paper using the finite element terminology in 1960 it attracted the attention of his friend, Professor O. C. Zienkiewicz, who was then on the faculty at Northwestern University. A few weeks after the presentation of the paper Zienkiewicz invited Clough to present a seminar on the finite element method to his students. Zienkiewicz was considered one of the world’s experts on the application of the finite difference method to the solution of continuum mechanics problems in Civil Engineering; therefore, Clough was prepared to debate the relative merits of the two methods. However, after a few penetrating questions about the finite element method, Zienkiewicz was almost an instant convert to the method.

Z. Bazant, then a visiting associate research engineer at UC Berkeley, recalled: …… The founding of FEM was paper in the ASCE Conf. on Electronic Computation Clough [14], which was the first to derive, by virtual work, the finite element stiffness matrix of an element of a continuum (a triangular constant strain element). ……. I recall Ray Clough showing to me his 1962 report to the U.S. Engineer District, Little Rock, Corps of Engineers on his analysis of a crack observed in Norfolk Dam [15], during my stay in Berkeley in 1969. I was mesmerized by seeing that 1962 report. It presented a 2D stress analysis of large crack observed in Norfork dam. The dam was subdivided into about 200 triangular elements and provided stress contours for a number of loading cases. …… Clough was at that time way ahead of anybody else.

To collaborate Z. Bazant’s recollection, we cite E.L. Wilson’s recounting of that part of FEM history:

In 1956, Ray, Shirley, and three small children spent a year in Norway at the Ship Research Institute in Trondheim. The engineers at the institute were calculating stresses due to ship vibrations to predict fatigue failures at the stress concentrations. This is when Ray realized his element research should be called the Finite Element Method which could solve many different types of problems in continuum mechanics. Ray realized the FEM was a direct competitor to the Finite Difference Method. At that time FDM was being used to solve many problems in continuum mechanics. His previous work at Boeing, the Direct Stiffness Method, was only used to calculate displacements not stresses.

In the fall semester of 1957, Ray returned from his sabbatical leave in Norway and immediately posted a page on the student bulletin board asking students to contact him if they were interested in conducting finite element research for the analysis of membrane, plate, shell, and solid structures. Although Ray did not have funding for finite element research, a few graduate students who had other sources of funds responded. At that time, the only digital computer in the College of Engineering was an IBM 701 that was produced in 1951 and was based on vacuum tube technology. The maximum number of linear equations that it could solve was 40. Consequently, when Ray presented his first FEM paper in September 1960, “The Finite Element Method in Plane Stress Analysis,” at the ASCE 2 nd  Conference on Electronic Computation in Pittsburgh, Pennsylvania, the coarse-mesh stress-distribution obtained was not very accurate. Therefore, most of the attendees at the conference were not impressed. After the improvement of the speed and capacity of the computers on the Berkeley campus, Professor Clough’s paper was a very fine mesh analysis of an existing concrete dam. The paper was first presented in September 1962 at a NATO conference in Lisbon, Portugal. Within a few months, the paper was republished in an international Bulletin, which had a very large circulation, as “Stress Analysis of a Gravity Dam by the Finite Element Method”, (with E. Wilson), International Bulletin RILEM, No. 10, June 1963.

The Lisbon paper reported on the finite element analysis of the 250-foot-high Norfork Dam in Arkansas, which had developed a vertical crack during construction in 1942. The FEM analysis correctly predicted the location and size of the crack due to the temperature changes and produced realistic displacements and stresses within the dam and foundation for both gravity and several hydrostatic load conditions. Because of this publication, many international students and visiting scholars came to Berkeley to work with Professor Clough. Also, he freely gave the FORTRAN listing of their finite element analysis computer program to be used to evaluated displacement and stresses in other two-dimensional plane structures with different geometry, materials, and loading. Therefore, professional engineers could easily use the powerful new FEM to solve for the stress distributions in their structural engineering problems in continuum mechanics. However, he did not capitalize on his success in the development of the FEM. He returned to the task of building the earthquake engineering program at Berkeley – the task he given when he was first hired in 1949.

For their decisive contributions to the creation and developments of FEM, R. W. Clough was awarded the National Medal of Science in 1994 by the then vice-president of the United States Al Gore, while O. C. Zienkiewicz was honored as a Commander of the Order of the British Empire (CBE). Today, the consensus is that J.H. Argyris, R.W. Clough and O. C. Zienkiewicz made the most pivotal, critical, and significant contributions to the birth and early developments of finite element method following an early contribution to its mathematical foundation from R. Courant.

It is worth noting that E.L. Wilson of UC Berkeley was the first person to develop finite element open-source software. An excerpt from Clough and Wilson’s paper in 1999 stated 13: In 1958 Wilson, under the direction of Clough, initiated the development of an automated finite element program based on the rectangular plane stress finite element developed at Boeing. After several months of learning to program the IBM 701, Wilson produced a limited capacity, semiautomated program which was based on the force method. An MS research report was produced, which has long since been misplaced, with the approximate title of Computer Analysis of Plane Stress Structures. …… In 1959 the IBM 704 computer was installed on the Berkeley Campus. It had 32K of 32-bit memory and a floating-point arithmetic unit which was approximately 100 times faster than the IBM 701. This made it possible to solve practical structures using fine meshes.

It is also worth mentioning that, Oden [16], a senior structural engineer in the research and development division of General Dynamics Corporation at Fort Worth at the time, wrote a 163-page comprehensive technical report with G.C. Best, in which they developed a long list of solid and structural finite elements, including tetrahedral, hexahedral, thin plate, thick plate, plate element with stringers or stiffeners, composite sandwich plate elements, and shallow shell elements [17]. In fact, Oden and Best wrote one of the first general purpose FEM computer codes at the time. The Fortran FEM code developed by Oden and Best had an element library that includes elements for 3D elasticity, 2D plane elasticity, 3D beam and rod elements, composite layered plate and shell elements, and elements for general composite materials. Their work also included hybrid methods and stress based FEMs, which may be even earlier than those of Pian [18]. Moreover, Oden and Best’s FEM code was also able to handle FEM eigenvalue modal analysis, and numerical integration over triangle and tetrahedra elements, and it had linear system solvers for general FEM static analysis that were among the most effective at that time (see [19]). This FEM computer code was used for many years in aircraft analysis and design in the aerospace and defense industry.

It appears that, I.C. Taig [20] of English Electric Aviation first introduced the concept of the isoparametric element when he used the matrix-displacement method to conduct stress analysis of aerospace structures, which was later formally dubbed ``isoparametric element’’ by Ergatoudis et al. [21].

It should be mentioned that there are some other pioneers who made some significant contributions in the early developments of FEM, such as Levy [22], Comer [23], Langefors [24], Denke [25], Wehle and Lansing [26], Hoff et al. [27], and Archer [28], among others. These individuals came together made remarkable and historic contributions to the creation of finite element method. Among them, some notable contributions were made by J. S. Przemieniecki (Janusz Stanisław Przemieniecki), who was a Polish engineer and a professor and then dean at the Air Force Institute of Technology in Ohio in the United States from 1961 to 1995. Przemieniecki conducted a series pioneering research works on using FEMs to perform stress and buckling analyses of aerospace structures such as plates, shells, and columns (see Przemieniecki [29], Przemieniecki and Denke [30], Przemieniecki [28,29,30,31,32,33]).

In terms of worldwide research interest, by 1965, FEM research had become a highly active field, with the total number of papers published in the literature exceeding 1000. During this period, there were two seemingly unrelated events for FEM development, which significantly affected future FEM developments. These events were the discovery of mixed variational principles in elasticity. In 1950, E. Reissner [34] rediscovered E. Hellinger’s mixed variational principle from 1914 [35], in which both the displacement field and the stress field are the primary unknowns. This variational principle is called the Hellinger–Reissner variational principle. Shortly after, Hu [36] and Washizu [37] proposed a three-field mixed variational principle in elasticity, which was called the Hu-Washizu variational principle. As early as 1964, Pian [18] recognized the potential of using these variational principles to formulate Galerkin weak form-based FE formulations and proposed the assumed stress FEM. This began the use of mixed variational principles to formulate Galerkin FEMs, which was followed by the assumed strain FEM developed by J. Simo and his co-workers in the later period of FEM developments.

2 (1966–1991) The Golden Age of the Finite Element Method

The mid 1960s saw rapid developments in finite element method research and applications. As T.J.R. Hughes recalled, “I first heard the words ‘the Finite Element Method” in 1967—which changed my life. I started to read everything that was available and convinced my boss to start the Finite Element Method Development Group, which he did. Dr. Henno Allik was Group Leader, and I was the Group, then we added programmers. In one year, we had a 57,000-line code, GENSAM (General Structural Analysis and Matrix Program, or something like that). That was 1969, and the code was continually developed thereafter and may still be in development and use at GD/Electric Boat and General Atomics, originally a division of General Dynamics.”

Starting from the end of 1960s, the rigorous approximation theory that underpins the FEM started to be developed (e.g.  Aubin [38], Zlamal [39], Birkhoff [40], Nitsche [41], Aziz [42], Bubuska [43], Bubuska and Aziz [44], Dupont [45], Douglas and Dupont [46]. Nitsche and Schatz [47], and Bubuska [48]). This movement was first highlighted by the proof of optimal and superconvergence of FEMs. This attracted the interest of some distinguished mathematicians all over the world, including G. Birkhoff, M.H. Schultz, R.S. Varaga, J. Bramble, M. Zlamal, J. Cea, J.P. Aubin, J. Douglas, T. Dupont, L.C. Goldstein, LR. Scott, J. Nitsche, A.H. Schatz, P.G. Cialet, G. Strang, G. Fix, JL. Lions, M. Crouzeix, P.A. Raviart, and I. Babuska, A.K. Aziz, and J.T. Oden (see Bramble, Notsche and Schatz [49], Bubuska, Oden, and Lee [50]) . Some notable results developed for the proof of FEM convergence are the Cea lemma and the Bramble-Hubert lemma. The fundamental work of Nitsche [41] on L∞ estimates for general classes of linear elliptic problems also stands out as one of the most important contributions for mathematical foundation of FEM in 1970s. It may be noted that unlike other mathematics movements, the convergence study of FEMs was an engineering-oriented movement. The mathematicians soon found that, in practice, engineers were using either non-conforming FEM interpolants or numerical quadrature that violates variational principles or the standard bilinear form in Hilbert space. G. Strang referred to these numerical techniques as “variational crimes”. To circumvent complicated convergence proofs, the early FEM patch tests were invented by B. Iron and R. Melosh, which were proven to be instrumental for ensuring convergence to the correct solution.

Following T.H.H. Pian’s invention of the assumed stress element, attention shifted to the mixed variational principle based FEM [51, 52]). In 1965, L.R. Hermann [53] proposed a mixed variational principle for incompressible solids. However, most mixed variational principles are not extreme variational principles, and thus suffer from numerical instability. In early 1970s, I. Babuska and F. Brezzi discovered their groundbreaking results, known today as the Babuska-Brezzi condition, or the LBB condition, giving tribute to O. Ladyzhenskya—a Russian mathematician, who provided the early insight of this problem. The so-called LBB condition, or the inf–sup condition, provides a sufficient condition for a saddle point problem to have a unique solution that depends continuously on the input data; thus, it provides a guideline to construct shape functions for the mixed variational principle-based FEM (see: [54]).

Entering the 1970s, FEM development began to focus on using FEM to simulate the dynamic behavior of structures, including crashworthiness in the automotive industry. Various time integration methods had been developed, including the Newmark-beta method, Wilson-theta method [55], the Hilbert-Hughes-Taylor alpha method [56], the Houbolt integration algorithm, and the explicit time integration algorithm [57].

Since the early 1970s, FEM explicit time integration methods had been used to solve various engineering problems, e.g., dynamic contact problems [58], dynamics problems of solid structure [59], fluid dynamics problems [60], and it was extensively used in large scale FEM Lagrangian hydrocodes developed in Lawrence Livermore National Laboratory [61]. However, an important development came in the late 1970s when T. Belytschko, K. C. Park, and later TJR Hughes proposed using explicit or implicit-explicit, or implicit time integration with damping control to solve nonlinear structural dynamics problems. Other notable contributors to this key development are A. Combescure, A. Gravouil among others (see e.g., [62]). It turned out that the explicit time integration was a game changer for the automotive industry, establishing FEM technology as the main tool of passenger vehicle design and crashworthiness analysis. By the end of 1980s, there were thousands of workstations running explicit time integration-based FEM codes in the three major automakers in the United States.

One of the main FEM research topics in the 1980s was using FEM techniques to solve the Navier–Stokes equation as an alternative to finite difference and finite volume methods. Starting the early 1970s, J.T. Oden and his co-workers had begun working on FEM solutions for fluid dynamics [63]. The main challenge in using FEM solving fluid dynamics problems is that the Navier–Stokes equation is not an elliptical partial different equation, and the minimization or variational principle-based Petrov–Galerkin procedure may suffer both stability as well as convergence issues. To resolve these issues, T.J.R. Hughes, and his co-workers such as A.N. Brooks and T.E. Tezduyar developed the streamline upwind/Petrov Galerkin method and later Stabilized Galerkin FEM to solve Navier–Stokes equations under various initial and boundary conditions [64]. Furthermore, Hughes and co-workers later developed space–time FEMs and variational multiscale FEMs (see Hughes and Tezdyuar [65] Brooks and Hughes [66], Mizukami and Hughes [67], and Hughes [68]). For this work, from 1986 to 1991, Hughes and his co-workers such as L.P. Franca and others wrote a ten-part series on FEM formulations for computational fluid dynamics Hughes et al. [69,70,71,71] and Shakib et al. [72]. Several years later, Hughes developed a Green function-based subgrid model acting as a stabilized finite element method [68].

Towards the mid-1980s, advanced FEM mesh generation techniques have been developed, which incorporate various solid modeling techniques by using interactive computer graphics and adaptive mesh generator [73], Bennett and Botkin [74] as well as improved quadtree approach to generate FEM meshes for complex geometric shapes and objects [75]. Today, FEM mesh generation has become an integrated part of solid modeling and engineering design, which has the capabilities of automatic node insertion and refinement Wordenweber [76] and Ho-Le [77].

Among the many advances in FEM technologies in 1980s, the most notable may belong to J. Simo at the University of California, Berkeley and later at Stanford University. Simo and Taylor [78] developed the consistent tangent operator for computational plasticity, which was a milestone after the original concept of consistent linearization proposed by Hughes and Pister [79]. Moreover, after the Hughes-Liu 3D degenerated continuum shell and beam elements and the Belytchko-Tsay single-point element, Simo and his co-workers, such as L. Vu-Quoc and D.D. Fox, developed geometrically exact beam and shell theories and their FEM formulations (see Simo and Vu-Quoc [80], Simo and Fox [81], and Simo et al. [82]. Moreover, Simo and his co-workers such as MS Rifai and F. Armero also developed various assumed strain or enhanced strain methods for mixed variational formulations (see Simo and Rifai [83] and Simo and Armero [84]). It should be noted that E. Ramm and his colleagues at the University of Stuttgart have also made significant contributions on geometrically nonlinear shell element formulations over a span of more than thirty years e.g., Andelfinger and Ramm [85] and Bischoff and Ramm [86].

Another highlight of FEM technology is the development of FEM solvers for fluid–structure interaction. From the mid-1970s to early 1990s, there was an urgent need to develop techniques to solve large-scale fluid–structure interaction problems in the aerospace and civil engineering industries. A class of FEM fluid–structure interaction solvers were developed, and some early contributors include J. Donea, A. Huerta (see [87,88,89,89]), and the Hughes-Liu-Zimmermann Arbitrary Lagrangian–Eulerian (ALE) fluid–structure FEM formulation (see [90, 91]), which describes the moving boundary problem. ALE-based FEM simulations were used due to their ability to alleviate many of the drawbacks of traditional Lagrangian-based and Eulerian-based FEM formulations.

When using the ALE technique in engineering modeling and simulations, the computational mesh inside the domains can move arbitrarily to optimize the shapes of elements, while the mesh on the boundaries and interfaces of the domains can move along with materials to precisely track the boundaries and interfaces of a multi-material system. The invention of the ALE FEM may be credited to Hirt, Amsden, and Cook [92]. C. Farhat was the first person to use a large-scale parallel ALE-FEM solver to compute fluid–structure interaction problems [93]. He and his group systematically applied FEM-based computational fluid dynamics (CFD) solvers for aircraft structure design and analysis. They developed the finite element tearing and interconnecting (FETI) method for the scalable solution of large-scale systems of equations on massively parallel processors. Today, some fundamental ALE concepts have also been applied to numerical modeling engineering fields other than FEM, such as meshfree modeling. For example, to alleviate tensile instability and distorted particle distributions in smoothed particle hydrodynamic (SPH) simulations, a so-called shifting technique has been adopted in many of today’s SPH simulations (e.g., Oger et al. [94]).

The FEM fluid–structure interaction research had a major impact on many practical applications, such as providing the foundation for the patient specific modeling of vascular disease and the FEM-based predictive medicine later developed by C.A. Taylor and T.J.R. Hughes and their co-workers in the mid-1990s (see Taylor et al. [95]). Holzapfel, Eberlein, Wriggers, and Weizsäcker also developed large strain FEM formulations for soft biological membranes (see Holzapfel et al. [96]).

Another major milestone in the development of FEMs was the invention and the development of nonlinear probabilistic or random field FEMs, which was first developed by W. K. Liu and T. Belytschko in the late 1980s. (e.g., Liu et al. [97]). By considering uncertainty in loading conditions, material behavior, geometric configuration, and support or boundary conditions, the probabilistic FEM provided a stochastic approach in computational mechanics to account for all these uncertain aspects, which could then be applied in structure reliability analysis. The random field FEM research has become crucial in civil and aerospace engineering and the field of uncertainty quantification.

In the early 1980s, M.E. Botkin at General Motors research Lab [98], and N. Kikuchi and his group at the University of Michigan developed structural shape optimization FEM for the automotive industry (see [99], [100], [101]). Other contributors include M. H. Imam from Uman Al-Qura University [102]). This preceded the seminal 1988 paper of Bendsoe and Kikuchi, who developed a homogenization approach to finding the optimal shape of a structure under prescribed loading. Later developments in topology optimization were driven by Gengdong Cheng, Martin Bendsoe, and his student, Ole Sigmund.

One of the driving forces in the FEM development during in the early decades was the safety analysis of big dams, at first, then of concrete nuclear reactor vessels for gas-cooled reactors, [103], nuclear containments, hypothetical nuclear reactor accidents [104, 105] and of tunnels and of foundations for reinforced concrete structures. To simulate concrete failure, the vertical stress drops in FEMs, and progressive softening technique were introduced already in 1968. However, the spurious mesh sensitivity and its impact on strain localization was generally overlooked until demonstrated mathematically by Bazant [105], because strain softening states of small enough test specimens in stiff enough testing frames are stable. Numerically this was demonstrated by crack band FEM calculations in Bazant and Cedolin [106], where it was also shown that spurious mesh sensitivity causing the sudden stress drop can be avoided by adjusting the material strength to ensure the correct energy release rate. Hillerborg et al. [107] avoided mesh sensitivity by using an interelement cohesive softening called the fictitious crack model. Despite the success of the early FEM calculations, the concept of progressive strain-softening damage was not generally accepted by mechanicians until its validity and limitations were demonstrated by Bazant and Belytschko [108] and Bazant and Chang [109]. They showed that the existence of elastic unloading stiffness (previously ignored) makes waves propagation in a strain-softening state possible. In 1989, Lubliner, Oliver, Oller, and Onate ([110] developed a plastic-damage theory-based FEM formulation to model concrete materials by introducing internal variables, which has the capacity of modeling concrete material degradation and cracking. Today, the multiscale based homogenization and damage analysis method is the state-of-the-art FEM modeling for concrete materials [111].

In the mid-1980s, the mesh sensitivity issue in calculating strain softening or strain localization problems in computational plasticity became a challenging topic. It was eventually accepted that the partial differential equations that are associated with the classical plasticity become ill-posed after the material passes the yield point and enters the softening stage. This especially became a dire situation when civil engineers applied FEM to solve complex structural and geotechnical engineering problems, which involved complex plastic deformations of concrete, rock, soil, clay, and granular material in general. Bazant and others realized that this is because the classical continuum plasticity theory lacks an internal length scale. To remedy this problem, starting from the middle 1980s, many efforts were devoted to establishing FEM formulations of nonlocal (Bazant et al. [112]), strain-gradient, strain-Laplacian media, micropolar or Cosserat continua, because they provided an internal length scale, allowing FEM simulations to capture, in a mesh-independent manner, the strain softening, strain localization or shear band formation. Pijaudier-Cabot and Bazant [113] developed an effective nonlocal FEM in which nonlocality is applied only to the damage strain. Other influential and representative works in this topic are from De Borst [114], Peerlings, De Borest and others [115]. Steinmann and Willam [116], Dietsche et al. [117], Steinmann [118], and Iordache and Willam [119].

In 1976, TJR Hughes, R.L. Taylor, J.L. Sackman, A. Curnier, and W. Kanoknukulchai published a paper entitled “A finite element method for a class of contact-impact problems.” [120]. This is one of the earliest FEM analyses in computational contact mechanics. It is the very first work on FEM modeling of dynamic contact and impact problems, and it plays an important role in the simulation accuracy for engineering problems involving interaction between different continuous objects. Examples include sheet metal forming, target impact and penetration, and interaction between pavement and tires. Developing accurate FEM contact algorithms has been a focal point since the 1980s. Various FEM contact algorithms have been developed, and some main contributors are N. Kikuchi, J.T. Oden, J. Simo, P, Wriggers, R.L. Taylor, P. Papadopoulos, and T.A. Laursen. FEM contact algorithm research remained an active research topic until late 1990s (see Simo et al. [121], Simo and Lausen [122], Kikuchi and Oden [123], and Papadopoulos and Taylor [124]).

In this period, one interesting emerging area was the development of FEM exterior calculus by Arnold et al. [125]. The FEM exterior calculus uses tools from differential geometry, algebraic topology, and homological algebra to develop FEM discretizations that are compatible with the underlying geometric, topological, and algebraic structures of the problems that are under consideration.

3 (1992–2017) Broad Industrial Applications and Materials Modeling

The first major event in FEM development this period was the formulation of the Zienkiewicz-Zhu error estimator [126], which was a major contribution to the mathematical approximation theory of FEMs in the 1990s. The Zienkiewicz-Zhu posteriori error estimators provide the quality control of a FEM solution with an optimal use of computational resources by refining the mesh adaptively. The idea and spirit of Zienkiewicz-Zhu was further carried out by Ainsworth and Oden (see [127, 128, 129]), and today using posteriori error estimation to improve the quality has been elevated to the height of Bayesian inference and Bayesian update. This research topic is now intimately related with what is now called validation and verification (V&V).

Since the late 1970s, Szabo [130] and Babuska [131, 132] started to develop hp versions of FEMs based on piecewise-polynomial approximations that employ elements of variable size (h) and polynomial degree (p). They discovered that the FEM converges exponentially when the mesh is refined using a suitable combination of h-refinements (dividing elements into smaller ones) and p-refinements (increasing their polynomial degree). The exponential convergence makes the method a very attractive choice compared to most other FEMs, which only converge with an algebraic rate. This work continued until the late 1990s, spearheaded by M. Ainsworth, L. Demkowicz, J.T. Oden, C.A. Duarte, O.C. Zienkiewicz, and C.E. Baumann (see Demkowicz et al. [133], Demkowicz et al. [134], Demkowicz et al. [135], Oden et al. [136], Oden et al. [137], Baumann and Oden [138]. J. Fish. [139] also proposed a s-version FEM by superposing additional mesh(es) of higher-order hierarchical elements on top of the original mesh of C0 FEM discretization, so that it increases the resolution of the FEM solution.

To solve material and structural failure problems, research work in the 1990s focused on variational principle based discretized methods to solve fracture mechanics problems or strain localization problems. In 1994, Xu and Needleman [140] developed a FEM cohesive zone model (CZM) that can simulate crack growth without remeshing, which was later further improved by M. Ortiz and his co-workers, e.g. Camacho and Ortiz [141] and Ortiz and Pandolfi [142], who later used CZM FEM to solve fragmentation and material fatigue problems [143]. It should be noted that long before the invention of cohesive zone model, Pietruszczak and Mroz [144] developed the first cohesive FEM for shear fracture in soil. Later, Bazant’s group at Northwestern University developed various interface FEMs, such as the microplane model, to study size effects of concrete and other brittle composite materials (Brocca and Bazant [145], Caner and Bazant [146]). These models became a standard tool for simulating missile impact and explosions at, e.g., ERDA Vicksburg. A microplane FEM model for fiber composites has been developed for Chrysler and Ford Co. to compare various designs of automobile crush-cans (Smilauer et al. [147]). An anisotropic poromechanical microplane model has been formulated and used for FEM analysis of hydraulic fracturing (Rahimi et al. [148]).

To alleviate mesh bias issues in modeling material fracture and damage problems, T. Belytschko and WK. Liu developed meshfree particle methods, namely the element-free Galerkin (EFG) method [149] and the reproducing kernel particle method (RKPM) [150, 151], which are based on the moving least square method and the wavelet multiresolution analysis, respectively. RKPM provides consistency and thus convergence enhancements as compared to the popular smoothed particle hydrodynamics (SPH) method. Li et al. [152, 153] successfully employed meshfree Galerkin methods to accurately simulate adiabatic shear band formation and propagation with minimum mesh adaptation. At the same time, Duarte and Oden developed the so-called hp-Cloud method, Onate et al. developed a finite point method [154], and Atluri and Zhu [155] proposed a meshless local Petrov–Galerkin (MLPG method, among many other meshfree methods. Fleming and Belytschko also showed that singularity functions could be included in the approximation functions to greatly improve simulations involving fracture mechanics (see Fleming et al. [156]).

One of the most challenging problems in the development of meshfree Galerkin methods is how to integrate the weak form, because the meshfree interpolants are highly irregular and it is difficult to make them variationally consistent. In 2001, J. S. Chen and his co-workers proposed a stabilized conforming nodal integration method for meshfree RKPM method, which is not only simple and stable, but also variationally consistent with the Galerkin weak formulation (see Chen et al. [157]).

Shortly after the meshfree method developments, I. Babuska and his co-workers developed the partition of unity finite element method (PUFEM), which was later coined generalized finite element method (GFEM) (see Melenk and Babuska [158] and Babuska and Melenk [159]). PUFEM is a powerful method because it can be used to construct FE spaces of any given regularity, which is a generalization of the h, p, and hp versions of the FEM, as well as providing the ability to embed an analytic solution into the FEM discretization instead of relying upon a generic polynomial interpolant.

A significant breakthrough in computational fracture mechanics and FEM refinement technology came in the late 1990s, when Belytschko and his co-workers, including Black, Moes, and Dolbow, developed the eXtended finite element (X-FEM) (see [160, 161], which uses various enriched discontinuous shape functions to accurately capture the morphology of a cracked body without remeshing. Because the adaptive enrichment process is governed by the crack tip energy release rate, X-FEM provides an accurate solution for linear elastic fracture mechanics (LEFM). In developing X-FEM, T. Belytschko brilliantly utilized the PUFEM concept to solve fracture mechanics problems without remeshing. Entering the new millennium, Bourdin, Francfort, Marigo developed a phase-field approach for modeling material fracture [162]. Almost simultaneously, Karma and his co-workers [163], [164]) also proposed and developed the phase-field method to solve crack growth and crack propagation problems, as the phase field method can accurately predict material damage for brittle fracture without remeshing. The main advantage of the phase field approach is that by using the Galerkin FEM to solve the continuum equations of motion as well as a phase equation, one can find the crack solution in continuum modeling without encountering stress singularity as well as remeshing, and the crack may be viewed as the sharp interface limit of the phase field solution. Some of the leading contributors for this research are Bourdin, Borden, Hughes, Kuhn, Muller, Miehe, Landis, among others (see Bourdin and Chambolle [165], Kuhn and Müller [166]. Miehe et al. [167], and Borden et al. [168], Wilson et al. [169], and Pham et al. [170]).

As mentioned before, the main reason for the huge success of FEMs is their broad applicability to engineering analysis and design across scientific disciplines. On the other hand, most mechanical engineering designs are performed by using various computer-aided design (CAD) tools, such as solid modeling. To directly blend the FEM into CAD design tools, T.J.R. Hughes, and his co-workers such as J.A. Cottrell and Y. Bazilevs [171], [172] developed the isogeometric analysis (IGA) FEM, which established the Galerkin variational weak formulation in the control mesh and uses the non-uniform rational basis spline (NURBS) functions as the FEM shape function to solve the problem at design stage. IGA-FEM method successfully integrates FEM into conventional NURBS-based CAD tools, without converting data between CAD and FEA packages in analyzing new designs during development stage.

Due to the emergence of nanotechnology, various multiscale methods have been developed to couple atomistic methods such as molecular dynamics and density functional theory (DFT) and other ab initio methods with continuum scale FEMs. The most notable contributions in this area are hand-shake method [173] quasi-continuum FEM [174], and the bridging scale method developed by Wagner and Liu [175]. In 2007, Gavini, Bhattacharya, and Ortiz developed quasi-continuum orbital-free density-functional theory (DFT) FEM for multi-million atom DFT calculations (see Gavini et al. [176]).

Broadly speaking, there has been a rapid development in computational homogenization methods since the late 1990s, which is an alternate approach to obtaining continuum-scale properties based on smaller scale microstructures. The computational homogenization method or the FEM homogenization methods for composite materials may be divided into two main categories: (1) Computational asymptotic homogenization method or multiscale computational homogenization, which is aimed for modeling composite materials with periodic microstructure. The pioneer of multiscale homogenization FEM may be credited to T.Y. Hou and his co-workers [177], and other earlier contributors are N. Kikuchi [101], Ghosh et al. [178], Fish et al. [179, 180]) and M. S. Shephard, (2) Computational micromechanics method, which is mainly aimed for composite materials with random microstructure, even though it may also be applied to materials with periodic microstructure, where the main contributors are: P. Suquet et al. at the French National Centre for Scientific Research (CNRS), Dvorak et al.[181] and J. Fish at Rensselaer Polytechnic Institute and C. Miehe and his co-workers (e.g. [182]) at the University of Stuttgart. During the same period, many novel FEM fluid–structure interaction solvers have been developed, for instance, the immersed FEM developed by Zhang et al. [183], which was motivated by the immersed boundary method pioneered by C. Peskin of Courant Institute of Mathematical Sciences, New York University.

An important FEM application area emerged with the developments of computational plasticity. The early finite element computational plasticity formulation is based on hypoelastic–plastic rate formulation. To satisfy the objectivity requirement, T.J, R. Hughes and J. Winget first proposed the so-called incremental objectivity (see Hughes and Winget [184]), which was probably in the first time the esoteric continuum mechanics theory was applied to practical FEM computational formulations and computations, and it in turn promoted the development of nonlinear continuum mechanics in 1980s and 1990s. Soon afterwards. Simo and Hughes then extended Hughes-Winget incremental objectivity algorithm to the finite deformation case in computational plasticity (see Simo and Hughes [185]). The notion of consistency between the tangent stiffness matrix and the integration algorithm employed in the solution of the incremental problem was introduced by Nagtegaal [186] and Simo and Taylor [78]. Consistent formulations have been subsequently developed for finite deformation plasticity by Simo and Oritz [187] (1988) within the framework of multiplicative decomposition of the deformation gradient and hyperelasticity. Also in 1980s, based on the Gurson model, Tvergaard and Needleman [188] developed the FEM based Gurson-Tvergaard-Needleman model, which is probably the most widely used FEM computational plasticity constitutive model used in material modeling, though recently attention has turned to developing machine-learning based or data-driven computational plasticity models e.g., F. Chinesta, et al. [189], and the unsupervised machine learning data-driven finite element methods, called as the Self-consistent Clustering Analysis (SCA) (see Z. Liu et al. [190, 191]).

An important advance of the FEM is the development of the crystal plasticity finite element method (CPFEM), which was first introduced in a landmark paper by Pierce et al. [192]. In the past almost four decades, there are numerous researchers who have made significant contributions to the subject, for example, A. Arsenlis and DM. Parks from MIT and Lawrence Livermore National Laboratory [193], [194], [195], Dawson et al. at Cornell University (Quey et al. [196]; Mathur and Dawson [197]; Raabe et al. at Max-Planck-Institute fur Eisenforschung [196,197,198,201], among others. Based on crystal slip, CPFEM can calculate dislocation, crystal orientation and other texture information to consider crystal anisotropy during computations, and it has been applied to simulate crystal plasticity deformation, surface roughness, fractures and so on. Recently, S. Li and his co-workers developed a FEM-based multiscale dislocation pattern dynamics to model crystal plasticity in single crystal [202, 203]. Yu et al., [204] reformulated the self-consistent clustering analysis (SCA) for general elasto-viscoplastic materials under finite deformation. The accuracy and efficiency for predicting overall mechanical response of polycrystalline materials are demonstrated with a comparison to traditional full-field FEMs.

In 2013, a group of Italian scientists and engineers led by L. Beirão da Veiga and F. Brezzi proposed a so-called virtual element method (VEM) (see Beirao et al. [205, 206]). The virtual element method is an extension of the conventional FEM for arbitrary element geometries. It allows the polytopal discretizations (polygons in 2-D or polyhedra in 3-D), which may be even highly irregular and non-convex element domains. The name virtual derives from the fact that knowledge of the local shape function basis is not required, and it is in fact never explicitly calculated. VEM possesses features that make it superior to the conventional FEM for some special problems such as the problems with complex geometries for which a good quality mesh is difficult to obtain, solutions that require very local refinements, and among others. In these special cases, VEM demonstrates robustness and accuracy in numerical calculations, when the mesh is distorted.

As early as 1957, R. Clough introduced the first graduate FEM course at UC-Berkeley, and since then FEM courses at both graduate and undergraduate levels have been added into engineering higher education curriculums in all the major engineering schools and Universities all over the world. As J.T. Oden recalled in his 1963 paper, “I went on to return to academia in 1964 and among my first chores was to develop a graduate course on finite element methods. At the same time, I taught mathematics and continuum mechanics, and it became clear to me that finite elements and digital computing offered hope of transforming nonlinear continuum mechanics from a qualitative and academic subject into something useful in modern scientific computing and engineering.”.

By the end of 2015, there have been more than several hundred monographs and textbooks on the FEM published in dozens of languages worldwide. An exposition of FEM mathematical theory by Strang and Fix [207] was among the earliest of those FEM books. J.T. Oden and his collaborators such as GF Carey and JN Reddy and others wrote a five-volume FEM monograph in 1980s (see [167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,213, 214]). Other influential finite element books or monographs are those by Zienkiewicz and Cheung [215], and Zienkiewicz and Taylor [216], ora later Zienkiewicz, Taylor and Zhu [217], TJR. Hughes [218], Cook et al. [219] or later Cook [220]), Bathe [221], and the Nonlinear FEM monographs by Belytschko et al. [222] and by De Borst et al. [223], among others, all have made major impacts on FEM educations and applications. Among all these FEM monographs and textbooks, the book by Zienkiewicz and Taylor [216] or Zienkiewicz et al. [217] probably have had most impacts on FEM technology popularization, which may be because Taylor wrote a FEM research computer program code named FEAP, which was placed in the appendix of that book, providing an immediate guidance and instruction on how to implement FEM in computers.

The ready availability of FEM textbooks and tutorials, along with the large FEM software market, have made the FEM accessible to users across academia and industry. The development of FEM software technology started in the early1960s. In 1963, E.L. Wilson and R. Clough developed a structural mechanics computing code called Symbolic Matrix Interpretive System, SMIS, which was initially intended to fill the gap between the matrix method and hand calculating in structure mechanics. It turns out the development of SMIS led to the birth of FEM software. Then, based on SMIS, Wilson initiated and developed a general-purpose static and dynamic Structural Analysis Program, SAP. In late 1960s and the early 1970s, K, J. Bathe developed the nonlinear FEM code ADINA based on SAP IV and NONSAP. Today, the brand name SAP2000 has become synonymous with the state-of-the-art FEM structural analysis and design methods since its introduction over 55 years ago.

At the same period, to compete with the Soviet Union’s space program, NASA developed its own FEM code called NASTRAN (NASA STRuctural ANalysis Program). The first version of NASTRAN was called COSMIC Nastran, which debuted in 1969, with a key figure in its development being R. H. MacNeal. As early as 1963, R. H. MacNeal founded the MacNeal-Schwendler Software Corporation (MSC) along with R. Schwendler. Under his leadership, MSC developed its first structural analysis software called SADSAM (Structural Analysis by Digital Simulation of Analog Methods), which showed the early form of FEM analysis software technology. In response to NASA’s request for proposals in 1965 for a general purpose structural analysis program, Dr. MacNeal contributed significantly to the early efforts of the aerospace industry by successfully simulating on-the-ground physical testing through computing to deliver the right answers and physics needed to take humans to the moon. In 1971, MSC Software released a commercial version of Nastran, named MSC Nastran.

About the same time in 1960s, J. A. Swanson worked at Westinghouse Astronuclear Laboratory in Pittsburgh, and he was responsible for stress analysis of the components in NERVA nuclear reactor rockets. While there he then developed 3D FEM model and computer codes to analyze and predict transient stresses and displacements of the reactor system. To integrate different computer codes and streamline the processing, Swanson asked his employer Westinghouse to develop a general-purpose FEM computer code, but his suggestion was rejected, and then he left the company and developed the initial ANSYS FEM code. Today, ANSYS has become one of the major FEM commercial software worldwide.

Several years later, J.O. Hallquist at Lawrence Livermore National Laboratory also developed a 3D nonlinear FEM code called DYNA3D, which was extensively used impact, dynamic contact, and failure analysis of structures, which later evolved to LS-DYNA. LS-DYNA is the major FEM software used in automobile design and crashworthiness analyses. In 2018, Livermore Software Technology was purchased by ANSYS, and LS-DYNA became a part of ANSYS as ANSYS LS-DYNA.

In the early development of nonlinear FEM software, there were two early pioneers: P. V. Marcal and D. Hibbitt, who was Marcal’s student at Brown University. In the early 1970s, Marcal founded the MARC corporation to develop the first general purpose nonlinear FEM program, which is still widely used today in industry and academia for analysis of complex structures, such as nuclear reactors, car crashworthiness and manufacturing processes. While Hibbitt together with B. Karlsson and P. Sorenson formed a FEM company called HKS in 1978, and they released a large-scale commercial FEM software called ABAQUS. One of the major features of ABAQUS is that it allows the user defined subroutines, which greatly facilitates the researchers to conduct their researches by using the standard FEM solvers with reliability and efficiency. Together with ANSYS, ABAQUS is one of the two commercial FEM software that dominate the market.

By the late 1990s to early 2000s, the FEM software industry had become a multi-billion-dollar business. There were several household FEM software company names such as ANSYS, ABAQUS, ADINA, LS-DYNA, NASTRAN, COMSOL Multiphysics, CSI, among others. Today, there are also a plethora of open-source FEM software available online, such as FreeFEM, OpenSees, Elmer, FEBio, FEniCS Project, DUNE, among some others.

4 Present and Future

The modern form of the FEM can routinely solve many large and complex industrial problems. It enables developing a fundamental understanding and allows for the predictive analysis for product design. For new scientific discoveries and engineering innovations, the development of new scientific principles often trails the pace of new inventions with the sheer volume of data that are generated across multiple spatial, temporal, and design parameter (spatial–temporal-parameter) spaces. For this reason, FE researchers are studying various forms of machine and deep learning methods, of which this class of methods covers the largest class of interpolations. According to the universal approximation theorem, a neural network (NN) can be designed and trained to approximate any given continuous function with any desired accuracy [224, 225, 226, 227] which is believed to drive new discoveries and enable future computational discretization technologies. In this context, Mechanistic Data Science (MDS) FEMs, which combine known scientific principles with newly collected data, will provide the critically needed research that can be a boon for new inventions.

Scientific and engineering problems typically fall under three categories: (1) problems with abundant data but undeveloped or unavailable scientific principles, (2) problems that have limited data and limited scientific knowledge, and (3) problems that have known scientific principles with uncertain parameters, with possible high computational load [228]. In essence, mechanistic data science (MDS) mimics the way human civilization has discovered solutions to difficult and unsolvable problems from the beginning of time. Instead of heuristics, MDS uses machine learning methods like active deep learning and hierarchical neural network(s) to process input data, extract mechanistic features, reduce dimensions, learn hidden relationships through regression and classification, and provide a knowledge database. The resulting reduced order form can be utilized for design and optimization of new scientific and engineering systems (see Liu et. al. [229]). Thus, the new focus of the FEM research has shifted towards the development of machine learning based FEMs and reduced order models.

With the recent development of machine learning and deep learning methods, solving FEM by constructing a deep neural network has become a state-of-the-art technology. Earlier research focused on building up a shallow neural network following the FEM structure to solve boundary value problems. Takeuchi and Kosugi [230] proposed a neural network representation of the FEM to solve Poisson equation problems. Yagawa and Aoki [231] replaced the FEM functional with the network energy of interconnected neural networks (NNs) to solve a heat conduction problem. Due to the limitation of computationliual power and slow convergency rate in shallow neural networks, earlier applications could only solve simple PDE problems. After the 2010s, neural networks for solving computational mechanics problems have become increasingly popular with the rapid growth of deep learning techniques and the development of more sophisticated neural network structures, such as convolutional neural networks (CNN), Generative Adversarial Networks (GAN) and residual neural networks (ResNet). For its high dimensional regression ability, some researchers, for example, Ghavamiana and Simone [232] used deep neural networks as a regression model to learn the material behavior or microstructure response. Other works focus on solving PDEs using deep learning neural networks. G. Karniadakis and his coworkers (see Raissi et al. [233, 234], Karniadakis et al. [235]) proposed a Physics-Informed Neural Networks (PINNs) to solve high dimensional PDEs in the strong form with constraints to accommodate both natural and essential boundary conditions. The idea of constructing deep neural networks following the FEM structure is investigated again with advanced neural network methodologies. Weinan E and his co-workers [236] and B. Yu [237] proposed a Deep Ritz Method for solving variational problems. Sirignano and Spiliopoulos [238] proposed the so-called Deep Galerkin Method (DGM) to solve high-dimensional PDEs. Zabaras and his co-worders proposed a CNN-based physics-constrained deep learning framework for high-dimensional surrogate modeling and uncertainty quantification (see Zhu et al. [239]). Rabczuk [240] systematically explore the potential to use NNs for computational mechanics by solving energetic format of the PDE (see Samaniego et al. [241]). Lee [242] proposed a partition of unity network for deep hp approximation of PDEs and extensively the training and initialization strategy to accelerate the convergence of the solution process (see Lee et al. [242]). The constructing of element shape function by activation functions has been studied by J. Opschoor and his coworkers (See [243, 244]).

Inspired by the universal approximation of deep neural networks (DNN), Zhang et al. [245] published the first paper on the construction of the FEM shape functions-based on the hierarchical nature of the DNN, called Hierarchical Deep-learning Neural Networks (HiDeNN). Specifically, the authors demonstrated the construction of a few classes of deep learning interpolation functions, such as the reproducing kernel particle method (RKPM), non-uniform rational B-spline (NURBS), and isogeometric analysis (IGA), among other approximation techniques. Saha et al. [228] generalized HiDeNN to a unified Artificial intelligence (AI)-framework, called HiDeNN-AI. HiDeNN-AI can assimilate many data-driven tools in an appropriate way, which provides a general approach to solve challenging science and engineering problems with little or no available physics as well as with extreme computational demand.

To reduce the FE computational cost, the so-called two-stage data-driven methods have been proposed of which during the offline stage, where a database generated by the FEM is first developed, and the final solutions are computed during the online stage. C. Farhat and his group at Stanford University have developed several dimensional reductions of nonlinear FEM dynamic models, including mesh sampling and weighting for the hyper reduction of nonlinear Petrov‐Galerkin reduced‐order models (see Farhat et al. [246], and Grimberg et al. [247]). To further reduce the FEM computational burden in multiscale analysis, Liu et al. [190] applied the unsupervised machine learning techniques, such as the k-mean clustering method to group the material points during the offline stage and obtain the final solutions by solving the reduced-order Lippmann–Schwinger micromechanics equations. This class of data-driven approaches circumvents the computational burden of the well-established FE square method, the offline-online database approach to solve the concurrent FEM problems. They named the method the Finite Element-self-consistent clustering analysis (FE-SCA) of which the computational cost of the microscale analysis is reduced tremendously in multiple orders of magnitude speed up (see Li et al. [248]). Gao et al. proposed an alternative (FE-SCAxSCA…xSCA) clustering analysis, of which the continuum FEM scale is concurrently solved with the (n-1) coupled-scale Lippmann–Schwinger micromechanics equations [249, 250]. It has been recently extended to (FE-SCA**2) by He et al. [251, 252].

Ortiz and his co-workers at Caltech developed data driven FEMs for dynamics and noisy data (see [253, 254]). Chen and his co-workers  have developed a physics-constrained data-driven RKPM method based on locally convex reconstruction for noisy databases (see He and Chen [255]). S. Li and his group at UC-Berkeley utilized FEM solution generated data to develop a machine learning based inverse solution to predict pre-crash data of car collision (see Chen et al. [256]). Bessa et al. [257] proposed a data-driven framework to address the longstanding challenge of a two-scale analysis and design of materials under uncertainty applicable to problems that involve unacceptable computational expense when solved by standard FEM analysis of representative volume elements. The paper defined a framework that incorporates the SCA method to build large databases suitable for machine learning. The authors believe that this will open new avenues to finding innovative materials with new capabilities in an era of high-throughput computing (“big-data”).

Reduced order modeling has been an active research field over the last decades. Early research works focused on the proper orthogonal decomposition (POD) method (also known as Karhunen-Loève transform, or principal component analysis) with the purpose of reducing the degrees of freedom of the discretized equations. The POD based model reduction has shown great success in computational fluid dynamics, see e.g., the works of Berkooz et al. [258]. For further accelerating the simulations, K. Willcox and her coworkers (see [259, 260]) proposed a missing point estimation method, which is known later as a hyper reduction method. Other notable works related to POD and hyper reduction methods are the Gauss–Newton with approximated tensors (GNAT) method, Grassmann manifold based reduced basis adaptation, thanks to C. Farhat and his coworkers (see Carlberg and Farhat [261], Carlberg et al. [262], Grimberg et al. [247], Amsallem and Farhat [263], and Farhat et al. [246]). For solid mechanics, Ryckelynck et al. [264] proposed a hyper reduction method based on FEM for dealing with nonlinear problems. Lu et al. [265] proposed an adaptive hyper reduction for coupled thermal-fluid analysis. Another type of model reduction method, which is based on mathematics and has a rigorous error bound estimate, is called reduced basis method, as proposed by Maday and Rønquist [266]. The proper generalized decomposition (PGD) based model reduction, as an extension of POD, can be

dated back to 1980s, and it was introduced by P. Ladevèze (See Ladeveze [267], Ladeveze and Rougee [268]) under the name of radial time–space approximation. F. Chinesta et al. [269, 270] developed a PGD method to account for the parameter space, aiming at building offline computational vademecum for fast online predictions. It is noted that PGD methods are based on the idea of separation of variables and in particular a canonical tensor decompostion (TD).

Recent works have been conducted to combine deep machine learning methods with reduced order modeling methods. Zhang et al. [271] consolidated the various attributes of TD and PGD methods with HiDeNN and proposed the so-called HiDeNN-TD and HiDeNN-PGD methods. The comparison of FEM, TD/PGD, HiDeNN-TD/PGD, HiDeNN, and DNN has been conducted in terms of accuracy and speed. It is shown that the HiDeNN-TD/PGD outperforms other methods with a good balance between accuracy and speed. The proposed HiDeNN-TD/PGD method is expected to provide novel powful tools for solving large scale high dimensional problems while maintaining a high accuracy. Various applications, including multiphysics coupled additive manufacturing, multiscale composite modeling, and structural topology optimization, have been discussed in a generalized reduced order machine learning finite element framework [272]. In particular, this reduced order machine learning framework is expected to enable ultra large-scale high resolution topology design that is currently challenging for the FEM based topology optimization. In this regard, Lu et al. [273] recently developed a convolution HiDeNN-TD formulation with a built-in density filter for high resolution topology design. This convolution formulation incorporates the concept of meshfree approximation into the finite element function approximation and allows smoother solutions and automatic length-scale control in topology design. It is shown that the convolution HiDeNN-TD leads to better design with smoother and fine structures. This general convolution formulation opens new perspectives to resolve the length scale effect and can be applied to many orther problems, such as additive manufacturing and microstructure modeling.

The extra-ordinary interpolating capability of neural network has resulted in numerous research in approximating and solving the ordinary and partial differential equations e.g. Chen et al. [274]. Currently, the researchers are looking into approximating the mathematical operators directly from universal approximation of operator theory (see Chen and Chen [275]). The goal of this endeavor is to solve the integral equations such as the Lippmann–Schwinger micromechanics equations (Z. Liu et. al. [190]). A major step forward towards this target is made by Li et al. [276] in their proposed Graph Kernel Network, in which they developed neural operators for solving partial differential equations. In this work, a kernel-based graph neural network is shown to be able to mimic the Green’s function method for solving partial differential equation. One drawback of this method is the associated computational cost and storage requirement increases with the size of the problem. A more general version of approximating the Green’s function operator is proposed by Lu et al. in their DeepONet [277] and by Gin et al. in their DeepGreen methods (see [278]). These approximation works are vital in computational mechanics as these can directly solve the micromechanics equation for multiscale analysis like FE-SCA. Moreover, these networks have suggested that the analytical calculus method such as differentiation and integration, and solution of differential and integration equations can be directly expressed as an approximation of neural network. In this regard, recurrent neural networks have shown promised to be identical in structure of wave equations [279]. These researches are directed towards a future when the research and education in the science, technology, engineering, and mathematics (STEM) sector with discrete calculus will be transformed with the aid of deep learning. We are envisioning this field as deep learning discrete calculus, a new perspective in teaching calculus by the integration of calculus definition, numerical analysis, and deep learning (see Liu et al. [280]).

The development of model reduction methods meets the urgent demand in the industry for fast and nearly real time simulations of engineering problems, such as online dynamic system control, structural health monitoring, vehicle health monitoring, on-line advanced manufacturing feedback control, automated driving controls and decisions, etc. Such applications usually require an intensive interaction between sensors, control algorithms, and simulation tools. Practical optimal control may require a reliable prediction within the range of milli- or sub-milli-second. Reducing the computation cost of simulations has been one of the major motivations for developing model reduction methods. Other reduced order modeling related topics are the feature engineering and data analytics, which constitute an extensive literature in field of machine learning. Thus, the reduced order modeling and the machine learning have intrinsic connections. Developing reduced order machine learning methods may enable physics-data combined models that can overcome the current bottleneck in model reduction methods and purely data-driven machine learning approaches.