FreeFEM Documentation on GitHub

stars - forks

Developers

File formats

Mesh file data structure

The mesh data structure, output of a mesh generation algorithm, refers to the geometric data structure and in some case to another mesh data structure.

In this case, the fields are

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
MeshVersionFormatted 0

Dimension [DIM](int)

Vertices
[Number of vertices](int)
X_1(double) Y_1(double) (Z_1(double)) Ref_1(int)
...
X_nv(double) Y_nv(double) (Z_nv(double)) Ref_nv(int)

Edges
[Number of edges](int)
Vertex1_1(int) Vertex2_1(int) Ref_1(int)
...
Vertex1_ne(int) Vertex2_ne(int) Ref_ne(int)

Triangles
[Number of triangles](int)
Vertex1_1(int) Vertex2_1(int) Vertex3_1(int) Ref_1(int)
...
Vertex1_nt(int) Vertex2_nt(int) Vertex3_nt(int) Ref_nt(int)

Quadrilaterals
[Number of Quadrilaterals](int)
Vertex1_1(int) Vertex2_1(int) Vertex3_1(int) Vertex4_1(int) Ref_1(int)
...
Vertex1_nq(int) Vertex2_nq(int) Vertex3_nq(int) Vertex4_nq(int) Ref_nq(int)

Geometry
[File name of geometric support](char*)

   VertexOnGeometricVertex
   [Number of vertex on geometric vertex](int)
   Vertex_1(int) VertexGeometry_1(int)
   ...
   Vertex_nvg(int) VertexGeometry_nvg(int)

   EdgeOnGeometricEdge
   [Number of geometric edge](int)
   Edge_1(int) EdgeGeometry_1(int)
   ...
   Edge_neg(int) EdgeGeometry_neg(int)

CrackedEdges
[Number of cracked edges](int)
Edge1_1(int) Edge2_1(int)
...
Edge1_nce(int) Edge2_nce(int)

When the current mesh refers to a previous mesh, we have in addition

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
MeshSupportOfVertices
[File name of mesh support](char*)

   VertexOnSupportVertex
   [Number of vertex on support vertex](int)
   Vertex_1(int) VertexSupport_1(int)
   ...
   Vertex_nvsv(int) VertexSupport_nvsv(int)

   VertexOnSupportEdge
   [Number of vertex on support edge](int)
   Vertex_1(int) EdgeSupport_1(int) USupport_1(double)
   ...
   Vertex_nvse(int) EdgeSupport_nvse(int) USupport_nvse(double)

   VertexOnSupportTriangle
   [Number of vertex on support triangle](int)
   Vertex_1(int) TriangleSupport_1(int) USupport_1(double) VSupport_1(double)
   ...
   Vertex_nvst(int) TriangleSupport_nvst(int) USupport_nvst(double) VSupport_nvst(double)

   VertexOnSupportQuadrilaterals
   [Number of vertex on support quadrilaterals]
   Vertex_1(int) TriangleSupport_1(int) USupport_1(double) VSupport_1(double)
   ...
   Vertex_nvsq(int) TriangleSupport_nvsq(int) USupport_nvsq(double) VSupport_nvsq(double)
  • nv means the number of vertices
  • ne means the number of edges
  • nt means the number of triangles
  • nq means the number of quadrilaterals
  • nvg means the number of vertex on geometric vertex
  • neg means the number of edges on geometric edge
  • nce means the number of cracked edges

bb file type to Store Solutions

The file is formatted such that:

1
2
3
4
5
2 [Number of solutions](int) [Number of vertices](int) 2

U_1_1(double) ... U_ns_1(double)
...
U_1_nv(double) ... U_ns_nv(double)
  • ns means the number of solutions
  • nv means the number of vertices
  • U_i_j is the solution component i at the vertex j on the associated mesh.

BB file type to store solutions

The file is formatted such that:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
2 [Number of solutions](int) [Type 1](int) ... [Type ns](int) [Number of vertices](int) 2

