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