Skip to content

Step 2: Assembly and boundary conditions

This step assembles stiffness matrix and rhs and applies homogeneous Dirichlet constraints.

Learning goals

After this page, you should be able to:

  • use dune-localfunctions Q1 basis evaluation (evaluateFunction, evaluateJacobian),
  • assemble element contributions into BCRSMatrix<double> / BlockVector<double>,
  • enforce zero Dirichlet data by matrix/rhs modification.

Local basis from dune-localfunctions

We use a predefined Q1 Lagrange element from dune-localfunctions:

using LocalFiniteElement = Dune::LagrangeCubeLocalFiniteElement<double,double,2,1>;
LocalFiniteElement lfe;
const auto& basis = lfe.localBasis();

For this tutorial we only need localBasis() to evaluate shape functions and derivatives.

Useful trait types:

  • RangeType = FieldVector<T,1> for scalar basis values,
  • JacobianType = FieldMatrix<T,1,m> for derivatives in local coordinates (m = local dimension, so in 2D: FieldMatrix<T,1,2>).

At each quadrature point \(\hat x\):

  • basis.evaluateFunction(local, phi) gives \(\{\hat\varphi_i(\hat x)\}_i\),
  • basis.evaluateJacobian(local, jacobian) gives \(\{\nabla_{\hat x}\hat\varphi_i(\hat x)\}_i\).

Let \(F_K : \hat K \to K\) be the element map with Jacobian \(J_K\). For each local shape function, the gradient transformation is

\[ \nabla_x \varphi_i(x) = J_K(\hat x)^{-T}\,\nabla_{\hat x}\hat\varphi_i(\hat x). \]

In code this is exactly the mv(...) call with the inverse-transposed Jacobian:

using Traits = typename LocalFiniteElement::Traits::LocalBasisType::Traits;
using JacobianType = typename Traits::JacobianType;
std::vector<JacobianType> jacobian;
basis.evaluateJacobian(local, jacobian); // reference gradients

const auto jacInvT = geo.jacobianInverseTransposed(local);
std::array<Dune::FieldVector<double,2>,numElementDofs> grad;
for (int i = 0; i < numElementDofs; ++i)
  jacInvT.mv(jacobian[i][0], grad[i]);  // grad[i] = J^{-T} * gradHat[i]

The resulting grad[i] vectors are the physical gradients used in the stiffness term.

Assembly loop

Element and quadrature traversal

Assembly iterates over all elements and, inside each element, over all quadrature points:

for (const auto& e : Dune::elements(gv)) {
  const auto v = elementDofs(gv, e); // local vertex dofs -> global dofs
  const auto geo = e.geometry();
  const auto quadrature = Dune::QuadratureRules<double,2>::rule(geo.type(), 2);

  for (const auto& qp : quadrature) {
    const auto local = qp.position();
    basis.evaluateFunction(local, phi);
    basis.evaluateJacobian(local, jacobian);

    const auto jacInvT = geo.jacobianInverseTransposed(local);
    const double dx = qp.weight() * geo.integrationElement(local);

    for (int i = 0; i < numElementDofs; ++i)
      jacInvT.mv(jacobian[i][0], grad[i]);

    // local-to-global accumulation comes here
  }
}

dx is the quadrature contribution to the integration measure on the physical element.

Local-to-global accumulation

Using the global DoF indices v, local contributions are scattered into global A and b:

for (int i = 0; i < numElementDofs; ++i) {
  for (int j = 0; j < numElementDofs; ++j)
    A[v[i]][v[j]] += dx * (grad[i] * grad[j]);
  b[v[i]] += dx * phi[i][0];
}

Here grad[i] * grad[j] is the Euclidean dot product \(\nabla\varphi_i \cdot \nabla\varphi_j\).

Dirichlet boundary treatment

We mark boundary vertices (x=0/1 or y=0/1) and then set:

  • all entries touching boundary rows/columns to 0,
  • boundary diagonal entries to 1,
  • boundary rhs entries to 0.

Boundary DoF index extraction

std::vector<unsigned char> isDirichlet(numDofs, false);
const auto& indexSet = gv.indexSet();
const double eps = 1e-12;

for (const auto& vtx : Dune::vertices(gv)) {
  const auto x = vtx.geometry().center();
  const auto i = indexSet.index(vtx); // global vertex / DoF index
  if (x[0] < eps || x[0] > 1.0 - eps || x[1] < eps || x[1] > 1.0 - eps)
    isDirichlet[i] = true;
}

Matrix and rhs modification

for (std::size_t i = 0; i < numDofs; ++i) {
  for (auto entry = A[i].begin(); entry != A[i].end(); ++entry) {
    const auto j = entry.index();
    if (isDirichlet[i] || isDirichlet[j])
      *entry = (isDirichlet[i] && i == j) ? 1.0 : 0.0;
  }
  if (isDirichlet[i])
    b[i] = 0.0;
}

This yields a linear system with built-in homogeneous Dirichlet constraints.

Summary and next steps

Assembly and boundary conditions are complete. Continue with Step 3: Solve and VTK output.