U_1_1_1(double) ... U_(type_k)_1_1(double)
...
U_1_1_1(double) ... U_(type_k)_nbv_1(double)

...

U_1_1_ns(double) ... U_(type_k)_1_ns(double)
...
U_1_nbv_ns(double) ... U_(type_k)_nbv_ns(double)
  • ns means the number of solutions
  • type_k mean the type of solution k:
    • 1: the solution is scalar (1 value per vertex)
    • 2: the solution is vectorial (2 values per vertex)
    • 3: the solution is a \(2\times 2\) symmetric matrix (3 values per vertex)
    • 4: the solution is a \(2\times 2\) matrix (4 values per vertex)
  • nbv means the number of vertices
  • U_i_j_k is the value of the component i of the solution k at vertex j on the associated mesh

Metric file

A metric file can be of two types, isotropic or anisotropic.

The isotropic file is such that

1
2
3
4
[Number of vertices](int) 1
h_0(double)
...
h_nv(double)
  • nv is the number of vertices
  • h_i is the wanted mesh size near the vertex i on associated mesh.

The metric is \(\mathcal{M}_i = h_i^{-2}I\) where \(I\) is the identity matrix.

The anisotropic file is such that

1
2
3
4
[Number of vertices](int) 3
a11_0(double) a21_0(double) a22_0(double)
...
a11_nv(double) a21_nv(double) a22_nv(double)
  • nv is the number of vertices
  • a11_i, a21_i and a22_i represent metric \(\mathcal{M}_i = \left(\begin{array}{cc}a_{11,i} & a_{12,i}\\a{12}_i & a_{22,i}\end{array}\right)\) which define the wanted size in a vicinity of the vertex i such that \(h\) in direction \(u \in \mathbb{R}^2\) is equal to \(|u|/\sqrt{u\cdot\mathcal{M}_i\, u}\), where \(\cdot\) is the dot product in \(\mathbb{R}^2\), and \(|\cdot|\) is the classical norm.

List of AM_FMT, AMDBA Meshes

The mesh is only composed of triangles and can be defined with the help of the following two integers and four arrays:

  • nbt the number of triangles
  • nbv the number of vertices
  • nu(1:3, 1:nbt) an integer array giving the three vertex numbers counterclockwise for each triangle
  • c(1:2, 1:nbv) a real array giving tje two coordinates of each vertex
  • refs(1:nbv) an integer array giving the reference numbers of the vertices
  • reft(1:nbt) an integer array giving the reference numbers of the triangles

AM_FMT Files

In Fortran the am_fmt files are read as follows:

1
2
3
4
5
6
7
open (1, file='xxx.am_fmt', form='formatted', status='old')
read (1, *) nbv, nbt
read (1, *) ((nu(i, j), i=1, 3), j=1, nbt)
read (1, *) ((c(i, j), i=1, 2), j=1, nbv)
read (1, *) ( reft(i), i=1, nbt)
read (1, *) ( refs(i), i=1, nbv)
close(1)

AM Files

In Fortran the am files are read as follows:

1
2
3
4
5
6
7
open (1, file='xxx.am', form='unformatted', status='old')
read (1, *) nbv, nbt
read (1) ((nu(i, j), i=1, 3), j=1, nbt),
& ((c(i, j), i=1, 2), j=1, nbv),
& (reft(i), i=1, nbt),
& (refs(i), i=1, nbv)
close(1)

AMDBA Files

In Fortran the amdba files are read as follows:

1
2
3
4
5
open (1, file='xxx.amdba', form='formatted', status='old')
read (1, *) nbv, nbt
read (1, *) (k, (c(i, k), i=1, 2), refs(k), j=1, nbv)
read (1, *) (k, (nu(i, k), i=1, 3), reft(k), j=1, nbt)
close(1)

msh Files

First, we add the notions of boundary edges

  • nbbe the number of boundary edge
  • nube(1:2, 1:nbbe) an integer array giving the two vertex numbers of boundary edges
  • refbe(1:nbbe) an integer array giving the reference numbers of boundary edges

