代码之家 › 专栏 › 技术社区 › DJames

Freefem++:用数值函数求解泊松方程

freefem++ finite-element-analysis numerical-methods

2

DJames · 技术社区 · 10 年前

我使用Freefem++来求解泊松方程

梯度^2 u(x,y,z)=-f(x,y,z)

当我有一个f的解析表达式时,它工作得很好,但现在我有了一个数值定义的f(即在网格上定义的一组数据),我想知道我是否仍然可以使用Freefem++。

例如,典型代码(在本例中为2D问题)如下

mesh Sh= square(10,10); // mesh generation of a square
fespace Vh(Sh,P1); // space of P1 Finite Elements
Vh u,v; // u and v belongs to Vh

func f=cos(x)*y; // analytical function

problem Poisson(u,v)= // Definition of the problem
    int2d(Sh)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear form
    -int2d(Sh)(f*v) // linear form
    +on(1,2,3,4,u=0); // Dirichlet Conditions
Poisson; // Solve Poisson Equation
plot(u); // Plot the result

我想知道我是否可以用数字来定义f,而不是用分析来定义f。

2 回复 | 直到 10 年前

1

4

JasonMArcher TWE 10 年前

网格&空间定义

我们定义了一个Nx=10网格和Ny=10的正方形单元,这在x轴上提供了11个节点,在y轴上也提供了同样的节点。

int Nx=10,Ny=10;
int Lx=1,Ly=1;
mesh Sh= square(Nx,Ny,[Lx*x,Ly*y]); //this is the same as square(10,10)
fespace Vh(Sh,P1); // a space of P1 Finite Elements to use for u definition

条件和问题陈述

我们不打算使用求解,但我们将处理矩阵(使用FreeFem求解的更复杂的方法)。

首先,我们为我们的问题(Dirichlet问题)定义CL。

varf CL(u,psi)=on(1,2,3,4,u=0); //you can eliminate border according to your problem state
    Vh u=0;u[]=CL(0,Vh);
    matrix GD=CL(Vh,Vh);

然后我们定义问题。而不是写作 dx(u)*dx(v)+dy(u)*dy(v) 我建议使用宏,所以我们将div定义如下,但要注意宏的完成方式 // 不是 ; .

 macro div(u) (dx(u[0])+dy(u[1])) //

因此泊松双线性形式变为:

varf Poisson(u,v)= int2d(Sh)(div(u)*div(v));

在我们提取Stifness矩阵之后

matrix K=Poisson(Vh,Vh);
matrix KD=K+GD; //we add CL defined above

我们继续进行求解,UMFPACK是FreeFem中的一个求解器,对此没有太多关注。

set(KD,solver=UMFPACK);

这里是你需要的。您需要在某些特定节点上定义函数f的值。我要给你一个秘密,泊松线性形式。

real[int] b=Poisson(0,Vh);

在要执行的任何节点上定义函数f的值。

b[100]+=20; //for example at node 100 we want that f equals to 20
b[50]+=50; //and at node 50 , f equals to 50

我们解决我们的系统。

u[]=KD^-1*b;

最后我们得到了剧情。

plot(u,wait=1);

我希望这会对你有所帮助,感谢我的实习主管Olivier,他总是给我一些关于FreeFem的技巧。我测试过了,效果很好。祝你好运

2

0

freude 9 年前

该方法由 美国空军 在函数 f 是独立的。对于以下术语 int2d(Sh)(f*u*v) ,需要另一种解决方案。我建议(实际上我在赫克特的手册中的某个地方写过)一种涵盖这两种情况的方法。然而,它只适用于 P1 有限元,其自由度与网格节点一致。

fespace Vh(Th,P1);
Vh f;

real[int] pot(Vh.ndof);

for(int i=0;i<Vh.ndof;i++){
    pot[i]=something;   //assign values or read them from a file 
}

f[]=pot;