fenics入门笔记

fenics名称释义：

fe：finite element的简写

cs：computational software的简写

ni：有了fe和cs后，由于最初fenics软件是在芝加哥大学(简称为phoenix）编译的，故而在其间加入ni就很自然而然了。其实就是取个谐音。

 

fenics有如下几个核心部件：

1. dolfin：

dynamic object-oriented library for finite element computation的简写。

意思就是用于有限元计算的动态面向对象库。它本质上是fenics的C++后端，配了python的接口，它实现了mesh、function space、finite element assembly、 refinement等数据结构及算法。

 

2. UFL：

The Unified Form Language的简写，这是一种用于声明有限元离散变分格式(variational form)的专用语言，它是建立在python语言上的，所以可以看作是一种推广了的python语言。具体地讲，利用这种语言我们可以以非常贴近数学记法的方式来编写计算程序。

 

dolfin和UFL都是FEniCS项目的一部分，前者可看作是用于求解问题的环境，它提供了积分域、边界条件等信息；后者则可视为前者的工具，也就是说：dolfin使用UFL的语法来指定（表出）具体的变分形式。

 

3. FFC：

The Fenics Form Compiler，也就是fenics中有限元变分型的编译器，可以把用UFL语言写的多线性型编译为更高效的C++代码。例如，把定义了一个变分问题的双线性型（左端）和线性型（右端）编译成能被用于组装对应的线性系统的代码。

 

4.FIAT：

Finite element automatic tabulator：在线段、三角形、四面体上生成任意阶的有限元

 

5. mshr：

fenics的网格生成器。

 

不同组件间的交互关系：

--------------------------------------------------------------------------

应用层          fenics apps

--------------------------------------------------------------------------

接口              dolfin

--------------------------------------------------------------------------

核心组件       ufl    syfi     ffc     ufc                                         ufl：unified form language；ffc：form compiler；

                      fiat   instant     ferari                                         ufc：unified form-assembly code

--------------------------------------------------------------------------

外部库           pets  numpy ...

--------------------------------------------------------------------------

fenics特点：在python/C++语法基础上，进一步针对有限元的数学语言进行了语法抽象，从而可以以一种非常贴近数学语言的方式来编写有限元程序。其实就是fenics花了大功夫自己写了一套有限元语言的编译器，从而可以自动生成底层代码。由于用户在非常高的抽象层次上进行编程，所以不需要考虑很多底层的处理。例如矩阵组装，数值积分，线性代数处理等等。

如果采用python接口，只需提供一个.py文件；如果采用C++接口，需提供.ufl和.cc文件。大多数fenics提供的C++库都可以通过python来调用，故而使用python和C++在效率上可能不会差太多。

先采用UFL语言来编写程序。写好之后，使用ffc编译器可将其编译成低级的按UFC（Unified Form-assembly Code）格式的C++代码。生成的这些代码可以被基于UFC的组装器，例如DOLFIN，组装成离散的方程组。

使用C++利用fenics的方式如图示：

demo.ufl——(FFC)——>demo.h

                                        demo.cpp——(C++ compiler)——>demo

在这里，form compiler FFC生成了头文件demo.h，它定义了一些适用于特定问题的类，这些类可以被用在main程序中，也可以被传递给dolfin库。它们提供了有限元组装过程的最内层循环所需要的函数功能，即：

* 基函数及其导数值的求取

* 局部自由度到全局自由度的映射关系

* 局部单元上刚度矩阵的求取

实现这些函数功能的的代码是由.ufl文件中所提供的顶层描述经form compiler解释后自动生成的。相反地，如果使用python写程序的话， 则form compiler会在需要的时候被自动调用（just in time）。

一个例子（.ufl+.cpp）：

解Poisson方程。

1. 第一步，写ufl文件。

先使用ufl语言来表述变分形式，下面这些是写在文件example_poisson.ufl中的：

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)

a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

2. 第二步，生成当前问题所用的C++对象。

对上一步编写的ufl文件进行编译(解释)，使用FFC编译得到一个头文件example_poisson.h，命令为：

ffc -l dolfin example_poisson.ufl

此时，FFC会自动生成如下类：FunctionSpace，BilinearForm和LinearForm。

3. 第三步，编写Poisson C++ code。

此时进行main程序的编写，在其中会用到dolfin下的类Function，SubDomain，UnitSquare，Constant，DirichletBC，VariationalProblem和File。此外，也会用到example_poisson.h头文件中我们自己定义的FunctionSpace，BilinearForm和LinearForm类。

#include <dolfin.h>
#include "example_poisson.h"

// For convenience
using namespace dolfin;

然后进行源项、边界法向导数和dirichlet边界条件的定义。这些都采用dolfin类的派生类来定义，我们只需要自己定义这些项的函数功能即可。

// Source term (right-hand side)
class Source : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    // The indexes are the dimension (x-dim := x[0], y-dim := x[1])
    // dx and dy are not derivatives, just differences.
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;

    // values is passed by reference, so the output to the function is
    // is not void as one may suspect.
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

// Normal derivative (Neumann boundary condition)
class dUdN : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = sin(5*x[0]);
  }
};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    // if x = 0 or x = 1 within epsilon (machine precision)
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
  }
};

然后是main函数，我们将其分为两个部分：建立问题和求解问题两部分。

