deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-93 tutorial program

This tutorial depends on step-4, step-8.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

This program was contributed by Sam Scheuerman and Wolfgang Bangerth.

Introduction

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

\begin{align*} -\Delta u =& f \text{ in }\Omega\\ u =&0 \text{ on } \partial\Omega, \end{align*}

which is as close as possible to a desired solution \(\overline{u}\). We can formulate this as an optimization problem

\begin{align*} \min_{u, f} & \frac 12 ||u-\overline{u}||^2\\ \text{s.t.} & -\Delta u = f. \end{align*}

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:

\begin{align*} \min_{u, C^k} & \frac 12 ||u-\overline{u}||^2\\ \text{s.t.} & -\Delta u = \sum_k C^k f_k. \end{align*}

To write this optimization problem in a more compact form, we introduce a Lagrangian

\begin{align*} L = \frac{1}{2}||u - \overline{u}||^2 - \left(\lambda, \left(-\Delta u - \sum_k C^k f_k\right)\right), \end{align*}

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

\begin{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]. \end{align*}

where

\begin{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. \end{align*}

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:

Visualization of four characteristic functions together on the same plot

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:

\begin{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. \end{align*}

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\):

\begin{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.\\ \end{align*}

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

\begin{align*} \mathcal{M}U - \mathcal{U} - \mathcal{N}^T\Lambda =& 0\\ -\mathcal{N}U + \mathcal{F}^TC =& 0\\ \mathcal{F}\Lambda =& 0. \end{align*}

Moving known quantities to the right-hand-side and combinging the matrix equations gives us an idealized system to solve:

\begin{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). \end{align*}

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.

What is the problem?

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 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.

How do we address this?

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.

What concrete steps do we need to take?

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.

The commented program

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_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;
 
 
Definition hp.h:117

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;
  };
 
Definition point.h:113

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;
  }
 
 
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471

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;
 
  DoFHandler<dim> dof_handler;
 
  hp::FECollection<dim> fe_collection;
  hp::QCollection<dim> quadrature_collection;
 
 
  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;
  };
 
 
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

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(
  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(2), 2, FE_DGQ<dim>(0), 1));
 
Definition fe_q.h:554

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:
  break;
  }
  }
 
#define DEAL_II_ASSERT_UNREACHABLE()

Step93::make_grid()

The make_grid() function makes a hypercube grid, see step-4:

  template <int dim>
  void Step93<dim>::make_grid()
  {
  triangulation.refine_global(7);
 
  std::cout << "Number of active cells: " << triangulation.n_active_cells()
  << std::endl;
  }
 
 
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)

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);
 
std::size_t size
Definition mpi.cc:749

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 =
  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;
 
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})

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;
 
 
IndexSet extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask)
Definition dof_tools.cc:416

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);
  0,
  constraints);
  constraints.close();
 
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

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);
 
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)

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,
 
  const FEValuesExtractors::Scalar lambda(1);
 
@ update_quadrature_points
Transformed quadrature points.

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]);
  }
  }
  }
  }
  }
 
const FEValues< dim, spacedim > & get_present_fe_values() const
Point< 2 > first
Definition grid_out.cc:4632

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());
  }
 
 

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,
 
  FullMatrix<double> cell_matrix;
  Vector<double> cell_rhs;
 
  std::vector<types::global_dof_index> local_dof_indices;
 
  const FEValuesExtractors::Scalar lambda(1);
 
  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);
 
 
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.

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);
  }
  }
 
 

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);
void initialize(const SparsityPattern &sparsity_pattern)
void vmult(Vector< double > &dst, const Vector< double > &src) const
  template <int dim>
  void Step93<dim>::solve()
  {
  Timer timer;
  timer.start();
 
  std::cout << "Beginning solve..." << std::endl;
 
  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;
  }
 
 
