### Getting the system matrix for a Poisson problem in C++

187

views

0

Hello,

I would like to use FEniCS as the initial interface to get my system matrix from a physical problem. Let's say I have a Poisson problem. The form file would be the same as here and the the main file would be here . with the difference that I use a quadrilateral element instead of a triangle element. But how do I get the values of the matrix into three arrays (the row pointer, column index and the value array) for the CSR format. I understand that I can use

which is, I would expect that if I have a square mesh with quadrilateral elements instead of triangle elements, I should get the system matrix same as that of the finite difference with a 5 point stencil. Am I wrong here ? Sorry, I am not that adept at Finite element method.

Thank you.

I would like to use FEniCS as the initial interface to get my system matrix from a physical problem. Let's say I have a Poisson problem. The form file would be the same as here and the the main file would be here . with the difference that I use a quadrilateral element instead of a triangle element. But how do I get the values of the matrix into three arrays (the row pointer, column index and the value array) for the CSR format. I understand that I can use

`A.get(double *block,size_t m,const int *row,size_t n,const int*col)`

to get the block of the matrix after assembling it from the bilinear form`assemble(A,a);`

but that doesnt seem to give me the values that I want, which leads to my next question.which is, I would expect that if I have a square mesh with quadrilateral elements instead of triangle elements, I should get the system matrix same as that of the finite difference with a 5 point stencil. Am I wrong here ? Sorry, I am not that adept at Finite element method.

Thank you.

Community: FEniCS Project

### 1 Answer

2

Here's python code to create assemble the stiffness matrix and save it in CSR format to a file.

```
import numpy as np
from fenics import *
parameters['linear_algebra_backend'] = 'Eigen'
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = (inner(grad(u), grad(v)) + u*v)*dx
A = assemble(a)
Aspa = as_backend_type(A).sparray()
np.savez_compressed('stiffness-matrix.npz', data=Aspa.data, indices=Aspa.indices, indptr=Aspa.indptr)
```

Later you can read the file with

```
stiffmat = np.load(file1)
data = stiffmat['data']
indices = stiffmat['indices']
indptr = stiffmat['indptr']
```

written
5 months ago by
dommath

Please login to add an answer/comment or follow this question.

https://www.allanswered.com/post/nqzoo/