作为偏微分方程求解器的有限元方法实现

FEM as a tool of PDEs solver and its open source implementation

(原文首发于http://hillyuan.blogspot.com )

1. Finite element function space

The construction of the finite element space Vh begins with subdividing the domain Ω into a set of non-overlapping or overlapping, such as those implemented in GFEM, cells T1,T2,...Tm . The domain Ω can now be approximated with a mesh domain,

Th(Ω)=∪Ti

The cells are typically simple polygonal shapes like intervals, triangles, quadrilaterals, tetrahedra or hexahedra but also could be described by B-splines as in isogeometric analysis.

Once a domain Ω has been partitioned into cells, one may define a local function space V on each cell Ti and use these local function spaces to build the global function space Vh . A cell Ti together with a local function space V and a set of rules for describing the functions in V is called a finite element. This definition reads as follows: a finite element is a triple (T,V,L) , where
• the domain T is a bounded, closed subset of Rd(for d=1,2,3...) with nonempty interior and piecewise smooth boundary;
• the space V=V(T) is a finite dimensional function space on T of dimension n;
• the set of degrees of freedom (nodes) L=1,2,...,n is a basis for the dual space V^ ; that is, the space of bounded linear functionals on V .

If finite element spaces for grad, div, and curl are implemented, most types of PDEs are solvable.

• H1Space - the most common FE space of continuous, piecewise-polynomial function:
H1(Ω)={v∈L2(Ω);∇u∈[L2(Ω)]2}

• HCurlSpace - FE space of vector-valued functions discontinuous along mesh edges, with continunous tangential component on the edges

H(curl,Ω)={E∈[L2(Ω)]2;∇×E∈L2(Ω)}

• HDivSpace - FE space of vector-valued functions discontinuous along mesh edges, with continunous normal component on the edges

H(div,Ω)={v∈[L2(Ω)]2;∇⋅v∈L2(Ω)}

2. Weak form of PDE

We consider a general linear variational problem written in the following canonical form: find u ∈ V such that

a(u,v)=L(v):∀v∈V^

where the test space V^ is defined by

V^={v∈H1(Ω):v=0 on ΓΩ}

and the trial space V contains members of V^ shifted by the Dirichlet condition:

V={v∈H1(Ω):v=u0 on ΓΩ}

We thus express the variational problem in terms of a bilinear form a and a linear form (functional) L :

a:V×V^→R

L:V^→R

The basic idea of FEM consists of approximating the big Hilbert space H01 , e.g., by some appropriate finite-dimensional space Hh , where h is a positive parameter, satisfying the following conditions:

• Hh⊆H01

• Solve a(uh,vh)=L(vh) obtain a solution vh and

• Hh→H01 as h→0

Example: Poisson problem

Find u such that

−Δu=f,x∈Ω⊂Rd

u=0 on ∂Ω

After integrating by part, we obtains its weak form

Find u∈V=H01(Ω):{v=0 on ∂Ω}

such that∫Ω∇u⋅∇vdΩ=∫ΩfvdΩ,∀v∈V

where a(u,v)=∫Ω∇u⋅∇vdΩ, L(v)=∫ΩfvdΩ

3. Implementation on some mathmatically oriented open source FEM software

The general FEM system could be thought of as having form a(uh,vh)=L(vh) and solved by following steps

• Build a mesh Th. It could be also read from existing mesh files.
• Define the finite element space Vh

• Define the variational problem P: a(uh,vh)=L(vh)

• Solve the problem P

• Postprocess

Lets see how it is really implemented in some open source FEM softwares

3.1 freeFEM

We need prepare following input file:

mesh Th=square(10,10) // generate a mesh Th

fespace Vh(Th,P1) // define FE space Vh on Th by element type P1

Vh u,v // declear funtion in Vh

func f = -1 // constant function

func g =0 // constant function

// a(u,v)=∫Ω∇u⋅∇vdΩ

problem Poisson(u,v) = int2d(Th)(dx(uh)*dx(vh) dy(uh)*dy(vh))

-int2d(f*v) // L(v)=∫ΩfvdΩ $

on(C, g) // u=g on boundary C

Poisson // solve it

plot(u) // visualize solution

3.2 FEniCS

Following python script is needed:

from fenics import *

Th = UnitSquareMesh(10, 10) // generate a mesh Th
Vh = FunctionSpace(mesh, ’P’, 1) // define FE space Vh on Th by element type P1

u = TrialFunction(Vh) // declear funtion in Vh
v = TestFunction(Vh) // declear funtion in Vh

f = Constant(-1.0) // constant function

g = Constant(0.0) // constant

a = dot(grad(u), grad(v))*dx // a(u,v)=∫Ω∇u⋅∇vdΩ
L = f*v*dx // L(v)=∫ΩfvdΩ

def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, g, boundary)