In Fortran the msh files are read as follows:

1
2
3
4
5
6
open (1, file='xxx.msh', form='formatted', status='old')
read (1, *) nbv, nbt, nbbe
read (1, *) ((c(i, k), i=1, 2), refs(k), j=1, nbv)
read (1, *) ((nu(i, k), i=1, 3), reft(k), j=1, nbt)
read (1, *) ((ne(i, k), i=1, 2), refbe(k), j=1, nbbe)
close(1)

ftq Files

In Fortran the ftq files are read as follows:

1
2
3
4
5
open(1,file='xxx.ftq',form='formatted',status='old')
read (1,*) nbv,nbe,nbt,nbq
read (1,*) (k(j),(nu(i,j),i=1,k(j)),reft(j),j=1,nbe)
read (1,*) ((c(i,k),i=1,2),refs(k),j=1,nbv)
close(1)

where if k(j) = 3 when the element j is a triangle and k(j) = 4 when the the element j is a quadrilateral.

sol and solb files

With the keyword savesol, we can store a scalar functions, a scalar finite element functions, a vector fields, a vector finite element fields, a symmetric tensor and a symmetric finite element tensor.

Such format is used in medit.

Extension file .sol

The first two lines of the file are :

  • MeshVersionFormatted 0
  • Dimension [DIM](int)

The following fields begin with one of the following keyword: SolAtVertices, SolAtEdges, SolAtTriangles, SolAtQuadrilaterals, SolAtTetrahedra, SolAtPentahedra, SolAtHexahedra.

In each field, we give then in the next line the number of elements in the solutions (SolAtVertices: number of vertices, SolAtTriangles: number of triangles, …). In other lines, we give the number of solutions, the type of solution (1: scalar, 2: vector, 3: symmetric tensor). And finally, we give the values of the solutions on the elements.

The file must be ended with the keyword End.

The real element of symmetric tensor :

(40) \[\begin{split}ST^{3d}=\left( \begin{array}{ccc} ST_{xx}^{3d} & ST_{xy}^{3d} & ST_{xz}^{3d}\\ ST_{yx}^{3d} & ST_{yy}^{3d} & ST_{yz}^{3d} \\ ST_{zx}^{3d} & ST_{zy}^{3d} & ST_{zz}^{3d} \end{array} \right) \quad ST^{2d}= \left( \begin{array}{cc} ST_{xx}^{2d} & ST_{xy}^{2d} \\ ST_{yx}^{2d} & ST_{yy}^{2d} \end{array} \right)\end{split}\]

stored in the extension .sol are respectively \(ST_{xx}^{3d}, ST_{yx}^{3d}, ST_{yy}^{3d}, ST_{zx}^{3d}, ST_{zy}^{3d}, ST_{zz}^{3d}\) and \(ST_{xx}^{2d}, ST_{yx}^{2d}, ST_{yy}^{2d}\)

An example of field with the keyword SolAtTetrahedra:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
SolAtTetrahedra
[Number of tetrahedra](int)
[Number of solutions](int) [Type of solution 1](int) ... [Type of soution nt](int)

U_1_1_1(double) ... U_nrs_1_1(double)
...
U_1_ns_1(double) ... U_(nrs_k)_ns_1(double)

...

U_1_1_nt(double) ... U_nrs_1_nt(double)
...
U_1_ns_nt(double) ... U_(nrs_k)_ns_nt(double)
  • ns is the number of solutions
  • typesol_k, type of the solution number k
    • typesol_k = 1 the solution k is scalar
    • typesol_k = 2 the solution k is vectorial
    • typesol_k = 3 the solution k is a symmetric tensor or symmetric matrix
  • nrs_k is the number of real to describe solution k
    • nrs_k = 1 if the solution k is scalar
    • nrs_k = dim if the solution k is vectorial (dim is the dimension of the solution)
    • nrs_k = dim*(dim+1)/2 if the solution k is a symmetric tensor or symmetric matrix
  • U_i_j_^k is a real equal to the value of the component i of the solution k at tetrahedron j on the associated mesh

