Hey all,

I am solving a linear system (A*U=F) and need to gather information about the A matrix to select an appropriate iterative solver for my simulation. Specifically, I want to know about the sparsity pattern of the matrix, including:

- the number of zeros in the matrix,
- their locations, and
- the condition number of the matrix.

Below is a minimal code example for elasticity where Matrix A is the one I want to analyze. However, I’m not sure how to extract this information. What is the best way to do that?

```
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import dolfinx
# -------------------------Material Properties
E = 1 * 130e9
nu = 0.28
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# ......................................... Reading mesh file
h = 0.5
mesh= create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[15, 15, 15], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200.0
Area2 = 0.5**2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)) # Displacement Element
Space = functionspace(mesh, el)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], h, atol=tol))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# # ...................................... BCs
# Clamp BC for left side
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
clamp_dofs = locate_dofs_topological(Space, mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(np.array([0, 0, 0], dtype=default_scalar_type), clamp_dofs, Space)
bc = [bc0] # BC 4
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l)) # Elasticity tensor
def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
# ......................................... Define Variational Formulation
VarForm = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()
with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [Result]) as vtx:
vtx.write(0.0)
```

I found a previous post here that seemed helpful, but I encountered the follwing error when I ran it. I am using DOLFINx version 0.8.0.