基于FreeFEM++的有限元编程--1
1 FreeFEM++简介
FreeFEM是开源的有限元模拟系统,有法国利翁斯实验室、埃尔和玛丽居里大学共同开发,在世界范围内广泛使用[i]。VS Code插件商店中有专门针对FreeFEM++的插件方便代码编辑,其lanuch.json的配置如下:
{
// Use IntelliSense to learn about possible attributes.
// Hover to view descriptions of existing attributes.
// For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
"version": "0.2.0",
"configurations": [
{
"name": "(Windows) Launch",
"type": "cppvsdbg",
"request": "launch",
"program": "D:/Program Files (x86)/FreeFem++/FreeFem++.exe",
"args": ["-f", "${file}"],
"stopAtEntry": false,
"cwd": "D:/Program Files (x86)/FreeFem++",
"environment": [],
"externalConsole": false
}
]
}
使用Free FEM++进行有限元模拟包括以下几个步骤。
建模:从DFE的强形式到弱形式,必须知道如何用FreeFEM表示变分形式。
编码:在文本编辑器中用FreeFEM语言编写程序。
运行:通过命令行执行程序代码,比如FreeFem++ mycode.edp,其中mycode.edp为文本文件。
可视化:在mycode.edp文件中直接用plot关键词展示函数,在程序运行过程中显示图像,利用图形参数wait=1可以在每个显示图形的计算步暂停一下。在plot图形显示界面,切换“f”键显示图像的颜色充填或等值线,切换“v”控制是否显示色标。
调试:全局变量debug可以帮助调整程序,使wait=true或wait=false
比如下列程序代码,把debug改为true可以连续的显示图像
//除了内置的plot关键词显示结果,也可以保存问vtk格式,采用paraview等软件查看结果
load "iovtk" //加载vtk模块,
bool debug = true;
border a(t=0, 2.*pi){x=cos(t); y=sin(t); label=1;};
border b(t=0, 2.*pi){x=0.3+0.3*cos(t);y=0.3*sin(t);label=2;};
plot(a(50)+b(-30), wait=debug);//plot the borders to see the intersection
//if debug == true, press Enter to continue
mesh Th = buildmesh(a(50)+b(-30));
plot(Th, wait=debug);//plot Th then press Enter
fespace Vh(Th, P2);
Vh f=sin(pi*x)*cos(pi*y);
Vh g=sin(pi*x+cos(pi*y));
plot(f, wait=debug);//plot the function f
plot(g, wait=debug); //plot the function g
string vtkFileName="./test1.vtk";
// 需要添加场名称后处理时方可选择变量
savevtk(vtkFileName,Th,f,g,dataname ="Th f g");
程序运行结果如下:
当border b(t=0, 2.*pi){x=0.8+0.3*cos(t);y=0.3*sin(t);label=2;};时边界a和边界b相交,无法根据两个边界剖分网格,后续执行失败。
border b(t=0, 2.*pi){x=0.3+0.3*cos(t);y=0.3*sin(t);label=2;};时边界a与边界b不相交,可以执行后续网格剖分和计算。
显示Th网格划分
函数f的计算结果
下图是用paraview显示的结果
函数g的计算结果
下图是用paraview的显示结果
2 基于实例的学习
2.1 从这里开始
对于给定的函数f(x,y),找到一个函数u(x,y)满足如下2个条件:
-Δu(x,y)=f(x,y) 对所有(x,y)在空间域Ω上成立
u(x,y)=0 对所有(x,y)在边界域上
内部函数u表示为
load "iovtk" //载入模块
border C(t=0,2*pi) {x=cos(t);y=sin(t);};
int nmsh=60;
mesh Th = buildmesh(C(nmsh));
plot(Th, wait=true);
fespace Vh(Th,P1);
Vh u,v,fh;
func f= x*y;
fh =f;
macro Grad(u) [dx(u),dy(u)] //
solve Poisson(u,v)=
int2d(Th)(Grad(u)'*Grad(v))
+int2d(Th)( f*v)
+on(C,u=0);
plot(u,fill=1);
string outDir = "./";
string vtkFileName = outDir + "Poisson_mesh_" + nmsh + ".vtk";
// 需要添加场名称后处理时方可选择变量
savevtk(vtkFileName,Th,u,f,dataname ="u source");
[i] https://freefem.org/