int main()
{
    // Create mesh and function space

    
    UnitSquare mesh(32, 32);          // Unit Square Mesh [0x1], [0x1] with 32 rectangles in each direction
    Poisson::FunctionSpace V(mesh);   // Define the FEM space relative to the mesh

    // Direchlet BC setup
    Constant u0(0.0); // The constant value
    DirichletBoundary boundary; // WHERE the BC exists (defined in class above)
    DirichletBC bc(V, u0, boundary); // Apply the constant to the boundary on the mesh

    // Define the variational problem
    // The coefficient matrices
    example_poisson::BilinearForm a(V,V);
    example_poisson::LinearForm L(V);

    Source f;   // the source term
    dUdN g;     // the boundary flux

    // attach the source and boundary terms to the linear term (sum matrices)
    L.f = f;
    L.g = g;
=========================================================================================
    // Setup the problem
    VariationalProblem problem(a, L, bc);
    Function u(V);

    // Solve the variational problem
    problem.parameters["linear_solver"] = "iterative";
    problem.solve(u);
    
    File file("poisson.pvd");  // Save solution in VTK format
    file << u;
    plot(u);                   // Plot the solution

    return 0;
}

Note：

注意用C++和python来写同样一套程序时在语法上的区别。例如创建对象的时候，在C++中是采用Class obj;的形式，而在python中是obj = Class()；的形式。

一个例子（全python）：

打开python交互窗口，输入如下代码：

>>> from dolfin import *

>>> mesh = UnitSquareMesh(16, 16)

>>> V=FunctionSpace(mesh, 'Lagrange', 3)

此时会调用即时编译器，在交互窗口会显示：Calling FFC just-in-time (JIT) compiler, this may take some time。

这里会使用FFC生成C++代码再利用C++编译器来编译。这只执行一次，结果缓存在~/.instant

>>>f=Expression("sin(6.0*pi*x[0])*sin(2.0*pi*x[1])", degree=2)

此时会调用即时编译器

生成的C++表达式用于计算会很省时

>>>def boundary(x, on_boundary):

...     return on_boundary        #注意这里的缩进

... 

>>>bc=DirichletBC(V, 0.0, boundary)

函数boundary定义了计算域的边界，bc代表了空间V上的狄利克雷边界条件

>>>u=TrialFunction(V)

>>>v=TestFunction(V)

>>>a=inner(grad(u), grad(v))*dx

>>>L=f*v*dx

上面这段代码使用了UFL语言定义了双线性型a和线性型L

>>>u=Function(V)

>>>solve(a==L, u, bc)

调用即时编译器

调用即时编译器

求解线性变分问题

最后我们在空间V上创建了有限元函数u（TRialFunction和TESTFunction仅仅被看作是多线性型的参数而已，并不是在内存中有实际值的函数）。然后要求DOLFIN组装出线性系统并使用后端的线性代数工具进行求解。对于前者，也是由FFC来生成C++代码，再使用C++编译器来编译。这也是与网格无关的，因此即便网格被加密，这也不会再次执行。

>>>plot(u, interactive=True)

上述过程也可以写到一个.py脚本中，再在shell中执行。

 

FFC组件的介绍

我对于fenics中的FFC组件比较感兴趣，它能将专门的ufl语言翻译并自动生成C++代码，用到一些较高级的编程技术。所以我在这里专门做一些介绍。

在上面介绍到，使用命令ffc -l dolfin example_poisson.ufl可自动生成一个C++头文件，它里面将包含了在当前问题中可能用到的有限元类。这一看上去挺高级的技术是怎么实现的呢？在fenics book中对此进行了介绍。FFC把编译过程分为了如下几个阶段，每一阶段的输出结果会作为下一阶段的输入。

compile stage 0：Parsing：在这个阶段，用户指定的form被解释并保存为UFL abstract syntax tree(AST)。实际的parsing过程被交给python，从.ufl文件到UFL form对象的转化则通过UFL的操作符重载来实现。

input：python code或.ufl文件，output：UFL form

compile stage 1：Analysis：在这个阶段对UFL form进行预处理，抽取出form元数据(FormData)，例如：用于定义form的是哪种有限元，系数的个数，单元类型(间隔，三角形或四面体)。

input：UFL form，output：预处理的UFL form和form元数据(metadata)

compile stage 2：Code representation：这个阶段检查输入并生成所有代码生成所需的数据。包括有限元基函数的生成，自由度映射的数据，可能的积分预计算。绝大多数复杂的翻译过程都发生在这一阶段。中间数据被保持为一个字典dictionary，把UFC函数名映射为代码生成所需的数据。最简单的情况，例如ufc::form::rank，这个数据可能是个简单数字，如2。更复杂的情况，如ufc::cell_tensor::tabulate_tensor，这个数据可能是个复杂的数据结构。

input：预处理的UFL form和form元数据，output：中间代表(intermediate representation)

compile stage 3：Optimization：这一阶段检查intermediate representation并实施优化。保存在intermediate representation字典中的数据被替换为优化后的数据，它编码了问题所需函数的优化版本。

input：intermediate representation，output：optimized intermediate representation

compile stage 4：Code generation：这个阶段检查优化后的IR并生成实际对应了每个UFC函数的C++代码。这个代码被保存为一个字典，它把UFC函数名映射为包含了C++代码的string。例如，ufc::form::rank生成的数据可能是字符串“return 2;”。把2/3/4三个阶段分开是有重要意义的，它使得阶段2和3可以将重心放在有关有限元和变分形式的算法方面，而阶段4则只集中于生成C++代码。

input：优化后的intermediate representation(OIR)，output：C++代码

compile stage 5：Code formatting：这个阶段检查生成的C++代码，并将它按照UFC的代码格式进行格式化，生成一个或多个.h/.cpp文件。这里是实际写出C++代码文件的地方。这一阶段依赖于UFC代码的模板。

input：C++代码，output：C++代码文件。

 

 

卸载fenics：

1. 仅仅卸载fenics

sudo apt-get remove fenics

2. 卸载fenics及其依赖文件

sudo apt-get remove --auto-remove fenics

3. 清除本地的config/data

sudo apt-get purge fenics

或：sudo apt-get purge --auto-remove fenics