ValueError: too many values to unpack (expected 1), while solving 2D elasticity problem

I have the following dolfinx code, where I am trying to solve a 2D linear elasticity problem on a 2D plate with 2 materials

import dolfinx
import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, io
import meshio
import gmsh
import numpy as np
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace,
                         form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from ufl import (TestFunction, TrialFunction,
                 dx, grad, inner)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import gmsh
from dolfinx.io import gmshio
rve_mesh, cell_tags, facet_tags = gmshio.read_from_msh("/mnt/c/Users/ABC/sfepy/simplerve.msh", MPI.COMM_WORLD,gdim=2)
ds = ufl.Measure("ds", domain=rve_mesh, subdomain_data=facet_tags)
dx = ufl.Measure("dx", domain=rve_mesh, subdomain_data=cell_tags)
from dolfinx.fem import VectorFunctionSpace
V = VectorFunctionSpace(rve_mesh, ("Lagrange", 1))

def clamped_boundary(x):
    return np.isclose(x[0], 0)

lam_values = [0.005,121]
mu_values = [0.007, 81.9]
mu  = Function(V)
lam = Function(V)
for cell_no in range(len(cell_tags.values)):
    subdomain_no = cell_tags.values[cell_no]
    subdomain_no = int(subdomain_no-1)
    mu.x.array[cell_no] = mu_values[subdomain_no]
    lam.x.array[cell_no] = lam_values[subdomain_no]
u, v = TrialFunction(V), TestFunction(V)
fdim = rve_mesh.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(rve_mesh, fdim, clamped_boundary)

u_D = np.array([0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma_matrix(u):
    return lam_values[0] * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu_values[0]*epsilon(u)

def sigma_fiber(u):
    return lam_values[1] * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu_values[1]*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(rve_mesh, ScalarType((0, 0.01)))
T = fem.Constant(rve_mesh, ScalarType((0, 0)))
a = (ufl.inner(sigma_matrix(u), epsilon(v)) * ufl.dx(1) + ufl.inner(sigma_fiber(u), epsilon(v)) * ufl.dx(2))
L = (ufl.dot(f, v) * ufl.dx(1) + ufl.dot(f, v) * ufl.dx(2) + ufl.dot(T, v) * ds)
problem = fem.petsc.LinearProblem(a, L,[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

The 2nd to last line throws the error.

ValueError: too many values to unpack (expected 1)

I am following the following tutorial
https://jorgensd.github.io/dolfinx-tutorial/chapter2/linearelasticity_code.html

You would need to make the problem reproducible for others, by:

  1. Using a mesh that is accessible for everyone
  2. Making sure that all import statements are included.

As a side note, I do not understand why you are creating two different sigmas

as you could just use mu and lam from

and have a single sigma.

Thank You for the side note
I have linked my mesh (a .msh file) below. Code has been edited in the post to import modules
GDrive link for .msh file

I tried running the again, and now I get a different error

Error: error code 73
[0] KSPSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/ksp/interface/itfunc.c:407
[0] PCSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/pc/interface/precon.c:990
[0] PCSetUp_ILU() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/ksp/pc/impls/factor/ilu/ilu.c:141
[0] MatILUFactorSymbolic() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/mat/interface/matrix.c:6697
[0] MatILUFactorSymbolic_SeqAIJ() at /home/conda/feedstock_root/build_artifacts/petsc_1656827302664/work/src/mat/impls/aij/seq/aijfact.c:1675
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 0

I looked it up and found this, where you have instructed to create a ‘submesh’

I wish to get the solution over the whole domain, which is divided into 2 subdomains. Is it required for me to create a submesh . (Since this step was not required in dolfin).

Unfortunately I am unable to figure out and implement your tip regarding the mu and lambda values. Since the sigma function depends on lambda and mu values, how will I pass the input which will indicate which value of mu/lambda to take, give that the mu.x.array has around 560 entries?

Also regarding the error with missing diagonal, I believe it is happening because I am not implementing the subdomains properly. I do not get any error if I use 1 equation with only ‘dx’, instead of using dx(1) and dx(2) as given in my code. Please refer to my message above

EDIT

I have redefined sigma to accept 2 inputs, one is u, the other being an integer which will either be one or two.

def sigma(u, v):
    return lam_values[v]* ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu_values[v]*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(rve_mesh, ScalarType((0, 0)))
T = fem.Constant(rve_mesh, ScalarType((0, 0)))
a = (ufl.inner(sigma(u, 0), epsilon(v))) * ufl.dx(1) + (ufl.inner(sigma(u, 1), epsilon(v))) * ufl.dx(2) 
L = (ufl.dot(f, v)) * ufl.dx(1) + (ufl.dot(f, v)) * ufl.dx(2) + ufl.dot(T, v) * ds

If you instead of defining mu and lambda in V

And use W = FunctionSpace(mesh, DG,0), the relationship between the cell index and the degree of freedom is 1 to 1, and the relationship

is correct, as explained in:
https://jorgensd.github.io/dolfinx-tutorial/chapter3/em.html#meshing-a-complex-structure-with-subdomains
then any integral with lam or mu will have the correct value in the correct subdomain.