FEniCS是一个由University of Chicago和Chalmers University of Technology主导开发的开源C++/Python的FEM代码。它可在Win/Mac/Unix下运行，但官方Win平台编译版只跟进到1.0，想玩最新的还是在Lin下吧。

它的结构为：

对于初级用户来说，用Python上手学习怎么使用FEniCS还是很容易的，也不必它在意你的PDEs或ODEs怎么去求解的，用户只需要懂得一点点有限元变分原理，对自己的方程要求解什么有概念即可。

以Poisson问题为例：

本代码定义了有限元网格、网格和单元类型离散函数空间，边界条件，以及线性弱形式方程，从而计算得到未知u的试函数。

变分方法：PDE中乘以一个试函数v，在计算域上积分，如遇2阶项便采用分部积分。

试函数v 即为 test function

而未知函数 u 则为 trial function

在FEniCS中，需要分别给这两个函数分配离散空间

有限元法作为变分法和分片插值法结合的产物，它在单元内使用低次多项式分片插值，然后在把它们组合起来，形成全域内的函数，从而逼近真解。

而变分问题最终可转化为     a(u,v) = L(v) 的形式，它由一个双线性项a(u,v)和一个线性项L(v)组成。

对于Poisson问题，则：

代码：

from dolfin import *

# Python语言导入库 dolfin

# Create mesh and define function space
mesh = UnitSquare(6, 4)

# 划分网格，将正方形先划成6×4的方格，然后在分两半，因此共48个三角形单元
#mesh = UnitCube(6, 4, 5)

# 正方体
V = FunctionSpace(mesh, 'Lagrange', 1)
# 单元类型定义，lagrange单元，自由度1个，既线性lagrange单元

#  Finite element        Short name         Sobolev space       Conforming
#  (Quintic) Argyris        ARG                  H2                 Yes
#  Arnold–Winther            AW               H(div; S)             Yes
#  Brezzi–Douglas–Marini    BDM                H(div)               Yes
#  Crouzeix–Raviart          CR                  H1                 No
#  Discontinuous Lagrange    DG                  L2                 Yes
#  (Cubic) Hermite          HER                  H2                 No
#  Lagrange                  CG                  H1                 Yes
#  Mardal–Tai–Winther       MTW                H1/H(div)            No/Yes
#  (Quadratic) Morley       MOR                  H2                 No
#  Nédélec first kind       NED1               H(curl)              Yes
#  Nédélec second kind      NED2               H(curl)              Yes
#  Raviart–Thomas            RT                 H(div)              Yes

# Define boundary conditions
u0 = expression_r('1 + x[0]*x[0] + 2*x[1]*x[1]')

# 坐标轴 x y z 分别表示为 x[0] x[1] x[2]，定义u0的表达式

def u0_boundary(x, on_boundary):
   return on_boundary

# 是否在边界上的布尔运算 也可以定义为：

#  def u0_boundary(x):
#      return x[0] == 0 or x[1] == 0 or x[0] == 1 or x[1] == 1

# 但是 == 这种强制的算符，并不推荐，推荐写法：

#  def u0_boundary(x):
#      tol = 1E-15
#      return  abs(x[0]) < tol or
#              abs(x[1]) < tol or
#              abs(x[0] - 1) < tol or
#              abs(x[1] - 1) < tol

bc = DirichletBC(V, u0, u0_boundary)

# 在函数空间V上，寻找 u0_boundary上的点，把u0的函数值赋给它们，边界类型为Dirichlet

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

# FEniCS区分test 和trial spaces，它们的区别只在是否引入边界条件

f = Constant(-6.0)

# 如果是常数推荐用Constant()
a = inner(nabla_grad(u), nabla_grad(v))*dx

# inner()表示积分符号， nabla_grad()梯度算符，它和grad()的区别在于，是否含有tensor算法
L = f*v*dx

# 组装得到 双线性项 a 和 线性项 L

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# 求解得到u

# 也可写为：

#  equation = inner(nabla_grad(u), nabla_grad(v))*dx == f*v*dx
#  u = Function(V)
#  solve(equation, u, bc)

# 求解器控制：

#  solver_parameters={"linear_solver": "cg","preconditioner": "ilu"})

# 默认使用 Krylov求解器，此外还有"PETSc", "uBLAS", "Epetra", or "MTL4".

#  info(parameters, True)

#  prm = parameters["krylov_solver"] # short form
#  prm["absolute_tolerance"] = 1E-10
#  prm["relative_tolerance"] = 1E-6
#  prm["maximum_iterations"] = 1000

# 求解也可以这样写

#  u = Function(V)
#  problem = LinearVariationalProblem(a, L, u, bc)
#  solver = LinearVariationalSolver(problem)
#  solver.solve()

# Plot solution and mesh
plot(u)
plot(mesh)

# Dump solution to file in VTK format
file = File('poisson.pvd')
file << u

# Hold plot
interactive()

# 获取节点值和坐标

#  u_nodal_values = u.vector()

#  u_array = u_nodal_values.array()

#  coor = mesh.coordinates()
#  if mesh.num_vertices() == len(u_array):
#  for i in range(mesh.num_vertices()):
#  print ’u(%8g,%8g) = %g’ % (coor[i][0], coor[i][1], u_array[i])

# 获取节点插值

# u_e = interpolate(u0, V)

# 获取某点指

center = (0.5, 0.5)

u(center)

mxio

2013.6.18