The format .solb is the same as format .sol but in binary (read/write is faster, storage is less).

A real scalar functions \(f1\), a vector fields \(\mathbf{\Phi} = [\Phi1, \Phi2, \Phi3]\) and a symmetric tensor \(ST^{3d}\) (40) at the vertices of the three dimensional mesh Th3 is stored in the file f1PhiTh3.sol using :

1
savesol("f1PhiST3dTh3.sol", Th3, f1, [Phi(1), Phi(2), Phi(3)], VV3, order=1);

where \(VV3 = [ST_{xx}^{3d}, ST_{yx}^{3d}, ST_{yy}^{3d}, ST_{zx}^{3d}, ST_{zy}^{3d}, ST_{zz}^{3d}]\).

For a two dimensional mesh Th, A real scalar functions \(f2\), a vector fields \(\mathbf{\Psi} = [\Psi1, \Psi2]\) and a symmetric tensor \(ST^{2d}\) (40) at triangles is stored in the file f2PsiST2dTh3.solb using :

1
savesol("f2PsiST2dTh3.solb", Th, f2, [Psi(1), Psi(2)], VV2, order=0);

where \(VV2 = [ST_{xx}^{2d}, ST_{yx}^{2d}, ST_{yy}^{2d}]\)

The arguments of savesol functions are the name of a file, a mesh and solutions. These arguments must be given in this order.

The parameters of this keyword are :

  • order = 0 is the solution is given at the center of gravity of elements. 1 is the solution is given at the vertices of elements.

In the file, solutions are stored in this order : scalar solutions, vector solutions and finally symmetric tensor solutions.

Adding a new finite element

Some notations

For a function \(\boldsymbol{f}\) taking value in \(\mathbb{R}^{N},\, N=1,2,\cdots\), we define the finite element approximation \(\Pi_h\boldsymbol{f}\) of \(\boldsymbol{f}\).

Let us denote the number of the degrees of freedom of the finite element by \(NbDoF\). Then the \(i\)-th base \(\boldsymbol{\omega}^{K}_{i}\) (\(i=0,\cdots,NbDoF-1\)) of the finite element space has the \(j\)-th component \(\mathbf{\omega}^{K}_{ij}\) for \(j=0,\cdots,N-1\).

The operator \(\Pi_{h}\) is called the interpolator of the finite element.

We have the identity \(\boldsymbol{\omega}^{K}_{i} = \Pi_{h} \boldsymbol{\omega}^{K}_{i}\).

Formally, the interpolator \(\Pi_{h}\) is constructed by the following formula:

(41) \[\Pi_{h} \boldsymbol{f} = \sum_{k=0}^{\mathtt{kPi}-1} \alpha_k \boldsymbol{f}_{j_{k}}(P_{p_{k}}) \boldsymbol{\omega}^{K}_{i_{k}}\]

where \(P_{p}\) is a set of \(npPi\) points,

In the formula (41), the list \(p_{k},\, j_{k},\, i_{k}\) depend just on the type of finite element (not on the element), but the coefficient \(\alpha_{k}\) can be depending on the element.

Tip

Classical scalar Lagrange finite element

With the classical scalar Lagrange finite element, we have \(\mathtt{kPi}=\mathtt{npPi}=\mathtt{NbOfNode}\) and

  • \(P_{p}\) is the point of the nodal points.
  • the \(\alpha_k=1\), because we take the value of the function at the point \(P_{k}\).
  • \(p_{k}=k\) , \(j_{k}=k\) because we have one node per function.
  • \(j_{k}=0\) because \(N=1\).

Tip

The Raviart-Thomas finite element

