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

More
1 year 6 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.

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

More
1 year 6 months ago - 1 year 6 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: 1 year 6 months ago by lkogler.

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

More
1 year 6 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.

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

Time to create page: 0.122 seconds