Automated computational modeling (ACM) is a fairly recent technique that exploits a hybrid approach by combining symbolic computation with numerical computation. These potentialities can be exploited in finite element formulations so it is possible to automate some procedures, including derivation ones, by performing the calculations that until recently were carried out manually in a more efficient and more precise way. This opens up new frontiers in automation and in the degree of precision and complexity of FEM simulations.

Finite Element Method

Finite Element Analysis (FEM) is one of the most powerful tools for the numerical simulation of complex industrial problems. There are commercial software for specific or more general uses and many research software where new formulations and discretization schemes still need to be developed. There are continuous evolutions in order to take into account the needs of advanced engineering applications. Finite element simulations are a general method for the numerical solution of partial differential equations. They therefore allow to describe new materials, structures, plants and biological processes and are a key element for predictive science in many complex processes in engineering, medicine or physics.

The FEM was born in the 1960s, but following the development of IT tools, it has an exponential evolution and development. It is established as one of the best tools for investigating those complex systems. The generality of the method, initially developed by engineers and later demonstrated by mathematicians, has allowed many studies and applications. This has paved the way for new lines of research which currently address issues of considerable theoretical and practical interest.

Shape functions

The idea behind the FEM is to see the domain of interest as the union of many sub-domains of elementary form. Therefore on each elementary subdomain the single partial differential equations (PDE) are considered. In a continuum problem, whatever the variable of interest, temperature, speed, density, displacement, etc., the problem presents an infinite number of unknowns. This is due to the fact that the field of interest is variable point by point. The discretization introduced with the FEM reduces it to a finite problem of unknowns. Going to expresses the unknown field in terms of approximating functions that identify the variable of interest in some points of the domain called nodes.

Sixteen triangular basis functions used to reconstruct J0 FEM
Summation of basis functions fen

The basic idea of the approximation used in the FEM is to approximate the true trend of the unknown function with that of some particular functions with a known trend, typically polynomial functions. There is therefore a reduction in the number of degrees of freedom which are infinite in the continuous medium, while, considering only some points (nodes) of the structure, they are finite in number.

We then speak of discretization of the structure as that operation that allows you to pass from the real structure to the discretized (or approximate) one for which it is possible to apply the finite element method in order to obtain an engineering solution to the problem.

Weak formulation

Also knowing that the solution through the use of numerical methods takes place by means of electronic computers, the idea of discretization is linked to the physical limit that these machines have in terms of data storage. We pass from a differential problem, which the computer cannot solve, to an algebraic problem that can be solved by inverting the system matrix. To arrive at this algebraic system, after considering the approximation of the unknown function as the sum of the shape functions for the values of the unknowns in the nodes.

fem mesh
fem solution
FEM Solution

We accept that the approximate function has a certain residue (that is, it is not exact) and we must minimize this residue in the sense of reducing its weighted average (I integrate it and I ask that it be null). We introduce some test functions on which we pass the derivative problem. That is, with this last statement, we pass from the problem in strong form to the problem in weak form, as a series of mathematical requirements for the unknown function are missing.

Then we introduce the function of forms by transforming the continuous function to a sum of nodal values and weights. This allows to take the terms out of the integral and to have a series of coefficients of the single element that multiply the nodal unknowns. By choosing then the test functions in a number congruent to the unknowns, an algebraic system is obtained which is more easily solvable.

\bold K\bold u=\bold f\Longrightarrow \bold u=\bold f\bold K^{-1}

Symbolic approach

At the base of every computation automation process there is the need to write the basic equations that describe the problem in an input form that allows to manipulate the variables in a symbolic way. Although computer algebra can be considered a subfield of scientific computation, they are generally considered to be distinct fields because scientific computation is usually based on numerical computations with approximate floating-point numbers, while symbolic computation emphasizes exact computation with expressions containing variables which have no given value and are manipulated as symbols.

Software applications that perform symbolic computations are called algebraic computer systems, with the term system alluding to the complexity of major applications that include at least:

  • Methods to represents numerical data inside a computer
  • Coding languages for the user (usually different from the language used for the implementation)
  • Dedicated memory manager
  • User interface for the input / output of mathematical expressions
  • Large set of routines to perform usual operations, such as simplification of expressions, differentiation by chain rule, polynomial factorization, indefinite integration, etc.

It is widely used for research and to design formulas used in number programs. Used for full scientific calculations, when purely numerical methods fail, such as in public key cryptography, or for some non-linear problems. a very slow process to implement. It is necessary to consider the theory of derivation of complex tensor fields in order to obtain tangent and residual matrices necessary to implement the formulations. The revolution came with the introduction of symbolic computation and automated steps.

Automated Computational Modelling

There are infinite ways to set a FEM resolution but from the point of view of automation we can summarize it in some elementary steps:

There are infinite ways to set a FEM resolution but from the point of view of automation we can summarize it in some elementary steps:

  1. Strong formulation and of the boundary problem
  2. Switch to weak formulation
  3. Definition of the discretization of the domain and approximation of the variables and introduction of the test functions
  4. Derivation and solutions of the equations at the element level
  5. Derivation of algebraic equations describing the contribution for the stiffness matrix and the tangent matrix
  6. Coding of derivatives in the language appropriate for the FEM solver
  7. Generating a Mesh
  8. Solution of the global problem
  9. Viewing and analyzing the results

Clearly there is a clear separation between a part of the problem relating to the single elements and independent of the physical problem and a part related to the laws governing the physical phenomenon of interest.

