基于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// Meshreal 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);   //显示网格// Fespacefespace femp1(disk, P1);     //在disk网格上创建由P1组成的有限元空间femp1 u, v;// Problem定义问题problem laplace(u, v)= int2d(disk)(  // bilinear formdx(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);// Displaycout << "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;

计算结果如下：

参考文献

基于MPI_FreeFem_的有限元并行计算_尚月强

基于Freefem_的有限元数值计算在计算物理中的应用_钟振

基于FreeFEM++的一阶有限元解法

基于FreeFEM++的有限元编程--2相关推荐

基于VTK的有限元计算结果可视化尹小刚，奚金荣
基于VTK的有限元计算结果可视化尹小刚,奚金荣河海大学工程力学系,南京 (210098) 摘要:有限元计算结果的可视化是分析有限元计算结果的重要手段.随着可视化的广泛应用,人们对于三维真实感图形的 ...
自动售货机编程_Rosmaro中基于视觉自动机的编程简介
自动售货机编程 by Łukasz Makuch 通过ŁukaszMakuch Rosmaro中基于视觉自动机的编程简介 (An introduction to visual automata-bas ...
arcgis下的python编程-基于ArcGIS的Python编程秘笈(第2版)
基于ArcGIS的Python编程秘笈(第2版) 译者简介审阅人简介 1 面向ArcGIS的Python语言基础 1.1 使用IDLE进行Python脚本开发 1.1.1 Python Shell窗 ...
基于MATLAB的面向对象编程（1）——类，属性，方法
这里写目录标题 1 入门 1.1 如何把面条抽象成Class 1.2 文件类 2 基于MATLAB的面向对象编程入门 2.1 如何定义一个类(Class) 2.2 创建一个对象 2.3 类的属性 2. ...
基于对象和面向对象编程范式辨析和主流编程语言中的应用
基于对象和面向对象编程范式辨析和主流编程语言中的应用前言本文的目的是想告诉大家,为什么C++的模板这么强大.为什么Ruby的Duck Typing(像鸭子那样编程)这么强大! 基于对象和面向对象编 ...
pyaudio：基于pyaudio利用Python编程从电脑端录制音频保存到指定文件夹+将录音上传服务器+录音进行识别并转为文本保存
pyaudio:基于pyaudio利用Python编程从电脑端录制音频保存到指定文件夹+将录音上传服务器+录音进行识别并转为文本保存目录输出结果代码实现输出结果代码实现 # -*- codi ...
基于Native Client的编程框架
上面是一个浏览器的架构图,JS引擎+渲染引擎+外壳Shell+内核 4部分构成了浏览器的主体,传统的插件(上图左部)通过NPAPI与浏览器通信,Native Client(上图右部)通过PPAPI与浏 ...
R开发(part12)--基于RC的面向对象编程
学习笔记,仅供参考,有错必纠参考自:<R的极客理想>-- 张丹文章目录 R开发基于RC的面向对象编程创建RC类和对象对象赋值定义对象的方法 RC对象内置方法和内置属性 RC类的 ...
R开发(part11)--基于S4的面向对象编程
学习笔记,仅供参考,有错必纠参考自:<R的极客理想>-- 张丹文章目录 R开发基于S4的面向对象编程创建S4对象访问S4对象的属性 S4的泛型函数查看S4对象的函数 R开发基 ...
R开发(part10)--基于S3的面向对象编程
学习笔记,仅供参考,有错必纠参考自:<R的极客理想>-- 张丹文章目录 R开发面向对象面向对象的R语言实现基于S3的面向对象编程泛型函数和方法调用查看S3对象的函数 S3对象 ...

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