基于FreeFEM++的有限元编程--2

基于FreeFem++的有限元模拟

1、引言

有限元方法是２０世纪５０年代伴随电子计算机的诞生，在计算数学和计算工程领域里诞生的一种高效而灵活的计算方法，它将古典变分法与分片多项式插值相结合，易于处理复杂的边值问题，具有有限差分法无可比拟的优越性，广泛应用于求解热传导、电磁场、流体力学等相关问题，已成为当今工程分析和计算中不可缺少的最重要的工具之一。有限元方法的“化整为零、积零为整”的基本思想与并行处理技术的基本原则“分而治之”基本一致，因而具有高度的内在并行性。

有限元基于分片插值函数，相比于有限差分方法，边界更接近于真实情况，误差相对较小。目前主流有限元软件众多，如 Ansys、MSC、Abaqus、COMSOL等,这些商业闭源软件在业界具有良好的口碑靠性较好。但这些软件由于相关模块被封装，难以理解有限元的计算原理,较难对相关模块作进一步的科研开发，限制了教学和科研的使用。除上述软件外,也有不少开源的有限元软件,如OpenFoam, Coder Aster、OOFEM、Emer、 ParaFEM和 Freefem这些软件不但对公众开放源代码,而且可靠性较高,相关模块的说明文档也较丰富,非常适合于科研。其中, Freefem++语句非常简洁,只需要将原方程转换成弱解形式,使用Freefem++很容易编写有限元程序,是最接近有限元的“计算语言”。

2、FreeFem++简介

FreeFem++是一款免费的、开放源码的二维、三维偏微分方程有限元计算软件，由巴黎第六大学研究人员开发，它集成网格生成器、线性方程组的求解器、后处理及计算结果可视化于一体，能快速而高效地实现复杂区域问题的有限元数值计算。它采用Delaunay算法生成离散偏微分方程所需的网格，具有网格自适应和移动网格生成的功能；提供了协调与非协调多种有限元离散方法，以及诸如LU分解、CG迭代法、GMRES迭代法等多种线性方程组的求解器，还提供了特征值问题的求解器等；采用稀疏矩阵存储格式，内存需求少，计算速度快；其编程的语义环境类似于C++ ，轻松易学，使用十分方便；同时还提供了多种计算结果的输出格式，可供可视化软件如Tecplot、Matlab等调用，受到了广大研究者与工程计算人员的青睐。

3、FreeFem++数据类型

在本质上，FreeFEM++是一个编译器:它的语言是类型化，多样化的，有可重用性。每个变量必须在说明语句中声明一个确定的类型;每个说明语句用标点符号“;”与下一个分离。它的语言与C++很相似，允许基本的类型定义，如:整数型、实数型、字符串、数组、二维有限元网格(mesh)、二维有限元空间(fespace)、解析函数func)线性和双线性运算符、稀疏矩阵、向量等。

在FreeFEM++里，类型声明是强制性的;这个特征是他的一个优点，因为一个拥有很多暗含类型的语言里很容易出现故障。变量的名字仅仅是一个字母数字串，下划线“_”是不允许使用的。

主要类型有：bool, int, string, real, complex, ofstream, ifstream, real[int], real[string], string[string], func, mesh, square, fespace, problem, solve, varf, matrix

bool 是布尔类型，用于逻辑运算和流程控制；

int 定义整型变量；

string 定义字符变量，一般用双引号括起来；

real 定义实型变量；

complex 定义复型变量；

ofstream 定义一个输出文档

ifstream 定义一个输入文档；

real[int] 定义一个变量，可以利用整型索引存储多重实型整数。比如：

real[string]定义一个变量,可以利用字符串索引储存多重实型数据。

string[string]定义一个变量,可以利用字符串索引储存字符串。

func 定义一个没有参数的函数,自变量为x,y,例如: func f=cos(x)+sin(y)

