Hi,

I am trying to solve a linear homogenization problem as described in here but using the total displacement as the unknown.

TLDR; I am solving a simple linear elasticity problem, looking to constrain the displacement on the right (u_R) and left (u_L) edges of a simple unit square as:

\hspace{15em} u_R - u_L = \widetilde{u} (*prescribed*)

I am thinking of using a lagrange multiplier (\lambda \in \mathbb{R}^2 (\Gamma)) defined on the left and right edge to impose the above constrain, and tweaking the bilinear and linear forms respectively as:

\hspace{1em} a((u, \lambda), (\delta u, \delta\lambda)) = {\large\int}_\Omega \sigma(u):\varepsilon(\delta u)dx + {\large \int}_\Gamma\ \ \delta\lambda\cdot (\delta u_R - \delta u_L) d\Gamma + {\large \int}_\Gamma\ \ \lambda \cdot (\delta u_R - \delta u_L) d\Gamma

\hspace{5em} L(\delta u, \delta\lambda) = {\large\int}_\Gamma \delta\lambda\cdot\widetilde{u} d\Gamma

I am using @cdaversin 's mixed-dimensional branch.

Here is a simple MWE that fails:

```
from dolfin import *
mesh = UnitSquareMesh(100, 100)
left = CompiledSubDomain("near(x[0], 0.) && on_boundary")
top = CompiledSubDomain("near(x[1], 1.) && on_boundary")
bottom = CompiledSubDomain("near(x[1], 0.) && on_boundary")
right = CompiledSubDomain("near(x[0], 1.) && on_boundary")
origin = CompiledSubDomain("near(x[0], 0.) && near(x[1], 0.) && on_boundary")
facetMarks= MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
funcMarks = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
E, nu = Constant(2.), Constant(0.25)
mu = E/2./(1+nu)
lmbda = 2*mu*nu/(1-2*nu)
eps = lambda x: sym(grad(x))
sigma = lambda x: lmbda*tr(eps(x)) * Identity(len(x)) + 2*mu*eps(x)
facetMarks.set_all(0)
funcMarks.set_all(0)
bottom.mark(facetMarks, 1)
top.mark(facetMarks, 2)
left.mark(facetMarks, 3)
right.mark(facetMarks, 4)
bottom.mark(funcMarks, 1)
top.mark(funcMarks, 1)
left.mark(funcMarks, 2)
right.mark(funcMarks, 2)
mesh_B = MeshView.create(facetMarks, 1)
mesh_T = MeshView.create(facetMarks, 2)
mesh_L = MeshView.create(facetMarks, 3)
mesh_R = MeshView.create(facetMarks, 4)
mesh_LR = MeshView.create(funcMarks, 2)
mesh_TB = MeshView.create(funcMarks, 1)
Wu = VectorElement("CG", mesh.ufl_cell(), 2)
Wlam_LR = VectorElement("R", mesh_LR.ufl_cell(), 0)
Wlam_TB = VectorElement("R", mesh_TB.ufl_cell(), 0)
V_u = FunctionSpace(mesh, Wu)
V_lam_LR = FunctionSpace(mesh_LR, Wlam_LR)
V_lam_TB = FunctionSpace(mesh_TB, Wlam_TB)
V = MixedFunctionSpace(V_u, V_lam_LR, V_lam_TB)
u, lam_R, lam_T = TrialFunctions(V)
uh, lam_Rh, lam_Th = TestFunctions(V)
# volume measures
ds_R = Measure("dx", domain=mesh_R)
ds_L = Measure("dx", domain=mesh_L)
ds_T = Measure("dx", domain=mesh_T)
ds_B = Measure("dx", domain=mesh_B)
# dirichelt BC's
bcs = [DirichletBC(V.sub_space(0), Constant((0,0)), origin, method="pointwise"),
DirichletBC(V.sub_space(0).sub(1), Constant(0.0), CompiledSubDomain("near(x[0], 1.) && near(x[1], 0.0) && on_boundary"), method="pointwise"),
DirichletBC(V.sub_space(0).sub(0), Constant(0.0), CompiledSubDomain("near(x[0], 0.) && near(x[1], 1.0) && on_boundary"), method="pointwise") ]
a = inner(sigma(u), eps(uh))*dx
a += inner(lam_Rh, u)* ds_R - inner(lam_Rh, u)*ds_L + inner(lam_R, uh)* ds_R - inner(lam_R, uh)* ds_L
a += inner(lam_Th, u)* ds_T - inner(lam_Th, u)* ds_B + inner(lam_T, uh)* ds_T - inner(lam_T, uh)* ds_B
L = inner(lam_Rh, Constant(("0.1", "0"))) * ds_R + inner(lam_Th, Constant(("0.", "-0.1")))* ds_T
usol = Function(V, name="u")
solve(a == L, usol, bcs, solver_parameters={"linear_solver": "lu"})
```

I get the following error:

```
PetscDolfinErrorHandler: line '503', function 'MatSetValues_SeqAIJ', file '/tmp/petsc-src/src/mat/impls/aij/seq/aij.c',
: error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------
New nonzero at (80782,0) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
```