understand how to implement such a setting, let us first recall how these
arise in the mathematical formulation and how they are implemented in a
matrix-based variant. In essence, an inhomogeneous Dirichlet condition sets
-some of the nodal values in the solution to given values rathern than
+some of the nodal values in the solution to given values rather than
determining them through the variational principles,
@f{eqnarray*}
-u_h(\mathbf{x}) = \sum_{i=1}^{n_\mathrm{dofs}} \varphi_i(\mathbf{x}) u_i =
-\sum_{i=1}^{n_\mathrm{dofs}-n_\mathrm{dofs,D}} \varphi_i(\mathbf{x}) u_i +
-\sum_{i=n_\mathrm{dofs}-n_\mathrm{dofs,D}+1}^{n_\mathrm{dofs}} \varphi_i(\mathbf{x}) g_i,
+u_h(\mathbf{x}) = \sum_{i\in \mathcal N} \varphi_i(\mathbf{x}) u_i =
+\sum_{i\in \mathcal N \setminus \mathcal N_D} \varphi_i(\mathbf{x}) u_i +
+\sum_{i\in \mathcal N_D} \varphi_i(\mathbf{x}) g_i,
@f}
-where $u_i$ denotes the nodal values of the solution and $g_i =
-g(\mathbf{x}_i)$ denotes the interpolation of boundary values on the
-Dirichlet-constrained node points <i>i</i>. In this formula, where we have
-assumed that the degrees of freedom are numbered such that the ones subject to
-a Dirichlet condition are all appended to the end of the list of degrees of
-freedom. We then insert this solution representation into the weak form,
-e.g. the Laplacian shown above, and move the known quantities to the right
-hand side:
+where $u_i$ denotes the nodal values of the solution and $\mathcal N$ denotes
+the set of all nodes. The set $\mathcal N_D\subset \mathcal N$ is the subset
+of the nodes that are subject to Dirichlet boundary conditions where the
+solution is forced to equal $u_i = g_i = g(\mathbf{x}_i)$ as the interpolation
+of boundary values on the Dirichlet-constrained node points $i\in \mathcal
+N_D$. In this formula, where we have assumed that the degrees of freedom are
+numbered such that the ones subject to a Dirichlet condition are all appended
+to the end of the list of degrees of freedom. We then insert this solution
+representation into the weak form, e.g. the Laplacian shown above, and move
+the known quantities to the right hand side:
@f{eqnarray*}
-(\varphi_i, u_h)_\Omega &=& (\varphi_i, f)_\Omega \quad \Rightarrow \\
-\sum_{j=1}^{n_\mathrm{dofs}-n_\mathrm{dofs,D}}(\varphi_i,\varphi_j)_\Omega u_j &=&
+(\nabla \varphi_i, \nabla u_h)_\Omega &=& (\varphi_i, f)_\Omega \quad \Rightarrow \\
+\sum_{j\in \mathcal N \setminus \mathcal N_D}(\nabla \varphi_i,\nabla \varphi_j)_\Omega \, u_j &=&
(\varphi_i, f)_\Omega
--\sum_{j=n_\mathrm{dofs}-n_\mathrm{dofs,D}+1}^{n_\mathrm{dofs}} (\varphi_i,\varphi_j)_\Omega g_j.
+-\sum_{j\in \mathcal N_D} (\nabla \varphi_i,\nabla\varphi_j)_\Omega\, g_j.
@f}
In this formula, the equations are tested for all basis functions $\varphi_i$
-with $i=1,\ldots,n_\mathrm{dofs}-n_\mathrm{dofs,D}$.
+with $i\in N \setminus \mathcal N_D$ that are not related to the nodes
+constrained by Dirichlet conditions.
In the implementation in deal.II, the integrals $(\varphi_i,\varphi_j)_\Omega$
on the right hand side are already contained in the local matrix contributions
we assemble on each cell. When using
-ConstraintMatrix::distributed_local_to_global as first described in the step-6
-and step-7 tutorial programs, we can account for the contribution of
+ConstraintMatrix::distributed_local_to_global() as first described in the
+step-6 and step-7 tutorial programs, we can account for the contribution of
inhomogeneous constraints <i>j</i> by multiplying the columns <i>j</i> and
rows <i>i</i> of the local matrix according to the integrals $(\varphi_i,
\varphi_j)_\Omega$ by the inhomogeneities and subtracting the resulting from
-the position <i>i</i> in the global vector, see also the @ref constraints
-module. In essence, we use some of the integrals that get eliminated from the
-left hand side of the equation to finalize the right hand side
-contribution. Similar mathematics are also involved when first writing all
-entries into a left hand side matrix and then eliminating matrix rows and
+the position <i>i</i> in the global right-hand-side vector, see also the @ref
+constraints module. In essence, we use some of the integrals that get
+eliminated from the left hand side of the equation to finalize the right hand
+side contribution. Similar mathematics are also involved when first writing
+all entries into a left hand side matrix and then eliminating matrix rows and
columns by MatrixTools::apply_boundary_values().
In principle, the components that belong to the constrained degrees of freedom
constraints.distribute(solution);
@endcode
-We could then pass the vector @p solution to the @p LaplaceOperator::vmult_add()
-function and add the new contribution to the @p system_rhs vector that gets
-filled in the @p LaplaceProblem::assemble_rhs() function. However, this idea does
-not work because the vmult() functions eliminates the entries filled into
-Dirichlet-constrained degrees of freedom on the fly by construction. Instead,
-we could select one of two following strategies.
+We could then pass the vector @p solution to the @p
+LaplaceOperator::vmult_add() function and add the new contribution to the @p
+system_rhs vector that gets filled in the @p LaplaceProblem::assemble_rhs()
+function. However, this idea does not work because the
+FEEvaluation::read_dof_values() call used inside the vmult() functions assumes
+homogeneous values on all constraints (otherwise the operator would not be a
+linear operator but an affine one). To also retrieve the values of the
+inhomogeneities, we could select one of two following strategies.
<h5> Use FEEvaluation::read_dof_values_plain() to avoid resolving constraints </h5>
-The class FEEvaluation has a facility that is addresses precisely this
+The class FEEvaluation has a facility that addresses precisely this
requirement: For non-homogeneous Dirichlet values, we do want to skip the
-resolution of the (Dirichlet) constraints upon reading the data from the
-vector @p solution. For example, we could extend the
-@p LaplaceProblem::assemble_rhs() function to deal with inhomogeneous Dirichlet
+implicit imposition of homogeneous (Dirichlet) constraints upon reading the
+data from the vector @p solution. For example, we could extend the @p
+LaplaceProblem::assemble_rhs() function to deal with inhomogeneous Dirichlet
values as follows, assuming the Dirichlet values have been interpolated into
the object @p constraints:
@code
ignores all constraints. Due to this setup, we must make sure that other
constraints, e.g. by hanging nodes, are correctly distributed to the input
vector already as they are not resolved as in
-FEEvaluation::read_dof_values(). Inside the loop, we then evaluate the
+FEEvaluation::read_dof_values_plain(). Inside the loop, we then evaluate the
Laplacian and repeat the second derivative call with
FEEvaluation::submit_gradient() from the @p LaplaceOperator class, but with the
sign switched since we wanted to subtract the contribution of Dirichlet
// proceed as usual with integration of right hand side function...
}
@endcode
+
+Instead of adding a second DoFHandler / ConstraintMatrix pair to the same
+MatrixFree::reinit() call, one could of course also construct an independent
+MatrixFree object that feeds the second @p LaplaceOperator instance, see also
+the discussion in MatrixFreeOperators::Base.