### Create a FEniCS matrix from a numpy matrix

104
views
0
6 weeks ago by
For purposes of exploring the way that linear solvers work, I want to create FEniCS Matrix's and Vector's with specified numeric values. I figured out how to do Vectors:

%matplotlib inline
import numpy as np
import fenics as fe
import petsc4py

n = 8

vec1 = np.random.normal(size=(8))
vec1
Out[]: array([-0.07560725,  0.24136354,  0.33426323, -0.48336659, -1.06575918,
-0.260127  ,  1.07919618, -0.11751131])

fvec1 = fe.Vector()
fvec1.init(n)
fvec1[:] = vec1
[ fvec1.size(), fvec1.array() ]

Out[]: [8, array([-0.07560725,  0.24136354,  0.33426323, -0.48336659, -1.06575918,
-0.260127  ,  1.07919618, -0.11751131])]
​

But I can't make it work for Matrix's. Here's my best attempt:

mat1 = np.random.normal(size=(n,n))
np.linalg.det(mat1)

Out[]: 19.882799874591505

fmat1 = fe.Matrix()

fmat1.init(8)
TypeError                                 Traceback (most recent call last)
<ipython-input-10-008e63fd4c54> in <module>()
----> 1 fmat1.init(8)

TypeError: in method 'GenericTensor_init', argument 2 of type 'dolfin::TensorLayout const &'

I looked at the documentation for this TensorLayout object and utterly failed to understand how to create one of them. If I try to go ahead with the Matrixobject I just created, e.g.

rows = cols = np.arange(n, dtype=np.intc)
fmat1.set(mat1.flatten(), rows, cols)

Then, (in a Jupyter notebook) I get a popup window saying "The kernel appears to have died. It will restart automatically." and I'm back to square one.

This can't really be so difficult.

How can I create an 8x8 Matrix and fill in particular numerical values?

Thanks.
Community: FEniCS Project

### 1 Answer

2
6 weeks ago by
OK, I figured it out, with some help from this question. Here's how you can do it:
n = 8

mat1 = np.random.normal(size=(n,n))
np.linalg.det(mat1)

Out[]: 6.6521476157180395

Mat = petsc4py.PETSc.Mat

pmat1 = Mat()
pmat1.create()
pmat1.setSizes((n, n))
pmat1.setType('seqaij')
pmat1.setUp()

Out[]: <petsc4py.PETSc.Mat at 0x7f246e0bd830>

pmat1[:,:] = mat1
pmat1.assemble()

pmat1[:, :]

Out[]:
array([[-0.13598332, -0.9572488 , -0.15780098, -0.1075275 , -1.38674517,
1.52496236, -0.85228079,  0.7654958 ],
[ 0.1593528 ,  0.31010832,  0.47407576,  0.43375233, -0.85509651,
0.71740213,  0.7670229 , -0.73699382],
[-0.95699298,  0.21158925,  0.07723483, -0.82385875,  0.43406038,
-0.49720557,  0.10707859, -0.39464126],
[ 0.25099417, -0.61912983, -0.89437139, -0.14872791,  0.82352802,
-1.69781   , -1.41689846,  0.63905837],
[ 2.04737697, -1.04603641, -0.86743069,  0.92215806,  0.59132639,
0.48615139, -0.74210908, -0.66687863],
[-0.55613246, -0.60924786,  0.25821336, -0.30575691, -0.11805725,
0.31867436, -0.73913848,  0.76352502],
[-0.03747193,  0.8156354 , -0.64619599, -0.23006657, -1.03477191,
-1.67011824, -0.40921348,  1.33355617],
[-0.92012716,  2.0877233 , -0.64820358, -0.12697107, -0.15428069,
0.31432339, -1.68075949, -1.78545492]])

fmat1 = fe.PETScMatrix(pmat1)
fmat1

Out[]: <dolfin.cpp.la.PETScMatrix; proxy of <Swig Object of type 'std::shared_ptr< dolfin::PETScMatrix > *' at 0x7f269c099660> >

fmat1.array()

Out[]:
array([[-0.13598332, -0.9572488 , -0.15780098, -0.1075275 , -1.38674517,
1.52496236, -0.85228079,  0.7654958 ],
[ 0.1593528 ,  0.31010832,  0.47407576,  0.43375233, -0.85509651,
0.71740213,  0.7670229 , -0.73699382],
[-0.95699298,  0.21158925,  0.07723483, -0.82385875,  0.43406038,
-0.49720557,  0.10707859, -0.39464126],
[ 0.25099417, -0.61912983, -0.89437139, -0.14872791,  0.82352802,
-1.69781   , -1.41689846,  0.63905837],
[ 2.04737697, -1.04603641, -0.86743069,  0.92215806,  0.59132639,
0.48615139, -0.74210908, -0.66687863],
[-0.55613246, -0.60924786,  0.25821336, -0.30575691, -0.11805725,
0.31867436, -0.73913848,  0.76352502],
[-0.03747193,  0.8156354 , -0.64619599, -0.23006657, -1.03477191,
-1.67011824, -0.40921348,  1.33355617],
[-0.92012716,  2.0877233 , -0.64820358, -0.12697107, -0.15428069,
0.31432339, -1.68075949, -1.78545492]])

isinstance(fmat1, fe.GenericMatrix)

Out[]: True​
Here's the notebook in nbviewer
Please login to add an answer/comment or follow this question.

Similar posts:
Search »