(42) \[RT0_{h} = \{ \mathbf{v} \in H(div) / \forall K \in \mathcal{T}_{h} \quad \mathbf{v}_{|K}(x,y) = \vecttwo{\alpha_{K}}{\beta_{K}} + \gamma_{K}\vecttwo{x}{y} \}\]

The degrees of freedom are the flux through an edge \(e\) of the mesh, where the flux of the function \(\mathbf{f} : \mathbb{R}^2 \longrightarrow \mathbb{R}^2\) is \(\int_{e} \mathbf{f}.n_{e}\), \(n_{e}\) is the unit normal of edge \(e\) (this implies a orientation of all the edges of the mesh, for example we can use the global numbering of the edge vertices and we just go to small to large number).

To compute this flux, we use a quadrature formula with one point, the middle point of the edge. Consider a triangle \(T\) with three vertices \((\mathbf{a},\mathbf{b},\mathbf{c})\).

Let denote the vertices numbers by \(i_{a},i_{b},i_{c}\), and define the three edge vectors \(\mathbf{e}^{0},\mathbf{e}^{1},\mathbf{e}^{2}\) by \(sgn(i_{b}-i_{c})(\mathbf{b}-\mathbf{c})\), \(sgn(i_{c}-i_{a})(\mathbf{c}-\mathbf{a})\), \(sgn(i_{a}-i_{b})(\mathbf{a}-\mathbf{b})\).

The three basis functions are:

\[\boldsymbol{\omega}^{K}_{0}= \frac{sgn(i_{b}-i_{c})}{2|T|}(x-a),\quad \boldsymbol{\omega}^{K}_{1}= \frac{sgn(i_{c}-i_{a})}{2|T|}(x-b),\quad \boldsymbol{\omega}^{K}_{2}= \frac{sgn(i_{a}-i_{b})}{2|T|}(x-c),\]

where \(|T|\) is the area of the triangle \(T\).

So we have \(N=2\), \(\mathtt{kPi}=6; \mathtt{npPi}=3;\) and:

  • \(P_{p} = \left\{\frac{\mathbf{b}+\mathbf{c}}{2}, \frac{\mathbf{a}+\mathbf{c}}{2}, \frac{\mathbf{b}+\mathbf{a}}{2} \right\}\)
  • \(\alpha_{0}= - \mathbf{e}^{0}_{2}, \alpha_{1}= \mathbf{e}^{0}_{1}\),
    \(\alpha_{2}= - \mathbf{e}^{1}_{2}, \alpha_{3}= \mathbf{e}^{1}_{1}\), \(\alpha_{4}= - \mathbf{e}^{2}_{2}, \alpha_{5}= \mathbf{e}^{2}_{1}\) (effectively, the vector \((-\mathbf{e}^{m}_{2}, \mathbf{e}^{m}_{1})\) is orthogonal to the edge \(\mathbf{e}^{m}= (e^m_{1},e^m_{2})\) with a length equal to the side of the edge or equal to \(\int_{e^m} 1\)).
  • \(i_{k}=\{0,0,1,1,2,2\}\),
  • \(p_{k}=\{0,0,1,1,2,2\}\) , \(j_{k}=\{0,1,0,1,0,1,0,1\}\).

Which class to add?

Add file FE_ADD.cpp in directory FreeFem-sources/src/femlib for example first to initialize :

1
2
3
4
5
6
7
8
9
#include "error.hpp"
#include "rgraph.hpp"
using namespace std;
#include "RNM.hpp"
#include "fem.hpp"
#include "FESpace.hpp"
#include "AddNewFE.h"

namespace Fem2D { ... }

