For some reasons I need to solve the resulting linear system of FEM discretization manually. What I do is I compute the stiffness matrix A and the right hand side vector L using FeniCS and then solve the resulting system of equations Ax = L with numpy. Here is my code for this for a simple poisson problem.
from fenics import *
import numpy as np
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x):
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",degree=2)
a = inner(grad(u), grad(v))dx
L = fv*dx
M = assemble(a)
anp = np.matrix( M.array() )
M = assemble(L)
Lnp = np.matrix( M.array() )
Lnp = np.transpose(Lnp);
unp = np.linalg.solve(anp,Lnp)
Somehow the stiffness matrix A (or anp in the code) is singular. I do not understand why is that, probably some stupid mistake...
Thanks for the help.
