Overview

LehrFEM++ provides a number of convenient functions for working with essential boundary conditions.

To fix degrees of freedom (DoFs) on the boundary, we can use the functions lf::assemble::FixFlaggedSolutionComponents or lf::assemble::FixFlaggedSolutionCompAlt (see also Quick Reference - Assembly). Both functions take a function as an argument that returns a pair of a boolean and a double: std::pair<bool, double> selector(unsigned int dof_idx). The selector function returns whether the DoF is to be fixed and, if so, the value it should be fixed to.

We can obtain the boundary flags as a MeshDataSet using:

lf::mesh::utils::CodimMeshDataSet<bool> bd_flags =

lf::mesh::utils::flagEntitiesOnBoundary(dofh.Mesh(), 2);

The second argument is the co-dimension of the entities we are interested in. To fix all DoFs on the boundary (including those associated with edges and cells), omit the second argument.

One function on the whole boundary

Let \(g\) be a function defining the values on the boundary. To get the function values of \(g\) at the boundary nodes, wrap it in a MeshFunction:

auto mf_g = lf::mesh::utils::MeshFunctionGlobal(g);
 
std::vector<std::pair<bool, double>> boundary_val =
    lf::fe::InitEssentialConditionFromFunction(*fe_space, bd_flags, mf_g);

This returns a std::vector<std::pair<bool, scalar_t>>. To create our selector:

auto selector = [&](unsigned int dof_idx) -> std::pair<bool, double> {
  return boundary_val[dof_idx];
};

If \(g\) has different definitions on different parts of the boundary (e.g., \(g = 0\) on \(\gamma_0\) and \(g = 1\) on \(\gamma_1\)), try to express \(g\) as a lambda function with an if-else statement.

Constant value on the boundary

The selector becomes:

auto selector = [&](unsigned int dof_idx) -> std::pair<bool, double> {
  if (bd_flags(dofh.Entity(dof_idx))) {
    return std::make_pair(true, 1.0);  // value 1.0 on the boundary
  } else {
    return std::make_pair(false, 0.0);  // value irrelevant
  }
};

BC only on part of the boundary

We need to create our own bd_flags. Initialize a MeshDataSet with default value false:

lf::mesh::utils::AllCodimMeshDataSet<bool> bd_flags_part(mesh_p, false);

Then loop over nodes and edges separately and set bd_flags to true for the nodes/edges where the boundary condition should be applied.

for (const auto& edge : fe_space->Mesh()->Entities(1)) {
  if (condition) bd_flags_part(*edge) = true;
}
 
for (const auto& node : fe_space->Mesh()->Entities(2)) {
  //...
}

Example (Essential boundary conditions)

auto mesh_factory = std::make_unique<mesh::hybrid2d::MeshFactory>(2);
auto gmsh_reader = io::GmshReader(std::move(mesh_factory), "mesh.msh");
auto mesh = gmsh_reader.mesh();
 
auto fe_space =
    std::make_shared<lf::uscalfe::FeSpaceLagrangeO1<double>>(mesh);
 
// We want to solve the pde
// - \laplace u = 0
// u = cos(x)*sin(y) on the boundary
 
// 1) Setup a mesh function that represents the prescribed values of u on the
// boundary
auto mf_boundary =
    mesh::utils::MeshFunctionGlobal([&](const Eigen::Vector2d& x) {
      return std::cos(x[0]) * std::sin(x[1]);
    });
 
// 2) determine the dofs on the boundary and to what value the should be set
auto boundary_cond = InitEssentialConditionFromFunction(
    *fe_space, base::PredicateTrue{}, mf_boundary);
 
// 3) Assemble the stiffness matrix:
assemble::COOMatrix<double> lhs(fe_space->LocGlobMap().NumDofs(),
                                fe_space->LocGlobMap().NumDofs());
auto mf_one = mesh::utils::MeshFunctionConstant(1.);
auto mf_zero = mesh::utils::MeshFunctionConstant(0.);
auto matrix_provider = lf::uscalfe::ReactionDiffusionElementMatrixProvider(
    fe_space, mf_one, mf_zero);
assemble::AssembleMatrixLocally(0, fe_space->LocGlobMap(),
                                fe_space->LocGlobMap(), matrix_provider, lhs);
 
// 4) Modify the system of equations such that the boundary values are
// enforced
Eigen::VectorXd rhs(fe_space->LocGlobMap().NumDofs());
assemble::FixFlaggedSolutionComponents(
    [&](unsigned int idx) { return boundary_cond[idx]; }, lhs, rhs);
 
// 5) solve the problem:
auto lhs_sparse = lhs.makeSparse();
Eigen::SimplicialLDLT<decltype(lhs_sparse)> solver;
solver.compute(lhs_sparse);
auto x = solver.solve(rhs);
 

Table of Contents

Overview

One function on the whole boundary

Constant value on the boundary

BC only on part of the boundary

Example (Essential boundary conditions)