记录 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)

3. 将形成的刚度矩阵转成 numpy.array() 格式

a = inner(grad(u), grad(v)) * dx
aa = assemble(a)
aM = aa.array()

L = inner(f, q) * dx
ll = assemble(L)
l1[:]  # 这样可以列出右端项所有的值

4. 得到 Constant 的值

aa = Constant(1e-3)
print(aa) # Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 49)
aa_ = float(aa)
print(aa_)

5. 去掉求解器中 “Solving linear variational problem” 提示

参考: How to disable message “Solving linear variational problem.”

set_log_level(level)
CRITICAL = 50, // errors that may lead to data corruption and suchlike
ERROR = 40, // things that go boom
WARNING = 30, // things that may go boom later
INFO = 20, // information of general interest
PROGRESS = 16, // what’s happening (broadly)
TRACE = 13, // what’s happening (in detail)
DBG = 10 // sundry

添加 set_log_active(False) 可以关掉所有提示 (包括 ‘Error’, ‘Warning’ …)
自己测试的是可以直接 set_log_active(50) 是可以关掉 “Solving linear variational problem” 提示.

6. 得到插值函数在某点处的值

例如定义了 fenics 空间 Vh, 以及给定 真解 Phi = Expression("x[0] + x[1]", degree=3) (即 Phi = x + y), 然后 Phi 插值到 Vh 空间中: Phi_h = interpolate(Phi, Vh), 我们想求在 点 P=(0.3, 1.2) 处的函数值

Phi_h = interpolate(Phi, Vh)
P = (0.3, 1.2)
val = Phi_h(P)