Next: 2 Conversational tests Up: 1 Introduction Prev: 1 Introduction Contents

1.1 Solving a concrete problem using MODULEF

The following simple example enables us to illustrate how to solve a problem with the help of MODULEF. This also gives us the opportunity to introduce the finite element notation used.

We start by formulating the boundary value problem, from which the variational formulation is derived. The finite element approximations and computer implementation are then discussed briefly, describing the various steps to be executed to solve the problem numerically using the modules in the MODULEF library. To achieve this, let us consider the Dirichlet problem described below:

1.1.1 The mathematical formulation

Consider the following Dirichlet boundary value problem:



, and
, i.e.,
, where is a regular open subset of with a regular boundary .

The variational (weak) formulation derived from the boundary value problem is given by:

Find , the space of admissible functions, such that


for all .

The variational formulation of the problem is the point of departure for the Finite Element Method.

Remark: A detailed analysis of the solutions to problems (1.1) and (1.2) is not in the scope of this report and the reader should consult additional literature if necessary (for example "The Finite Element Method for Elliptic Problems" by P. G. Ciarlet (1978)).

At first sight the above type of problem cannot be solved by MODULEF. However, closer analysis will show that we can treat it as a thermal type problem.

1.1.2 The finite element approximation


The Finite Element Method consists of finding an approximate solution in a finite dimensional subspace. In order to solve the variational problem (1.2), we must therefore pose the problem in a finite dimensional (discrete) subspace, so that equation (1.2) becomes:

Find the approximate solution such that


for all , and where the finite dimensional (discrete) subspaces and are defined by:



denotes the space of polynomials of degree less or equal to k, where k depends on the choice of finite element.

Finite element notation

Before writing equation (1.3) in terms of the finite element approximations, let us define the following notation:

Equation (1.3) in terms of the above finite element notation can now be written as follows:

Find which satisfies , for , such that


for any satisfying , .

At this stage, the derivation from formulation (1.3) to formulation (1.6) is quite general (independent of MODULEF). It remains to define the corresponding mathematical operators, after which we must find their equivalents in MODULEF.

Let us first consider the classical thermal problem, described in [100], which is fully solved by the MODULEF system, where the element stiffness and mass matrices and force vector is given by:


We notice that if we set:
[k] = [I] : the thermal conductivity ([I] : unit matrix),
g=0 : the coefficient of heat transfer at ,
: the mass density, and
: the boundary force.

in the classical thermal system, we can define the element stiffness and mass matrices and force vector for the Dirichlet problem in terms of the thermal problem. The resulting simplified system is then given by:


By combining these simplified element calculations, equation (1.6) can now be rewritten, so that our problem can be completely solved by MODULEF, as a thermal type problem.

Remark: A list and the description of the thermal finite elements available in the MODULEF library are given in [100].

In terms of the above notation, (1.6) becomes:

Find which satisfies , for , such that


for all , such that , for , which results in a linear system with a positive definite symmetric matrix.

Assembly of matrices

Let be an element matrix of element T, with corresponding matrix

which constitutes the assembly operation, independent of the structure of .

There are two possible interpretations of the above example: If is the system matrix to be solved, it may be written as follows:

  1. or equivalently
  2. so that

From a mathematical point of view, these formulas are identical. However, they may lead to two different computer problems. In both cases it is necessary to calculate the element mass and stiffness matrices, after which we either:

  1. or else
Either of the above possibilities may be adopted, depending on the nature of the problem under consideration.

The finite element calculations

After assembly of the system matrices, equation (1.9) becomes:

Find satisfying , for all , such that


for all satisfying , for all .

Recall that, for the discrete space span, is given by:

Let us consider the decomposition of and , such that


where is a matrix.

With this decomposition, and by using the properties of and in subspace into account, equation (1.10) is equivalent to:


Thus, in order to calculate , it is necessary to solve (1.12), or equivalently:


where is a -dimensional array, whose components are the first components of .


  1. From a mathematical point of view, the relation (1.12) implies (1.13), therefore, to obtain u, it is possible to solve (1.12) or (1.13). However, from a practical (computer) point of view, it is preferable to solve (1.12) rather than (1.13), as (1.13) requires a second enumeration of the degrees of freedom, or otherwise the imposed boundary conditions will not be taken into account.
  2. The direct computation of (equation (1.12)), by assembling the elements , necessitates a distinction between the elements for and , so that the condition is taken into account.

    It is therefore simpler to first calculate and then compute . This operation is called "imposing the boundary conditions".

