FEniCS中如何提取出刚度矩阵

我在https://fenicsproject.org/qa/12488/jacobian-matrix-extracted-array-sparsity-pattern-visualized/中找到了类似的解答，整理下来方便以后查阅。

网页中所提问题是：

Hello. In my code, I compute a Jacobian using the following command:
J = derivative(F_sys,current_sol)
where F_sys is the system of equations to solve and current_sol is the solution at the current time step.
How can I:
Extract the Jacobian matrix as a numpy array, and
Visualize its sparsity pattern (preferably with a built-in FEniCS command if one exists)?
Thank you for helping. Let me know if you need any additional info.

其他人给出的回答：
代码1，代码1给出的是一个非稀疏的矩阵

from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)
u, v = Function(V), TestFunction(V)

F = dot(grad(u), grad(v))*dx - Constant(1.0)*v*dx
J = derivative(F, u)

J_mat = assemble(J)
J_array = J_mat.array()
print J_array
plt.spy(J_array)
plt.show()

代码2，代码2给出了如何输出稀疏矩阵的方法

from dolfin import *
import ufl
import scipy.sparse as sp
import matplotlib.pyplot as plt


def spy_form_matrix(A, l=None, bcs=None, marker_size=2.0, show=True, pattern_only=True, **kwargs):
    assert isinstance(A, ufl.form.Form)
    if l: assert isinstance(l, ufl.form.Form)
    eigen_matrix = EigenMatrix()
    if not l:
        assemble(A, tensor=eigen_matrix)
        if not bcs is None:
            for bc in bcs: bc.apply(eigen_matrix)
    else:
        SystemAssembler(A, l, bcs).assemble(eigen_matrix)
    A = eigen_matrix

    row, col, data = A.data()
    if pattern_only:
        data[:] = 1.0

    sp_mat = sp.csr_matrix((data, col, row), dtype='float')
    plt.spy(sp_mat, markersize=marker_size, precision=0, **kwargs)
    if show: plt.show()


mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)
u, v = Function(V), TestFunction(V)

F = dot(grad(u), grad(v))*dx - Constant(1.0)*v*dx
J = derivative(F, u)

spy_form_matrix(J)

下面的代码可以直接从assemble后得到的矩阵变成python里的稀疏矩阵：

fe.parameters['linear_algebra_backend'] = 'Eigen'
def tran2SparseMatrix(A):
    row,col,val = fe.as_backend_type(A).data()
    return sps.csr_matrix((val,col,row))

下面的代码给出了导出numpy矩阵后和直接用FEniCS计算的对比

import numpy as np 
import matplotlib.pyplot as plt 
import fenics as fe 
import scipy.sparse as sps
import scipy.sparse.linalg as spsl
import time

fe.parameters['linear_algebra_backend'] = 'Eigen'
def tran2SparseMatrix(A):
    row,col,val = fe.as_backend_type(A).data()
    return sps.csr_matrix((val,col,row))

mesh = fe.UnitSquareMesh(100, 100)
V = fe.FunctionSpace(mesh, 'P', 2)

def boundary(x, on_boundary):
    return on_boundary

du, v = fe.TrialFunction(V), fe.TestFunction(V)
f = fe.Constant(1.0)

af = fe.inner(fe.grad(du), fe.grad(v))*fe.dx
bf = fe.inner(f, v)*fe.dx
A_c = fe.assemble(af)
b_c = fe.assemble(bf)
bc = fe.DirichletBC(V, fe.Constant(0.0), boundary)
bc.apply(A_c, b_c)

u = fe.Function(V)
start = time.time()
fe.solve(A_c, u.vector(), b_c)
end = time.time()
print(end-start)

uu1 = fe.Function(V)
uu2 = fe.Function(V)
start = time.time()
A = tran2SparseMatrix(A_c)
b = np.array([b_c[:], b_c[:]]).T
temp = spsl.spsolve(A, b)
end = time.time()
print(end-start)
uu1.vector()[:] = temp[:, 0]
uu2.vector()[:] = temp[:, 1]

en = fe.assemble(fe.inner(u-uu1, u-uu1)*fe.dx)
print('Error is ', en)
en = fe.assemble(fe.inner(u-uu2, u-uu2)*fe.dx)
print('Error is ', en)