Hello everyone,

I’m trying to implement the mixed method from this paper for the Linear Elasticity problem. When it comes to solving the resulting linear system, I’m using

```
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None}
problem = LinearProblem(a, L, bcs=[], petsc_options=petsc_options)
wh = problem.solve()
```

which, from my understanding, employs the standard LU solver.

By doing this, I get a NaN solution. However, if I change the solver package to superlu

```
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "superlu", "ksp_monitor": None}
```

the method now works as expected.

Can someone help me figure out what is happening? I’m using the 0.7.0 version of dolfinx installed through conda. Here is the full code for context

```
import dolfinx.mesh as mesh
from dolfinx.fem import (FunctionSpace, assemble_scalar, form)
from dolfinx.fem.petsc import LinearProblem
import ufl
from ufl import (inner, outer, sym, div, grad, tr, dx, Identity, TrialFunction,
TestFunction, FiniteElement, VectorElement, MixedElement,
split, as_vector)
from mpi4py import MPI
import numpy as np
# Problem definition ==========================================================
# Domain (rectangle)
cornerA = [0, 0]
cornerB = [1, 1]
# Mesh
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array(cornerA), np.array(cornerB)],
[30, 30], cell_type=mesh.CellType.quadrilateral)
# Lamé parameters
mu = 0.3
lamb = 1
# Elasticity tensor
def tensorC(u):
return 2*mu*sym(grad(u)) + lamb*div(u)*Identity(2)
# Compliance tensor
def tensorA(s):
return (1/(2*mu))*(s - (lamb/(2*lamb+2*mu))*tr(s)*Identity(2))
# Exact solutions
def u_aux(mod):
return [ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1]),
ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1])]
x = ufl.SpatialCoordinate(msh)
u_ufl = ufl.as_vector(u_aux(x))
sig_ufl = tensorC(u_ufl)
r_ufl = (grad(u_ufl)[0,1] - grad(u_ufl)[1,0])/2
# Load function
load = div(sig_ufl)
# Mixed method with wak symmetry ==============================================
# Approximation spaces
degree_s = 2
degree_u = 1
degree_r = 1
el_sig = FiniteElement('RTCF', msh.ufl_cell(), degree_s)
el_u = VectorElement('DQ', msh.ufl_cell(), degree_u)
el_rot = FiniteElement('DPC',msh.ufl_cell(), degree_r)
W = FunctionSpace(msh,MixedElement([el_sig, el_sig, el_u, el_rot]))
# Test and trial functions
e1 = as_vector([1, 0])
e2 = as_vector([0, 1])
TrialF = TrialFunction(W)
TestF = TestFunction(W)
(sig1, sig2, u, r) = split(TrialF)
(tau1, tau2, v, q) = split(TestF)
sig = outer(e1, sig1) + outer(e2,sig2)
tau = outer(e1, tau1) + outer(e2, tau2)
# Variational problem
a = inner(tensorA(sig), tau)*dx
a += inner(u,div(tau))*dx
a += r*(tau[0,1]-tau[1,0])*dx
a += inner(v,div(sig))*dx
a += q*(sig[0,1]-sig[1,0])*dx
L = inner(load,v)*dx
# Linear system solver
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None}
#petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "superlu", "ksp_monitor": None}
problem = LinearProblem(a, L, bcs=[], petsc_options=petsc_options)
wh = problem.solve()
(sig1h, sig2h, uh, rh) = wh.split()
sigh = outer(e1, sig1h) + outer(e2, sig2h)
# Approximation errors ========================================================
esig = sig_ufl - sigh
eu = u_ufl - uh
er = r_ufl - rh
eL2sig = np.sqrt(assemble_scalar(form(inner(esig,esig)*dx)))
eL2u = np.sqrt(assemble_scalar(form(inner(eu,eu)*dx)))
eL2r = np.sqrt(assemble_scalar(form(inner(er,er)*dx)))
print(eL2sig, eL2u, eL2r)
# =============================================================================
```