Then add a class which derive for public TypeOfFE like:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
class TypeOfFE_RTortho : public TypeOfFE { public:
    static int Data[]; //some numbers
    TypeOfFE_RTortho():
    TypeOfFE(
        0+3+0,  //nb degree of freedom on element
        2,      //dimension N of vectorial FE (1 if scalar FE)
        Data,   //the array data
        1,      //nb of subdivision for plotting
        1,      //nb of sub finite element (generaly 1)
        6,      //number kPi of coef to build the interpolator
        3,      //number npPi of integration point to build interpolator
        0       //an array to store the coef \alpha_k to build interpolator
        //here this array is no constant so we have
        //to rebuilt for each element
    )
    {
        const R2 Pt[] = {R2(0.5, 0.5), R2(0.0, 0.5), R2(0.5, 0.0) };
        // the set of Point in hat{K}
        for (int p = 0, kk = 0; p < 3; p++){
            P_Pi_h[p] = Pt[p];
            for (int j = 0; j < 2; j++)
                pij_alpha[kk++] = IPJ(p, p, j);
        }
    } //definition of i_k, p_k, j_k in interpolator

    void FB(const bool *watdd, const Mesh &Th, const Triangle &K,
        const R2 &PHat, RNMK_ &val) const;

    void Pi_h_alpha(const baseFElement &K, KN_<double> &v) const;
} ;

where the array data is formed with the concatenation of five array of size NbDoF and one array of size N.

This array is:

1
2
3
4
5
6
7
8
9
int TypeOfFE_RTortho::Data[] = {
    //for each df 0, 1, 3:
    3, 4, 5, //the support of the node of the df
    0, 0, 0, //the number of the df on the node
    0, 1, 2, //the node of the df
    0, 0, 0, //the df come from which FE (generally 0)
    0, 1, 2, //which are the df on sub FE
    0, 0
}; //for each component j=0, N-1 it give the sub FE associated

where the support is a number \(0,1,2\) for vertex support, \(3,4,5\) for edge support, and finally \(6\) for element support.

The function to defined the function \(\boldsymbol{\omega}^{K}_{i}\), this function return the value of all the basics function or this derivatives in array val, computed at point Phat on the reference triangle corresponding to point R2 P=K(Phat); on the current triangle K.

