Applying Multiple Dirichlet & Boundary Conditions

General
1

I am trying to follow the approach suggested in this dolfinx tutorial, to set multiple Dirichlet and Neumann Boundary conditions.

I do fail to see how to apply a zero Neumann condition. I would have thought that having a null function as

u_ex

in the above demo would do the trick, but I get an error

after the line

f = -div(grad(u_ex(x)))

with the error message

Cannot determine geometric dimension from expression.
ERROR:UFL:Cannot determine geometric dimension from expression.

What is the problem with a zero function?

This is a minimum working example, I am using dolfinx 0.6.0

from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
from dolfinx.fem import (Constant,  Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import create_vtk_mesh

from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)

import numpy as np
import pyvista

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)


u_ex = lambda x: 0.0 + 0.0*x[0]**2 + 0.0*x[1]**2
x = SpatialCoordinate(mesh)
# Define physical parameters and boundary condtions
s = u_ex(x)
f = -div(grad(u_ex(x)))
n = FacetNormal(mesh)
g = -dot(n, grad(u_ex(x)))


kappa = Constant(mesh, default_scalar_type(1))
r = Constant(mesh, default_scalar_type(1000))
# Define function space and standard part of variational form
V = FunctionSpace(mesh, ("Lagrange", 1))
u, v = TrialFunction(V), TestFunction(V)
F = kappa * inner(grad(u), grad(v)) * dx - inner(f, v) * dx
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]
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)
class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
                       BoundaryCondition("Dirichlet", 2, u_ex),
                       BoundaryCondition("Robin", 3, (r, s)),
                       BoundaryCondition("Neumann", 4, g)]
2

If you set everything to zero, UFL simplifies it to be a single zero function, that does not contain the geometrical information you would like.
If you want to keep the current structure, you would have to adapt the u_ex function, to use

def u_ex(mod, x):
    if mod == np:
        zero = dolfinx.fem.Constant(mesh, 0.0)
    elif mod == dolfinx:
        zero = 0.0
    else:
        raise RuntimeError(f"Unsupported {mod=}")
    return zero + zero*x[0]**2 + zero*x[1]**2

s = u_ex(dolfinx, x)

and if you want to use the expression in an interpolation

s = Function(V)
s.interpolate(lambda x: u_ex(np, x))
3

@dokken thank for this much appreciated.
I am still not sure though I fully get To use the newly defined function I use the np argument, and I do not get anymore the previous error

f = -div(grad(u_ex(np, x)))
n = FacetNormal(domain)
g = -dot(n, grad(u_ex(np, x)))

Then I follow the tutorial and go

boundaries = [(1, lambda x: np.isclose(x[0], 0)), 
              (2, lambda x: np.isclose(x[0],  1)), 
               (3, lambda x: np.isclose(x[1],  0)),  
               (4, lambda x: np.isclose(x[1],  1)) 
              ]  

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

but the kernel crashes systematically (using a Jupyter notebook), am I missing anything?

I fail to see what your comment on the interpolation refers to, is the way to go to

s = Function(V)
s.interpolate(lambda x: u_ex(np, x))
f = -div(grad(s))
n = FacetNormal(domain)
g = -dot(n, grad(s))

I tried this one too, still got the kernel crashing

Also, I see in the Singular Poisson tutorial another mod argument is used, u_exact = u_ex(ufl, x), trying to make sense of it all

Thank you

4

You should send in dolfinx not numpy if you want to use it with ufl.SpatialCoordinate and ufl.div, ufl.grad etc.

5

I am sorry I am being likely very dense.

If I try

f = -div(grad(u_ex(dolfinx, x)))

I get the error

ERROR:UFL:Cannot determine geometric dimension from expression.
Cannot determine geometric dimension from expression.

If on the other hand

s = Function(V)
s.interpolate(lambda x: u_ex(dolfinx, x))
f = -div(grad(s))
n = FacetNormal(domain)
g = -dot(n, grad(s))

no such error, but then the kernel crashes after

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0],  1)),
               (3, lambda x: np.isclose(x[1],  0)),
                (4, lambda x: np.isclose(x[1],  1))
              ]


facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, 1),
                       BoundaryCondition("Dirichlet", 2, -1),
                       BoundaryCondition("Neumann", 3, g),
                       BoundaryCondition("Neumann", 4, g)]
6

Please provide a full reproducible example, rather than partial snippets of your code.

7

Yes, let me attach the exact example bugging me. I think I made some progress thanks to your hints, but it is still not solving with a type error

TypeError: assemble_matrix(): incompatible function arguments. The following argument types are supported:
    1. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], bcs: List[dolfinx::fem::DirichletBC<double>], unrolled: bool = False) -> None
    2. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], rows0: numpy.ndarray[numpy.int8], rows1: numpy.ndarray[numpy.int8], unrolled: bool = False) -> None

### Example


