### Create a FEniCS matrix from a numpy matrix

104

views

0

For purposes of exploring the way that linear solvers work, I want to create FEniCS

How can I create an 8x8

Thanks.

`Matrix`

's and `Vector`

's with specified numeric values. I figured out how to do `Vector`

s:```
%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 `Matrix`

object 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.

How can I create an 8x8

`Matrix`

and fill in particular numerical values?Thanks.

Community: FEniCS Project

### 1 Answer

2

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.