DirichletBC on part of a FEniCS boundary: facet markers and the whole-boundary mistake

# DirichletBC on Part of a FEniCS Boundary: Facet Markers and the Whole-Boundary Mistake

Scope: legacy FEniCS (DOLFIN, `from fenics import *`). In DOLFINx the API is
different (`fem.locate_dofs_topological` + `fem.dirichletbc`); do not mix the two.

## The concrete problem

You want `u = VALUE` on only *part* of the boundary (e.g. the inlet of a channel,
the left edge of a square) and different (or no) conditions elsewhere. The classic
mistake is marking the whole boundary and pinning every boundary DOF:

```python
# WRONG: this pins u on the ENTIRE boundary, not just the left edge
def whole_boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, Constant(1.0), whole_boundary)
```

The correct approach is to label boundary *facets* with integer markers in a
`MeshFunction` and pass the marker id to `DirichletBC`.

## The recipe

```python
from fenics import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "P", 1)

# 1. Define the boundary pieces as SubDomain objects
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

# 2. Build a facet MeshFunction: dimension = mesh dim - 1
#    (edges in 2D, faces in 3D; cell-dim markers are a common bug)
facet_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facet_markers.set_all(0)          # 0 = "unmarked / not a Dirichlet boundary"

# 3. Stamp the marker ids onto the facets
left = Left()
right = Right()
left.mark(facet_markers, 1)
right.mark(facet_markers, 2)

# 4. Construct one DirichletBC per marked piece, keyed by marker id
#    DirichletBC(V, g, sub_domains, sub_domain, method="topological")
bc_left  = DirichletBC(V, Constant(0.0), facet_markers, 1)
bc_right = DirichletBC(V, Constant(1.0), facet_markers, 2)
bcs = [bc_left, bc_right]

# 5. Use in a solve as usual
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
uh = Function(V)
solve(a == L, uh, bcs)
```

The compact alternative using `CompiledSubDomain` (tolerance-safe `near`):

```python
left  = CompiledSubDomain("near(x[0], 0.0) && on_boundary")
right = CompiledSubDomain("near(x[0], 1.0) && on_boundary")
```

Reuse the same `facet_markers` MeshFunction for Neumann terms on marked pieces:

```python
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
L = f*v*dx + g_in*v*ds(1)   # traction/heat flux through the left edge only
```

## Classic mistakes and how to catch them

1. **BC applied to the whole boundary.** `inside()` returning bare `on_boundary`
   (or `DirichletBC(V, g, mesh)`-style whole-mesh markers) pins every boundary DOF.
   Symptom: the solution is flat along edges you meant to leave free.
2. **Wrong entity dimension.** A `MeshFunction` over cells (`dim()`) marks volumes,
   not boundary pieces. The `DirichletBC(V, g, markers, id)` constructor needs a
   *facet* MeshFunction (`dim() - 1`). Using a cell mesh function silently marks
   nothing (or the wrong thing).
3. **Forgot to mark / wrong id.** The marker id passed to `DirichletBC` must match
   the id used in `.mark()`. Facets left at `0` are excluded, so a typo gives a BC
   on zero facets: no error, just a missing constraint.
4. **Overlapping marks.** If two `SubDomain`s overlap, the later `.mark()` wins on
   shared facets. Order your marks intentionally; do not double-mark a facet with
   two different ids.
5. **DOF identification method.** The default `"topological"` method only finds DOFs
   lying on a fully-marked facet; for discontinuous elements use
   `method="geometric"`. `DirichletBC(V, g, markers, id, method="geometric")`.
6. **Vector-valued components.** For a `VectorFunctionSpace`, constrain one
   component with `DirichletBC(W.sub(0), value, facet_markers, 1)`.

## Verify the BC covers what you think it covers

```python
# Count constrained DOFs before trusting the solution
for bc, name in [(bc_left, "left"), (bc_right, "right")]:
    n_dofs = len(bc.get_boundary_values())
    print(name, "constrained DOFs:", n_dofs)

# Sanity check: the value actually appears on the marked boundary
import numpy as np
coords = mesh.coordinates()
left_dofs = np.array([d for d in bc_left.get_boundary_values().keys()])
```

Rule of thumb: if `get_boundary_values()` is empty for a BC you defined, the
marker id has no marked facets; if it is far larger than expected, you marked the
whole boundary.

## Why this works (grounding)

`DirichletBC` accepts three ways to specify boundary indicators: a `SubDomain`
object (facets found on first `apply()`), a `MeshFunction` over facets plus an
integer id selecting which facets to include, or boundary data attached to the
mesh. The facet-marker form is the reliable one for multi-part boundaries because
each piece is labeled once, reused by both `DirichletBC` and `Measure("ds",
subdomain_data=...)`, and checked by id.