import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
import dolfinx 
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import create_vtk_mesh
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
import pyvista
from mpi4py import MPI
#from petsc4py import PETSc
pyvista.set_jupyter_backend("trame")

print (dolfinx.__version__)

## GMSH meshing

import gmsh
gmsh.initialize()
gmsh.model.add("Cracked_plate_approx")
## SQ stands for size quarter
SQ = 1
eps = 0.01

p1 = gmsh.model.occ.addPoint(0, 0, 0, 0, 1)
p2 = gmsh.model.occ.addPoint(0+eps, SQ, 0, 0, 2)
p3 = gmsh.model.occ.addPoint(SQ, SQ, 0, 0, 3)
p4 = gmsh.model.occ.addPoint(SQ, - SQ, 0, 0, 4)
p5 = gmsh.model.occ.addPoint(-SQ, - SQ, 0, 0, 5)
p6 = gmsh.model.occ.addPoint(-SQ, SQ, 0, 0, 6)

p7 = gmsh.model.occ.addPoint(0-eps, SQ, 0, 0, 7)


# Creating connection lines
l1 = gmsh.model.occ.addLine(1, 2, 10)
l2 = gmsh.model.occ.addLine(2, 3, 11)
l3 = gmsh.model.occ.addLine(3, 4, 12)
l4 = gmsh.model.occ.addLine(4, 5, 13)

l5 = gmsh.model.occ.addLine(5, 6, 14)
l6 = gmsh.model.occ.addLine(6, 7, 15)
l7 = gmsh.model.occ.addLine(7, 1, 16)



loop_outer = gmsh.model.occ.addCurveLoop([l1,l2,l3,l4, l5, l6, l7], 30)
#loop_circle = gmsh.model.occ.addCurveLoop([l5,l6,l7,l8], 19)
gmsh.model.occ.addPlaneSurface([30], 100)

gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(2, [100], 400)

gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.05)
gmsh.model.mesh.generate(gdim)

from dolfinx.io import gmshio
from mpi4py import MPI

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

# Plotting mesh
ndim = domain.geometry.dim

pyvista.set_jupyter_backend('trame')
#pyvista.OFF_SCREEN 
topology, cell_types, geometry = dolfinx.plot.create_vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, domain.geometry.x)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
#plotter.set_scale(5,5)
#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
if not pyvista.OFF_SCREEN:
    plotter.show()

## Variational Problem

V = FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
# Original problem
a = dot(grad(u), grad(v)) * dx

u_dirich =  lambda x:  x[0]


def u_ex(mod, x):
    if mod == np:
        zero = dolfinx.fem.Constant(domain, 0.0)
    elif mod == dolfinx:
        zero = 0.0
    else:
        raise RuntimeError(f"Unsupported {mod=}")
    return zero + zero*x[0]**2 + zero*x[1]**2
x = SpatialCoordinate(domain)


s = Function(V)
s.interpolate(lambda x: u_ex(dolfinx, x))

n = FacetNormal(domain)
g = -dot(n, grad(s))

## Neumann boundaries on the crack faces
## Dirichlet boundarues on the left and right edges

boundaries = [(1, lambda x: (np.isclose(x[1], x[0]*SQ/eps))&(x[0] > 0)),
              (2, lambda x: (np.isclose(x[1],  -x[0]*SQ/eps))&(x[0] <0)),
               (3, lambda x: np.isclose(x[0],  -SQ)),
                (4, lambda x: np.isclose(x[0],  SQ))
              ]


facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
bcs = [BoundaryCondition("Dirichlet", 1, u_dirich),
                       BoundaryCondition("Dirichlet", 2, u_dirich),
                       BoundaryCondition("Neumann", 3, g),
                       BoundaryCondition("Neumann", 4, g)]

source = 6.
f = Constant(domain, default_scalar_type(source))
L = f * v * dx 

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
8

Ok, I have been working on it and made some progress I believe. It solces, yet, I get only nans in the solution making it ever so had to debug, please find.a woking MWE below, thanks

import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
import dolfinx 
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import create_vtk_mesh
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
import pyvista
from mpi4py import MPI
#from petsc4py import PETSc
pyvista.set_jupyter_backend("trame")
print (dolfinx.__version__)
# import gmsh
# gmsh.initialize()
# gmsh.model.add("Cracked_plate_approx")
domain = create_unit_square(MPI.COMM_WORLD, 10, 10)

ndim = domain.geometry.dim

pyvista.set_jupyter_backend('trame')
#pyvista.OFF_SCREEN 
topology, cell_types, geometry = dolfinx.plot.create_vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, domain.geometry.x)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
#plotter.set_scale(5,5)
#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
if not pyvista.OFF_SCREEN:
    plotter.show()

#from pathlib import Path
#Path("output").mkdir(parents=True, exist_ok=True)
#figure = plotter.screenshot("output/mesh.png")
## Variational Problem

V = FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
# Original problem
#
source = 0.0
f = Constant(domain, default_scalar_type(source))
#L = f * v * dx 

F =  inner(grad(u), grad(v)) * dx - f*inner(f, v) * dx

u_dirich =  lambda x:  0.0+1.0*x[0]+0.0*x[1] # applying constant displacement on the outer vertical faces

def u_ex(mod, x):
    if mod == np:
        zero = dolfinx.fem.Constant(domain, 0.0)
    elif mod == dolfinx:
        zero = 0.0
    else:
        raise RuntimeError(f"Unsupported {mod=}")
    return zero + zero*x[0]**2 + zero*x[1]**2
x = SpatialCoordinate(domain)


s = Function(V)
s.interpolate(lambda x: u_ex(dolfinx, x))

n = FacetNormal(domain)
g = -dot(n, grad(s))

## Neumann boundaries on the top/botim faces
## Dirichlet boundary conds on the left and right faces

boundaries = [(1, lambda x: (np.isclose(x[1], 0.0))),
                (2, lambda x: (np.isclose(x[1],  1.0))),
                (3, lambda x: np.isclose(x[0],  0.)),
                    (4, lambda x: np.isclose(x[0],  1.0))
                ]



facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
bcs = [BoundaryCondition("Neumann", 1, g),
        BoundaryCondition("Neumann", 2, g),
        BoundaryCondition("Dirichlet", 3, u_dirich),
        BoundaryCondition("Dirichlet", 4, u_dirich)]
bcs = []
for condition in bcs:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc



# #### SINGLE DIRICHLET BOUNDARY COND
# ### JUST TO CHECK IT WORKS IN THE BASIC SETTING

# from dolfinx import fem
# uD = fem.Function(V)
# #uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# uD.interpolate(lambda x: 0.0 + 1.0*x[0])



# import numpy as np
# # # Create facet to cell connectivity required to determine boundary facets
# tdim = domain.topology.dim
# fdim = tdim - 1

# dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))
# u_L = Function(V)
# u_L.interpolate(lambda x: 0.0 + 1.0 * x[0])
# bc_L = dirichletbc(u_L, dofs_L)

# dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1.0))
# u_R = Function(V)
# u_R.interpolate(lambda x: 0.0 + 1.0 * x[0])
# bc_R = dirichletbc(u_R, dofs_R)


# bcs_dirich_only = [ bc_R, bc_L]

## DIRICH_ONLY: BASIC CHECK
#problem = LinearProblem(a, L, bcs=bcs_dirich_only, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

### DIRICHLET PLUS NEUMANN BCs defined as in the tutorial
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uh = problem.solve()
#pyvista.start_xvfb()

pyvista_cells, cell_types, geometry = create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")

plotter = pyvista.Plotter()
#plotter.set_viewup([2,0,0])
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=False,cmap= 'bwr')

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("neumann_dirichlet.png")
9

There is so much noise in this code. Please remove all plotting and all code that is commented out.
Currently, you are defining the variable bcs multiple times, making it hard to know what you really want to use

10

Indeed a silly naming error on my end, for which I apologise, was hiding in the festering mess.

It seems to solve now fine

import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
import dolfinx 
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import create_vtk_mesh
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)

from mpi4py import MPI
print (dolfinx.__version__)
domain = create_unit_square(MPI.COMM_WORLD, 10, 10)

ndim = domain.geometry.dim

V = FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
source = 0.0
f = Constant(domain, default_scalar_type(source))
F =  inner(grad(u), grad(v)) * dx - f*inner(f, v) * dx

u_dirich =  lambda x:  0.0+1.0*x[0]+0.0*x[1] # applying constant displacement on the outer vertical faces

def u_ex(mod, x):
    if mod == np:
        zero = dolfinx.fem.Constant(domain, 0.0)
    elif mod == dolfinx:
        zero = 0.0
    else:
        raise RuntimeError(f"Unsupported {mod=}")
    return zero + zero*x[0]**2 + zero*x[1]**2
x = SpatialCoordinate(domain)


s = Function(V)
s.interpolate(lambda x: u_ex(dolfinx, x))

n = FacetNormal(domain)
g = -dot(n, grad(s))

## Neumann boundary conds on the top/bottom faces
## Dirichlet boundary conds on the left and right faces

boundaries = [(1, lambda x: (np.isclose(x[1], 0.0))),
                (2, lambda x: (np.isclose(x[1],  1.0))),
                (3, lambda x: np.isclose(x[0],  0.)),
                    (4, lambda x: np.isclose(x[0],  1.0))
                ]



facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Neumann", 1, g),
        BoundaryCondition("Neumann", 2, g),
        BoundaryCondition("Dirichlet", 3, u_dirich),
        BoundaryCondition("Dirichlet", 4, u_dirich)]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc


a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uh = problem.solve()