The above analysis enables us to extract the different mathematical operators, so that the problem can be decomposed into the following steps:

  1. mesh of the domain,  
  2. choice of finite elements,  
  3. physical data and calculation of the element matrices and arrays,
  4. solution of the linear problem, and
  5. display of the results.

Remark: Step 1) and 2) result from the chosen variational formulation.

Each step can be decomposed into sub-steps. The decomposition of each step, as well as the choice of modules which compute the mathematical operators, are briefly discussed in the remainder of this section.

MODULEF terminology

Let us first define the following terminology:

Variational unknown:
unknown of the continuous problem;
Degree of freedom (d.o.f.):
unknown of a discrete system (here a linear system);
associated with a d.o.f. in order to define it clearly, for example: VN is the value at a node, DX is the derivative with respect to x, etc. The complete list is given in array MAI8 of the data structure (D.S.) MAIL (consult [MODULEF User Guide - 2] for a description of data structures);
support for the degrees of freedom;
defines the geometry of the element;
usual geometric sense;
Sub-domain number:
number associated with each material, permitting the distinction of the sub-domains with different characteristics or different treatments of the same material;
Reference number:
similar notion to that of the sub-domains but this time applied to the boundaries (surfaces in , lines and points in or ). These numbers allows for the treatment "by block" of the boundary conditions, the loads, or to define the curved lines (or surfaces).
  1. Reference numbers:

    In the figure 1.1, the body is clamped at boundary AF, and pulled along boundary CD. It is composed of two materials, with a circular hole and a half-circle.

    In this configuration there are two sub-domains:

    1. the copper domain ABEF, and
    2. the iron domain BCDE,

    and at least four reference numbers:

    1. to define the circles which form the boundary of the hole and half-circle, and the line which partitions the domain into two sub-domains,
    2. to define the fixed portion of the boundary, and
    3. to define the edge DC where the load is applied.

    Remark: Lines AC and FD do not have reference numbers assigned as no particular conditions applies there.

    Figure 1.1: Clamped body with applied load 

  2. Nodes and points:

    Consider the two elements illustrated in figure 1.2: a straight P2 triangle and a curved P2 triangle. A straight element is an element whose points are the vertices. A curved element is an element whose points are the vertices plus some other characteristics. In the case of a P2 triangle, the straight element is defined by 3, and the curved element by 6 points, both having 6 nodes.

    Figure 1.2: P2 triangular elements 

1.1.3 Relation between the mathematical and computer operators


So far we have defined various mathematical operators to solve a typical problem. With each of these mathematical operators we will now associate a computer module, to solve the problem numerically. The data structures (D.S.) provides the interfaces between the modules. In this section we will mention the modules and data structures associated with each step.

Step 1:
Mesh of the domain

This step is the simplest with respect to the mathematical operator. The discretization must respect the general rules of the finite element mesh. Let and be elements of then:

The corresponding implementation is not trivial (see for example [George-1986]) and, even if all the necessary operators are present in the library, it proves labour costly. A well adapted mesh is an essential condition for good results.

The mesh generation technique adopted in MODULEF has a "tree" (top-down, bottom-up) structure, where the domain is partitioned into geometrically simple sub-domains. These sub-domains are then manipulated by simple operators (symmetry, rotation, translation, etc.) and then glued together to obtain the final mesh. This procedure is implemented by the module APNOPO for a two-dimensional domain. The mesh generation technique, as well as the different modules available, are described in detail in [MODULEF User Guide - 3].

The resulting mesh is stored in the D.S. NOPO. This D.S. contains the mesh description: the elements, the nodes (with or without vertices), as well as the coordinates at the vertices, that have been generated.

At this stage the nodes have been numbered. However, the interpolation has not been introduced, which will be discussed in the next step.



A plot of the mesh is generated by module TRNOPO
[4][MODULEF User Guide - 6], which will be discussed in the next section.

Step 2:
Choice of finite element

In this step the interpolation of functions are taken into account. This procedure is implemented by module COMACO [13]. Two cases may occur:

  1. The chosen finite element exists in one of the following libraries:
    Thermal elements [100]
    Linear elastic elements [101]
    Nonlinear elements [97]

  2. The chosen element does not exist, in which case
    • it can either be incorporated in one of the following libraries: THER, ELAS, ELNL, FLUI, or PERS [95] (this operation is described in [95]), after which we return to the first case, or
    • it can be directly described in COMACO using data cards.

