基于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/