AFEPack 使用 Tutorial（一）：解二维泊松方程

申明

这只是阅读了 AFEPack 自带的文档，并结合简单的 Poisson Equation 的理解。

今后需要使用 AFEPack，就将这个做一个学习笔记。关于 AFEPack 的编译和安装方面笔记，以后再写。

关于 AFEPack

AFEPack 是北大理学院李若教授写的一个开源 PDE 求解软件包。功能非常强大。

路径设置

AFEPack 编译完成后，需要设置以下的路径信息。

AFEPack 基础路径

该路径用于定义 AFEPack 源码所在。假设我们将 AFEPack 安装在 /usr/local/src/AFEPack 这个路径，我们可以通过配置 ~/.bashrc 文件，增加如下的定义。

export AFEPACK_PATH="/usr/local/src/AFEPack"
           

请务必和自己机器上的路径保持一致。

这个路径在 AFEPack 提供的标准例程中，是在对应的 examples 的 run 文件中。

AFEPack 模板单元路径

用于定义网格文件中使用的模板单元所在路径。

这个需要根据自己使用模板进行添加。一般常用的模板是三角形（triangle），那么我们的配置如下：

export AFEPACK_TEMPLATE_PATH="$AFEPACK_PATH/template/triangle"
           

如果使用其他模板，可以自己添加对应的模板单元。AFEPack 目前提供的模板单元比较多，如下所示：

对隐私信息进行删除处理，请见谅。

AFEPack 库文件

用于定义 AFEPack so 文件所在的路径。

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib
           

AFEPack Step by Step Tutorial 1: solve Poisson Equation

Poisson Equation 的源码在 example/possion_equation 目录下。

Poisson Equation 使用

网格文件定义为 D.d。

编译 EasyMesh 网格文件

~/AFEPack/example/possion_equation$ easymesh D
Working.  Please wait !
Improving the grid quality.  Please wait !
Renumerating nodes, elements and sides. Please wait !
Processing material marks. Please wait !
           

这样，我们就可以获得 D.e、D.n 和 D.s。

网格如上图所示。

编译程序

使用 make 进行编译，具体编译信息可以自己查看 Makefile。

~/AFEPack/example/possion_equation$ make
g++ -c -o possion_equation.o possion_equation.cpp -I. -I/usr/include/deal.II/ -std=c++11 -g -O2 -fPIC -D__SERIALIZATION__ -DMULTITHREAD -pthread -fte
mplate-depth-256 -I/root/Pkg/AFEPack/library/include
g++ -o main  possion_equation.o  -ldeal_II -L/usr/lib/x86_64-linux-gnu/ -ltbb -ldl -lm  -pthread -L../../library/lib -lAFEPack -ldeal_II -L/usr/lib/x
86_64-linux-gnu/ -ltbb -ldl -lm  -pthread
           

运行结果

~/AFEPack/example/possion_equation$ ./run D
Reading easymesh data file ...
        reading node data ... OK!
        reading side data ... OK!
        reading element data ... OK!
AFEPack library file found: /root/Pkg/AFEPack/template/triangle/triangle.tmp_geo
AFEPack library file opened: /root/Pkg/AFEPack/template/triangle/triangle.tmp_geo
Opening shared library /root/Pkg/AFEPack/template/triangle/./triangle.geometry.so ...
        tried /root/Pkg/AFEPack/template/triangle/./triangle.geometry.so: success
AFEPack library file found: /root/Pkg/AFEPack/template/triangle/triangle.crd_trs
AFEPack library file opened: /root/Pkg/AFEPack/template/triangle/triangle.crd_trs
Opening shared library /root/Pkg/AFEPack/template/triangle/./triangle.geometry.so ...
        tried /root/Pkg/AFEPack/template/triangle/./triangle.geometry.so: success
AFEPack library file found: /root/Pkg/AFEPack/template/triangle/triangle.1.tmp_dof
AFEPack library file opened: /root/Pkg/AFEPack/template/triangle/triangle.1.tmp_dof
AFEPack library file found: /root/Pkg/AFEPack/template/triangle/triangle.1.bas_fun
AFEPack library file opened: /root/Pkg/AFEPack/template/triangle/triangle.1.bas_fun
Opening shared library /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Opening shared library /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Opening shared library /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /root/Pkg/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Building degree of freedom for the FEM space ...
        total 506 degree of freedom found.
