Setting Dirichlet (essential) boundary conditions

Dirichlet boundary conditions in cpppetsc are set up in three steps: tag the relevant mesh points with a label value, describe which dofs on each label value are constrained and to what value, then resolve everything into the locally-owned (global dof index, value) pairs the solver expects. This mirrors how examples/ForceElement.cc sets up the boundary conditions for its example problem.

1. Tag boundary points with `LabelConstraints`

A LabelConstraints pairs a predicate on a mesh point's coordinates with a label value:

struct LabelConstraints {
  std::function<bool(std::span<const PetscReal>)> condition;
  PetscInt label_value{-1};
};

For every mesh point whose coordinates satisfy condition, label_value is written into the mesh's "Face Sets" DMLabel.

examples/ForceElement.cc

  constexpr auto left_bc_tag = 5;
  constexpr auto lower_bc_tag = 6;
  const auto functional_bc = std::vector<LabelConstraints<dimension>>{
      {.condition =
           [](std::span<const PetscReal> coords) {
             const auto tolerance = .1;
             return std::abs(coords[0] - 0.) < tolerance; // x == 0
           },
       .label_value = left_bc_tag},
      {.condition =
           [](std::span<const PetscReal> coords) {
             const auto tolerance = .1;
             return std::abs(coords[1] - 0.) < tolerance; // y == 0
           },
       .label_value = lower_bc_tag},
  };

  mesh.markBoundaryCells(functional_bc);

markBoundaryCells must be called before extracting boundary dofs (step 3).

If you import the mesh via Mesh::fromGmsh(), the "Face Sets" label is already populated from the gmsh Physical Groups, so you can skip this step and reuse the gmsh physical-group ids as label values.

2. Describe constrained dofs with `EssentialBoundaryCondition`

An EssentialBoundaryCondition says, for a given label_value, which local dof indices (e.g. 0/1/2 for x/y/z) are constrained and to what value:

struct EssentialBoundaryCondition {
  inline static std::string_view label_name{"Face Sets"};
  PetscInt label_value{-1};
  std::unordered_map<PetscInt, PetscScalar> constrained_dofs;
};

Continuing the example — fixing both x and y displacement (dofs 0 and 1) to 0.0 on both tagged edges:

const std::vector<EssentialBoundaryCondition> essential_bcs{
    {.label_value = left_bc_tag, .constrained_dofs{{0, 0.0}, {1, 0.0}}},
    {.label_value = lower_bc_tag, .constrained_dofs{{0, 0.0}, {1, 0.0}}},
};

3. Resolve to local dofs with `extractBoundaryDofs`

extractBoundaryDofs turns the label-based conditions into a std::vector<ConstrainedDOF> ({PetscInt dof; PetscScalar value;}) — the concrete global dof indices, restricted to the ones owned by the calling rank:

const auto boundary_conditions =
    extractBoundaryDofs<Mesh, dof_per_vertex>(mesh, essential_bcs);

This is MPI-aware: it gathers the resolved dofs across all ranks and keeps only those in the calling rank's ownership range, so it produces correct results on distributed meshes without any extra work.

4. Pass the result to the solver

solver.computeSolution() takes the resolved boundary conditions directly, alongside an initial guess vector:

auto initial_guess = Vector::fromGlobalMesh(mesh);
auto result = solver.computeSolution(
    boundary_conditions, std::move(initial_guess), time, &assembler);

Minimal end-to-end example

constexpr int left_bc_tag = 5;
constexpr int lower_bc_tag = 6;

const std::vector<LabelConstraints<dimension>> functional_bc{
    {.condition = [](std::span<const PetscReal> c) {
       return std::abs(c[0] - 0.) < .1; // x == 0
     },
     .label_value = left_bc_tag},
    {.condition = [](std::span<const PetscReal> c) {
       return std::abs(c[1] - 0.) < .1; // y == 0
     },
     .label_value = lower_bc_tag},
};
const std::vector<EssentialBoundaryCondition> essential_bcs{
    {.label_value = left_bc_tag, .constrained_dofs{{0, 0.0}, {1, 0.0}}},
    {.label_value = lower_bc_tag, .constrained_dofs{{0, 0.0}, {1, 0.0}}},
};

mesh.markBoundaryCells(functional_bc);
const auto boundary_conditions =
    extractBoundaryDofs<Mesh, dof_per_vertex>(mesh, essential_bcs);

auto initial_guess = Vector::fromGlobalMesh(mesh);
auto result = solver.computeSolution(
    boundary_conditions, std::move(initial_guess), time, &assembler);

Beyond single-point constraints: `GeneralizedVertexBoundaryCondition`

EssentialBoundaryCondition can only fix a dof to a constant. If a constraint instead needs to combine multiple source vertices (e.g. periodic or linear-combination boundary conditions, or coupling one vertex's dof to another's), use GeneralizedVertexBoundaryCondition and its own extractBoundaryDofs overload from ae108/cpppetsc/GeneralizedVertexBoundaryCondition.h instead.

Rather than a label value, each condition targets a vertex directly by its global_vertex_id, and constrains it per dof component to either a plain value (empty sources) or a linear combination of other vertices' dofs plus an offset — value = Σ(factor · source_dof) + value_offset:

struct GeneralizedVertexConstraint {
  struct SourceItem {
    PetscScalar factor{0.};
    PetscInt global_vertex_id{-1};
    PetscInt dof_component{-1};
  };
  PetscInt dof_component{-1};
  std::vector<SourceItem> sources;
  PetscScalar value_offset{0.};
};

struct GeneralizedVertexBoundaryCondition {
  PetscInt global_vertex_id{-1};
  std::unordered_map<PetscInt, GeneralizedVertexConstraint> constrained_dofs;
};

Minimal example, adapted from examples/GeneralizedNonlinearSolver.cc: fix vertex 0 fully, prescribe a displacement of 0.5 on dof 0 of vertex 1, and couple vertex 2's dofs to vertex 3's dofs plus an offset:

const std::vector<GeneralizedVertexBoundaryCondition> essential_bcs{
    {.global_vertex_id = 0,
     .constrained_dofs = {{0, {.sources = {}, .value_offset = 0.0}},
                          {1, {.sources = {}, .value_offset = 0.0}}}},
    {.global_vertex_id = 1,
     .constrained_dofs = {{0, {.sources = {}, .value_offset = 0.5}},
                          {1, {.sources = {}, .value_offset = 0.0}}}},
    {.global_vertex_id = 2,
     .constrained_dofs =
         {{0,
           {.sources = {{.factor = 1., .global_vertex_id = 3,
                         .dof_component = 0}},
            .value_offset = 0.5}},
          {1,
           {.sources = {{.factor = 1., .global_vertex_id = 3,
                         .dof_component = 1}},
            .value_offset = 0.0}}}},
};

const auto constrained_dofs =
    extractBoundaryDofs<Mesh, dof_per_vertex>(mesh, essential_bcs);

const auto result = solver.computeSolution(
    constrained_dofs, Vector::fromGlobalMesh(mesh), time, &assembler);

There is no markBoundaryCells() step here since vertices are addressed directly rather than via a label. Note that global_vertex_id refers to the mesh's current (post-distribution) DMPlex vertex numbering, which can differ from the "natural" numbering once a mesh has been redistributed across ranks; use naturalToGlobalVertexIds(mesh) to translate natural vertex ids into that numbering before building essential_bcs if needed.

Setting Dirichlet (essential) boundary conditions

1. Tag boundary points with LabelConstraints

2. Describe constrained dofs with EssentialBoundaryCondition

3. Resolve to local dofs with extractBoundaryDofs