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

More
3 years 10 months 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:
Code:
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.
More
3 years 10 months ago - 3 years 10 months ago #3483 by lkogler
Instead of
Code:
V = ngs.H1(mesh, order=2, dirichlet='left', dim=mesh.dim)
use
Code:
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:
Last edit: 3 years 10 months ago by lkogler.
More
3 years 10 months 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.
Time to create page: 0.100 seconds