The results of this step are found in two data structures: MAIL and COOR. The D.S. MAIL contains the description of the mesh (nodes, points, etc.) and of the interpolation (number and name of the variational unknowns, mnemonics, etc.). The D.S. COOR contains the coordinates of the points (vertices or not).

Note: The D.S. COOR does not contain the coordinates of the nodes on exit from COMACO. Consider again the examples of the elements "straight P2 triangle" and "curved P2 triangle": In the first example the D.S. COOR contains the coordinates of the vertices. In the second example, the D.S. COOR contains the coordinates of the vertices and those of the mid-side points.

Step 3:
Entering of physical data, and creation of element matrices and arrays

The purpose of this step is to construct the element matrices and arrays, , , and . This calculation requires the knowledge of the physical data which describes the materials, loads, temperatures, etc. In the present example, this corresponds to f and a which appear on the right- and left-hand-side of the equation. All physical data must be known at the time of activation of the module calculating the element matrices and arrays.

There are two possibilities:

  1. For the general case (data not constant), the data is provided by a subroutine or array. The module COFORC [14] constructs the D.S. FORC, which indicates how the necessary data will be provided in order to calculate the element force vectors. Module COMILI [14] constructs the D.S. MILI, which indicates how the necessary information (physical properties) will be provided in order to be able to calculate the element matrices. The element matrices and arrays are constructed by module THENEW [14].
  2. When the material characteristics, loads or heat sources are constant by sub-domain or by boundary section characterized by reference number, module THERCT [14] is used for the thermal case and ELASCT [14] in elasticity, instead of COFORC, COMILI, and THENEW.

In all the cases, it is possible to provide these characteristics by sub-domain and/or reference number, which minimizes the amount of data. The results of this step are found in the D.S. TAE, containing the element matrices and arrays corresponding to all the elements in the mesh.


  1. This step changes according to the nature of the problem. The modules THENEW, THERCT and ELASCT performs this step for the cases of thermal and linear elasticity problems. These modules can be adapted for the solution of other problems.

    Other modules implement this step for some specific problems (ex:
    [2] COTAE in nonlinear elasticity with large deformations [97]; BIHAP1,
    [4] BIHAP2 for the biharmonic problem [93], etc.).

  2. In order to treat a new problem, it is often sufficient to construct the D.S. TAE from the D.S. MAIL and D.S. COOR, and from the material characteristics.

Step 4:
Solution of the linear problem

This step can be decomposed into several sub-steps. The sub-steps vary according to the desired method of solution, however the basic principles remain the same. The MODULEF library contains a large number of solution methods for linear systems, amongst which are:

In all these cases the matrix corresponding to the finite element discretization, is sparse. In order to minimize the space occupied in memory, an appropriate method of matrix storage is sought:

The non-linear problem can often be solved by regarding it as that of a series of linear problems, each treated as described above.

The choice of solution method for linear systems is not obvious, and depends strongly on the nature of the problem to be solved. In the case of simple systems, all methods produce good results. When there are many systems to solve, for example an iterative process, there are some points in favour of both types of solution methods:

