Hi, I’m having trouble inverting the sparse matrix using the script below. For any order p, I think the number of freedofs of the P space should equal the number of dofs of the A space, namely N*(p+1), where N is the number of 1D elements, resulting in a square matrix which can be shown to be invertible.
In the case p=0, both umfpack and pardiso give a warning ‘sizes don’t match’, but the inverse is computed anyway yielding the correct solution. When p>0, both umfpack and pardiso throw an exception ‘Symbolic factorization failed’ when the inverted matrix is multiplied by a vector of length N*(p+1). It behaves as if the non-freedof of the P space due to the Dirichlet boundary condition is not being accounted for. This issue does not occur if I use a product space, e.g. AP = FESpace([A, P])
[code]from ngsolve import x, dx, grad
import ngsolve as ng
from ngsolve.meshes import Make1DMesh
N = 10
p = 1
mesh = Make1DMesh(N)
A = ng.L2(mesh, order=p, all_dofs_together=True)
P = ng.H1(mesh, order=p+1, dirichlet=‘right’)
theta = A.TestFunction()
phi = P.TrialFunction()
a2p = ng.BilinearForm(trialspace=P, testspace=A)
a2p += grad(phi) * theta * dx
a2p.Assemble()
a2pinv = a2p.mat.Inverse(freedofs=P.FreeDofs(), inverse=‘umfpack’)
m = ng.BilinearForm(A)
u, v = A.TnT()
m += u * v * dx
m.Assemble()
gfalpha = ng.GridFunction(A)
gfphi = ng.GridFunction(P)
gfalpha.Set(x)
tvec = gfalpha.vec.CreateVector()
tvec.data = m.mat * gfalpha.vec
gfphi.vec.data = a2pinv * tvec
print(“solution:”, gfphi.vec)
[/code]
Thanks for taking a look!
Best,
Dow
Attachment: a9b4f8838241382bb8d112a8feea5edf