2

I am using Freefem++ to solve the poisson equation

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

It works well when I have an analytical expression for f, but now I have an f numerically defined (i.e. a set of data defined on a mesh) and I am wondering if I can still use Freefem++.

I.e. typical code (for a 2D problem in this case), looks like the following

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

I am wondering if I can define f numerically, rather than analytically.

4

2 回答 2

3

网格和空间定义

我们定义了一个 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

条件和问题说明

我们不会使用solve，但我们将处理矩阵（使用FreeFem 求解的更复杂的方法）。

首先，我们为我们的问题定义 CL（Dirichlet 问题）。

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 如下，但注意宏以//NOT 结尾;。

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

所以泊松双线性形式变为：

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

在我们提取刚度矩阵之后

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 上专门给我一些技巧。我测试了它，它工作得很好。祝你好运。

于 2014-03-25T16:59:24.353 回答
0

afaf的方法适用于函数f是独立函数的情况。对于类似的术语int2d(Sh)(f*u*v)，需要另一种解决方案。我建议（实际上我在 Hecht 的手册中某处将其标记为红色）一种涵盖这两种情况的方法。但是，它仅适用于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;
于 2015-03-28T17:04:04.327 回答