The index \(i,j,k\) of the array \(val(i,j,k)\) correspond to:

  • \(i\) is the basic function number on finite element \(i \in [0,NoF[\)

  • \(j\) is the value of component \(j \in [0,N[\)

  • \(k\) is the type of computed value \(f(P),dx(f)(P), dy(f)(P), ...\ i \in [0,\mathtt{last\_operatortype}[\).

    Note

    For optimization, this value is computed only if whatd[k] is true, and the numbering is defined with

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    enum operatortype {
    op_id = 0,
    op_dx = 1, op_dy = 2,
    op_dxx = 3,op_dyy = 4,
    op_dyx = 5,op_dxy = 5,
    op_dz = 6,
    op_dzz = 7,
    op_dzx = 8, op_dxz = 8,
    op_dzy = 9, op_dyz = 9
    };
    const int last_operatortype = 10;
    

The shape function:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
void TypeOfFE_RTortho::FB(const bool *whatd, const Mesh &Th, const Triangle & K,
    const R2 &PHat,RNMK_ &val) const
{
    R2 P(K(PHat));
    R2 A(K[0]), B(K[1]), C(K[2]);
    R l0 = 1 - P.x-P.y;
    R l1 = P.x, l2 = P.y;
    assert(val.N() >= 3);
    assert(val.M() == 2);
    val = 0;
    R a = 1./(2*K.area);
    R a0 = K.EdgeOrientation(0) * a;
    R a1 = K.EdgeOrientation(1) * a;
    R a2 = K.EdgeOrientation(2) * a;

    if (whatd[op_id]){ //value of the function
        assert(val.K() > op_id);
        RN_ f0(val('.', 0,0)); //value first component
        RN_ f1(val('.', 1,0)); //value second component
        f1[0] = (P.x - A.x)*a0;
        f0[0] = -(P.y - A.y)*a0;

        f1[1] = (P.x - B.x)*a1;
        f0[1] = -(P.y - B.y)*a1;

        f1[2] = (P.x - C.x)*a2;
        f0[2] = -(P.y - C.y)*a2;
    }

    if (whatd[op_dx]){ //value of the dx of function
        assert(val.K() > op_dx);
        val(0,1,op_dx) = a0;
        val(1,1,op_dx) = a1;
        val(2,1,op_dx) = a2;
    }
    if (whatd[op_dy]){
        assert(val.K() > op_dy);
        val(0,0,op_dy) = -a0;
        val(1,0,op_dy) = -a1;
        val(2,0,op_dy) = -a2;
    }

    for (int i = op_dy; i < last_operatortype; i++)
        if (whatd[op_dx])
            assert(op_dy);
}

The function to defined the coefficient \(\alpha_{k}\):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
void TypeOfFE_RT::Pi_h_alpha(const baseFElement &K, KN_<double> &v) const
{
    const Triangle &T(K.T);

    for (int i = 0, k = 0; i < 3; i++){
        R2 E(T.Edge(i));
        R signe = T.EdgeOrientation(i) ;
        v[k++] = signe*E.y;
        v[k++] = -signe*E.x;
    }
}

Now , we just need to add a new key work in FreeFEM.

Two way, with static or dynamic link so at the end of the file, we add:

With dynamic link it is very simple (see section Dynamical link), just add before the end of FEM2d namespace:

1
2
3
   static TypeOfFE_RTortho The_TypeOfFE_RTortho;
   static AddNewFE("RT0Ortho", The_TypeOfFE_RTortho);
} //FEM2d namespace

Try with ./load.link command in examples++-load/ and see BernardiRaugel.cpp or Morley.cpp new finite element examples.

Otherwise with static link (for expert only), add

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
//let the 2 globals variables
static TypeOfFE_RTortho The_TypeOfFE_RTortho;
//the name in freefem
static ListOfTFE typefemRTOrtho("RT0Ortho", &The_TypeOfFE_RTortho);

//link with FreeFEM do not work with static library .a
//so add a extern name to call in init_static_FE
//(see end of FESpace.cpp)
void init_FE_ADD() { };
//end
} //FEM2d namespace

To inforce in loading of this new finite element, we have to add the two new lines close to the end of files src/femlib/FESpace.cpp like:

1
2
3
4
5
6
7
8
//correct problem of static library link with new make file
void init_static_FE()
{ //list of other FE file.o
    extern void init_FE_P2h() ;
    init_FE_P2h() ;
    extern void init_FE_ADD(); //new line 1
    init_FE_ADD(); //new line 2
}

and now you have to change the makefile.

First, create a file FE_ADD.cpp contening all this code, like in file src/femlib/Element_P2h.cpp, after modify the Makefile.am by adding the name of your file to the variable EXTRA_DIST like:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# Makefile using Automake + Autoconf
# ----------------------------------
# Id

# This is not compiled as a separate library because its
# interconnections with other libraries have not been solved.

EXTRA_DIST=BamgFreeFem.cpp BamgFreeFem.hpp CGNL.hpp CheckPtr.cpp        \
ConjuguedGradrientNL.cpp DOperator.hpp Drawing.cpp Element_P2h.cpp      \
Element_P3.cpp Element_RT.cpp fem3.hpp fem.cpp fem.hpp FESpace.cpp      \
FESpace.hpp FESpace-v0.cpp FQuadTree.cpp FQuadTree.hpp gibbs.cpp        \
glutdraw.cpp gmres.hpp MatriceCreuse.hpp MatriceCreuse_tpl.hpp          \
MeshPoint.hpp mortar.cpp mshptg.cpp QuadratureFormular.cpp              \
QuadratureFormular.hpp RefCounter.hpp RNM.hpp RNM_opc.hpp RNM_op.hpp    \
RNM_tpl.hpp     FE_ADD.cpp

and do in the FreeFEM root directory

1
2
3
autoreconf
./reconfigure
make

For codewarrior compilation add the file in the project an remove the flag in panal PPC linker FreeFm++ Setting Dead-strip Static Initializition Code Flag.

Table of content