Sparse inverse issue in 1D

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

Hi Dow,

Inverse(freedofs=...) 

does not distinguish between rows and cols, so this does not work.

You can eliminate the right-most dof completely from the system, by setting it as unused-dof, and then compress the space to kick-out the unused dofs completely (P.ndof changes):

https://ngsolve.org/media/kunena/attachments/744/blftest_post.py
https://ngsolve.org/media/kunena/attachments/744/blftest_post.py

P = ng.H1(mesh, order=p+1, dirichlet='right')
P.couplingtype[N] = COUPLING_TYPE.UNUSED_DOF
P = ng.Compress(P)
print (P.ndof, A.ndof)

best,
Joachim

Attachment: 638a32b3c4b373604cba265c221f9340