Geo
January 30, 2024, 3:33pm
1
Hi,
I was wondering when I assemble a matrix, say in the Poisson tutorial in the step
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
can I print in the screen A and L or transform them to np matrix or array without resorting to the source code?
dokken
January 30, 2024, 3:48pm
2
For small matrices you can first assemble the matrix with A=dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a))
and call print(A[:,:])
.
For larger matrices one should use the scipy csr matrix. There are plenty of posts with this on the forum.
The vector is easy to print, simply assemble it
b = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L))
print(b.vector.array)
or use
b = dolfinx.fem.assemble_vector(dolfinx.fem.form(L))
print(b.array)
See for instance: Solving PDEs with different scalar (float) types — DOLFINx 0.7.3 documentation
1 Like
Geo
January 30, 2024, 4:09pm
3
Thanks, worked perfectly. I will copy-paste a complete code below in case someone needs it in the future.
from mpi4py import MPI
import dolfinx
import matplotlib.pyplot as plt
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import dolfinx.fem.petsc
import scipy
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim-1, boundary_facets)
print(boundary_dofs)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(domain, ScalarType(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
uh = problem.solve()
# Assemble the sparce matrix
A=dolfinx.fem.assemble_matrix(dolfinx.fem.form(a))
As = scipy.sparse.csr_matrix((A.data, A.indices, A.indptr))
print(As)
# Convert the scipy sparse matrix to a dense numpy array
A_dense = As.toarray()
# Convert the numpy array to a numpy matrix
A_matrix = np.matrix(A_dense)
print(A_matrix)
b = dolfinx.fem.assemble_vector(dolfinx.fem.form(L))
print(b.array)