记录 fenics 的使用过程

偶尔要用到 fenics 来计算, 但是各种不熟悉, 每次网上找, 资料很多很乱, 而且 fenics 的版本更新后, 很多函数的用法都报错, 所以这里还是记录一下自己的使用过程.

同时需要参考 fenics 中的基本类型介绍: Docs » User manual » Form language).

1. 定义网格和空间后, 得到空间的*度编号、函数在*度处的值、

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)

# 网格节点处的坐标值
mesh_coor = mesh.coordinates()

# 建立标量和向量空间
space_S1 = FunctionSpace(mesh, 'CG', 1)  # 1 表示 Lagrange 1 次元
space_S3 = VectorFunctionSpace(mesh, 'CG', 1, dim=3)

# 得到*度编号
dof_f1 = np.array(df.dof_to_vertex_map(space_S1), dtype=int)
dof_f3 = np.array(df.dof_to_vertex_map(space_S3), dtype=int)

# 建立 space_S1 上的插值 f1, 并得到函数在*度处的值
f1 = interpolate(df.Expression("1.0+0*x[1]", degree=1), space_S1)
val_dof_f1 = f1.vector().get_local()  

# 计算在网格节点处的函数值
v_vertex = f1.compute_vertex_values(mesh)

# 取梯度的分量 
Df0 = grad(f1)[0]
Df1 = grad(f1)[1]
# # 在 Docs » User manual » Form language 中搜索 grad 可以看到:
# # In UFL, the following pairs of declarations are equivalent:
# Dfi = grad(f)[i]
# Dfi = f.dx(i)
# Dvi = grad(v)[i, j]
# Dvi = v[i].dx(j)
# DAi = grad(A)[..., i]
# DAi = A.dx(i)
# # for a scalar expression f, a vector expression v, and a tensor expression A of arbitrary rank.

上面 并得到函数在*度处的值 需要注意的是 fenics 开始用的是 f1.vector().array(), 但这最新版本中已经不可用了 (参见 Replace vector().array() to vector().get_local()).

2. 将两个函数作为分量形成一个向量函数

from dolfin import *

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)

u = Function(V)
u.interpolate(Expression("1+0*x[0]", degree=1))
v = Function(V)
v.interpolate(Expression("2+0*x[0]", degree=1))
W = VectorFunctionSpace(mesh, 'CG', 1, dim=2)
w = Function(W)
ua = u.vector().get_local()
va = v.vector().get_local()
wa = w.vector().get_local()
wa[::2] = ua
wa[1::2] = va
w.vector().set_local(wa)