AMGSolver initializing ... OK! grid levels: 3
AMGSolver begin with initial residual 19.7721 ...
        converge with residual 3.02858e-08 at step 6.

L2 error = 0.00396625
           

从上面的反馈可以看出，最终的 L 2 L^2 L2 Norm 误差。

泊松方程解如上图所示。

下面，我们开始进一步了解 AFEPack 如何实现 Poisson Equation 求解。

Poisson Equation 分析

泊松方程定义为 − Δ u = f -\Delta u = f −Δu=f。我们假设其真解为 u = S i n ( π x ) S i n ( 2 π y ) u=Sin(\pi x)Sin(2\pi y) u=Sin(πx)Sin(2πy)，这样带入到泊松方程中。

我们可以先计算 u x x u_{xx} uxx​， u x = π C o s ( π x ) S i n ( 2 π y ) u_{x}=\pi Cos(\pi x)Sin(2\pi y) ux​=πCos(πx)Sin(2πy)， u x x = − π 2 S i n ( π x ) S i n ( 2 π y ) u_{xx}=-\pi^2 Sin(\pi x)Sin(2\pi y) uxx​=−π2Sin(πx)Sin(2πy)。

同理 u y = 2 π S i n ( π x ) C o s ( 2 π y ) u_{y}=2\pi Sin(\pi x)Cos(2\pi y) uy​=2πSin(πx)Cos(2πy)， u y y = − 4 π 2 S i n ( π x ) S i n ( 2 π y ) u_{yy}=-4\pi^2Sin(\pi x)Sin(2\pi y) uyy​=−4π2Sin(πx)Sin(2πy)。

这样， f = − Δ u = − ( u x x + u y y ) = π 2 S i n ( π x ) S i n ( 2 π y ) + 4 π 2 S i n ( π x ) S i n ( 2 π y ) = 5 π 2 S i n ( π x ) S i n ( 2 π y ) f=-\Delta u=-(u_{xx}+u_{yy})=\pi^2 Sin(\pi x)Sin(2\pi y)+4\pi^2Sin(\pi x)Sin(2\pi y)=5\pi^2Sin(\pi x)Sin(2\pi y) f=−Δu=−(uxx​+uyy​)=π2Sin(πx)Sin(2πy)+4π2Sin(πx)Sin(2πy)=5π2Sin(πx)Sin(2πy)

现在，一个完整的泊松问题定义如下

