Error while printing stiffness matrix and saving it to an external file

Question

Hi, This is a trivial question but any help is highly appreciated. I want to view my stiffness matrix and save it to an external file? Is it possible?

import numpy as np

a = inner(sigma, grad(v))*dx 
L = -(dot(r, v)*dx + dot(f,v)*ds(1))

A,b = assemble_system(a,L,bc)

# View entire array ( not truncated form) 
np.set_printoptions(threshold='nan')
# print entire stiffness matrix and it's largest entries
print A.array(), np.sort(A.array()[abs(A.array())>1e2])     
Y = A.array()

The above is not giving me any output and when I use command File("stiffness.pvd") << Y, that doesn't work either.

Comments

It will be easier for someone to help you if you post the full code (as long as you manage to write a short complete code). You are missing definitions of at least mesh, function space, test and trial functions above, so it is not easy for others to understand what the real problem is.

---

Ok this is the code. It just keeps running, never gives an output. Thanks for your time and help. Highly appreciate it.

from dolfin import *
import numpy as np

# Load mesh and define function space
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d0*6)
n1 = int(d1*6)
n2 = int(d2*6)
o = Point(0, 0, 0)
l = Point(8, 1, 1)

# Mesh for cantilever beam
Mesh = BoxMesh(o, l, n0, n1, n2)

Mesh.init()

#Note that for efficiency, only entities of dimension zero (vertices) and entities of the maximal dimension (cells) exist when creating a Mesh. Other entities must be explicitly created by calling init(). For example, all edges in a mesh may be created by a call to mesh.init(1). 

V  = VectorFunctionSpace(Mesh, "CG", 1)

tol = 1E-14 # tolerance http://www.stata.com/manuals14/m-1tolerance.pdf

# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

# Right boundary
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)      

# Top face at which load is added
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 1.0)

#Boundary segments
left = Left()
top = Top()
right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", Mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)      

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, -1.9, 0.0))

# Elasticity parameters
E = 200*(10**9)
nu = 0.24

mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Strain
def eps(v):
    return sym(grad(v))

# Stress
def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d) 

# Define new measures associated with the interior domains and exterior boundaries
ds = Measure("ds")[boundaries]

# Define variational problem
g_T= Constant((0.0, -8500000.0, 0.0))

a = inner(grad(v),sigma(u))*dx 

L= dot(v,g_T)*ds(2)+dot(v,g_T)*dx

# Boundary condition
left_end_displacement = Constant((0.0, 0.0, 0.0))

bc = [DirichletBC(V, left_end_displacement, left)] 

# Compute solution
u = Function(V)
A,b = assemble_system(a,L,bc)
u = Function(V)
U = u.vector()
solve(A,U,b)

np.set_printoptions(threshold='nan')
# print entire stiffness matrix and it's largest entries
print A.array(), np.sort(A.array()[abs(A.array())>1e2])     
Y = A.array()
File("stiffnessmatrix.pvd", "compressed") << Y

# Save solution to VTK format
File("Output_elas3D_PointLoad.pvd", "compressed") << u
plot(u, mode="displacement",interactive=True)

Answer 1

Maybe you should to do this to export:

np.savetxt("stiffnessmatrix.txt", A.array())

and this to print:

print("stiffness_matrix: "), A.array()

(*) It is perhaps a good idea to reduce (considerably) the number of nodes in the mesh to see if what you are printing it's o.k., e.g, run the code with this mesh:

Mesh = BoxMesh(o, l, 1, 1, 1)

Comments

Thanks a lot hernan_mella! You are awesome. Even though I still don't understand why my version did not work and yours did.

Now I want to make sure that what I am printing is correct. I think it can try to view the stiffness matrix for the same problem in Fenics and Matlab Puffin and compare them. Do you have any personal suggestions?

Thank you!

---

In your code the line that cause the problem is:

np.set_printoptions(threshold='nan')

honestly i don't know why this happens (i don't use python language to write my codes, but maybe is a bug), but if you eliminate this line your code works perfectly.

Answering at your last question, i don't have any suggestion for you (I
I just found out about of the existence of puffin =( ).

regards.