Conversion of NGSolve matrices to numpy/scipy for 3D function space

Main Forum Old-Forum - do not post here
1

Dear Forum,

I would like to inspect the stiffness and the mass matrix in numpy/scipy for a problem in 3D elasticity. To convert the ngsolve matrices to scipy I use the code from the how-to in the documentation. The code works well for an example with 1D Poisson equation. However, when going to a 3D space there are some problems which I cannot sort out and for which I cannot find documentation. A minimum example which leads to the problem is this:

import ngsolve as ngs
import netgen.csg as csg
import scipy.sparse as sparse

# Construction of geometry and mesh
left  = csg.Plane(csg.Pnt(0,0,0), csg.Vec(-1,0,0) )
right = csg.Plane(csg.Pnt(1,1,1), csg.Vec( 1,0,0) )
front = csg.Plane(csg.Pnt(0,0,0), csg.Vec(0,-1,0) )
back  = csg.Plane(csg.Pnt(1,1,1), csg.Vec(0, 1,0) )
bot   = csg.Plane(csg.Pnt(0,0,0), csg.Vec(0,0,-1) )
top   = csg.Plane(csg.Pnt(1,1,1), csg.Vec(0,0, 1) )
left.bc('left')

cube = left * right * front * back * bot * top
geo = csg.CSGeometry()
geo.Add(cube)

mesh = ngs.Mesh(geo.GenerateMesh(maxh=0.1))

# Definition of function space
V = ngs.H1(mesh, order=2, dirichlet='left', dim=mesh.dim)

# Definition of bilinear forms
mu = 1  # all constants set to zero
lmbda = 1
rho = 1

def epsilon(u):
    return 0.5 * (ngs.Grad(u) + ngs.Grad(u).trans)

def sigma(u):
    return 2 * mu * epsilon(u) + lmbda * ngs.Trace(epsilon(u)) * ngs.Id(3)

u, v = V.TnT()
k = ngs.BilinearForm(V, symmetric=True)
k += ngs.InnerProduct(sigma(u), epsilon(v)) * ngs.dx
k.Assemble()

m = ngs.BilinearForm(V, symmetric=True)
m += rho * ngs.InnerProduct(u, v) * ngs.dx
m.Assemble()

# Conversion of matrices to sparse scipy matrices

#rows,cols,vals = k.mat.COO()
rows,cols,vals = k.mat.CSR()
#k_np = sparse.csr_matrix((vals,(rows,cols)))

print(type(rows))
print(type(cols))
print(type(vals))

# Output
# <class 'ngsolve.bla.FlatVectorD'>
# <class 'pyngcore.FlatArray_I_S'>
# <class 'pyngcore.FlatArray_S_S'>

The method COO() does not seem to exist for the form matrix k. Instead I was using the CSR method. However, the CSR method returns a FlatVectorD, a FlatArray_I_S and a FlatArray_S_S object instead of the raw row and column indices and the raw values. I could not find any information how to constrcut from these objects a sparse matrix in scipy. I would be very happy if you could help me with this issue. Thanks in advance.

2

Instead of

V = ngs.H1(mesh, order=2, dirichlet='left', dim=mesh.dim)

use

V = ngs.VectorH1(mesh, order=2, dirichlet='left')

The difference is that the first space generates a matrix where every entry is a small 2x2 block, whereas the second one generates a normal sparse matrix.

Edit:
You can also get the block-sparse-matrix into scipy, but it requires some loops in python (see attached file)

https://ngsolve.org/media/kunena/attachments/779/tosparse.py

Best,
Lukas

Attachment: tosparse.py

3

Dear Lukas,

Thanks for your quick reply. Both of your solutions run on my computer and are a great help.

Just one other question about the conversion to scipy matrices which I was actually wondering about. In the how-to in the documentation, the matrix data is extracted from ngsolve using the COO() method (rows,cols,vals = a.mat.COO()). I assume that one should create a spcipy matrix in COO format. However, in the next line a matrix in csr format is created (A = sp.csr_matrix((vals,(rows,cols)))). What is the thought behind that? Should this be used for symmetric bilinear froms? Thanks in advance.