mesh 进行三角剖分,值得注意的是,如果是简单的矩形平面,可以应用命令" square"将单元划分成规则的三角形,例如mesh Th= square(10, 10);可以在单位正方形上产生10×10的网格;

fespace定义一个新的有限元空间;

problem定义一个偏微分问题的弱形式;

solve定义一个问题并且解出它;

varf 定义一个完整的变分形式;

matrix定义一个稀疏矩阵;

4、FreeFem++里的全局变量

x, y, z, Label, P, N, lenEdge, hTriangle, nuTriangle, nuEdge, nTonEdge, area, Volume, cout, cin, endl, true, false, pi

x 现有点的x坐标(实数值)

y 现有点的y坐标(实数值)

z 现有点的z坐标(实数值)

Label 指明包含该点的边界编号,如果该点在边界上,否则显示为0(整数)

P 给出现有点信息,通过输入px,p.y可以得到该点的坐标。

N 给出外部单元在该点的法向量,条件是该点必须在通过 border命令定义的曲线上。N.x和N.y分别输出在x和y方向的法向量分量。

lenEdge给出现有边的长度

hTriangle给出现有三角形的尺寸。

nuTriangle给出现有三角形的索引编号(整数值)。

nuEdge给出三角形里现有边界的索引编号(整数值)

nTonEdge给出现有边界相邻三角形的编号(整数值)。

area 给出现有三角形的面积值(实数值)

Volume给出现有四面体的体积值(实数值)

cout 标准的输出命令

cin 标准的输入命令

endl 输出或输入命令行的结束语

true 指布尔值为真

false 指布尔值为假

pi 指π的近似数值(实数值)

plot 可视化

5、实例应用

下面的示例展示定义在圆形区域上的拉普拉斯方程的求解，体现了Freefem++定义求解边界、创建网格，创建有限元空间，求解，显示结果的基本过程。而且采用了一般的网格划分方法，和自适应网格划分，展示了不同划分方法的效果、求解结果，及误差结果。

// This test shows some powerful features of FreeFEM on a

// simple example: -\Delta(u) = 1 in the unit cercle with u=0 on the

// border of the unit cercle. This problem has an analytical solution

// u = (1-x^2-y^2)/4



// Mesh

real pi = 4*atan(1);

border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}  //定义边界



mesh disk = buildmesh(a(50));    //定义网格

plot(disk, wait=true);   //显示网格



// Fespace

fespace femp1(disk, P1);     //在disk网格上创建由P1组成的有限元空间

femp1 u, v;



// Problem定义问题

problem laplace(u, v)

    = int2d(disk)(  // bilinear form

          dx(u)*dx(v)

        + dy(u)*dy(v)

    )

    + int2d(disk)(  // linear form

        - 1*v

    )

    + on(1, u=1)    // boundary condition

    ;



// Solve求解

laplace;



// Error定义误差

femp1 err = u - (1-x^2-y^2)/4;



// Plot显示求解结果和误差

plot(u, value=true, wait=true);

plot(err, value=true, wait=true);



// Display (on terminal)

cout << "error L2 = " << sqrt(int2d(disk)( (u-(1-x^2-y^2)/4)^2 )) << endl;

cout << "error H10 = " << sqrt(int2d(disk)((dx(u)+x/2)^2)

                             + int2d(disk)((dy(u)+y/2)^2) )<< endl;



// Mesh adaptation 创建自适应网格

disk = adaptmesh(disk, u, err=0.01);

plot(disk, wait=1);  



// Solve求解

laplace;



// Error计算误差

err =u-(1-x^2-y^2)/4;



// Plot 显示求解结果和误差

plot(u, value=true, wait=true);

plot(err, value=true, wait=true);



// Display

cout << "error L2 = " << sqrt(int2d(disk)( (u-(1-x^2-y^2)/4) ^2 )) << endl;

cout << "error H10 = " << sqrt(int2d(disk)((dx(u)+x/2)^2)

                             + int2d(disk)((dy(u)+y/2)^2)) << endl;

计算结果如下：