--- /dev/null
+<i>This program was contributed by Sam Scheuerman and Wolfgang
+Bangerth.</i>
+<br>
+
+<a name="step-93-Intro"></a>
+<h1>Introduction</h1>
+
+In all previous tutorial programs, we have only ever dealt with
+discretizations of partial differential equations posed on a
+Triangulation where every degree of freedom was associated with
+a local geometric element (vertices, lines, quadrangles, etc.) which
+is clearly part of one cell or another (or, in the worst case, some
+small number of adjacent cells). This was true
+whether we were solving a single differential equation (say, in step-4
+or step-6), or a system of coupled equations (step-8, step-22, and
+many others). In other words, the *unknowns* we were considering were
+always functions of space (and sometimes time) that we discretized on
+a Triangulation.
+
+There are also situations where we have unknowns in the equations
+we consider that are *not* functions but are *scalars* instead. The
+example we will consider here to illustrate this kind of situation is
+an optimization problem: Let's assume we have a body to which we have
+attached four heaters. How should we set the power levels $C_1,\ldots,
+C_4$ of these four heaters so that the (steady state) temperature $u(\mathbf x)$ of
+the body is as close as possible to a desired state? The issue here is
+that we have one unknown that is a function (namely, the temperature
+$u(\mathbf x)$) for which we can apply the usual finite element
+approximation, but we also have four scalars $C_i$ that are *not*
+functions and that, after discretization of the whole problem, are not
+tied to one of the cells of the mesh. As a consequence, we call these
+unknowns *non-local degrees of freedom* to represent that they are not
+tied to a specific locality: They represent something that does not
+depend to the spatial variable $\mathbf x$ at all, and is not the
+result of a discretization process.
+
+Before going into how we solve this issue, let us first state the concrete set of equations.
+
+We want to find a solution $u$ to the Poisson equation
+
+@f{align*}{
+ -\Delta u =& f \text{ in }\Omega\\
+ u =&0 \text{ on } \partial\Omega,
+@f}
+
+which is as close as
+possible to a desired solution $\overline{u}$. We can formulate this as an optimization problem
+
+@f{align*}{
+ \min_{u, f} & \frac 12 ||u-\overline{u}||^2\\
+ \text{s.t.} & -\Delta u = f.
+@f}
+
+Instead of allowing *any* source function $f$, we will assume that the heat source is described by four
+heating pads $f_k(\mathbf x)$ that are characteristic functions of a subset of the domain, times
+a heater setting $C^k$ (perhaps the power level you set each heating pad to).
+In other words, we have that $f(\mathbf x) = \sum_{k=0}^m C^k f_k(\mathbf x)$.
+The optimization problem then reads as follows:
+@f{align*}{
+ \min_{u, C^k} & \frac 12 ||u-\overline{u}||^2\\
+ \text{s.t.} & -\Delta u = \sum_k C^k f_k.
+@f}
+
+To write this optimization problem in a more compact form, we introduce a Lagrangian
+
+@f{align*}{
+ L = \frac{1}{2}||u - \overline{u}||^2 - \left(\lambda, \left(-\Delta u - \sum_k C^k f_k\right)\right),
+@f}
+
+where $\lambda = \lambda(\mathbf{x})$ is a Lagrange multiplier and $(\cdot, \cdot)$ denotes the usual $L^2$ scalar
+product (note that $\lambda$ is a function that lives in the space of test functions, not just a scalar). We do this because finding a critical point of $L$
+is equivalent to solving the original optimization problem, and we can find a critical point
+of $L$ by finding roots of its derivatives. At this point, we have a choice: we can either
+derive the optimality conditions first, then discretize the problem (in the literature this
+is called the optimize then discretize approach), or we can first discretize $L$ and then
+determine the optimality conditions for the discrete system (discretize then optimize). We
+will use the second approach here, because the formulation is a more straightforward.
+The discretized Lagrangian is
+
+@f{align*}{
+ L =& \frac{1}{2}||u_h - \overline{u}||^2 - (\lambda_h, (-\Delta u_h - f))\\
+ =& \frac{1}{2}\left[\sum_{i=0}^n\sum_{j=0}^n U^iU^j\int_{\Omega}\varphi_i\varphi_j
+ - 2\sum_{i=0}^n U^i\int_\Omega \varphi_i \overline{u} + \int_\Omega\overline{u}^2\right]
+ - \sum_{i=0}^n\sum_{j=0}^n\left[\Lambda^iU^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right]
+ + \sum_{i=0}^n \sum_{k=0}^m \left[\Lambda^i C^k \int_{\Omega}\varphi_i f_k\right].
+@f}
+
+where
+
+@f{align*}{
+ u_h(x) =& \sum_i^n U^i\varphi_i(x),\\
+ \lambda_h =& \sum_j^n \Lambda^j\varphi_j,\text{ and}\\
+ f =& \sum_{k=0}^m C^k f_k.
+@f}
+
+Note the application of integration by parts above, which makes use of the Dirichlet boundary conditions to eliminate the boundary integral. In practice, we assume that each $f_k$ is something simple. In this example step, each $f_k$ is a characteristic function as mentioned, and $m=4$. These functions, plotted together, look like this:
+
+<center><img src="https://www.dealii.org/images/steps/developer/step-93.shape_of_characteristic_functions.png"
+ alt="Visualization of four characteristic functions together on the same plot"
+ width="50%"></center>
+
+
+_Note:_ At this point, we can clearly see where the nonlocal degrees of freedom are coming from. The term $\int_\Omega \varphi_i f_k$ is an integral which cannot a priori be confined to a single cell or small collection of cells, but rather must be evaluated on all cells that lie in the domain of $f$.
+
+To find a critical point, we take a partial derivative with respect to each $U^i$,
+ $\Lambda^j$, and $C^k$. The derivatives have the following form:
+
+@f{align*}{
+ \frac{\partial L}{\partial U^i} =& \sum_{j=0}^nU^j\int_{\Omega}\varphi_i\varphi_j
+ - \int_\Omega \varphi_i\overline{u}
+ - \sum_{j=0}^n\left[\Lambda^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right]\\
+ \frac{\partial L}{\partial \Lambda^j} =& -\sum_{i=0}^n\left[U^i\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right]
+ + \sum_{k=0}^m C^k\int_{\Omega}\varphi_j f_k\\
+ \frac{\partial L}{\partial C^k} =& \sum_{j=0}^n \Lambda^j\int_{\Omega}\varphi_j f_k.
+@f}
+
+Note the paired indices on the left- and right-hand sides. As with many other minimization problems, we can observe that for a critical point to occur,
+these derivatives must all be 0. This gives us $n + n + m$ linear equations to solve
+for the coefficient vectors $U$, $\Lambda$, and $C$:
+
+@f{align*}{
+ \frac{\partial L}{\partial U^i} =& \sum_{j=0}^nU^j\int_{\Omega}\varphi_i\varphi_j
+ - \int_\Omega\varphi_i\overline{u}
+ - \sum_{j=0}^n\left[\Lambda^j\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right] = 0\\
+ \frac{\partial L}{\partial \Lambda^j} =& -\sum_{i=0}^n\left[U^i\int_{\Omega}\nabla\varphi_i\cdot\nabla\varphi_j\right]
+ + \sum_{k=0}^m C^k\int_{\Omega}\varphi_j f_k = 0\\
+ \frac{\partial L}{\partial C^k} =& \sum_{j=0}^n \Lambda^j\int_{\Omega}\varphi_j f_k = 0.\\
+@f}
+
+The summations are suggestive of matrix-vector multiplication, and in fact we can write the equations
+in such a form. Take $(\mathcal{M}_{i, j}) = \int_\Omega\varphi_i\varphi_j$, $(\mathcal{N}_{i, j})
+ = \int_\Omega \nabla\varphi_i\cdot\nabla\varphi_j$, $(\overline{\mathcal{U}}_i) = \int_\Omega\varphi_i\overline{u}$,
+ and $(\mathcal{F}_{j, k}) = \int_\Omega\varphi_j
+ f_k$, and note that $\mathcal{F}$ is an $n\times m$ matrix. Then we get the matrix system
+
+ @f{align*}{
+ \mathcal{M}U - \mathcal{U} - \mathcal{N}^T\Lambda =& 0\\
+ -\mathcal{N}U + \mathcal{F}^TC =& 0\\
+ \mathcal{F}\Lambda =& 0.
+ @f}
+
+Moving known quantities to the right-hand-side and combinging the matrix equations gives us an
+idealized system to solve:
+
+@f{align*}{
+ \left(\begin{array}{c c c}
+ \mathcal{M} & -\mathcal{N}^T & 0\\
+ -\mathcal{N} & 0 & \mathcal{F}^T\\
+ 0 & \mathcal{F} & 0
+ \end{array}\right)
+ \left(\begin{array}{c}
+ U\\\Lambda\\ C
+ \end{array}\right) =
+ \left(\begin{array}{c}
+ \overline{\mathcal{U}}\\0\\0
+ \end{array}\right).
+@f}
+
+Unfortunately, setting up this type of block matrix, where some of the dofs couple in a
+nonlocal way, does not work well with deal.II's internal structure. So, we will need to
+modify how this problem is constructed to work well with deal.II, which is discussed below.
+
+<h3>What is the problem?</h3>
+
+The straightforward way to deal with the equations above is to define
+an FESystem object with two components, one for the temperature
+$u(\mathbf x)$ and one for the Lagrange multiplier $\lambda(\mathbf
+x)$. This would work in much the same way as you define system
+elements in all other @ref vector_valued "vector valued problems".
+You will then end up with $N_U$ variables for the finite element
+approximation of the temperature, and $N_\Lambda$ variables for the
+finite element approximation of the Lagrange multiplier, and all of
+these $N_U+N_\Lambda$ unknowns live on the triangulation in the form
+of a DoFHandler object. While in principle $N_U$ and $N_\Lambda$ can
+be different (if you chose different element types, or different polynomial
+degrees to discretize $u$ and $\lambda$), in practice it is best to keep them the same, so we have
+that $N_U = N_\Lambda = n$.
+
+Now the problem is that what you really want is a linear system of
+$N=n + n + m$ unknowns that includes the $m = 4$ degrees of
+freedom for the source strengths we are trying to optimize for: We
+need a matrix of size $N\times N$, and solution and right hand side
+vectors of size $N$. Moreover, if you try to build a block linear
+system because that's how your linear solver is structured, you
+probably want a $3\times 3$ block structure for the matrix, and three
+blocks for the vectors, even though the finite element and the
+DoFHandler only know of two solution components.
+
+Conceptually, of course, there is nothing wrong with creating matrices
+and vectors that have sizes that don't match the number of unknowns in
+a DoFHandler. You just need a way to translate between the index of a
+degree of freedom in a DoFHandler and the corresponding row and column
+in a linear system. Here, perhaps one would make the choice that the
+mesh-based degrees of freedom (which the DoFHandler numbers from zero
+to $n+n-1$) come first, and the four $C_k$ additional scalar
+unknowns come last. However, the fact that we need to spell out such
+a convention should give you a hint that something with this kind of
+design is wrong. Most codes which would use a convention like this
+will, almost certainly, leave this convention implicit (something
+for the programer to know about) rather than enforcing it
+explicitly (say, in the form of a function `mesh_dof_to_linear_system_index()`).
+It is then very easy for a number of bugs to emerge.
+For example, we can not renumber all degrees of freedom with the
+functions in DoFRenumbering in arbitrary ways if we want the $m$
+additional variables always at the end.
+
+Precisely because it is so easy to make mistakes and because it is so
+restricting to have implicit conventions like "Yes, your matrix can be
+larger than the number of DoFs in the DoFHandler; but it can't be
+smaller", deal.II a long time ago made the decision that a lot of
+functions test that the sizes of vectors and matrices are equal to the
+number of unknowns in a DoFHandler. For example,
+DoFTools::make_sparsity_pattern() takes a DoFHandler and expects that
+the sparsity pattern has as many rows and columns as there are DoFs in
+the DoFHandler. So does
+DoFTools::make_hanging_node_constraints(). VectorTools::interpolate_boundary_values()
+uses the same convention. Many many other functions do too, and they
+all in fact check that this is true and will terminate the program if
+it is not. This is because one can find many many bugs in codes this
+way, and it also allows for an assumption that is so obvious that we
+have never stated it so far: The degree of freedom with index $i$
+corresponds to row $i$ of the linear system and column $i$ of the
+matrix. If we allowed sizes to differ, we would need to have a way to
+make this connection. Further, that way would need to be be explicit (by way of a
+function that translates from DoF index to row and column, as already mentioned above) or there
+would be no end to the possibility of introducing bugs.
+
+In other words, the fact that we want to have $2n$ degrees
+of freedom in the DoFHandler, but matrices and vectors with
+$2n + m$ rows and columns does not work with deal.II's
+underlying design decisions. Creating matrices and vectors that happen
+to be larger is going to trigger lots of assertions, and for good
+reasons.
+
+
+
+<h3>How do we address this?</h3>
+
+So then how are we going to address this? At its core, what needs to
+happen is that the DoFHandler knows about the additional degrees of
+freedom. The question is how we are going to teach it about these
+"nonlocal" degrees of freedom -- they are, after all, not ones that
+are associated with vertices, edges, faces, or cells of the mesh as
+the usual finite element degrees of freedom are. The approach we are
+going to use is not what one would generally call "obvious" or
+"natural", but it does fit into the infrastructure deal.II provides
+(even though this infrastructure was originally built for other
+purposes). In fact, at the time this program is written in 2024, there
+is substantial debate on whether we can implement *better* ways to achieve
+the same result that do not require this kind of "abuse" of features intended
+for other purposes; if consensus is ever reached on these alternative
+approaches, this program may also be converted.
+
+Whatever the future, let us turn to how our goal can be achieved today.
+The key to the following is that deal.II (i) allows using different
+finite elements on different cells, originally intended to support
+$hp$ finite elements; and (ii) has elements that have no degrees of
+freedom. We have used these features in other tutorial programs:
+step-27 and step-75 make use of feature (i) to implement the $hp$
+finite element method, and step-46 makes use of feature (ii) for
+multiphysics problems where some
+variables live only on some cells.
+
+We will (ab)use the general approach used in step-46 for our purposes
+here. Specifically, on the first four cells of the mesh, we will
+define a finite element that has three components: a temperature component and a
+Lagrange multiplier component (both of which are the usual FE_Q
+element), plus a third vector component that is piecewise constant
+(FE_DGQ(0)). On all other
+cells, the finite element used has the same temperature and Lagrange
+multiplier components, but the third component is of type FE_Nothing
+and so has no degrees of freedom at all. Note that with this setup,
+the first four cells will have one additional degree of freedom.
+The end result of this
+approach is that nearly every function you have ever used will work as
+one would expect. In particular:
+
+- We have exactly four additional degrees of freedom, located on the
+ first four cells.
+- We can use DoFRenumbering::component_wise() to sort degrees of freedom
+ explicitly into the order described above (rather than using an
+ implicit convention). In particular, if we really wanted to, we could
+ also ask for the nonlocal degrees of freedom to be numbered *first*
+ rather than last -- because they are on equal footing to all of the other
+ field-based variables.
+- We can use functions such as
+ DoFTools::count_dofs_per_fe_block() or DoFTools::count_dofs_per_fe_component()
+ to create a $3\times 3$ block structure
+ of the linear system, because the finite elements used on all cells have
+ three vector components.
+- The size of matrices and vectors corresponds to the number of degrees
+ of freedom in the DoFHandler.
+- The additional degrees of freedom end up in the output files created by
+ DataOut because they are finite element fields -- just fields of a very
+ specific kind.
+
+What we have with this setup is a good way to not break things in deal.II.
+However, it is not so natural that we can use the full extent of features
+in deal.II without issue. In particular, to the DoFHandler the
+nonlocal degrees of freedom only live on one cell each. On the other
+hand, from the perspective of the PDE, we need these nonlocal DoFs on
+all cells that intersect the areas in which the heating is applied. So,
+we cannot use `cell->get_dof_indices()` to find DoF
+indices for the third vector component on most cells, because the DoF handler
+does not see that these dofs reach beyond the cell they are defined on. As a consequence,
+we will have to do some extra work when building sparsity patterns, and later when building the
+system matrix/system right hand side.
+
+
+<h3>What concrete steps do we need to take?</h3>
+
+There are four places in the standard finite element code where we will
+need to take special care: assigning FE_DGQ elements, constructing the sparsity pattern, constructing
+the system matrix/system right hand side, and outputting results.
+
+1. First, we need to determine on which cells the FE_Nothing element becomes
+an FE_DGQ element. We do this in `setup_system()`, before dofs are distributed.
+Since we are using FE_DGQ elements, we don't need to worry about boundary
+conditions affecting the nonlocal dofs. We simply enumerate the cells
+until we have collected enough dofs, which in this program depends on the
+dimension of the problem. Each time we encounter a new cell, we
+change the hp finite element from an FESystem with a 3rd FE_Nothing
+component to an FESystem with a third FE_DGQ component. Calling
+DoFHandler::distribute_dofs() after this is all we need to do to have the system set
+up correctly. Next, we need to grab the indices of the FE_DGQ dofs. We
+do this by calling DoFTools::extract_dofs(), using a component mask to
+pick just the third component of the FESystem. We store these indices
+in a standard vector, and use them for the three parts below.
+
+2. When we construct the sparsity pattern, deal.II implicitly assumes that
+only dofs that share a cell could have a non-zero entry in the matrix.
+Consequently, when working with non-local dofs, the sparsity pattern that
+DoFTools::make_sparsity_pattern() creates will omit entries, on all
+but the first four cells, that in fact
+we want to fill. To fix this, after creating the sparsity pattern we
+need to loop over all cells in the mesh and check whether the dofs on
+that cell overlap with regions where each $f_k$ is nonzero. If they do,
+then we need to add an entry to the sparsity pattern at location
+(`dof_index`, `nonlocal_dof_index`) and (`nonlocal_dof_index`, `dof_index`).
+Note that we need to add two entries because it doesn't matter whether
+the nonlocal dof is the row or column index, the entry could be
+nonzero. In this example, we further refine this by only adding entries
+where the $\Lambda$ component of the FESystem is nonzero, because we only
+deal with the nonlocal functions in the equations stemming from the
+$\Lambda^j$ and $C^k$ derivatives.
+
+3. We take a similar approach when constructing the system matrix. We do
+the standard loop over cells, quadrature points, and dof indices. Normally,
+we would simply integrate the finite elements on the cell and add these
+values to a local matrix. However, because we are dealing with nonlocal dofs
+this approach won't work. Instead, we use the standard approach for
+constructing the $\mathcal{M}$ and $\mathcal{N}$ portions of the block system
+matrix, but for the $\mathcal{F}$ portions we need to numerically integrate
+each finite element against the nonlocal $f_k$ functions, and enter the
+resulting value in the correct location in the system matrix. To do the integration,
+we simply instantiate a Point object for the quadrature point and pass this
+to the $f_k$ functions, keeping all other aspects of integration the same. We
+then use the extracted dof indices to add the integrated value to the system
+matrix, using the loop's current dof index as the other index. The system right hand
+side is constructed as usual, integrating the right hand side function on each cell.
+
+4. Finally, to properly visualize the solution and constituent parts of the program,
+we need to correctly interpret the nonlocal dofs. We do this by interpolating the
+$f_k$ functions onto a smaller DoFHandler (one which has only one FE_Q element),
+then multiplying by the corresponding coefficient $c_k$. For this program, we also
+add together all the individual $C^kf_k$ to visualize $f$, and interpolate the
+target function as well. These steps (for the $f_k$) all occur in a loop, which makes
+scaling the number of nonlocal dofs relatively straightforward.
+
+To see how we actually do these things, let's look at the code.
--- /dev/null
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2024 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * Part of the source code is dual licensed under Apache-2.0 WITH
+ * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+ * governing the source code and code contributions can be found in
+ * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+ *
+ * ------------------------------------------------------------------------
+ *
+ * Authors : Sam Scheuerman, Wolfgang Bangerth, Colorado State University, 2024.
+ */
+
+
+// @sect3{Include files and other top matter}
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_enriched.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_minres.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+
+namespace Step93
+{
+ using namespace dealii;
+
+
+ // @sect3{Target and indicator functions}
+
+ // We start by defining a function class for $\bar u$; it represents the heat
+ // profile we want to match.
+ //
+ // This class has two member variables:
+ //
+ // - `center`: A Point object representing the center of the function
+ // - `radius`: A double representing the radius of the circular step, or the
+ // standard deviation of a Gaussian
+ template <int dim>
+ class TargetFunction : public Function<dim>
+ {
+ public:
+ TargetFunction()
+ : Function<dim>(1)
+ , center()
+ , radius(0.1)
+ {}
+
+ // Next, we define an overloaded constructor for TargetFunction, so we can
+ // choose to set the center and radius
+ //
+ // Parameters:
+ //
+ // - `center`: A constant Point pointer used to set the member variable
+ // center
+ // - `radius`: A constant double used to set the member variable radius
+ TargetFunction(const unsigned int n_components,
+ const Point<dim> ¢er,
+ const double radius = .3)
+ : Function<dim>(n_components)
+ , center(center)
+ , radius(radius)
+ {}
+
+ virtual double value(const Point<dim> &p,
+ const unsigned int component = 0) const override;
+
+ private:
+ const Point<dim> center;
+ const double radius;
+ };
+
+ // The value() function returns the value of the function at point `p` and
+ // component index `component`. In this case, if the component corresponds to
+ // the solution $u$, then the function returns a value based on either a step
+ // function or a Gaussian. If the component corresponds to any of the other
+ // variables, the function always returns 0.
+ //
+ // Note: In the documentation we discuss using a Gaussian target function, but
+ // this is not included in the code for readability. To run the code with a
+ // Gaussian target function, simply replace the code in the `component == 0`
+ // braces by
+ //
+ // `return std::exp(-((p-center).norm()*(p-center).norm())/(radius*radius));`
+ template <int dim>
+ double TargetFunction<dim>::value(const Point<dim> &p,
+ const unsigned int component) const
+ {
+ if (component == 0)
+ {
+ if ((p - center).norm() <= radius)
+ return 1;
+ else
+ return 0;
+ }
+ else
+ return 0;
+ }
+
+ // The next class we define is one for a circular indicator function. Unlike
+ // the target function, this does not need a component argument because we
+ // have to manually address where it gets used in the code. Objects of this
+ // function type correspond to the nonlocal dofs.
+ //
+ // Parameters:
+ //
+ // - `center`: A constant Point object giving the center of the indicator
+ // region.
+ // - `radius`: The radius of the region.
+ template <int dim>
+ class CircularIndicatorFunction : public Function<dim>
+ {
+ public:
+ CircularIndicatorFunction();
+
+ CircularIndicatorFunction(const Point<dim> ¢er, const double radius);
+
+ virtual double value(const Point<dim> &p,
+ const unsigned int component = 0) const override;
+
+ private:
+ const Point<dim> center;
+ const double radius;
+ };
+
+ template <int dim>
+ CircularIndicatorFunction<dim>::CircularIndicatorFunction()
+ : center(Point<dim>())
+ , radius(20)
+ {}
+
+
+ template <int dim>
+ CircularIndicatorFunction<dim>::CircularIndicatorFunction(
+ const Point<dim> ¢er,
+ const double radius)
+ : center(center)
+ , radius(radius)
+ {}
+
+ template <int dim>
+ double CircularIndicatorFunction<dim>::value(
+ const Point<dim> &p,
+ const unsigned int /* component */) const
+ {
+ if ((center - p).norm() <= radius)
+ return 1;
+ else
+ return 0;
+ }
+
+
+ // @sect3{The principal class}
+
+ // The main class is very similar to step-4 in structure, given that this is
+ // a relatively simple program that does not use adaptive mesh refinement.
+ // However, there are three new member variables:
+ //
+ // - `nonlocal_dofs`: A `std::vector` of dof indices that stores the
+ // dof index for the nonlocal dofs.
+ //
+ // - `heat_functions`: A `std::vector` of CircularIndicatorFunction objects;
+ // these are the heat sources.
+ //
+ // - `target_function`: This is the function we want to match. We store it
+ // as a class variable because it is used both in assemble_system() and
+ // output_results().
+ template <int dim>
+ class Step93
+ {
+ public:
+ Step93();
+
+ void run();
+
+ private:
+ void make_grid();
+ void setup_system();
+ void assemble_system();
+ void solve();
+ void output_results() const;
+
+ Triangulation<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_collection;
+
+ AffineConstraints<double> constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ std::vector<types::global_dof_index> nonlocal_dofs;
+
+ std::vector<CircularIndicatorFunction<dim>> heat_functions;
+
+ const TargetFunction<dim> target_function;
+ };
+
+
+ // @sect4{The Step93 constructor}
+ // In this constructor, we set up several of the fundamental data
+ // structures this program needs. Specifically, we generate the
+ // finite element collection, which is basically a list of all the
+ // possible finite elements we could use on each cell. This
+ // collection has two elements: one FESystem that has two degree 2
+ // FE_Q elements and one FE_Nothing element (to be used on all cells
+ // that are not "special"), and one FESystem that has two degree 2
+ // FE_Q elements and one degree 0 FE_DGQ element (to be used on
+ // those "special" cells we use to anchor the non-local degrees of
+ // freedom, as discussed in the introduction).
+ //
+ // Where we have a collection of elements, we then also need a
+ // collection of quadratures -- which here has only a single element
+ // because we can use the same quadrature on all cells.
+
+ // The default constructor also initializes
+ // dof_handler and target_function,
+ // and it generates a vector of CircularIndicatorFunctions objects which
+ // represent the heat functions.
+ template <int dim>
+ Step93<dim>::Step93()
+ : dof_handler(triangulation)
+ , target_function(3,
+ (dim == 1 ? Point<dim>(0.5) :
+ dim == 2 ? Point<dim>(0.5, 0.5) :
+ Point<dim>(0.5, 0.5, 0.5)))
+ {
+ // Here, we generate the finite element collection, which is basically a
+ // list of all the possible finite elements we could use on each cell. This
+ // collection has two elements: one FESystem that has two degree 2 FE_Q
+ // elements and one FE_Nothing element, and one FESystem that has two
+ // degree 2 FE_Q elements and one degree 0 FE_DGQ element.
+ fe_collection.push_back(
+ FESystem<dim>(FE_Q<dim>(2), 2, FE_Nothing<dim>(), 1));
+ fe_collection.push_back(FESystem<dim>(FE_Q<dim>(2), 2, FE_DGQ<dim>(0), 1));
+
+ // The quadrature collection is just one degree 3 QGauss element.
+ quadrature_collection.push_back(QGauss<dim>(3));
+
+ // Here, we choose the center points for the heat functions to be the
+ // vertices of a lattice, in this case the corners of a hypercube
+ // centered at 0, with side length 1. The block of code below
+ // then creates the CircularIndicatorFunction
+ // objects for the main class. To make these, we check the dimension of the
+ // problem and create 2, 4, or 8 function objects.
+ switch (dim)
+ {
+ case (1):
+ heat_functions.emplace_back(Point<dim>({0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5}), 0.2);
+ break;
+ case (2):
+ heat_functions.emplace_back(Point<dim>({0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({0.5, -0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, -0.5}), 0.2);
+ break;
+ case (3):
+ heat_functions.emplace_back(Point<dim>({0.5, 0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({0.5, 0.5, -0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({0.5, -0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({0.5, -0.5, -0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, 0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, 0.5, -0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, -0.5, 0.5}), 0.2);
+ heat_functions.emplace_back(Point<dim>({-0.5, -0.5, -0.5}), 0.2);
+ break;
+ default:
+ DEAL_II_ASSERT_UNREACHABLE();
+ break;
+ }
+ }
+
+ // @sect4{Step93::make_grid()}
+
+ // The make_grid() function makes a hypercube grid, see step-4:
+ template <int dim>
+ void Step93<dim>::make_grid()
+ {
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(7);
+
+ std::cout << "Number of active cells: " << triangulation.n_active_cells()
+ << std::endl;
+ }
+
+
+ // @sect4{Step93::setup_system()}
+
+ // The `setup_system()` function is similar to step-4, except we have to add a
+ // few steps to prepare for the nonlocal dofs.
+ template <int dim>
+ void Step93<dim>::setup_system()
+ {
+ // Here, we loop over the cells and set the FESystem index to 1, which
+ // corresponds to the system with 2 FE_Q elements and one FE_DGQ element. We
+ // do this until we have enough dofs for each heat function. Note that we
+ // use the global active cell index to measure when to stop the
+ // loop. This allows the loop to be run in parallel with no alteration.
+ // Then, we call DoFHandler::distribute_dofs() to actually enumerate all
+ // degrees of freedom.
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ if (cell->global_active_cell_index() < heat_functions.size())
+ {
+ cell->set_active_fe_index(1);
+ }
+ else
+ {
+ break;
+ }
+ }
+ dof_handler.distribute_dofs(fe_collection);
+
+ // Once we've assigned dofs, the code block below counts the
+ // number of dofs in the system, and outputs to the console. In
+ // other contexts, we might want to use *block* matrices (see, for
+ // example, step-20 or step-22) to build more efficient linear
+ // solvers; here, we will just put everything into one big matrix
+ // and so knowing the number of unknowns for each of the variables
+ // is purely for informational purposes.
+ const std::vector<types::global_dof_index> dofs_per_component =
+ DoFTools::count_dofs_per_fe_component(dof_handler);
+ const unsigned int n_dofs_u = dofs_per_component[0],
+ n_dofs_l = dofs_per_component[1],
+ n_dofs_c = dofs_per_component[2];
+ std::cout << "Number of degrees of freedom: " << n_dofs_u << "+" << n_dofs_l
+ << "+" << n_dofs_c << " = " << n_dofs_u + n_dofs_l + n_dofs_c
+ << std::endl;
+
+ // Finally, we need to extract the indices of the finite elements which
+ // correspond to the non-local dofs.
+ //
+ // First, we make a component mask, which is `false` except for
+ // the third component. This will extract only the dofs from the
+ // third component of the FE system. Next, we actually extract the
+ // dofs, and store them in an IndexSet variable. Finally, we add
+ // each extracted index to the member array `nonlocal_dofs`.
+ const ComponentMask component_mask_c({false, false, true});
+ const IndexSet indices_c =
+ DoFTools::extract_dofs(dof_handler, component_mask_c);
+
+ for (const types::global_dof_index non_local_index : indices_c)
+ nonlocal_dofs.push_back(non_local_index);
+
+ std::cout << "Number of nonlocal dofs: " << nonlocal_dofs.size()
+ << std::endl;
+
+
+ // The mesh we created above is not locally refined, and so there
+ // are no hanging node constraints to keep track of. But it does
+ // not hurt to just use the same setup we have used starting in
+ // step-6 of building a constraints object that contains hanging
+ // node constraints (anticipating that perhaps we'd want to do
+ // adaptive mesh refinement in a later step) into which we then
+ // also put the constraints for boundary values on both the $u$
+ // and $\lambda$ variables.
+ //
+ // Because the nonlocal degrees of freedom use discontinuous
+ // elements, they do not contribute to boundary values
+ // (discontinuous elements do not have degrees of freedom
+ // logically located on the boundary that could be interpolated)
+ // and we do not need to exclude these solution components
+ // explicitly when calling
+ // VectorTools::interpolate_boundary_values()
+ constraints.clear();
+ DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 0,
+ Functions::ZeroFunction<dim>(3),
+ constraints);
+ constraints.close();
+
+ // The remainder of the function deals with building the sparsity
+ // pattern. It consists of two parts: The entries that result from
+ // the usual integration of products of shape functions, and then
+ // the entries that result from integrals that contain nonlocal
+ // degrees of freedom (which one can think of as associated with
+ // shape functions that are constant across the entire
+ // domain). The first part is easily built using standard tools:
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+
+ // The other part is more awkward. We have matrix entries that
+ // result from terms such as $\int_{\Omega}\varphi_j f_k$ where
+ // each $\varphi_j$ is a shape function associated to $\lambda$
+ // and each $f_k$ is a
+ // characteristic function of a part of the domain. As a
+ // consequence, we end up with a nonzero matrix entry $A_{jk}$
+ // (along with its transpose $A_{kj}$) if there is overlap between
+ // $\varphi_j$ and $f_k$. In practice, because we will use
+ // quadrature, this means that we end up with a quadrature point
+ // on a cell on which $\varphi_j$ lives and at which $f_k$ is not
+ // zero. (We will implicitly assume that a shape function that
+ // lives on the current cell is nonzero at all quadrature points
+ // -- an assumption that is generally true unless one chooses
+ // specific quadrature formulas.) Determining which sparsity
+ // pattern entries we need to add then essentially comes down to
+ // "simulating" what would happen if we actually computed
+ // integrals and which matrix entries would end up being non-zero
+ // in the process. As a consequence, the following code's general
+ // structure looks very similar to what we will do for the
+ // nonlocal contributions in the `assemble_system()` function
+ // below.
+ //
+ // To get this started, we create an `hp_fe_values` that we
+ // will only use to query the quadrature point locations. The
+ // non-local dofs will need to interact with the second component
+ // of the fe system (namely, $\lambda$), so we also declare a
+ // variable that will help us extract this scalar field for use
+ // below.
+ hp::FEValues<dim> hp_fe_values(fe_collection,
+ quadrature_collection,
+ update_quadrature_points);
+
+ const FEValuesExtractors::Scalar lambda(1);
+
+ // Then, we loop over the cells, then over the quadrature points, and
+ // finally over the indices, as if we were constructing a mass matrix.
+ // However, what we instead do here is check two things. First, we check if
+ // the quadrature point is within the radius of a circular indicator
+ // function that represents our non-local dof. If so
+ // then we add an entry to the sparse matrix at the
+ // (nonlocal dof index, lambda dof index) entry and the (lambda dof index,
+ // nonlocal dof index) entry for all lambda degrees of freedom. (Because the
+ // matrix we solve with has both the lambda-nonlocal interacting block and
+ // its transpose, we need to add two entries each time.)
+ std::vector<types::global_dof_index> local_dof_indices;
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ hp_fe_values.reinit(cell);
+
+ const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+ local_dof_indices.resize(fe_values.dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+
+ for (const unsigned int q_index : fe_values.quadrature_point_indices())
+ {
+ const Point<dim> q_point = fe_values.quadrature_point(q_index);
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ if (fe_values.get_fe().system_to_component_index(i).first ==
+ 1) // 'i' is a lambda shape function
+ {
+ for (unsigned int j = 0; j < heat_functions.size(); ++j)
+ if (heat_functions[j].value(q_point) != 0)
+ {
+ dsp.add(local_dof_indices[i], nonlocal_dofs[j]);
+ dsp.add(nonlocal_dofs[j], local_dof_indices[i]);
+ }
+ }
+ }
+ }
+ }
+
+ // The rest (below) is standard setup code, see step-4:
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+ }
+
+
+ // @sect4{Step93::assemble_system()}
+
+ // The `assemble_system()` function works very similar to how is
+ // does in other tutorial programs (cf. step-4, step-6, step-8, and
+ // for the vector-valued case see step-22). However, there is an
+ // additional component to constructing the system matrix, because
+ // we need to handle the nonlocal dofs manually.
+ template <int dim>
+ void Step93<dim>::assemble_system()
+ {
+ // First, we do a standard loop setup for constructing the system matrix.
+ hp::FEValues<dim> hp_fe_values(fe_collection,
+ quadrature_collection,
+ update_values | update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
+
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ const FEValuesExtractors::Scalar u(0);
+ const FEValuesExtractors::Scalar lambda(1);
+ const FEValuesExtractors::Scalar c(2);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
+
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+ hp_fe_values.reinit(cell);
+
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+ local_dof_indices.resize(fe_values.dofs_per_cell);
+
+ cell->get_dof_indices(local_dof_indices);
+
+
+ // In the loop over quadrature points, we start by building
+ // all of the usual terms that are bilinear in shape functions
+ // corresponding to the $u$ and $\lambda$ variables:
+ for (const unsigned int q_index : fe_values.quadrature_point_indices())
+ {
+ const double JxW = fe_values.JxW(q_index);
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ const double phi_i_u = fe_values[u].value(i, q_index),
+ phi_i_l = fe_values[lambda].value(i, q_index);
+
+ const Tensor<1, dim> grad_i_u =
+ fe_values[u].gradient(i, q_index),
+ grad_i_l =
+ fe_values[lambda].gradient(i, q_index);
+
+ for (const unsigned int j : fe_values.dof_indices())
+ {
+ const double phi_j_u = fe_values[u].value(j, q_index);
+
+ const Tensor<1, dim> grad_j_u =
+ fe_values[u].gradient(j, q_index),
+ grad_j_l =
+ fe_values[lambda].gradient(j,
+ q_index);
+
+ cell_matrix(i, j) += phi_i_u * phi_j_u * JxW;
+ cell_matrix(i, j) += -grad_i_u * grad_j_l * JxW;
+ cell_matrix(i, j) += -grad_i_l * grad_j_u * JxW;
+ }
+
+ const Point<dim> q_point = fe_values.quadrature_point(q_index);
+ cell_rhs(i) += (phi_i_u * target_function.value(q_point) * JxW);
+
+
+ // For the integrals that involve the nonlocal dofs,
+ // we make use of the quadrature point again. To
+ // compute the integrals, we loop over the
+ // heat functions, adding the numeric integral of each
+ // heat equation with each $\lambda$ shape function,
+ // at the appropriate indices (which we found in
+ // `setup_system()`). Note that if we try to add 0 to
+ // a matrix entry we have not previously indicated should
+ // be nonzero, there will not be a problem; but if we
+ // try to add a nonzero value to an entry not
+ // previously added to the sparsity pattern, we will
+ // get an error. In other words, the following lines
+ // of the code check that we adjusted the sparsity
+ // pattern correctly in the previous function.
+ for (unsigned int j = 0; j < heat_functions.size(); ++j)
+ {
+ system_matrix.add(local_dof_indices[i],
+ nonlocal_dofs[j],
+ heat_functions[j].value(q_point) *
+ phi_i_l * JxW);
+ system_matrix.add(nonlocal_dofs[j],
+ local_dof_indices[i],
+ heat_functions[j].value(q_point) *
+ phi_i_l * JxW);
+ }
+ }
+ }
+
+ // Finally, we copy the local contributions to the linear
+ // system into the global matrix and right hand side vector,
+ // taking into account hanging node and boundary values
+ // constraints:
+ constraints.distribute_local_to_global(
+ cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+ }
+ }
+
+
+ // @sect4{Step93::solve()}
+
+ // The solve() function works similar to how it is done in step-6
+ // and step-8, except we need to use a different solver because the
+ // linear problem we are trying to solve is a saddle point problem
+ // for which the Conjugate Gradient algorithm is not
+ // applicable. But, because the matrix is symmetric, we can use
+ // SolverMinRes, an iterative solver specialized for symmetric
+ // indefinite problems. This solver could be improved with the use
+ // of preconditioners, but we don't do that here for simplicity (see
+ // the Possibilities for Extensions section below).
+ //
+ // As you will see in the output, given that we are not using a
+ // preconditioner, we need a *lot* of iterations to solve this
+ // linear system. We set the maximum to one million, more than we
+ // need of course, but an indication that this is not an efficient
+ // solver. For smaller problems, one can also use a direct solver
+ // (see step-29) for which you would just replace the main part of
+ // this function by the following three lines of code:
+ //
+ // `SparseDirectUMFPACK direct_solver;`
+ //
+ // `direct_solver.initialize(system_matrix);`
+ //
+ // `direct_solver.vmult(solution, system_rhs);`
+ template <int dim>
+ void Step93<dim>::solve()
+ {
+ std::cout << "Beginning solve..." << std::endl;
+ Timer timer;
+
+ SolverControl solver_control(1'000'000, 1e-6 * system_rhs.l2_norm());
+ SolverMinRes<Vector<double>> solver(solver_control);
+
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ timer.stop();
+ std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl;
+ std::cout << "Solved in " << solver_control.last_step()
+ << " MINRES iterations." << std::endl;
+ }
+
+
+ // @sect4{Step93::output_results()}
+
+ // The `output_results()` function is a bit more robust for this program than
+ // is typical. This is because, in order to visualize the heat sources we have
+ // optimized, we need to do extra work and interpolate them onto a mesh. We do
+ // this by instantiating a new DoFHandler object and then using the helper
+ // function VectorTools::interpolate().
+ //
+ // The top of the function is as always when using vector-valued
+ // elements (see, for example, step-22) and simply outputs all of
+ // the solution variables on the mesh cells they are defined on:
+ template <int dim>
+ void Step93<dim>::output_results() const
+ {
+ const std::vector<std::string> solution_names = {"u", "lambda", "c"};
+
+ const std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ interpretation = {DataComponentInterpretation::component_is_scalar,
+ DataComponentInterpretation::component_is_scalar,
+ DataComponentInterpretation::component_is_scalar};
+
+ DataOut<dim> data_out;
+ data_out.add_data_vector(dof_handler,
+ solution,
+ solution_names,
+ interpretation);
+
+ // The non-local degrees of freedom are of course defined on the
+ // first several cells, but at least logically are considered to
+ // live everywhere (or, if you prefer, nowhere at all since they
+ // do not represent spatial functions). But, conceptually, we use
+ // them as multipliers for the heat sources, and so while the
+ // coefficients $C^k$ are non-local, the heat source $\sum_k C^k
+ // f_k(\mathbf x)$ *is* a spatially variable function. It would be
+ // nice if we could visualize that as well. The same is true for
+ // the target heat distribution $\bar u$ we are trying to match.
+ //
+ // To do so, we
+ // create a new dof handler to output the target function
+ // and heat plate values, and associate it with
+ // a finite element with a degree that matches what we used to solve for $u$
+ // and $\lambda$, although in reality this is an arbitrary choice:
+ DoFHandler<dim> new_dof_handler(triangulation);
+
+ const FE_Q<dim> new_fe(2);
+ new_dof_handler.distribute_dofs(new_fe);
+
+ // To get started with the visualization, we need a vector which
+ // stores the interpolated target function. We create the vector,
+ // interpolate the target function $\bar u$ onto the mesh, then
+ // add the data to our data_out object.
+ Vector<double> target(new_dof_handler.n_dofs());
+ VectorTools::interpolate(new_dof_handler, target_function, target);
+ data_out.add_data_vector(new_dof_handler, target, "u_bar");
+
+ // In order to visualize the sum of the heat sources $\sum_k C^k
+ // f_k(\mathbf x)$, we create a vector which will store the
+ // interpolated values of this function. Then, we loop through
+ // the heat functions, create a vector to store the interpolated
+ // data, call the VectorTools::interpolate() function to fill the
+ // vector, multiply the interpolated data by the nonlocal dof
+ // value $C^k$ (so that the heat plate is set to the correct
+ // temperature), and then add this data to the sum of heat
+ // sources. Because we can, we also add the vector for each source
+ // individually to the DataOut object, so that they can be
+ // visualized individually.
+ Vector<double> full_heat_profile(new_dof_handler.n_dofs());
+
+ for (unsigned int i = 0; i < heat_functions.size(); ++i)
+ {
+ Vector<double> hot_plate_i(new_dof_handler.n_dofs());
+
+ VectorTools::interpolate(new_dof_handler,
+ heat_functions[i],
+ hot_plate_i);
+
+ hot_plate_i *= solution[nonlocal_dofs[i]];
+ full_heat_profile += hot_plate_i;
+
+ const std::string data_name =
+ "Heat_Source_" + Utilities::int_to_string(i);
+ data_out.add_data_vector(new_dof_handler, hot_plate_i, data_name);
+ }
+
+ // Once all the heat functions have been combined, we add them to the
+ // data_out object, and output everything into a file:
+ data_out.add_data_vector(new_dof_handler,
+ full_heat_profile,
+ "Full_Heat_Profile");
+
+ data_out.build_patches();
+
+ std::ofstream output("solution.vtu");
+ data_out.write_vtu(output);
+
+ // Finally, we output the nonlocal coefficient values to the console:
+ std::cout << "The c coefficients are " << std::endl;
+ for (long unsigned int i = 0; i < nonlocal_dofs.size(); ++i)
+ {
+ std::cout << "\tc" << i + 1 << ": " << solution[nonlocal_dofs[i]]
+ << std::endl;
+ }
+ }
+
+
+ // @sect4{Step93::run()}
+
+ // The run() function runs through each step of the program, nothing new here:
+ template <int dim>
+ void Step93<dim>::run()
+ {
+ make_grid();
+ setup_system();
+ assemble_system();
+ solve();
+ output_results();
+ }
+} // namespace Step93
+
+
+// @sect3{The main() function}
+
+// The `main()` function looks essentially like that of most other tutorial
+// programs.
+int main()
+{
+ try
+ {
+ Step93::Step93<2> heat_optimization_problem;
+ heat_optimization_problem.run();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}