]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add the step-93 tutorial program. 17848/head
authorSam Scheuerman <sam.scheuerman.051@gmail.com>
Thu, 12 Sep 2024 16:57:52 +0000 (10:57 -0600)
committerSam Scheuerman <sam.scheuerman.051@gmail.com>
Tue, 15 Apr 2025 16:18:26 +0000 (10:18 -0600)
doc/doxygen/tutorial/tutorial.h.in
doc/news/changes/major/20241106Scheuerman [new file with mode: 0644]
examples/step-93/CMakeLists.txt [new file with mode: 0644]
examples/step-93/doc/builds-on [new file with mode: 0644]
examples/step-93/doc/intro.dox [new file with mode: 0644]
examples/step-93/doc/kind [new file with mode: 0644]
examples/step-93/doc/results.dox [new file with mode: 0644]
examples/step-93/doc/tooltip [new file with mode: 0644]
examples/step-93/step-93.cc [new file with mode: 0644]
tests/examples/step-93.diff [new file with mode: 0644]

index ee6a84be73169caed03a2832b6174db151a29a2f..9acd8df45724a260b46f08c4c74975ec47685fb8 100644 (file)
  *       <td> Solving the Laplace-Beltrami equation on a surface using the trace finite element method.
  *       <br/> Keywords: MeshWorker::mesh_loop(), NonMatching::FEImmersedSurfaceValues
  *       </td></tr>
+ * 
+ *   <tr valign="top">
+ *       <td>step-93</td>
+ *       <td> Using nonlocal dofs to solve a simple optimization problem with the Poisson equation
+ *       <br/> Keywords: nonlocal dofs, optimization, Poisson
+ *       </td></tr>
  * </table>
  *
  *
