Applying loading on interior facets in L-form

variational formulation
1

I am facing difficulty applying traction loading on interior facets in the L-form of the equation. I get the error “Discontinuous type Argument must be restricted”. From other users, I have come to this understanding that I need to use some kind of discontinous functionspace and then interpolate to apply loading on internal facets. Is my understanding correct.

Here is my code for reference.

cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {19};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;


import ufl

gdim = 2

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure
#rho = 7850
#g = 9.81
T = fem.Constant(domain, 1000000.0)

# Discontinous loading for Internal Facet
#VV = fem.functionspace(domain, ('CG', 1))
#Td = interpolate(Expression('std::max(0., x[0]-0.5)', degree=1), VV)

#fg = fem.Constant(domain, np.array([0, -rho*g]))
#fg = fem.Constant(domain, 78500.0)

#Self-weight on the surface
n = FacetNormal(domain)
#nt= FacetTangent(domain)
#L_form = dot(-fp*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_markers)

#L_form = dot(T*n,u_) * ds(3)
#L_form = dot(-T*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)
L_form = dot(T*n,u_) * dS(7)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]



#------------------------------------------------------------------------------------------------------------------------------------------
#periodic_boundary = facet_markers.find(2)

#Cyclic Symmetry Condition
#def periodic_relation(x):
 #   out_x = np.zeros(x.shape)
  #  out_x[0] = 1 - x[0]
   # out_x[1] = x[1]
   # out_x[2] = x[2]
   # return out_x

# Periodic boundary - topological
#facets = mesh.locate_entities_boundary(domain, gdim - 1, periodic_boundary)
#arg_sort = np.argsort(facet_markers)
#mt = mesh.meshtags(domain, gdim - 1, facet_markers[arg_sort], np.full(len(facet_markers), 2, dtype=np.int32))

#with common.Timer("~~Periodic: Compute mpc condition"):
 #   mpc = dolfinx_mpc.MultiPointConstraint(V)
  #  mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, bcs)
   # mpc.finalize()

#-------------------------------------------------------------------------------------------------------------------------------------------                      

#WPsolve = LinearProblem(a_form, L_form, mpc, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
#WPsolve.solve()

WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()


V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)

import pyvista
from dolfinx import plot

pyvista.set_jupyter_backend("static")

topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()

vtk = io.VTKFile(domain.comm, "wpex5.pvd", "u")
vtk.write_function(u, 0)
vtk.write_function(sig, 0)
vtk.close()
2

Please note that your example is not reproducible, and it is unclear where in this quite extensive code the loading condition is, and what the error is.

First of all, try making an example illustrating your problem (i.e. adding an internal load) on a built in mesh such as a unit cube or square.

Secondly, ensure that the code you supply can reproduce the error message when executed by others.

Finally, please remove all unrelated code to the loading. If the error with respect to loading happens in Line 53, we only need the first 53 lines of your code. Then you can also remove all unused variables in those 53 lines of code that isn’t needed to reproduce the error message.

3

Please see below the clean code below for your reference.

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx import default_scalar_type

from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

This is mesh

cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {19};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;

Here is the code:

# Mesh
mesh_path = "Mesh_fenics.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)

# Check for Measures

import pyvista
from dolfinx import plot

pyvista.set_jupyter_backend("static")

topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()

from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile

with XDMFFile(MPI.COMM_WORLD, "mtgaurav.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(cell_markers, domain.geometry)

import ufl

gdim = 2

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure

T = fem.Constant(domain, 1000000.0)

n = FacetNormal(domain)


ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_markers)


L_form = dot(T*n,u_) * dS(7)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]

WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()


V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)

import pyvista
from dolfinx import plot

pyvista.set_jupyter_backend("static")

topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()

vtk = io.VTKFile(domain.comm, "wpex5.pvd", "u")
vtk.write_function(u, 0)
vtk.write_function(sig, 0)
vtk.close()

Here is the error outcome

--------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[42], line 68
     64 #bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
     65 #bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
     66 bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]
---> 68 WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
     69 WPsolve.solve()
     72 V0 = fem.functionspace(domain, ("DG", 0, (3,)))

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/petsc.py:590, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
    588 self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
    589 self._A = create_matrix(self._a)
--> 590 self._L = _create_form(L, form_compiler_options=form_compiler_options, jit_options=jit_options)
    591 self._b = create_vector(self._L)
    593 if u is None:
    594     # Extract function space from TrialFunction (which is at the
    595     # end of the argument list as it is numbered as 1, while the
    596     # Test function is numbered as 0)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
    185         return list(map(lambda sub_form: _create_form(sub_form), form))
    186     return form
