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

1 week 3 days ago #3482 by Daniel
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.

Please Log in or Create an account to join the conversation.

1 week 2 days ago - 1 week 2 days ago #3483 by lkogler
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)

File Attachment:

File Name: tosparse.py
File Size:3 KB


Best,
Lukas
Attachments:

Please Log in or Create an account to join the conversation.

1 week 2 days ago #3484 by Daniel
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.

Please Log in or Create an account to join the conversation.

© 2019 Netgen/NGSolve