How to get Sparse matrix


191
views
1
5 months ago by
Hi guys,

If I compare the matrix matA after assemble I get a banded matrix.
p = TrialFunction(V)
q = TestFunction(V)

w=0.0025
mu = 8.9E-4
f=Wells(degree=3)

F = (w**3/(12*mu))*inner(grad(q),grad(p))*dx-f*q*dx

a, L = lhs(F),rhs(F)

matA = assemble(a)
​

How could I get the sparse matrix related to this matrix? I would like to use this sparse matrix to benchmark a FEM software that is developed in matlab.

Thank you very much
Community: FEniCS Project

2 Answers


2
5 months ago by
Ruben,

You can get a sparse matrix constructed in FEniCS, and written to disk in Matlab format (.mat file), which you can load into Matlab by issuing the command

>> load A

as follows:


Step 1: Include the following two line at the top of your python file.
from scipy.sparse import csr_matrix
import scipy.io as io


Step 2: Write your python lines which ultimately make the matrix, i.e.,
.
.
.

matA = Assemble(a)


Step 3: Write the matrix matA as A.mat on disk (readable by Matlab) by including the following lines.
A = as_backend_type(matA).mat()
r, c = A.getSize()
rptr, cind, vals = A.getValuesCSR()
A = csr_matrix((vals,cind,rptr), shape=(r,c))
io.savemat('A.mat',{'A':A})


The first line

A = as_backend_type(matA).mat()

retrieves the matrix from FEniCS's backend, whichever is being used.

The next line

r,c = A.getSize()

sets the number of rows in r, and the number of cols in c

The next line
rptr, cind, vals = A.getValuesCSR()

fills the arrays rptr, cind, and vals with the Compressed Row Sparse data of your matrix. Row pointers go in rptr, column indices go in cind, nonzero values go in vals

The next line
A = csr_matrix((vals,cind,rptr), shape=(r,c))

makes a scipy matrix A (over-riding the previous assignment of A) with the compressed row sparse data vals, cind, rptr.

The next line
io.savemat('A.mat',{'A':A})

saves this matrix on disk as A.mat. In this matrix the data of the scipy matrix A is stored (second part, second argument). When you finally load it in Matlab using the "load" command, the data is loaded in Matlab's workspace as the matrix A (this is specified in single quotes in the first part, second argument).


Cheers,

Hisham
Works perfect, thank you very much!!!
written 5 months ago by Ruben Gonzalez  
1
5 months ago by
You can use PETSc viewers through petsc4py to get the matrix in Matlab sparse format:

from dolfin import *

import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh,"Lagrange",1)

p = TrialFunction(V)
q = TestFunction(V)

w=0.0025
mu = 8.9E-4
f=Constant(0.0) #Wells(degree=3)

F = (w**3/(12*mu))*inner(grad(q),grad(p))*dx-f*q*dx

a, L = lhs(F),rhs(F)

matA = assemble(a)

viewer = PETSc.Viewer().createASCII("a.dat","w")
viewer.pushFormat(PETSc.Viewer.Format.ASCII_MATLAB)
viewer(as_backend_type(matA).mat())​
Thank you for your answer, I got the sparse matrix
written 5 months ago by Ruben Gonzalez  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »