### How print out stiffness matrix as a text file

134

views

2

Hi All,

I am trying to figure out how can I print out stiffness matrix of a simple script as a text file.

For example, This simple vector Laplacian.

I am trying to figure out how can I print out stiffness matrix of a simple script as a text file.

For example, This simple vector Laplacian.

from dolfin import * # Create mesh mesh = UnitSquareMesh(32, 32) # Define Exact solution u_D = Expression(('0.5*(x[0] - pow(x[0], 2))','x[1] - pow(x[1], 2)'),degree=2) # Define function spaces NED = FiniteElement("N1curl", mesh.ufl_cell(), 1) CG = FiniteElement("CG", mesh.ufl_cell(), 1) W = CG * NED V = FunctionSpace(mesh,W) # Define trial and test functions (sigma, u) = TrialFunctions(V) (tau, v) = TestFunctions(V) # Define bilinear form and right hand side a = (sigma*tau - dot(u, grad(tau))+ dot(grad(sigma), v) + inner(curl(u),curl(v)))*dx f = Constant((-1., -2.)) L = dot(f, v)*dx # solve and plot w = Function(V) solve(a == L, w) sigma, u = w.split() plot(u) import matplotlib.pyplot as plt plt.show() # Save solution to file in VTK format vtkfile = File('veclapNED2/solution.pvd') vtkfile << u # Compute error in L2 norm error_L2 = errornorm(u_D, u, 'L2') # Print errors print('error_L2 =', error_L2)

Community: FEniCS Project

### 1 Answer

4

You can obtain human-readable forms of the LHS matrix in a few different ways. I've included two possibilities in the following modified version of your script:

```
from dolfin import *
from petsc4py import PETSc
from numpy import array2string
import sys
# Create mesh
# If you want to print the entire matrix, including zeros, it will be huge,
# even for small meshes.
#mesh = UnitSquareMesh(32, 32)
mesh = UnitSquareMesh(2, 2)
# Define Exact solution
u_D = Expression(('0.5*(x[0] - pow(x[0], 2))','x[1] - pow(x[1], 2)'),degree=2)
# Define function spaces
NED = FiniteElement("N1curl", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = CG * NED
V = FunctionSpace(mesh,W)
# Define trial and test functions
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
# Define bilinear form and right hand side
a = (sigma*tau - dot(u, grad(tau))+ dot(grad(sigma), v) + inner(curl(u),curl(v)))*dx
f = Constant((-1., -2.))
L = dot(f, v)*dx
###############################################################################
# Obtain the LHS matrix as a GenericTensor. (You can get a GenericVector
# of the RHS as well, with assemble(L).)
matA = assemble(a)
# Possiblity 1:
# Produce Matlab/Octave data file. This writes out an ASCII representation
# of the sparse matrix that can be read in via Matlab/Octave and manipulated
# as needed.
viewer = PETSc.Viewer().createASCII("a.m","w")
viewer.pushFormat(PETSc.Viewer.Format.ASCII_MATLAB)
viewer(as_backend_type(matA).mat())
# Possibility 2:
# Print out the entire matrix, including zeros. (This only works in serial.)
stringA = str(matA.array())
outf = open("a.txt","w")
outf.write(array2string(matA.array(),max_line_width=sys.maxint))
outf.close()
exit()
###############################################################################
# solve and plot
w = Function(V)
solve(a == L, w)
sigma, u = w.split()
plot(u)
import matplotlib.pyplot as plt
plt.show()
# Save solution to file in VTK format
vtkfile = File('veclapNED2/solution.pvd')
vtkfile << u
# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')
# Print errors
print('error_L2 =', error_L2)
```

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