How print out stiffness matrix as a text file


134
views
2
3 months ago by
Ali  
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.

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
3 months ago by
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.

Similar posts:
Search »