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.