{ − Δ u = f u ∣ ∂ Ω = u b u b ( x , y ) = s i n ( π ∗ x ) ∗ s i n ( 2 ∗ π ∗ y )   i f   ( x , y ) ∈ ∂ Ω f ( x , y ) = 5 π 2 S i n ( π x ) S i n ( 2 π y )   f o r   ( x , y ) ∈ Ω \left\{ \begin{aligned} -\Delta u & = & f &\\ u |_{\partial \Omega} & = & u_b & \\ u_b(x,y) &= &sin(\pi*x)*sin(2*\pi*y) & \ if \ (x,y) \in \partial \Omega\\ f(x,y) &= & 5\pi^2Sin(\pi x)Sin(2\pi y) & \ for \ (x,y) \in \Omega \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​−Δuu∣∂Ω​ub​(x,y)f(x,y)​====​fub​sin(π∗x)∗sin(2∗π∗y)5π2Sin(πx)Sin(2πy)​ if (x,y)∈∂Ω for (x,y)∈Ω​

其中 u b u_b ub​ 是边界条件， f f f 是右端项。

Step 1: 网格文件

AFEPack 文档中提到支持 4 4 4 种网格文件（Mesh data）：MIT 开源的 easymesh、Dr Bojan Niceno 的 2d mesh generator、bamg、gmesh。

当前用的比较多的是 easymesh。数据文件保存为三个文件：后缀为 n 的文件，该文件保存节点信息（the data file to store the node information）、后缀为 s 的文件，该文件保存边信息（the data file to store the side information）和后缀为 e 的文件，该文件保存元素信息（the data file to store the element information）。

使用网格文件

MasyMesh

头文件

需要包含下面的头文件。

#include <AFEPack/EasyMesh.h>
           

对应源文件在 library/include/EasyMesh.h 和 library/src/EasyMesh.cpp。

读取

假设我们的 Mesh 文件名字为 a.n、a.s 和 a.e。

EasyMesh easy_mesh; // construct a EasyMesh object named "easy_mesh"
easy_mesh.readData("a"); // read in the mesh data in files "a.n", "a.s", "a.e"
           

写入

easy_mesh.writeData("b"); // write the mesh data in files "b.n", "b.s", "b.e"
           

RegularMesh

The hierarchy tree class HGeometryTree<2> can read in easymesh data file as its data source. Another class RegularMesh<2> also support outputting the data into these kind of file format.

头文件

需要包含下面的头文件。

#include <AFEPack/Geometry.h>
#include <AFEPack/HGeometry.h>
           

对应源文件在 library/include/HGeometry.h 和 library/src/HGeometry.cpp。

读取

HGeometryTree<2> h_geometry_tree; // construct a hierarchy geometry tree object named "h_geometry_tree"
h_geometry_tree.readEasyMesh("a"); // read in easymesh data file as macro elements
IrregularMesh<2> irregular_mesh(hierarchy_tree); // construct a irregular mesh object "irregular_mesh"
irregular_mesh.globalRefine(2); // global refine the irregular mesh for 2 times
irregular_mesh.semiregularize(); // turn the mesh into a semiregular one
irregular_mesh.regularize(); // generate regular mesh from the irregular mesh
RegularMesh<2>& regular_mesh = irregular_mesh.regularMesh(); // this is the regular mesh from the irregular mesh
           

写入

regular_mesh.writeEasyMesh("b")// write the data into files "b.n", "b.s", "b.e"
           

bamg

a very powerful mesh generator by Fridiric Hecht(France). Class DBMesh is designed for reading such mesh data file.

头文件

需要包含下面的头文件。

#include <AFEPack/DBMesh.h>
           

读取

DBMesh dbmesh; // construct a DBMesh object named "dbmesh"
dbmesh.readData("a_g.msh"); // read in the data file generated by bamg
           

写入

Mesh<2,2> mesh; // construct a Mesh<2,2> object named "mesh"
dbmesh.generateMesh(mesh); // convert the dbmesh data file to the library internal data format
           

gmsh

3D Mesh。Class GmshMesh is also derivated from class SimplestMesh<3,3>.

头文件

需要包含下面的头文件。

#include <AFEPack/GmshMesh.h>
           

读取

GmshMesh gmshmesh; // construct a GmshMesh object named "gmshmesh"
gmshmesh.readData("a.msh"); // read in the data file generated by gmsh
           

写入

Mesh<3,3> mesh; // construct a Mesh<3,3> object named "mesh"
gmshmesh.generateMesh(mesh); // convert the gmsh data file to the library internal data format
           

Step 2: 构建模板单元

Now let’s construct a template element, which is a main task to construct a finite element space. We will construct a linear Lagrangian element on a triangle, the most often used one in coding. Most of the work needed in this operation should be achieved by hand. We divide it into several steps.

几何信息

The geometry information includes the shape description of the triangle, a function to calculate the volume of the triangle and the numerical quadrature information on the triangle. So the data file have three part. The following is an example to describe the geometry information on a triangle element. The texts following a “#” are comment(When you write such a file of yourself, please remove those comments because until now the library doesn’t support comments in this file).

# the information of the function to calculate the volume of the triangle
./triangle.geometry.so                        # the function to calculate the volume of the triangle is in this shared library
    triangle_volume                           # the name of the function to calculate the volume of the triangle

# the geometry shape information of the triangle
3                                             # there are 3 points in this geometry
    0.0    0.0                                # these three lines are the cooridinate of the 3 points
    1.0    0.0
    0.0    1.0

3                                             # there are 3 vertices in this geometry
    0                                         # the information of number 0 vertex is
        1    0                                # 1 point, that's the number 0 point
        1    0                                # 1 boundary, that's the number 0 point
    1
        1    1
        1    1
    2
        1    2
        1    2

3                                             # there are 3 sides in this geometry
    0                                         # the information of number 0 side is
        2    1    2                           # 2 vertices, saying the number 1 and number 2 vertices
        2    1    2                           # 2 boundarys, saying the number 1 and number 2 vertices
    1
        2    0    2
        2    0    2
    2
        2    0    1
        2    0    1

1                                             # there is 1 element in this geometry
    0                                         # the information of number 0 element is
        3    0    1    2                      # 3 vertices, saying the number 0, 1 and 2 vertices
        3    0    1    2                      # 3 boundarys, saying the number 0, 1 and 2 boundarys

# the numerical quadrature information on the triangle
2                                             # there are 2 quadrature information on this geometry
    2                                         # for the first one, algebraic accuracy is 2
    1                                         # 1 quadrature point
    0.333333    0.333333    1.0               # coordinate of the quadrature point is (0.333333, 0.333333), weight 1.0

    3                                         # for the second one, algebric accuracy is 3
    3                                         # 3 quadrature points
    0.500000    0.000000    0.333333          # coordinate of the quadrature point is (0.500000, 0.000000), weight 0.333333
    0.500000    0.500000    0.333333          # coordinate of the quadrature point is (0.500000, 0.500000), weight 0.333333
    0.000000    0.500000    0.333333          # coordinate of the quadrature point is (0.000000, 0.500000), weight 0.333333
           

Let’s save those information in a file named, such as “triangle.tmp_geo”.

AFEPack 内部实现

Obviously, we need a C source code file for the shared library “triangle.geometry.so”. That’s easy. We edit such a source file as following:

* triangle.geometry.c */

double triangle_volume(const double * vertex) {
   /* we know the volume of the triangle is 0.5 */
    return 0.5;
}
/* end of file */
           

AFEPack 将三角形的几何信息编译成为 triangle.geometry。

$ gcc -c triangle.geometry.c
$ gcc -shared -o triangle.geometry.so triangle.geometry.
           

使用

我们只需要按照下面步骤来构建三角形模板即可。

#include <AFEPack/Geometry.h>
... ...
TemplateGeometry<2>    triangle_template_geometry;
triangle_template_geometry.readData("triangle.tmp_geo");
... ...
           

坐标变换

The coordinate transformation is the information to map the point from the template geometry to the real element, and vice visa. It’s in fact four C function in a shared library.

AFEPack 内部实现

It’s in fact four C function in a shared library. We need prepare a text file as

./triangle.geometry.so                            # the shared library in which those functions included
    local_to_global                               # the function to map a point from template to real element
    global_to_local                               # the function to map a point from real to template element
    local_to_global_jacobian                      # the jacobian determinant of the function to map a point from template to real element
    global_to_local_jacobian                      # the jacobian determinant of the function to map a point from real to template element
           

Save this file, and give it a name as “triangle.crd_trs”. The corresponding C source file is as:

/* triangle.crd_trs.c */

void local_to_global(const double * point_on_template_element,
        const double ** vertices_of_template_element,
        const double ** vertixes_of_real_element,
        double * point_on_real_element)
{
    /* write code here to map the point on the template element
       to the point on the real element */
}

void global_to_local(const double * point_on_real_element,
        const double ** vertices_of_real_element,
        const double ** vertices_of_template_element,
        double * point_on_template_element)
{
    /* write code here to map the point on the real element
       to the point on the template element */
}

double local_to_global_jacobian(const double * point_on_template_element,
        const double ** vertices_of_template_element,
        const double ** vertices_of_real_element)
{
    /* calculate the jacobian determinant */
    return jacobian_determinant;
}

double global_to_local_jacobian(const double * point_on_real_element,
        const double ** vertices_of_template_element,
        const double ** vertices_of_real_element)
{
    /* calculate the jacobian determinant */
    return jacobian_determinant;
}

/* end of file */

           

As before, we compile the file into shared library as

$ gcc -c triangle.crd_trs.c
$ gcc -shared -o triangle.geometry.so 
           

使用

#include <AFEPack/Geometry.h>
... ...
CoordTransform<2> triangle_coord_transform;
triangle_coord_transform.readData("triangle.crd_trs");
... ...
           

自由度

This is a very simple but important information. It’s only a text file to tell how many degrees of freedom on certain geometries. Let’s see the following example:

3                                            # there are 3 lines of information in this file

0    0    1                                  # there are 1 dof on the number 0 of 0 dimension geometries
0    1    1                                  # there are 1 dof on the number 1 of 0 dimension geometries
0    2    1                                  # there are 1 dof on the number 2 of 0 dimension geometries
           

使用

Let’s save it into a file named as “triangle.1.tmp_dof” and construct the degree of freedom object

#include <AFEPack/Geometry.h>
#include <AFEPack/TemplateElement.h>
... ...
TemplateDOF<2> triangle_template_dof(triangle_template_geometry); // the information of template geometry is needed
triangle_template_dof.readData("triangle.1.tmp_dof");
           

基函数

Because the basis function are very flexible, we provide certain mechanism to identify different basis functions and then build basis functions from shape functions and then distribute degree of freedom on the whole triangulation. We implemented the identification of basis functions as a protocol, though not very tough.

AFEPack 内部实现

Basis Function Identification Protocol(BFIP)

In the impelmented basis functions in the library, we adopted the following format to identify different basis functions:

template <int DIM>
class BasisFunctionIdentity
{
private:
    unsigned int order;                // the order of the polynomial on the boundary of the element
    int alpha[DIM];                    // the interpolation operator description
    unsigned int flag;                 // additional flag, for special basis function
};
           

In which “order” is the order of the polynomial on the boundary of the element, such that the basis functions of different order can be differed. Generally, the basis function has certain interpolation meaning such as value, or value of certain derivative operator( ∂ ∂ x i \frac \partial {\partial x^i} ∂xi∂​ or ∂ ∂ n \frac \partial {\partial n} ∂n∂​, and so on). The array “alpha” gives such information, such as the index of the differential operator. The “flag” is kept as the further distinguish flag to some really like basis functions. Though we called it a “protocol”, you can freely choose the value in this struct to identify your own basis functions. We tried to make it like a protocol because we attend to consider the case when many contributors would like to provide their basis functions that to make all those basis functions from different source able to cooperate together smoothly.

After appointing the identity of the basis function, we need to further tell which geometry in the template element this basis function relys on and the interpolate point of the basis function. At last, the information needed is the shared library of the implementation of the basis function and the symbol name. The following is the piecewise linear basis function on the first vertex of the triangle template element:

0	0				 // the basis function rely on the 0-th 0-dim geometry
0.0     0.0				 // the interpolate point is (0.0, 0.0)
1	0 0	0			 // its identity is: order=1, alpha=[0,0], flag=0
./triangle.1.bas_fun.so			 // the filename of the shared library in which there are the functions needed
        lambda_1   gradient_lambda_1     // the name of the value and gradient function of this basis function
           

Then you can write the implementation of the basis function in a C source file, saying “triangle.1.bas_fun.c”, as

void lambda_1(const double * p, const double ** v, void * value)
{
        double * val = (double *)value;
        double area = (v[1][0] - v[0][0])*(v[2][1] - v[0][1])
                - (v[1][1] - v[0][1])*(v[2][0] - v[0][0]);
        val[0] = (v[1][0] - p[0])*(v[2][1] - p[1])
                - (v[1][1] - p[1])*(v[2][0] - p[0]);
        val[0] /= area;
};

void gradient_lambda_1(const double * p, const double ** v, void * value)
{
        double * val = (double *)value;
        double area = (v[1][0] - v[0][0])*(v[2][1] - v[0][1])
                - (v[1][1] - v[0][1])*(v[2][0] - v[0][0]);
        val[0] = (v[1][1] - v[2][1])/area;
        val[1] = (v[2][0] - v[1][0])/area;
};
           

Note the prototype of the value and gradient function of a basis function are as

void function_name(const double * p, const double ** v, void * value)
    p: the coordinate of the point
    v: the coordinate of the vertices of the element
    value: the returned value will be stored here
           

Of course, you should write other basis functions to construct a complete template element. All basis functions on a template element is administrated by a struct “class BasisFunctionAdmin<value_type, DIM, TDIM>”, it will read in information from a text file in which there are the information of all those basis function. The text file, assumed with name “triangle.1.bas_fun”, is as

3							// there are total 3 basis functions

0       0						// the information of the 1st basis function	
0.0     0.0
1       0 0     0
./triangle.1.bas_fun.so
        lambda_1        gradient_lambda_1

0       1						// the information of the 2nd basis function
1.0     0.0
1       0 0     0
./triangle.1.bas_fun.so
        lambda_2        gradient_lambda_2

0       2						// the information of the 3rd basis function
0.0     1.0
1       0 0     0
./triangle.1.bas_fun.so
        lambda_3        gradient_lambda_3
           

使用

#include <AFEPack/Geometry.h>
#include <AFEPack/TemplateElement.h>

...
BasisFunctionAdmin<double,2,2> triangle_basis_function(triangle_template_dof);
triangle_basis_function.readData("triangle.1.bas_fun");
           

小节

如果只是纯粹使用，可以直接看本小节。

With all those steps, we at last can construct a template element for a finite element space. Let’s collect those codes together as following:

#include <AFEPack/Geometry.h>
#include <AFEPack/TemplateElement.h>

int main()
{
        // construct template geometry
        TemplateGeometry<2> triangle_template_geometry;
        triangle_template_geometry.readData("triangle.tmp_geo");

        // construct coordinate transformation
        CoordTransform<2,2> triangle_coord_transform;
        triangle_coord_transform.readData("triangle.crd_trs");

        // construct template degree of freedom
        TemplateDOF<2> triangle_template_dof(triangle_template_geometry);
        triangle_template_dof.readData("triangle.1.tmp_dof");

        // construct basis functions
        BasisFunctionAdmin<double,2,2> triangle_basis_function(triangle_template_dof);
        triangle_basis_function.readData("triangle.1.bas_fun");

        // construct template element, 
        // because to construct a finite element space, we general need a list of
        // template elements that we store the template elements in a vector and
        // provide it to a finite element space
        std::vector<TemplateElement<double,2,2> > template_element(1);
        template_element.reinit(triangle_template_geometry,
                triangle_template_dof,
                triangle_coord_transform,
                triangle_basis_function);

        ... ...
}
           

Step 3: 构建有限元空间

Now we already have mesh data and template elements, and the following is to build a finite element space. This is really easy now - the only thing you need to do is: appoint a template for every geometry element. Let’s continue the codes in last example:

// read in a mesh at first, assume that there are only triangles in the mesh
        Mesh<2,2> mesh;
        mesh.readData("mesh_data_file");

        // construct the finite element space on the mesh and the template element
        FEMSpace<double,2> fem_space(mesh, template_element);

        // this is the number of the geometry elements
        int n_element = mesh.n_geometry(2);
        // and set the number of elements in the finite element space the same as geometry elements
        fem_space.element().resize(n_element);
        // appoint that for every element, its template element is the 0-th one in the vector template_element
        for (int i = 0;i < n_element;i ++)
                fem_space.element(i).reinit(fem_space,i,0);

        // build the element information in the finite element space
        fem_space.buildElement();
        fem_space.buildDof();
        fem_space.buildDofBoundaryMark();

        // done! now we get a finite element space but it looks that nothing happened though the library
        // in fact have done a lot of things for you, :-)
           

Step 4: 构建刚度矩阵

Now it’s time to show that the library have done something for you. We will construct a stiff matrix on the finite element space. It’s again very easy as a continuing of the last example:

#include <AFEPack/BilinearOperator.h>

        ... ...

        StiffMatrix<2,double> stiff_matrix(fem_space);
        stiff_matrix.build();

        // done! now we output the sparsity pattern of this matrix to have a look
        ofstream os("a_file_name");
        stiff_matrix.get_sparsity_pattern().print_gnuplot(os);
           

Look, it really built a sparse matrix for you! (the figure is created by gnuplot.)

Step 5: 解拉普拉斯方程

With all those preparation, we can now try to solve a Laplacian equation on the mesh.

We have calculated the stiff matrix already, and in the following we need only calculate the right hand side and apply the boundary condition. See the following code:

#include <AFEPack/AMGSolver.h>
#include <AFEPack/Operator.h>

          ... ...

          // calculate the right hand side
          Vector<double> right_hand_side;
          Operator::L2Discretize(&f, right_hand_side, 3);

	      // construct the solution which is a finite element function
          FEMFunction<double,2> solution(fem_space);

          // construct the boundary condition for boundary with boundary marker 1,
          // boundary type is Dirichlet
          // and the value is decided by function u_b
          BoundaryFunction<double,2> boundary(BoundaryConditionInfo::DIRICHLET, 1, &u_b);

          // construct a boundary condition administrator to store boundary condition
          BoundaryConditionAdmin<double,2> boundary_admin(fem_space);
          // add the boundary condition to the administrator
          boundary_admin.add(boundary);
          // apply the boundary condition on the linear system
          boundary_admin.apply(stiff_matrix, solution, right_hand_side);

          // construct a algebraic multigrid solver for the linear system
          AMGSolver solver(stiff_matrix);
          // and solver it
          solver.solve(solution, right_hand_side);
          
          // save the solution in a OpenDX data file
          solution.writeOpenDXData("a_file_name");
		
          // calculate the L2 error
          double error = Functional::L2Error(solution, FunctionFunction<double>(&u), 3);
          std::cerr << "\nL2 error = " << error << std::endl;

          return 0;
}
           

Of course, you should write down the expression of functions u b u_b ub​ and f f f as

double u_b(const double * p) {
           return sin(PI*p[0]) * sin(2*PI*p[1]);
}

//右端项
double f(const double * p) {
           return 5*PI*PI*u(p);
}
           

u b u_b ub​ and f f f 的定义请参考上面关于 Poisson Equation 定义。Let’s have a look at the solution

it’s really the solution except the mesh is too coarse.

Step 6: 自适应网格

其实到这里为止，我们已经解决了泊松方程。如果网格不够，我们可以通过自适应加密的方式，增加网格数量。

Mesh adaptation is only supported for simplex element in the package. The class Mesh<DOW,DIM> doesn’t support mesh adaptation automaticly because the mesh adaptation is based on the finite element base part. A derived class RegularMesh will support mesh adaptation. The mesh adaptation is designed to support multimesh machenism at the beginning that even you want to adapt one mesh, it will be taken as multiple mesh adaptation. Let’s see how to write such code.

#include <AFEPack/HGeometry.h>		// the head file declaring those class support mesh adaptation

    int main(void)
    {

        // construct a hierarchy geometry tree for 2 dimensional case
        HGeometryTree<2> h_geometry_tree;
        // read in the macro mesh from easymesh formatted data files
        h_geometry_tree.readEasyMesh("mesh_file_name");
        // construct a irregular mesh based on the hierarchy geometry tree
        IrregularMesh<2> irregular_mesh(h_geometry_tree);
        // semiregularize the irregular mesh to a so called semi-regular mesh
        irregular_mesh.semiregularize();	// because the macro mesh from easymesh formatted data is always 
        					// semi-regular, we in fact can omit this line. but for a general
        					// irregular mesh, this line is essential.
        // regularize the mesh to get the regular mesh from the irregular mesh
        irregular_mesh.regularize(false);	// if the argument is true, the elements of the mesh will be sorted
        					// according certain principle to try to minimize the band width of
        					// the obtained matrix.
        // construct an indicator on the regular mesh from the irregular mesh
        Indicator<2> indicator(irregular_mesh.regularMesh());
        ... ... 	// you can set the value of the indicator to control how to adapt the mesh by certain way. 
        		// see example in the following part to see how to do that
        // construct a mesh adaptor which will provide facilities to implement the mesh adaptation
        MeshAdaptor<2> mesh_adaptor(irregular_mesh);
        // set the convergence order of the mesh adaptor as 1.0
        mesh_adaptor.convergenceOrder() = 1.;
        // and the refine step as 1
        mesh_adaptor.refineStep() = 1;
        // and the indicator is the one we gave just now
        mesh_adaptor.setIndicator(indicator);
        // the tolerence of the mesh adaptor is set as 1.0e-8
        mesh_adaptor.tolerence() = 1.0e-8;
        // implement the mesh adaption
        mesh_adaptor.adapt();

        // now the mesh is adapted. let's get the regular mesh from it to see what a mesh we get
        // this mesh is really not semiregular after the adaptation now. we semiregularize it at first
        irregular_mesh.semiregularize();
        // and then construct the regular mesh from it
        irregular_mesh.regularize(false);
        // save the regular mesh to certain data files
        irregular_mesh.regularMesh().writeEasyMesh("mesh_file_name");

    };
           

This is a comparative complete example to show how to implement a mesh adaptation on a 2D mesh except there are no detail how to set the indicator. In the following figures there is a example we get to refine the meshes in two contacted circle. The first figure is the original macro mesh and the second one is the refined mesh and the third one is the local detail of the refined mesh.

从上面的三个图像中，我们可以看到，通过自适应加密，中心点网格越来越密。