Definition timer.h:117
void start()
Definition timer.cc:176

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>
 
  DataOut<dim> data_out;
  data_out.add_data_vector(dof_handler,
  solution,
  solution_names,
  interpretation);
 
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_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");
 
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)

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);
  }
 
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466

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;
  }
  }
 
 

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
 
 

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;
  }

Results

When you run the program with the step target function (in 2D), the output looks something like this:

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

When you run the program with the Gaussian target function (in 2D), the output should look like this:

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

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:

Step Target Function
Gaussian Target Function
Step target function Gaussian target function

After solving the Lagrangian system, we arrive at solutions \(U_\text{step}\) and \(U_\text{gauss}\) that look like this:

\(U_\text{step}\)
\(U_\text{gauss}\)
Solution for step shaped target function Solution for Gaussian target function

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.

Heat plate settings for matching step function
Heat plate settings for matching Gaussian
Heat plate settings for matching step function Heat plate settings for matching Gaussian

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:

\(\Lambda_\text{step}\)
\(\Lambda_\text{gauss}\)
Lagrange multiplier for step target function Lagrange multiplier for Gaussian target function

Possibilities for extensions

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. For block preconditioners, the key realizations is that the blocks of the system matrix

    \begin{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) \end{align*}

    can often individually be solved with quite efficiently; for example, \(\mathcal{M}\) is a mass matrix that is easily solved with using a CG iteration, and \(\mathcal N\) is a Laplace matrix for which CG with a geometric or algebraic multigrid preconditioner is very effective. Using the ideas of step-20 and step-22, we should then create a \(3\times 3\) block preconditioner in which some blocks correspond to the inverses of \(\mathcal{M}\) or \(\mathcal{N}\), or some kind of Schur complement. A starting point for this kind of consideration is [17] .
  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.

The plain program

/* ------------------------------------------------------------------------
*
* 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.
*/
#include <fstream>
#include <iostream>
namespace Step93
{
using namespace dealii;
template <int dim>
class TargetFunction : public Function<dim>
{
public:
TargetFunction()
: Function<dim>(1)
, center()
, radius(0.1)
{}
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;
};
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;
}
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;
}
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;
DoFHandler<dim> dof_handler;
hp::FECollection<dim> fe_collection;
hp::QCollection<dim> quadrature_collection;
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;
};
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)))
{
fe_collection.push_back(
fe_collection.push_back(FESystem<dim>(FE_Q<dim>(2), 2, FE_DGQ<dim>(0), 1));
quadrature_collection.push_back(QGauss<dim>(3));
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:
break;
}
}
template <int dim>
void Step93<dim>::make_grid()
{
triangulation.refine_global(7);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
template <int dim>
void Step93<dim>::setup_system()
{
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);
const std::vector<types::global_dof_index> dofs_per_component =
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;
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;
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
0,
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
hp::FEValues<dim> hp_fe_values(fe_collection,
quadrature_collection,
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]);
}
}
}
}
}
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void Step93<dim>::assemble_system()
{
hp::FEValues<dim> hp_fe_values(fe_collection,
quadrature_collection,
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
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_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);
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 (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);
}
}
}
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
template <int dim>
void Step93<dim>::solve()
{
Timer timer;
timer.start();
std::cout << "Beginning solve..." << std::endl;
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;
}
template <int dim>
void Step93<dim>::output_results() const
{
const std::vector<std::string> solution_names = {"u", "lambda", "c"};
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
DataOut<dim> data_out;
data_out.add_data_vector(dof_handler,
solution,
solution_names,
interpretation);
DoFHandler<dim> new_dof_handler(triangulation);
const FE_Q<dim> new_fe(2);
new_dof_handler.distribute_dofs(new_fe);
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");
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);
}
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);
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;
}
}
template <int dim>
void Step93<dim>::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
} // namespace Step93
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;
}
void write_vtu(std::ostream &out) const
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
const unsigned int dofs_per_cell
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
double wall_time() const
Definition timer.cc:262
double stop()
Definition timer.cc:193
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)