- Thank you received: 0
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
Conversion of NGSolve matrices to numpy/scipy for 3D function space
3 years 3 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.
3 years 3 months ago - 3 years 3 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: 3 years 3 months ago by lkogler.
3 years 3 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.
Time to create page: 0.152 seconds