--> 188 return _create_form(form)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
    180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    181 return form argument"""
    182 if isinstance(form, ufl.Form):
--> 183     return _form(form)
    184 elif isinstance(form, collections.abc.Iterable):
    185     return list(map(lambda sub_form: _create_form(sub_form), form))

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:141, in form.<locals>._form(form)
    139 if mesh is None:
    140     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 141 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
    142                                        form_compiler_options=form_compiler_options,
    143                                        jit_options=jit_options)
    145 # For each argument in form extract its function space
    146 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
     51 @functools.wraps(local_jit)
     52 def mpi_jit(comm, *args, **kwargs):
     53 
     54     # Just call JIT compiler when running in serial
     55     if comm.size == 1:
---> 56         return local_jit(*args, **kwargs)
     58     # Default status (0 == ok, 1 == fail)
     59     status = 0

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    202 # Switch on type and compile, returning cffi object
    203 if isinstance(ufl_object, ufl.Form):
--> 204     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    205 elif isinstance(ufl_object, ufl.FiniteElementBase):
    206     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    197     except Exception:
    198         pass
--> 199     raise e
    201 obj, module = _load_objects(cache_dir, module_name, form_names)
    202 return obj, module, (decl, impl)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    187     for name in form_names:
    188         decl += form_template.format(name=name)
--> 190     impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
    191                             cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    192 except Exception as e:
    193     try:
    194         # remove c file so that it will not timeout next time

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    256 import ffcx.compiler
    258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
    259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
    262 ffibuilder = cffi.FFI()
    263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
    264                       extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/compiler.py:97, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
     95 # Stage 1: analysis
     96 cpu_time = time()
---> 97 analysis = analyze_ufl_objects(ufl_objects, options)
     98 _print_timing(1, time() - cpu_time)
    100 # Stage 2: intermediate representation

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:88, in analyze_ufl_objects(ufl_objects, options)
     85     else:
     86         raise TypeError("UFL objects not recognised.")
---> 88 form_data = tuple(_analyze_form(form, options) for form in forms)
     89 for data in form_data:
     90     elements += [convert_element(e) for e in data.unique_sub_elements]

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:88, in <genexpr>(.0)
     85     else:
     86         raise TypeError("UFL objects not recognised.")
---> 88 form_data = tuple(_analyze_form(form, options) for form in forms)
     89 for data in form_data:
     90     elements += [convert_element(e) for e in data.unique_sub_elements]

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:163, in _analyze_form(form, options)
    160 complex_mode = "_Complex" in options["scalar_type"]
    162 # Compute form metadata
--> 163 form_data = ufl.algorithms.compute_form_data(
    164     form,
    165     do_apply_function_pullbacks=True,
    166     do_apply_integral_scaling=True,
    167     do_apply_geometry_lowering=True,
    168     preserve_geometry_types=(ufl.classes.Jacobian,),
    169     do_apply_restrictions=True,
    170     do_append_everywhere_integrals=False,  # do not add dx integrals to dx(i) in UFL
    171     complex_mode=complex_mode)
    173 # If form contains a quadrature element, use the custom quadrature scheme
    174 custom_q = None

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py:328, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
    326 # Propagate restrictions to terminals
    327 if do_apply_restrictions:
--> 328     form = apply_restrictions(form)
    330 # If in real mode, remove any complex nodes introduced during form processing.
    331 if not complex_mode:

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:180, in apply_restrictions(expression)
    177 integral_types = [k for k in integral_type_to_measure_name.keys()
    178                   if k.startswith("interior_facet")]
    179 rules = RestrictionPropagator()
--> 180 return map_integrand_dags(rules, expression,
    181                           only_integral_type=integral_types)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:73, in map_integrand_dags(function, form, only_integral_type, compress)
     71 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
     72     """Map integrand dags."""
---> 73     return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
     74                           form, only_integral_type)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:25, in map_integrands(function, form, only_integral_type)
     23 """Apply transform(expression) to each integrand expression in form, or to form if it is an Expr."""
     24 if isinstance(form, Form):
---> 25     mapped_integrals = [map_integrands(function, itg, only_integral_type)
     26                         for itg in form.integrals()]
     27     nonzero_integrals = [itg for itg in mapped_integrals
     28                          if not isinstance(itg.integrand(), Zero)]
     29     return Form(nonzero_integrals)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:25, in <listcomp>(.0)
     23 """Apply transform(expression) to each integrand expression in form, or to form if it is an Expr."""
     24 if isinstance(form, Form):
---> 25     mapped_integrals = [map_integrands(function, itg, only_integral_type)
     26                         for itg in form.integrals()]
     27     nonzero_integrals = [itg for itg in mapped_integrals
     28                          if not isinstance(itg.integrand(), Zero)]
     29     return Form(nonzero_integrals)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:33, in map_integrands(function, form, only_integral_type)
     31 itg = form
     32 if (only_integral_type is None) or (itg.integral_type() in only_integral_type):
---> 33     return itg.reconstruct(function(itg.integrand()))
     34 else:
     35     return itg

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:73, in map_integrand_dags.<locals>.<lambda>(expr)
     71 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
     72     """Map integrand dags."""
---> 73     return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
     74                           form, only_integral_type)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/map_dag.py:34, in map_expr_dag(function, expression, compress, vcache, rcache)
     15 def map_expr_dag(function, expression,  compress=True, vcache=None, rcache=None):
     16     """Apply a function to each subexpression node in an expression DAG.
     17 
     18     If the same funtion is called multiple times in a transformation
   (...)
     32         The result of the final function call
     33     """
---> 34     result, = map_expr_dags(function, [expression], compress=compress,
     35                             vcache=vcache,
     36                             rcache=rcache)
     37     return result

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/map_dag.py:98, in map_expr_dags(function, expressions, compress, vcache, rcache)
     96 # Cache miss: Get transformed operands, then apply transformation
     97 if cutoff_types[v._ufl_typecode_]:
---> 98     r = handlers[v._ufl_typecode_](v)
     99 else:
    100     r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:104, in RestrictionPropagator.reference_value(self, o)
    102 f, = o.ufl_operands
    103 assert f._ufl_is_terminal_
--> 104 g = self(f)
    105 if isinstance(g, Restricted):
    106     side = g.side()

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/multifunction.py:95, in MultiFunction.__call__(self, o, *args)
     93 def __call__(self, o, *args):
     94     """Delegate to handler function based on typecode of first argument."""
---> 95     return self._handlers[o._ufl_typecode_](o, *args)

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:59, in RestrictionPropagator._require_restriction(self, o)
     57 """Restrict a discontinuous quantity to current side, require a side to be set."""
     58 if self.current_restriction is None:
---> 59     raise ValueError(f"Discontinuous type {o._ufl_class_.__name__} must be restricted.")
     60 return o(self.current_restriction)

ValueError: Discontinuous type Argument must be restricted.
4

Here is the formatted code for your reference.

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx import default_scalar_type

from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {19};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;

# Mesh
mesh_path = "Mesh_fenics.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)


import ufl

gdim = 2

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure

T = fem.Constant(domain, 1000000.0)

n = FacetNormal(domain)


ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_markers)


L_form = dot(T*n,u_) * dS(7)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)

bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]


WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()


V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)
5

Dear @Gaurav_Verma,
Your code is far from cleaned up and is missing the error message.
I will give you the benefit of the doubt and illustrate to you how an MWE with error trace should look like:



from mpi4py import MPI
from dolfinx import fem
import dolfinx.fem.petsc
from dolfinx.io import gmshio

from ufl import (FacetNormal, as_matrix, TestFunction, TrialFunction, as_vector, dot, dx, inner, grad, sym)
import ufl

mesh_path = "Mesh_fenics.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)



gdim = 2

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure

T = fem.Constant(domain, 1000000.0)

n = FacetNormal(domain)


ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_markers)


L_form = dot(T*n,u_) * dS(7)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)

bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]


WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)

which yields:

Traceback (most recent call last):
  File "/root/shared/mwe_french.py", line 81, in <module>
    WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py", line 759, in __init__
    self._L = _create_form(
              ^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py", line 307, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py", line 302, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py", line 240, in _form
    ufcx_form, module, code = jit.ffcx_jit(
                              ^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py", line 62, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py", line 212, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms
    raise e
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 94, in <genexpr>
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 180, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py", line 337, in compute_form_data
    form = apply_restrictions(form)
           ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_restrictions.py", line 183, in apply_restrictions
    return map_integrand_dags(rules, expression, only_integral_type=integral_types)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags
    return map_integrands(
           ^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 30, in map_integrands
    map_integrands(function, itg, only_integral_type) for itg in form.integrals()
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
                           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>
    lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
                ^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 101, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_restrictions.py", line 106, in reference_value
    g = self(f)
        ^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/multifunction.py", line 95, in __call__
    return self._handlers[o._ufl_typecode_](o, *args)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_restrictions.py", line 61, in _require_restriction
    raise ValueError(f"Discontinuous type {o._ufl_class_.__name__} must be restricted.")
ValueError: Discontinuous type Argument must be restricted.
Exception ignored in: <function LinearProblem.__del__ at 0x73ec59653a60>
Traceback (most recent call last):
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py", line 801, in __del__
    self._solver.destroy()
    ^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'

Then one can inspect the message, seeing that it states

ValueError: Discontinuous type Argument must be restricted.

which takes me to the following lines:

You are here using the dS measure which is an interior facet measure. However, looking at your facet markers, the boundary marked with 7 is an exterior boundary, so why would you need it?
Secondly, if you want to apply a load on an interior facet, you need to decide in what direction the facet normal should point. How to do this is for instance illustrated in:

6

Dear Joergen,

Yes you were right. I had marked the exterior facet twice in the physical group. Now, I have made the correction.

I am now getting a solution when I try provide direction to the facet normal. Pls check if this is correct. I only made this correction in the L_form.

L_form = dot(T * n("+"),u_("+")) * dS(5)
7

This is not necessarily correct, as I tried to emphasize with