- Posts: 2
- Thank you received: 0

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:

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.

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 ofuse

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)

Best,

Lukas

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

`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)

Best,

Lukas

Last Edit: 1 week 2 days ago by lkogler.

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.

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.