diff --git a/doc/news/changes/major/20241106Scheuerman b/doc/news/changes/major/20241106Scheuerman
new file mode 100644 (file)
index 0000000..a89ab28
--- /dev/null
@@ -0,0 +1,4 @@
+New: The step-93 tutorial program shows how to use nonlocal dofs in the
+deal.II framework, in the context of a simple optimization problem.
+<br>
+(Sam Scheuerman, Wolfgang Bangerth, 2024/11/06)
diff --git a/examples/step-93/CMakeLists.txt b/examples/step-93/CMakeLists.txt
new file mode 100644 (file)
index 0000000..af81e34
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-93 tutorial program:
+##
+
+# Set the name of the project and target:
+set(TARGET "step-93")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#    file(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#    file(GLOB_RECURSE TARGET_INC  "include/*.h")
+#    set(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+set(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+cmake_minimum_required(VERSION 3.13.4)
+
+find_package(deal.II 9.7.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ ../../../ ../../build-dealii $ENV{DEAL_II_DIR}
+  )
+if(NOT ${deal.II_FOUND})
+  message(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+endif()
+
+deal_ii_initialize_cached_variables()
+project(${TARGET})
+deal_ii_invoke_autopilot()
diff --git a/examples/step-93/doc/builds-on b/examples/step-93/doc/builds-on
new file mode 100644 (file)
index 0000000..11f3e73
--- /dev/null
@@ -0,0 +1 @@
+step-4 step-8
diff --git a/examples/step-93/doc/intro.dox b/examples/step-93/doc/intro.dox
new file mode 100644 (file)
index 0000000..03a121c
--- /dev/null
@@ -0,0 +1,370 @@
+<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.
diff --git a/examples/step-93/doc/kind b/examples/step-93/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-93/doc/results.dox b/examples/step-93/doc/results.dox
new file mode 100644 (file)
index 0000000..0f3af8f
--- /dev/null
@@ -0,0 +1,188 @@
+<h1>Results</h1>
+
+When you run the program with the step target function (in 2D), the output looks something like this:
+
+@code
+Number of active cells: 16384
+Number of degrees of freedom: 66049+66049+4 = 132102
+Number of nonlocal dofs: 4
+Beginning solve...
+Wall time: 63.9265s
+Solved in 39973 MINRES iterations.
+The c coefficients are
+       c1: 28.7408
+       c2: -6.51604
+       c3: -6.51604
+       c4: -1.62044
+@endcode
+
+When you run the program with the Gaussian target function (in 2D), the output should look like this:
+
+@code
+Number of active cells: 16384
+Number of degrees of freedom: 66049+66049+4 = 132102
+Number of nonlocal dofs: 4
+Beginning solve...
+Wall time: 98.4858s
+Solved in 62131 MINRES iterations.
+The c coefficients are
+       c1: 23.553
+       c2: -4.86562
+       c3: -4.86562
+       c4: -1.42344
+@endcode
+
+The goal of this program is to determine which temperature settings best match the target function, so first
+let's see what these targets look like:
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+      <center><b>Step Target %Function</b></center>
+    </td>
+    <td>
+      <center><b>Gaussian Target %Function</b></center>
+    </td>
+  </tr>
+  <tr>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.target_step.png"
+         alt="Step target function"
+         width="90%">
+    </td>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.target_gauss.png"
+         alt="Gaussian target function"
+         width="90%">
+    </td>
+  </tr>
+</table>
+
+After solving the Lagrangian system, we arrive at solutions $U_\text{step}$ and $U_\text{gauss}$ that
+look like this:
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+    <center><b>$U_\text{step}$</b></center>
+    </td>
+    <td>
+    <center><b>$U_\text{gauss}$</b></center>
+    </td>
+  </tr>
+  <tr>
+    <td>
+    <img src="https://www.dealii.org/images/steps/developer/step-93.U_step.png"
+         alt="Solution for step shaped target function"
+         width="90%">
+    </td>
+    <td>
+    <img src="https://www.dealii.org/images/steps/developer/step-93.U_gauss.png"
+         alt="Solution for Gaussian target function"
+         width="90%">
+    </td>
+  </tr>
+</table>
+
+Notice that $U_\text{gauss}$ matches the target much better than
+$U_\text{step}$. Intuitively, this makes sense: in general, solutions
+to the heat equation look something like Gaussians, so the
+Gaussian target function is a much more "natural" thing to match than
+a sharp step function. We can also see this in the optimal heat
+profiles.
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+      <center><b>Heat plate settings for matching step function</b></center>
+    </td>
+    <td>
+      <center><b>Heat plate settings for matching Gaussian</b></center>
+    </td>
+  </tr>
+  <tr>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.heat_profile_step.png"
+         alt="Heat plate settings for matching step function"
+         width="90%">
+    </td>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.heat_profile_gauss.png"
+         alt="Heat plate settings for matching Gaussian"
+         width="90%">
+    </td>
+  </tr>
+</table>
+
+Notice that for the Gaussian target, the 4 plates are set to less extreme values. In contrast, to try to match the step function, higher and lower temperatures must be applied.
+
+While it does not contain much useful information, we can also plot the Lagrange multiplier $\Lambda$, which has an interesting shape:
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+      <center><b>$\Lambda_\text{step}$</b></center>
+    </td>
+    <td>
+      <center><b>$\Lambda_\text{gauss}$</b></center>
+    </td>
+  </tr>
+  <tr>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.L_step.png"
+         alt="Lagrange multiplier for step target function"
+         width="90%">
+    </td>
+    <td>
+      <img src="https://www.dealii.org/images/steps/developer/step-93.L_gauss.png"
+         alt="Lagrange multiplier for Gaussian target function"
+         width="90%">
+    </td>
+  </tr>
+</table>
+
+
+
+<h3>Possibilities for extensions</h3>
+
+There are a few ways that this program could be extended, which we list below.
+
+1. As mentioned in the code documentation, this program does not make
+use of any preconditioners before solving. This is because, for a 2D
+problem, the code runs fast enough that optimization is not
+necessary. However, as shown in the screen output above, the number of
+iterations required to solve the linear system is quite large. Thus,
+for larger problems, it would be good if the solver ran more
+quickly. See the "Possibilities for extensions" section of step-6 for
+a more detailed discussion on how to change preconditioners. We should
+note that since the block matrix we use has many zeros on the
+diagonal, preconditioners like PreconditionJacobi will not work
+because they divide by diagonal entries. Instead, block
+preconditioners such as those discussed in step-20 or step-22 (among
+many others) will likely be useful.
+
+2. To validate the optimization problem is working correctly, we could
+try to match a target function which is itself a solution to the
+Poisson equation with prescribed heat profile. If the optimization
+problem is being solved correctly, it should be able to perfectly
+match this solution. To create such a function, we would need to first
+solve the Poisson problem on a scalar field, with a RHS described by
+the chosen heat profile. See step-7 for more information on the method
+of manufactured solutions.
+
+3. The program at the moment has the number of nonlocal degrees of freedom
+hardcoded as $2^d$ (see the constructor). We then assign each of these
+degrees of freedom to one of the first $2^d$ cells. This is not going to
+be much of a problem because there are always enough cells for this as
+long as you start with a mesh that is at least once refined. What would
+we do if we had a number of nonlocal DoFs that is not easily predictable,
+and that may be larger than the number of cells? Perhaps a better approach
+would be to come up with a way to assign *all* of these to the first cell,
+because there is *always* at least one cell. The way to achieve this would
+be to replace the use of FE_DGQ(0) (an element with exactly one degree of
+freedom) by an element that has more than one -- and in particular exactly
+the right number of degrees of freedom -- and that can be used on the first
+cell; all other cells would then use FE_Nothing. At the time of writing
+this program, there is no element class that can easily be given a *specific*
+number of degrees of freedom, but it would not be very difficult to write
+such a class.
diff --git a/examples/step-93/doc/tooltip b/examples/step-93/doc/tooltip
new file mode 100644 (file)
index 0000000..9a059f3
--- /dev/null
@@ -0,0 +1 @@
+Working with nonlocal DoFs. Solving an optimal control problem with constrained input.
diff --git a/examples/step-93/step-93.cc b/examples/step-93/step-93.cc
new file mode 100644 (file)
index 0000000..5057f4f
--- /dev/null
@@ -0,0 +1,821 @@
+/* ------------------------------------------------------------------------
+ *
+ * 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>  &center,
+                   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> &center, 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> &center,
+    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;
+}
diff --git a/tests/examples/step-93.diff b/tests/examples/step-93.diff
new file mode 100644 (file)
index 0000000..a2d4bc2
--- /dev/null
@@ -0,0 +1,33 @@
+diff --git a/examples/step-93/step-93.cc b/examples/step-93/step-93.cc
+index 2e7944c2b7..a0aac6263b 100644
+--- a/examples/step-93/step-93.cc
++++ b/examples/step-93/step-93.cc
+@@ -374,7 +374,7 @@ namespace Step93
+   void Step93<dim>::make_grid()
+   {
+     GridGenerator::hyper_cube(triangulation, -1, 1);
+-    triangulation.refine_global(7);
++    triangulation.refine_global(3);
+     std::cout << "Number of active cells: " << triangulation.n_active_cells()
+               << std::endl;
+@@ -655,7 +655,7 @@ namespace Step93
+   {
+     // Notice that we also time how long this process takes.
+     std::cout << "Beginning solve" << std::endl;
+-    Timer timer;
++    //Timer timer;
+     SolverControl solver_control(5000000, 1e-6 * system_rhs.l2_norm());
+     SolverCG<Vector<double>> solver(solver_control);
+@@ -665,8 +665,8 @@ namespace Step93
+     /* SparseDirectUMFPACK direct_solver;
+     direct_solver.initialize(system_matrix);
+     direct_solver.vmult(solution, system_rhs); */
+-    timer.stop();
+-    std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl;
++    //timer.stop();
++    //std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl;
+   }
+   // The output_results() function is a bit more robust for this program than

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.