solve(a == L, u, bc) // solve it

plot(u) // visualize solution

A Julia wrapper of FEniCS given here.

3.3 Feel

Following c program construct a solver for Poisson equation

auto mesh = ... // build a mesh

auto Vh = Pch<a>(mesh) // define FE space

auto a = form2(_trial=Vh, _test=Vh) // a(u,v)=∫Ω∇u⋅∇vdΩ

a = integrate( _range=elements(mesh), _expr=gradt(u)*trans(grad(u));

auto l= form1(_test=Vh) //declare a one form

l = integrate( _range=elements(mesh), _expr=f*id(u)); // L(v)=∫ΩfvdΩ
a.solve(_solution=u, _rhs=l ) // solve it

3.4 GetDP

The input file to solve Poisson equation look like:

FunctionSpace {
{ Name H1; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef vn; Function BF_Node; Support D; Entity NodesOf[All]; }
}
}
}
Formulation{
{ Name Poisson; Type FemEquation;
Quantity {
{ Name v; Type Local; NameOfSpace H1; }
}
Equation {
Integral { [ a[] * Dof{d v}, {d v} ] ; In D; Jacobian V; Integration I; }
Integral { [ f[], {v} ] ; In D; Jacobian V; Integration I; }
}
}
}

3.5 Sundance

As a part of big Trilinos, it seems that it is not on develop right now. Compile error may happen due to its old-styled code. My fork with modification, which is confirmed to work well, is given here.

// read a mesh

Mesh mesh = meshSrc.getMesh();

// define FE space
BasisFamily basis = new Lagrange(1);
Expr u = new UnknownFunction(basis, "u");
Expr v = new TestFunction(basis, "v");
QuadratureFamily quad1 = new GaussianQuadrature(1);

/** Write the weak form */
Expr eqn
= Integral(interior, (grad*u)*(grad*v), quad1)
Integral(east, v, quad1); /* Set up linear problem */
LinearProblem prob(mesh, eqn, bc, v, u, vecType);

/* --- solve the problem --- */

/* Create the solver as specified by parameters in an XML file */

LinearSolver<double> solver
= LinearSolverBuilder::createSolver(solverFile);

/* Solve! The solution is returned as an Expr containing a
* DiscreteFunction */
Expr soln = prob.solve(solver);

3.6 MFEM

You cannot define above two form and one form yourself, MFEM provides prepared classes of one form and two form integrators instead. Something like, e.g.

Therefore, it is not that strictly mathematically oriented than above ones but would expect more computational efficiency because software intepretation of one form, two form definition are not needed anymore.

Its implementation on Poisson problem looks like:

// read a mesh

Mesh *mesh = new Mesh(mesh_file, 1, 1);

// define FE space

FiniteElementCollection *fec;

FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);

// define one form

LinearForm *b = new LinearForm(fespace);

ConstantCoefficient one(1.0);

b->AddDomainIntegrator(new DomainLFIntegrator(one));

b->Assemble();

// define two form

BilinearForm *a = new BilinearForm(fespace);

a->AddDomainIntegrator(new MixedGradGradIntegrator(one));

a->Assemble();

Reference:

1) A Logg, K-A Mardel. GN Wells: Automated Solution of Differential Equations by the Finite Element Method, Springer-Verlag. 2012

2) Eric Sonnendru¨cker, Ahmed Ratnani : Advanced Finite Element Methods