- Thank you received: 0
Conversion of NGSolve matrices to numpy/scipy for 3D function space
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:
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:
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.
1 year 6 months ago - 1 year 6 months ago #3483
by lkogler
Replied by lkogler on topic Conversion of NGSolve matrices to numpy/scipy for 3D function space
Instead of
use
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
Code:
V = ngs.H1(mesh, order=2, dirichlet='left', dim=mesh.dim)
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)
Best,
Lukas
Attachments:
Last edit: 1 year 6 months ago by lkogler.
Please Log in or Create an account to join the conversation.
1 year 6 months ago #3484
by Daniel
Replied by Daniel on topic Conversion of NGSolve matrices to numpy/scipy for 3D function space
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.
Time to create page: 0.122 seconds