The implementation of automated calculation brings with it several advantages:

  1. Symbolic formulation is more easily understood and less prone to error
  2. Algebraic operations, including differentiation, are performed automatically
  3. The codes are automatically generated, highly efficient and easily transportable
  4. Multi-language and multi-environment capability
  5. Symbolic elements and inputs prepared for a wide range of problems and easily adaptable to the specific problem
    Easy implementation of multi-field and multi-physics problems.

It is based on the variational formulations of the problem, but using weak formulations it could be adapted to arbitrary partial differential equations. The main advantage of using symbolic code development is that development time, especially for complex materials or elements, is reduced by orders of magnitude.

Ace Fem

Thus a first division between commercial software and business software is possible. In the company software we try to have an interface that allows to have a high reproducibility of the tests in different cases. Furthermore, we try to have an easier debugging. In the research software there is a high level of customization and a fine control of the resolving methods. This often leads to having to reset part of the formulation as the cases and boundary conditions change.

Widely used Symbolic and Algebraic Computational systems (SAC) such as Mathematica or Maple have become an integrated computing environment that covers all aspects of a computational process, including numerical analysis and graphical presentation of results. The FEM method is implemented in general SAC systems usually as an add-on package or toolbox like AceGen Korelc (2011) for Mathematica. The major limitation of symbolic systems, when applied to complex engineering problems, is an uncontrollable growth of expressions and consequently redundant operations and inefficient codes.

fem solution
Hybrid (symbolic-numerical) approach to the automation of the FEM method

AceFEM is a general finite element system for Mathematica that effectively combines symbolic and numerical approaches. It includes a general finite element environment designed to solve multiphysics and multifield problems, the powerful AceGen automatic code generation package for symbolic generation of new finite elements, and access to the AceShare finite element file sharing system.


Using Mathematica, the program provides both extensive symbolic capabilities and the numerical efficiency of a finite element trading environment. It allows to solve large-scale industrial problems with several million elements and use advanced Mathematica features.

The AceFEM system is equipped with a large library of finite elements (solid, thermal, contact, 2D, 3D, etc.). Additional finite elements can be symbolically derived and encoded using the specialized AceGen automatic code generation package for generating highly numerically efficient finite element codes. AceFEM also allows you to easily create standalone, customized and finite element-based applications for Mathematica.


A general requirement for efficient use of the automation system is that the basic equations that define the problem are written in an input form that allows symbolic manipulation. The advantage of automation in computational mechanics becomes evident only when the description of the problem, that is, the way the basic equations are written, is appropriate for the symbolic description.


It is known that the problem of elastic equilibrium consists in determining the stress field, the deformations and the displacements of a continuum made of material with a linear elastic behavior, with constraints and loads. It is also known that the equilibrium solution can be seen as the minimization of a total energy. So, writing its functional, we know that its solution is the one that minimizes it. Furthermore we can see the integral directly on the relative domain in ξ and η by applying the map and taking into account the transformation of areas (determinant of the Jacobian).

\Pi_{e}(\mathbf{u})=\int_{\Omega_{e}} \frac{1}{2} \varepsilon(\mathbf{u}) \cdot \mathbb{C} \varepsilon(\mathbf{u}) d v=\int_{\Omega_{\square}} \frac{1}{2} \varepsilon(\mathbf{u}) \cdot \mathbb{C} \varepsilon(\mathbf{u}) \operatorname{det}\left(\mathbf{J}_{e}\right) d V

So if the solution displacement is the minimum of the free energy, the resolution passes to impose its stationarity. This, having an incorrect solution, will give us a residue. We have therefore returned to the weak form of equilibrium. Deriving again we will obtain the tangent matrix, as well as the stiffness of the formulation.

\mathbf{r}_{e}=\frac{d \mathbf{\Pi }_{e}}{d \mathbf{u}}\\ \:\\
\mathbf{K}_{e}=\frac{d \mathbf{r}_{e}}{d \mathbf{u}}

These quantities, for each element, will be passed to the FEM solver which assembles them into a global matrix and solving the associated system will give us the solution, that is the displacements that satisfy the equilibrium. In particular, we have a discretized solution where the integral is reduced to a weighted sum of the functions in the nodes and therefore in matrices that represent a linear problem.

\left[\mathbf{r}_{e}\right]_{i}=\frac{d \Pi_{e}}{d\left[\mathbf{u}_{e}\right]_{i}}=\sum_{a} w_{g} J_{e}\left(\xi_{g}, \eta_{g}\right)\left(\frac{d}{d\left[\mathbf{u}_{e}\right]_{i}} \Psi_{e}\left(\xi_{g}, \eta_{g}, \mathbf{u}_{e}\right)\right)\\
\left[\mathbf{K}_{e}\right]_{i j}=\frac{d\left[\mathbf{r}_{e}\right]_{i}}{d\left[\mathbf{u}_{e}\right]_{j}}

And so we’re going to use an iterative algorithm to find the solution in terms of displacements:

\mathbf{u}_{\text {tot }}=\mathbf{K}_{\text {tot }}^{-1} \mathbf{f}_{\text {tot }}:\mathbf{R}\left(\mathbf{u}_{t o t}\right)=\mathbf{K}_{t o t} \mathbf{u}_{t o t}-\mathbf{f}_{t o t}=0\\\:\\\Longrightarrow \begin{aligned}
&\mathbf{u}_{t o t}^{(1)}=-\left.\mathbf{K}_{t o t}^{-1} \mathbf{R}\right|_{0} \\
&\mathbf{u}_{t o t}^{(2)}=\mathbf{u}_{t o t}^{(1)}-\left.\mathbf{K}_{t o t}^{-1} \mathbf{R}\right|_{\mathbf{u}_{t o t}^{(1)}}