In order to analyze the different sub-steps, consider the simple example of solving a linear system by a direct Cholesky method in main memory (m.m.).

  1. Construction of pointers
      This step constructs the pointers required for the skyline storage by looking at all the elements in the mesh, and calculates the amount of memory to be occupied by the matrix. Computation will therefore only be started if the memory space available is sufficient.

    This step is performed by module PREPAC [MODULEF User Guide - 5]. The result is found in the D.S. MUA, which does not contain the matrix but only the pointers.

  2. Assembly
    Using the D.S. MUA created in the previous step, and the D.S. TAE containing the element matrices and arrays, the modules ASSMUA and ASEMBV respectively assembles the element matrices and the right-hand-side vectors, i.e., computes the following equations (see equation (1.9)):

    The results of this step is found in the D.S. MUA, containing the system matrix, and the D.S. B, containing the right-hand-side vector.

  3. Implementation of imposed boundary conditions
    Boundary conditions fall into the following two categories:
    1. those which are implemented in the variational formulation, for example the Neumann or Fourier conditions in the classical thermal problem (see [100]: Introduction), and

    2. those which are not taken into account in the variational formulation, for example the Dirichlet conditions.

    The purpose of this step, which consists of two sub-steps, is the implementation of the second type of boundary conditions.

    1. Construction of the D.S. BDCL
      The D.S. BDCL contains a description of the imposed boundary conditions (numbers and values of restricted degrees of freedom). This data structure is created by the modules COBDC1 or COBDCL [18]. The module COBDCL constructs the D.S. BDCL, degree of freedom by degree of freedom, from the data cards. Module COBDC1 constructs the same data structure by utilizing global notions, for example reference numbers and the names of the variational unknowns, and is simpler to use than COBDCL. Both modules treat boundary conditions defined as linear relations, which is particularly useful when periodic conditions are encountered.


      1. If the problem under consideration contains boundary conditions defined as linear relations, this step must be performed before step 1, as the matrix profile is modified.
      2. Module COBDC1 utilizes the D.S. COOR, which contains the nodal coordinates. This data structure can be generated by the module CORNOE [63].

    2. Effective implementation of the boundary conditions
      This step is concerned with the derivation from to , considering equations (1.11) and (1.12).

      There are several techniques for implementing boundary conditions, as described in [MODULEF User Guide - 5]: module CLIMPC, for skyline matrices. Conceptually, we can distinguish two classes:

      1. The penalized method:

        The point of departure is the following "computer" equality:

        if (for example: and ).

        In practice we start with and of equation (1.11), and replace the diagonal coefficients of by VTG and the coefficients of by , where



        is a vector of values of imposed degrees of freedom equal to zero in this example (see (1.12)).

        A variation on the above method consists of replacing the diagonal coefficients of by:

        and the right-hand-side by:
        where .

        Advantages and disadvantages:

        • In the first variation, the initial values of the diagonal coefficients are lost.
        • In the second variation, it is necessary to know the values of the diagonal coefficients in order to take the boundary conditions into account in the right-hand-side vector. It is therefore necessary to do both incorporations of boundary conditions simultaneously, which could be inconvenient in an iterative method where the matrix is constant.

        To summarize:

        Variation 1
        incorporates the boundary conditions into the matrix and right-hand-side vector separately
        [4] ([MODULEF User Guide - 5]: PROFIL and CLIMPC), with


        Variation 2
        incorporates the boundary conditions simultaneously in MUA and B, with

      2. The direct method:

        The idea is to derive equation (1.12) from (1.11) effectively. Note that in equation (1.12) an important property, namely the symmetry of the matrix, is lost. The importance lies in the property that for a symmetric matrix only half the matrix needs to be stored. We therefore want to derive (1.15), which is equivalent to equation (1.12), from equation (1.11), where (1.15) is symmetric.

        In the general case (imposed d.o.f. has a value of , zero or nonzero (1.14)), if we consider the decomposition of :

        is the vector of free degrees of freedom, and
        is the vector of the blocked degrees of freedom,

        then equation (1.12) can be written as follows:

        which is equivalent to:


        In practice, it is equation (1.15) that is calculated (using the notation in [MODULEF User Guide - 5]: PROFIL and CLIMPC, with NIVO = 3 and VTG = 1).

      Remark: The different techniques above are independent of the matrix storage technique, and can therefore be used in other methods of storage.

  4. Cholesky factorization
    This step computes the triangular matrix L which satisfies the relationship:
    and is performed by module CHOLPC [MODULEF User Guide - 5].

    The result is stored in the D.S. MUA which contains the factorized matrix.

  5. Forward- and backward-substitution
    The forward- and back-substitution solves two linear systems of triangular matrices. It is performed by module DRCHPC
    [4][MODULEF User Guide - 5].

    The result is stored in the D.S. B.

The solution of the linear system is decomposed into several steps in order to simplify the computation as much as possible. For example:

Step 5:
Interpretation of results

The graphic visualization of results is particularly important. Indeed, it is difficult to interpret the results by merely printing the output data structure B. In order to aid you with the interpretation of results, the following tools are available:

In the post-treatment, we can also utilize the modification modules for the data structures NOPO [MODULEF User Guide - 3] and B. For example, if the problem is symmetric, the computation need only be done on part of the domain. However, it is desirable to present the results on the complete domain, by "gluing" together the sub-domains.

1.1.4 Conclusions

The purpose of this example was to illustrate how find the suitable modules in the libraries, capable of solving the problem under consideration. However, it does not show how to use a given module. The examples treated in the remainder of this report will illustrate this process in detail.

Next: 2 Conversational tests Up: 1 Introduction Prev: 1 Introduction Contents