step-79.h

Go to the documentation of this file.

 
 *   
 *         for (unsigned int i = 0; i < dim; ++i)
 *           {
 *             for (unsigned int k = 0; k < dim; ++k)
 *               {
 *                 coupling[displacement<dim> + i]
 *                         [displacement_multiplier<dim> + k] = DoFTools::always;
 *                 coupling[displacement_multiplier<dim> + k]
 *                         [displacement<dim> + i] = DoFTools::always;
 *               }
 *           }
 *   
 *         /* Coupling for slack variables */
 *         coupling[density_lower_slack<dim>][density_lower_slack<dim>] =
 *           DoFTools::always;
 *         coupling[density_lower_slack<dim>][density_upper_slack<dim>] =
 *           DoFTools::always;
 *         coupling[density_upper_slack<dim>][density_lower_slack<dim>] =
 *           DoFTools::always;
 *   
 *         coupling[density_lower_slack_multiplier<dim>]
 *                 [density_lower_slack_multiplier<dim>] = DoFTools::always;
 *         coupling[density_lower_slack_multiplier<dim>]
 *                 [density_upper_slack_multiplier<dim>] = DoFTools::always;
 *         coupling[density_upper_slack_multiplier<dim>]
 *                 [density_lower_slack_multiplier<dim>] = DoFTools::always;
 *       }
 *   
 * @endcode
 * 
 * Before we can create the sparsity pattern, we also have to
 * set up constraints. Since this program does not adaptively
 * refine the mesh, the only constraint we have is one that
 * couples all density variables to enforce the volume
 * constraint. This will ultimately lead to a dense sub-block
 * of the matrix, but there is little we can do about that.
 * 
 * @code
 *       const ComponentMask density_mask =
 *         fe.component_mask(ValueExtractors::densities<dim>);
 *       const IndexSet density_dofs =
 *         DoFTools::extract_dofs(dof_handler, density_mask);
 *   
 *       types::global_dof_index last_density_dof =
 *         density_dofs.nth_index_in_set(density_dofs.n_elements() - 1);
 *       constraints.clear();
 *       constraints.add_line(last_density_dof);
 *       for (unsigned int i = 0; i < density_dofs.n_elements() - 1; ++i)
 *         constraints.add_entry(last_density_dof,
 *                               density_dofs.nth_index_in_set(i),
 *                               -1);
 *       constraints.set_inhomogeneity(last_density_dof, 0);
 *   
 *       constraints.close();
 *   
 * @endcode
 * 
 * We can now finally create the sparsity pattern for the
 * matrix, taking into account which variables couple with
 * which other variables, and the constraints we have on the
 * density.
 * 
 * @code
 *       DoFTools::make_sparsity_pattern(dof_handler, coupling, dsp, constraints);
 *   
 * @endcode
 * 
 * The only part of the matrix we have not dealt with is the
 * filter matrix and its transpose. These are non-local
 * (integral) operators for which deal.II does not currently
 * have functions. What we will ultimately need to do is go
 * over all cells and couple the unfiltered density on this
 * cell to all filtered densities of neighboring cells that
 * are less than a threshold distance away, and the other way
 * around; for the moment, we are only concerned with building
 * the sparsity pattern that would correspond to this kind of
 * matrix, so we perform the equivalent loop and where later
 * on we would write into an entry of the matrix, we now
 * simply add an entry to the sparsity matrix:
 * 
 * @code
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           const unsigned int i = cell->active_cell_index();
 *           for (const auto &check_cell : find_relevant_neighbors(cell))
 *             {
 *               const double distance =
 *                 cell->center().distance(check_cell->center());
 *               if (distance < filter_r)
 *                 {
 *                   dsp
 *                     .block(SolutionBlocks::unfiltered_density,
 *                            SolutionBlocks::unfiltered_density_multiplier)
 *                     .add(i, check_cell->active_cell_index());
 *                   dsp
 *                     .block(SolutionBlocks::unfiltered_density_multiplier,
 *                            SolutionBlocks::unfiltered_density)
 *                     .add(i, check_cell->active_cell_index());
 *                 }
 *             }
 *         }
 *   
 * @endcode
 * 
 * Having so generated the "dynamic" sparsity pattern, we can
 * finally copy it to the structure that is used to associate
 * matrices with a sparsity pattern. Because the sparsity
 * pattern is large and complex, we also output it into a file
 * of its own for visualization purposes -- in other words,
 * for "visual debugging".
 * 
 * @code
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       std::ofstream out("sparsity.plt");
 *       sparsity_pattern.print_gnuplot(out);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *   
 *   
 * @endcode
 * 
 * What is left is to correctly size the various vectors and
 * their blocks, as well as setting initial guesses for some
 * of the components of the (nonlinear) solution vector. We
 * here use the symbolic component names for individual blocks
 * of the solution vector and, for brevity, use the same trick
 * with `using namespace` as above:
 * 
 * @code
 *       nonlinear_solution.reinit(block_sizes);
 *       system_rhs.reinit(block_sizes);
 *   
 *       {
 *         using namespace SolutionBlocks;
 *         nonlinear_solution.block(density).add(density_ratio);
 *         nonlinear_solution.block(unfiltered_density).add(density_ratio);
 *         nonlinear_solution.block(unfiltered_density_multiplier)
 *           .add(density_ratio);
 *         nonlinear_solution.block(density_lower_slack).add(density_ratio);
 *         nonlinear_solution.block(density_lower_slack_multiplier).add(50);
 *         nonlinear_solution.block(density_upper_slack).add(1 - density_ratio);
 *         nonlinear_solution.block(density_upper_slack_multiplier).add(50);
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Creatingthefiltermatrix"></a> 
 * <h3>Creating the filter matrix</h3>
 * 
 
 * 
 * Next up, a function that is used once at the beginning of the
 * program: It creates a matrix @f$H@f$ so that the filtered density
 * vector equals @f$H@f$ times the unfiltered density.  The creation
 * of this matrix is non-trivial, and it is used in every
 * iteration, and so rather than reforming it as we do with the
 * Newton matrix, it is made only once and stored separately.
 *   
 
 * 
 * The way this matrix is computed follows the outline used above
 * already to form its sparsity pattern. We repeat this process here
 * for the sparsity pattern of this separately formed matrix, and
 * then actually build the matrix itself. You may want to check the
 * definition of this matrix in the introduction to this program.
 * 
 * @code
 *     template <int dim>
 *     void SANDTopOpt<dim>::setup_filter_matrix()
 *     {
 * @endcode
 * 
 * The sparsity pattern of the filter has already been determined
 * and implemented in the setup_system() function. We copy the
 * structure from the appropriate block and use it again here.
 * 
 
 * 
 * 
 * @code
 *       filter_sparsity_pattern.copy_from(
 *         sparsity_pattern.block(SolutionBlocks::unfiltered_density,
 *                                SolutionBlocks::unfiltered_density_multiplier));
 *       filter_matrix.reinit(filter_sparsity_pattern);
 *   
 * @endcode
 * 
 * Having so built the sparsity pattern, now we re-do all of
 * these loops to actually compute the necessary values of the
 * matrix entries:
 * 
 
 * 
 * 
 * @code
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           const unsigned int i = cell->active_cell_index();
 *           for (const auto &check_cell : find_relevant_neighbors(cell))
 *             {
 *               const double distance =
 *                 cell->center().distance(check_cell->center());
 *               if (distance < filter_r)
 *                 {
 *                   filter_matrix.add(i,
 *                                     check_cell->active_cell_index(),
 *                                     filter_r - distance);
 *                   
 *                 }
 *             }
 *         }
 *   
 * @endcode
 * 
 * The final step is to normalize the matrix so that for each
 * row, the sum of entries equals one.
 * 
 * @code
 *       for (unsigned int i = 0; i < filter_matrix.m(); ++i)
 *         {
 *           double denominator = 0;
 *           for (SparseMatrix<double>::iterator iter = filter_matrix.begin(i);
 *                iter != filter_matrix.end(i);
 *                iter++)
 *             denominator = denominator + iter->value();
 *           for (SparseMatrix<double>::iterator iter = filter_matrix.begin(i);
 *                iter != filter_matrix.end(i);
 *                iter++)
 *             iter->value() = iter->value() / denominator;
 *         }
 *     }
 *   
 * @endcode
 * 
 * This function is used for building the filter matrix. We create a set of
 * all the cell iterators within a certain radius of the cell that is input.
 * These are the neighboring cells that will be relevant for the filter.
 * 
 * @code
 *     template <int dim>
 *     std::set<typename Triangulation<dim>::cell_iterator>
 *     SANDTopOpt<dim>::find_relevant_neighbors(
 *       typename Triangulation<dim>::cell_iterator cell) const
 *     {
 *       std::set<unsigned int>                               neighbor_ids;
 *       std::set<typename Triangulation<dim>::cell_iterator> cells_to_check;
 *   
 *       neighbor_ids.insert(cell->active_cell_index());
 *       cells_to_check.insert(cell);
 *   
 *       bool new_neighbors_found;
 *       do
 *         {
 *           new_neighbors_found = false;
 *           for (const auto &check_cell :
 *                std::vector<typename Triangulation<dim>::cell_iterator>(
 *                  cells_to_check.begin(), cells_to_check.end()))
 *             {
 *               for (const auto n : check_cell->face_indices())
 *                 {
 *                   if (!(check_cell->face(n)->at_boundary()))
 *                     {
 *                       const auto & neighbor = check_cell->neighbor(n);
 *                       const double distance =
 *                         cell->center().distance(neighbor->center());
 *                       if ((distance < filter_r) &&
 *                           !(neighbor_ids.count(neighbor->active_cell_index())))
 *                         {
 *                           cells_to_check.insert(neighbor);
 *                           neighbor_ids.insert(neighbor->active_cell_index());
 *                           new_neighbors_found = true;
 *                         }
 *                     }
 *                 }
 *             }
 *         }
 *       while (new_neighbors_found);
 *       return cells_to_check;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="AssemblingtheNewtonmatrix"></a> 
 * <h3>Assembling the Newton matrix</h3>
 * 
 
 * 
 * Whereas the setup_filter_matrix function built a matrix that is the same as
 * long as the mesh does not change (which we don't do anyway in
 * this program), the next function builds the matrix to be solved
 * in each iteration. This is where the magic happens. The components
 * of the system of linear equations describing Newton's method for
 * finding the solution of the KKT conditions are implemented here.
 *   
 
 * 
 * The top of the function is as in most of these functions and just
 * sets up all sorts of variables necessary for the actual assembly,
 * including a whole bunch of extractors. The entire set up should
 * look familiar, though somewhat lengthier, if you've previously
 * looked at @ref step_22 "step-22".
 * 
 * @code
 *     template <int dim>
 *     void SANDTopOpt<dim>::assemble_system()
 *     {
 *       TimerOutput::Scope t(timer, "assembly");
 *   
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *   
 *       MappingQ<dim>     mapping(1);
 *       QGauss<dim>       quadrature_formula(fe.degree + 1);
 *       QGauss<dim - 1>   face_quadrature_formula(fe.degree + 1);
 *       FEValues<dim>     fe_values(mapping,
 *                               fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       FEFaceValues<dim> fe_face_values(mapping,
 *                                        fe,
 *                                        face_quadrature_formula,
 *                                        update_values | update_quadrature_points |
 *                                          update_normal_vectors |
 *                                          update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     dummy_cell_rhs(dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       std::vector<double>                    lambda_values(n_q_points);
 *       std::vector<double>                    mu_values(n_q_points);
 *       const Functions::ConstantFunction<dim> lambda(1.);
 *       const Functions::ConstantFunction<dim> mu(1.);
 *       std::vector<Tensor<1, dim>>            rhs_values(n_q_points);
 *   
 * @endcode
 * 
 * At this point, we apply the filter to the unfiltered
 * density, and apply the adjoint (transpose) operation to the
 * unfiltered density multiplier, both to the current best
 * guess for the nonlinear solution. We use this later to tell
 * us how far off our filtered density is from the filter
 * applied to the unfiltered density. That is because while at
 * the solution of the nonlinear problem, we have
 * @f$\rho=H\varrho@f$, but at intermediate iterations, we in
 * general have @f$\rho^k\neq H\varrho^k@f$ and the "residual"
 * @f$\rho^k-H\varrho^k@f$ will then appear as the right hand side
 * of one of the Newton update equations that we compute
 * below.
 * 
 * @code
 *       BlockVector<double> filtered_unfiltered_density_solution =
 *         nonlinear_solution;
 *       BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
 *         nonlinear_solution;
 *   
 *       filter_matrix.vmult(filtered_unfiltered_density_solution.block(
 *                             SolutionBlocks::unfiltered_density),
 *                           nonlinear_solution.block(
 *                             SolutionBlocks::unfiltered_density));
 *       filter_matrix.Tvmult(
 *         filter_adjoint_unfiltered_density_multiplier_solution.block(
 *           SolutionBlocks::unfiltered_density_multiplier),
 *         nonlinear_solution.block(SolutionBlocks::unfiltered_density_multiplier));
 *   
 *   
 *       std::vector<double>                  old_density_values(n_q_points);
 *       std::vector<Tensor<1, dim>>          old_displacement_values(n_q_points);
 *       std::vector<double>                  old_displacement_divs(n_q_points);
 *       std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
 *       std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
 *       std::vector<double>         old_displacement_multiplier_divs(n_q_points);
 *       std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
 *         n_q_points);
 *       std::vector<double> old_lower_slack_multiplier_values(n_q_points);
 *       std::vector<double> old_upper_slack_multiplier_values(n_q_points);
 *       std::vector<double> old_lower_slack_values(n_q_points);
 *       std::vector<double> old_upper_slack_values(n_q_points);
 *       std::vector<double> old_unfiltered_density_values(n_q_points);
 *       std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
 *       std::vector<double> filtered_unfiltered_density_values(n_q_points);
 *       std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
 *         n_q_points);
 *   
 *       using namespace ValueExtractors;
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_matrix = 0;
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           fe_values.reinit(cell);
 *   
 *           lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
 *           mu.value_list(fe_values.get_quadrature_points(), mu_values);
 *   
 * @endcode
 * 
 * As part of the construction of our system matrix, we need to
 * retrieve values from our current guess at the solution.
 * The following lines of code retrieve the needed values.
 * 
 * @code
 *           fe_values[densities<dim>].get_function_values(nonlinear_solution,
 *                                                         old_density_values);
 *           fe_values[displacements<dim>].get_function_values(
 *             nonlinear_solution, old_displacement_values);
 *           fe_values[displacements<dim>].get_function_divergences(
 *             nonlinear_solution, old_displacement_divs);
 *           fe_values[displacements<dim>].get_function_symmetric_gradients(
 *             nonlinear_solution, old_displacement_symmgrads);
 *           fe_values[displacement_multipliers<dim>].get_function_values(
 *             nonlinear_solution, old_displacement_multiplier_values);
 *           fe_values[displacement_multipliers<dim>].get_function_divergences(
 *             nonlinear_solution, old_displacement_multiplier_divs);
 *           fe_values[displacement_multipliers<dim>]
 *             .get_function_symmetric_gradients(
 *               nonlinear_solution, old_displacement_multiplier_symmgrads);
 *           fe_values[density_lower_slacks<dim>].get_function_values(
 *             nonlinear_solution, old_lower_slack_values);
 *           fe_values[density_lower_slack_multipliers<dim>].get_function_values(
 *             nonlinear_solution, old_lower_slack_multiplier_values);
 *           fe_values[density_upper_slacks<dim>].get_function_values(
 *             nonlinear_solution, old_upper_slack_values);
 *           fe_values[density_upper_slack_multipliers<dim>].get_function_values(
 *             nonlinear_solution, old_upper_slack_multiplier_values);
 *           fe_values[unfiltered_densities<dim>].get_function_values(
 *             nonlinear_solution, old_unfiltered_density_values);
 *           fe_values[unfiltered_density_multipliers<dim>].get_function_values(
 *             nonlinear_solution, old_unfiltered_density_multiplier_values);
 *           fe_values[unfiltered_densities<dim>].get_function_values(
 *             filtered_unfiltered_density_solution,
 *             filtered_unfiltered_density_values);
 *           fe_values[unfiltered_density_multipliers<dim>].get_function_values(
 *             filter_adjoint_unfiltered_density_multiplier_solution,
 *             filter_adjoint_unfiltered_density_multiplier_values);
 *   
 *           for (const auto q_point : fe_values.quadrature_point_indices())
 *             {
 * @endcode
 * 
 * We need several more values corresponding to the test functions
 * coming from the first derivatives taken from the Lagrangian,
 * that is the @f$d_{\bullet}@f$ functions. These are calculated here:
 * 
 * @code
 *               for (const auto i : fe_values.dof_indices())
 *                 {
 *                   const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
 *                     fe_values[displacements<dim>].symmetric_gradient(i, q_point);
 *                   const double displacement_phi_i_div =
 *                     fe_values[displacements<dim>].divergence(i, q_point);
 *   
 *                   const SymmetricTensor<2, dim>
 *                     displacement_multiplier_phi_i_symmgrad =
 *                       fe_values[displacement_multipliers<dim>].symmetric_gradient(
 *                         i, q_point);
 *                   const double displacement_multiplier_phi_i_div =
 *                     fe_values[displacement_multipliers<dim>].divergence(i,
 *                                                                         q_point);
 *   
 *                   const double density_phi_i =
 *                     fe_values[densities<dim>].value(i, q_point);
 *                   const double unfiltered_density_phi_i =
 *                     fe_values[unfiltered_densities<dim>].value(i, q_point);
 *                   const double unfiltered_density_multiplier_phi_i =
 *                     fe_values[unfiltered_density_multipliers<dim>].value(i,
 *                                                                          q_point);
 *   
 *                   const double lower_slack_multiplier_phi_i =
 *                     fe_values[density_lower_slack_multipliers<dim>].value(
 *                       i, q_point);
 *   
 *                   const double lower_slack_phi_i =
 *                     fe_values[density_lower_slacks<dim>].value(i, q_point);
 *   
 *                   const double upper_slack_phi_i =
 *                     fe_values[density_upper_slacks<dim>].value(i, q_point);
 *   
 *                   const double upper_slack_multiplier_phi_i =
 *                     fe_values[density_upper_slack_multipliers<dim>].value(
 *                       i, q_point);
 *   
 *   
 *                   for (const auto j : fe_values.dof_indices())
 *                     {
 * @endcode
 * 
 * Finally, we need values that come from the second round
 * of derivatives taken from the Lagrangian,
 * the @f$c_{\bullet}@f$ functions. These are calculated here:
 * 
 * @code
 *                       const SymmetricTensor<2, dim> displacement_phi_j_symmgrad =
 *                         fe_values[displacements<dim>].symmetric_gradient(j,
 *                                                                          q_point);
 *                       const double displacement_phi_j_div =
 *                         fe_values[displacements<dim>].divergence(j, q_point);
 *   
 *                       const SymmetricTensor<2, dim>
 *                         displacement_multiplier_phi_j_symmgrad =
 *                           fe_values[displacement_multipliers<dim>]
 *                             .symmetric_gradient(j, q_point);
 *                       const double displacement_multiplier_phi_j_div =
 *                         fe_values[displacement_multipliers<dim>].divergence(
 *                           j, q_point);
 *   
 *                       const double density_phi_j =
 *                         fe_values[densities<dim>].value(j, q_point);
 *   
 *                       const double unfiltered_density_phi_j =
 *                         fe_values[unfiltered_densities<dim>].value(j, q_point);
 *                       const double unfiltered_density_multiplier_phi_j =
 *                         fe_values[unfiltered_density_multipliers<dim>].value(
 *                           j, q_point);
 *   
 *   
 *                       const double lower_slack_phi_j =
 *                         fe_values[density_lower_slacks<dim>].value(j, q_point);
 *   
 *                       const double upper_slack_phi_j =
 *                         fe_values[density_upper_slacks<dim>].value(j, q_point);
 *   
 *                       const double lower_slack_multiplier_phi_j =
 *                         fe_values[density_lower_slack_multipliers<dim>].value(
 *                           j, q_point);
 *   
 *                       const double upper_slack_multiplier_phi_j =
 *                         fe_values[density_upper_slack_multipliers<dim>].value(
 *                           j, q_point);
 *   
 * @endcode
 * 
 * This is where the actual work starts. In
 * the following, we will build all of the
 * terms of the matrix -- they are numerous
 * and not entirely self-explanatory, also
 * depending on the previous solutions and its
 * derivatives (which we have already
 * evaluated above and put into the variables
 * called `old_*`). To understand what each of
 * these terms corresponds to, you will want
 * to look at the explicit form of these terms
 * in the introduction above.
 *                     
 
 * 
 * The right hand sides of the equations being
 * driven to 0 give all the KKT conditions
 * for finding a local minimum -- the descriptions of what
 * each individual equation are given with the computations
 * of the right hand side.
 * 
 
 * 
 * 
 * @code
 *                       /* Equation 1 */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (
 *   
 *                           -density_phi_i * unfiltered_density_multiplier_phi_j
 *   
 *                           + density_penalty_exponent *
 *                               (density_penalty_exponent - 1) *
 *                               std::pow(old_density_values[q_point],
 *                                        density_penalty_exponent - 2) *
 *                               density_phi_i * density_phi_j *
 *                               (old_displacement_multiplier_divs[q_point] *
 *                                  old_displacement_divs[q_point] *
 *                                  lambda_values[q_point] +
 *                                2 * mu_values[q_point] *
 *                                  (old_displacement_symmgrads[q_point] *
 *                                   old_displacement_multiplier_symmgrads[q_point]))
 *   
 *                           + density_penalty_exponent *
 *                               std::pow(old_density_values[q_point],
 *                                        density_penalty_exponent - 1) *
 *                               density_phi_i *
 *                               (displacement_multiplier_phi_j_div *
 *                                  old_displacement_divs[q_point] *
 *                                  lambda_values[q_point] +
 *                                2 * mu_values[q_point] *
 *                                  (old_displacement_symmgrads[q_point] *
 *                                   displacement_multiplier_phi_j_symmgrad))
 *   
 *                           + density_penalty_exponent *
 *                               std::pow(old_density_values[q_point],
 *                                        density_penalty_exponent - 1) *
 *                               density_phi_i *
 *                               (displacement_phi_j_div *
 *                                  old_displacement_multiplier_divs[q_point] *
 *                                  lambda_values[q_point] +
 *                                2 * mu_values[q_point] *
 *                                  (old_displacement_multiplier_symmgrads[q_point] *
 *                                   displacement_phi_j_symmgrad)));
 *   
 *                       /* Equation 2 */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (density_penalty_exponent *
 *                            std::pow(old_density_values[q_point],
 *                                     density_penalty_exponent - 1) *
 *                            density_phi_j *
 *                            (old_displacement_multiplier_divs[q_point] *
 *                               displacement_phi_i_div * lambda_values[q_point] +
 *                             2 * mu_values[q_point] *
 *                               (old_displacement_multiplier_symmgrads[q_point] *
 *                                displacement_phi_i_symmgrad))
 *   
 *                          + std::pow(old_density_values[q_point],
 *                                     density_penalty_exponent) *
 *                              (displacement_multiplier_phi_j_div *
 *                                 displacement_phi_i_div * lambda_values[q_point] +
 *                               2 * mu_values[q_point] *
 *                                 (displacement_multiplier_phi_j_symmgrad *
 *                                  displacement_phi_i_symmgrad))
 *   
 *                         );
 *   
 *                       /* Equation 3, which has to do with the filter and which is
 *                        * calculated elsewhere. */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (-1 * unfiltered_density_phi_i *
 *                            lower_slack_multiplier_phi_j +
 *                          unfiltered_density_phi_i * upper_slack_multiplier_phi_j);
 *   
 *   
 *                       /* Equation 4: Primal feasibility */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (
 *   
 *                           density_penalty_exponent *
 *                             std::pow(old_density_values[q_point],
 *                                      density_penalty_exponent - 1) *
 *                             density_phi_j *
 *                             (old_displacement_divs[q_point] *
 *                                displacement_multiplier_phi_i_div *
 *                                lambda_values[q_point] +
 *                              2 * mu_values[q_point] *
 *                                (old_displacement_symmgrads[q_point] *
 *                                 displacement_multiplier_phi_i_symmgrad))
 *   
 *                           + std::pow(old_density_values[q_point],
 *                                      density_penalty_exponent) *
 *                               (displacement_phi_j_div *
 *                                  displacement_multiplier_phi_i_div *
 *                                  lambda_values[q_point] +
 *                                2 * mu_values[q_point] *
 *                                  (displacement_phi_j_symmgrad *
 *                                   displacement_multiplier_phi_i_symmgrad)));
 *   
 *                       /* Equation 5: Primal feasibility */
 *                       cell_matrix(i, j) +=
 *                         -1 * fe_values.JxW(q_point) *
 *                         lower_slack_multiplier_phi_i *
 *                         (unfiltered_density_phi_j - lower_slack_phi_j);
 *   
 *                       /* Equation 6: Primal feasibility */
 *                       cell_matrix(i, j) +=
 *                         -1 * fe_values.JxW(q_point) *
 *                         upper_slack_multiplier_phi_i *
 *                         (-1 * unfiltered_density_phi_j - upper_slack_phi_j);
 *   
 *                       /* Equation 7: Primal feasibility - the part with the filter
 *                        * is added later */
 *                       cell_matrix(i, j) += -1 * fe_values.JxW(q_point) *
 *                                            unfiltered_density_multiplier_phi_i *
 *                                            (density_phi_j);
 *   
 *                       /* Equation 8: Complementary slackness */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (lower_slack_phi_i * lower_slack_multiplier_phi_j
 *   
 *                          + lower_slack_phi_i * lower_slack_phi_j *
 *                              old_lower_slack_multiplier_values[q_point] /
 *                              old_lower_slack_values[q_point]);
 *   
 *                       /* Equation 9: Complementary slackness */
 *                       cell_matrix(i, j) +=
 *                         fe_values.JxW(q_point) *
 *                         (upper_slack_phi_i * upper_slack_multiplier_phi_j
 *   
 *   
 *                          + upper_slack_phi_i * upper_slack_phi_j *
 *                              old_upper_slack_multiplier_values[q_point] /
 *                              old_upper_slack_values[q_point]);
 *                     }
 *                 }
 *             }
 *   
 * @endcode
 * 
 * Now that we have everything assembled, all we have to
 * do is deal with the effect of (Dirichlet) boundary
 * conditions and other constraints. We incorporate the
 * former locally with just the contributions from the
 * current cell, and then let the AffineConstraint class
 * deal with the latter while copying contributions from
 * the current cell into the global linear system:
 * 
 * @code
 *           MatrixTools::local_apply_boundary_values(boundary_values,
 *                                                    local_dof_indices,
 *                                                    cell_matrix,
 *                                                    dummy_cell_rhs,
 *                                                    true);
 *   
 *           constraints.distribute_local_to_global(cell_matrix,
 *                                                  local_dof_indices,
 *                                                  system_matrix);
 *         }
 *   
 * @endcode
 * 
 * Having accumulated all of the terms that belong
 * into the Newton matrix, we now also have to
 * compute the terms for the right hand side
 * (i.e., the negative residual). We already do this
 * in another function, and so we call that here:
 * 
 * @code
 *       system_rhs = calculate_test_rhs(nonlinear_solution);
 *   
 * @endcode
 * 
 * Here we use the filter matrix we have already
 * constructed. We only need to integrate this filter applied
 * to test functions, which are piecewise constant, and so the
 * integration becomes a simple multiplication by the measure
 * of the cell.  Iterating over the pre-made filter matrix
 * allows us to use the information about which cells are in
 * or out of the filter without repeatedly checking neighbor
 * cells again.
 * 
 * @code
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           const unsigned int i = cell->active_cell_index();
 *           for (typename SparseMatrix<double>::iterator iter =
 *                  filter_matrix.begin(i);
 *                iter != filter_matrix.end(i);
 *                ++iter)
 *             {
 *               const unsigned int j     = iter->column();
 *               const double       value = iter->value() * cell->measure();
 *   
 *               system_matrix
 *                 .block(SolutionBlocks::unfiltered_density_multiplier,
 *                        SolutionBlocks::unfiltered_density)
 *                 .add(i, j, value);
 *               system_matrix
 *                 .block(SolutionBlocks::unfiltered_density,
 *                        SolutionBlocks::unfiltered_density_multiplier)
 *                 .add(j, i, value);
 *             }
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="SolvingtheNewtonlinearsystem"></a> 
 * <h3>Solving the Newton linear system</h3>
 * 
 
 * 
 * 
 
 * 
 * We will need to solve a linear system in each iteration. We use
 * a direct solver, for now -- this is clearly not an efficient
 * choice for a matrix that has so many non-zeroes, and it will
 * not scale to anything interesting. For "real" applications, we
 * will need an iterative solver but the complexity of the system
 * means that an iterative solver algorithm will take a good deal
 * of work. Because this is not the focus of the current program,
 * we simply stick with the direct solver we have here -- the
 * function follows the same structure as used in @ref step_29 "step-29".
 * 
 * @code
 *     template <int dim>
 *     BlockVector<double> SANDTopOpt<dim>::solve()
 *     {
 *       TimerOutput::Scope t(timer, "solver");
 *   
 *       BlockVector<double> linear_solution;
 *       linear_solution.reinit(nonlinear_solution);
 *   
 *       SparseDirectUMFPACK A_direct;
 *       A_direct.initialize(system_matrix);
 *       A_direct.vmult(linear_solution, system_rhs);
 *   
 *       constraints.distribute(linear_solution);
 *   
 *       return linear_solution;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Detailsoftheoptimizationalgorithm"></a> 
 * <h3>Details of the optimization algorithm</h3>
 * 
 
 * 
 * The next several functions deal with specific parts of the
 * optimization algorithm, most notably with deciding whether the
 * direction computed by solving the linearized (Newton) system is
 * viable and, if so, how far we want to go in this direction.
 * 
 
 * 
 * 
 * <a name="Computingsteplengths"></a> 
 * <h4>Computing step lengths</h4>
 * 
 
 * 
 * We start with a function that does a binary search to figure
 * out the maximum step that meets the dual feasibility -- that
 * is, how far can we go so that @f$s>0@f$ and @f$z>0@f$. The function
 * returns a pair of values, one each for the @f$s@f$ and @f$z@f$ slack
 * variables.
 * 
 * @code
 *     template <int dim>
 *     std::pair<double, double> SANDTopOpt<dim>::calculate_max_step_size(
 *       const BlockVector<double> &state,
 *       const BlockVector<double> &step) const
 *     {
 *       double       fraction_to_boundary;
 *       const double min_fraction_to_boundary = .8;
 *       const double max_fraction_to_boundary = 1. - 1e-5;
 *   
 *       if (min_fraction_to_boundary < 1 - barrier_size)
 *         {
 *           if (1 - barrier_size < max_fraction_to_boundary)
 *             fraction_to_boundary = 1 - barrier_size;
 *           else
 *             fraction_to_boundary = max_fraction_to_boundary;
 *         }
 *       else
 *         fraction_to_boundary = min_fraction_to_boundary;
 *   
 *       double step_size_s_low  = 0;
 *       double step_size_z_low  = 0;
 *       double step_size_s_high = 1;
 *       double step_size_z_high = 1;
 *       double step_size_s, step_size_z;
 *   
 *       const int max_bisection_method_steps = 50;
 *       for (unsigned int k = 0; k < max_bisection_method_steps; ++k)
 *         {
 *           step_size_s = (step_size_s_low + step_size_s_high) / 2;
 *           step_size_z = (step_size_z_low + step_size_z_high) / 2;
 *   
 *           const BlockVector<double> state_test_s =
 *             (fraction_to_boundary * state) + (step_size_s * step);
 *   
 *           const BlockVector<double> state_test_z =
 *             (fraction_to_boundary * state) + (step_size_z * step);
 *   
 *           const bool accept_s =
 *             (state_test_s.block(SolutionBlocks::density_lower_slack)
 *                .is_non_negative()) &&
 *             (state_test_s.block(SolutionBlocks::density_upper_slack)
 *                .is_non_negative());
 *           const bool accept_z =
 *             (state_test_z.block(SolutionBlocks::density_lower_slack_multiplier)
 *                .is_non_negative()) &&
 *             (state_test_z.block(SolutionBlocks::density_upper_slack_multiplier)
 *                .is_non_negative());
 *   
 *           if (accept_s)
 *             step_size_s_low = step_size_s;
 *           else
 *             step_size_s_high = step_size_s;
 *   
 *           if (accept_z)
 *             step_size_z_low = step_size_z;
 *           else
 *             step_size_z_high = step_size_z;
 *         }
 *   
 *       return {step_size_s_low, step_size_z_low};
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Computingresiduals"></a> 
 * <h4>Computing residuals</h4>
 * 
 
 * 
 * The next function computes a right hand side vector linearized
 * around a "test solution vector" that we can use to look at the
 * magnitude of the KKT conditions.  This is then used for testing
 * the convergence before shrinking the barrier size, as well as in the
 * calculation of the @f$l_1@f$ merit.
 *   
 
 * 
 * The function is lengthy and complicated, but it is really just a
 * copy of the right hand side part of what the `assemble_system()`
 * function above did.
 * 
 * @code
 *     template <int dim>
 *     BlockVector<double> SANDTopOpt<dim>::calculate_test_rhs(
 *       const BlockVector<double> &test_solution) const
 *     {
 * @endcode
 * 
 * We first create a zero vector with size and blocking of system_rhs
 * 
 * @code
 *       BlockVector<double> test_rhs;
 *       test_rhs.reinit(system_rhs);
 *   
 *       MappingQ<dim>         mapping(1);
 *       const QGauss<dim>     quadrature_formula(fe.degree + 1);
 *       const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 *       FEValues<dim>         fe_values(mapping,
 *                               fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *       FEFaceValues<dim>     fe_face_values(mapping,
 *                                        fe,
 *                                        face_quadrature_formula,
 *                                        update_values | update_quadrature_points |
 *                                          update_normal_vectors |
 *                                          update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *       FullMatrix<double> dummy_cell_matrix(dofs_per_cell, dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       std::vector<double> lambda_values(n_q_points);
 *       std::vector<double> mu_values(n_q_points);
 *   
 *       const Functions::ConstantFunction<dim> lambda(1.), mu(1.);
 *       std::vector<Tensor<1, dim>>            rhs_values(n_q_points);
 *   
 *   
 *       BlockVector<double> filtered_unfiltered_density_solution = test_solution;
 *       BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
 *         test_solution;
 *       filtered_unfiltered_density_solution.block(
 *         SolutionBlocks::unfiltered_density) = 0;
 *       filter_adjoint_unfiltered_density_multiplier_solution.block(
 *         SolutionBlocks::unfiltered_density_multiplier) = 0;
 *   
 *       filter_matrix.vmult(filtered_unfiltered_density_solution.block(
 *                             SolutionBlocks::unfiltered_density),
 *                           test_solution.block(
 *                             SolutionBlocks::unfiltered_density));
 *       filter_matrix.Tvmult(
 *         filter_adjoint_unfiltered_density_multiplier_solution.block(
 *           SolutionBlocks::unfiltered_density_multiplier),
 *         test_solution.block(SolutionBlocks::unfiltered_density_multiplier));
 *   
 *   
 *       std::vector<double>                  old_density_values(n_q_points);
 *       std::vector<Tensor<1, dim>>          old_displacement_values(n_q_points);
 *       std::vector<double>                  old_displacement_divs(n_q_points);
 *       std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
 *       std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
 *       std::vector<double>         old_displacement_multiplier_divs(n_q_points);
 *       std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
 *         n_q_points);
 *       std::vector<double> old_lower_slack_multiplier_values(n_q_points);
 *       std::vector<double> old_upper_slack_multiplier_values(n_q_points);
 *       std::vector<double> old_lower_slack_values(n_q_points);
 *       std::vector<double> old_upper_slack_values(n_q_points);
 *       std::vector<double> old_unfiltered_density_values(n_q_points);
 *       std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
 *       std::vector<double> filtered_unfiltered_density_values(n_q_points);
 *       std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
 *         n_q_points);
 *   
 *       using namespace ValueExtractors;
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_rhs = 0;
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           fe_values.reinit(cell);
 *   
 *           lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
 *           mu.value_list(fe_values.get_quadrature_points(), mu_values);
 *   
 *           fe_values[densities<dim>].get_function_values(test_solution,
 *                                                         old_density_values);
 *           fe_values[displacements<dim>].get_function_values(
 *             test_solution, old_displacement_values);
 *           fe_values[displacements<dim>].get_function_divergences(
 *             test_solution, old_displacement_divs);
 *           fe_values[displacements<dim>].get_function_symmetric_gradients(
 *             test_solution, old_displacement_symmgrads);
 *           fe_values[displacement_multipliers<dim>].get_function_values(
 *             test_solution, old_displacement_multiplier_values);
 *           fe_values[displacement_multipliers<dim>].get_function_divergences(
 *             test_solution, old_displacement_multiplier_divs);
 *           fe_values[displacement_multipliers<dim>]
 *             .get_function_symmetric_gradients(
 *               test_solution, old_displacement_multiplier_symmgrads);
 *           fe_values[density_lower_slacks<dim>].get_function_values(
 *             test_solution, old_lower_slack_values);
 *           fe_values[density_lower_slack_multipliers<dim>].get_function_values(
 *             test_solution, old_lower_slack_multiplier_values);
 *           fe_values[density_upper_slacks<dim>].get_function_values(
 *             test_solution, old_upper_slack_values);
 *           fe_values[density_upper_slack_multipliers<dim>].get_function_values(
 *             test_solution, old_upper_slack_multiplier_values);
 *           fe_values[unfiltered_densities<dim>].get_function_values(
 *             test_solution, old_unfiltered_density_values);
 *           fe_values[unfiltered_density_multipliers<dim>].get_function_values(
 *             test_solution, old_unfiltered_density_multiplier_values);
 *           fe_values[unfiltered_densities<dim>].get_function_values(
 *             filtered_unfiltered_density_solution,
 *             filtered_unfiltered_density_values);
 *           fe_values[unfiltered_density_multipliers<dim>].get_function_values(
 *             filter_adjoint_unfiltered_density_multiplier_solution,
 *             filter_adjoint_unfiltered_density_multiplier_values);
 *   
 *           for (const auto q_point : fe_values.quadrature_point_indices())
 *             {
 *               for (const auto i : fe_values.dof_indices())
 *                 {
 *                   const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
 *                     fe_values[displacements<dim>].symmetric_gradient(i, q_point);
 *                   const double displacement_phi_i_div =
 *                     fe_values[displacements<dim>].divergence(i, q_point);
 *   
 *                   const SymmetricTensor<2, dim>
 *                     displacement_multiplier_phi_i_symmgrad =
 *                       fe_values[displacement_multipliers<dim>].symmetric_gradient(
 *                         i, q_point);
 *                   const double displacement_multiplier_phi_i_div =
 *                     fe_values[displacement_multipliers<dim>].divergence(i,
 *                                                                         q_point);
 *   
 *   
 *                   const double density_phi_i =
 *                     fe_values[densities<dim>].value(i, q_point);
 *                   const double unfiltered_density_phi_i =
 *                     fe_values[unfiltered_densities<dim>].value(i, q_point);
 *                   const double unfiltered_density_multiplier_phi_i =
 *                     fe_values[unfiltered_density_multipliers<dim>].value(i,
 *                                                                          q_point);
 *   
 *                   const double lower_slack_multiplier_phi_i =
 *                     fe_values[density_lower_slack_multipliers<dim>].value(
 *                       i, q_point);
 *   
 *                   const double lower_slack_phi_i =
 *                     fe_values[density_lower_slacks<dim>].value(i, q_point);
 *   
 *                   const double upper_slack_phi_i =
 *                     fe_values[density_upper_slacks<dim>].value(i, q_point);
 *   
 *                   const double upper_slack_multiplier_phi_i =
 *                     fe_values[density_upper_slack_multipliers<dim>].value(
 *                       i, q_point);
 *   
 *                   /* Equation 1: This equation, along with equations
 *                    * 2 and 3, are the variational derivatives of the
 *                    * Lagrangian with respect to the decision
 *                    * variables - the density, displacement, and
 *                    * unfiltered density. */
 *                   cell_rhs(i) +=
 *                     -1 * fe_values.JxW(q_point) *
 *                     (density_penalty_exponent *
 *                        std::pow(old_density_values[q_point],
 *                                 density_penalty_exponent - 1) *
 *                        density_phi_i *
 *                        (old_displacement_multiplier_divs[q_point] *
 *                           old_displacement_divs[q_point] *
 *                           lambda_values[q_point] +
 *                         2 * mu_values[q_point] *
 *                           (old_displacement_symmgrads[q_point] *
 *                            old_displacement_multiplier_symmgrads[q_point])) -
 *                      density_phi_i *
 *                        old_unfiltered_density_multiplier_values[q_point]);
 *   
 *                   /* Equation 2; the boundary terms will be added further down
 *                    * below. */
 *                   cell_rhs(i) +=
 *                     -1 * fe_values.JxW(q_point) *
 *                     (std::pow(old_density_values[q_point],
 *                               density_penalty_exponent) *
 *                      (old_displacement_multiplier_divs[q_point] *
 *                         displacement_phi_i_div * lambda_values[q_point] +
 *                       2 * mu_values[q_point] *
 *                         (old_displacement_multiplier_symmgrads[q_point] *
 *                          displacement_phi_i_symmgrad)));
 *   
 *                   /* Equation 3 */
 *                   cell_rhs(i) +=
 *                     -1 * fe_values.JxW(q_point) *
 *                     (unfiltered_density_phi_i *
 *                        filter_adjoint_unfiltered_density_multiplier_values
 *                          [q_point] +
 *                      unfiltered_density_phi_i *
 *                        old_upper_slack_multiplier_values[q_point] +
 *                      -1 * unfiltered_density_phi_i *
 *                        old_lower_slack_multiplier_values[q_point]);
 *   
 *   
 *   
 *                   /* Equation 4; boundary term will again be dealt
 *                    * with below. This equation being driven to 0
 *                    * ensures that the elasticity equation is met as
 *                    * a constraint. */
 *                   cell_rhs(i) += -1 * fe_values.JxW(q_point) *
 *                                  (std::pow(old_density_values[q_point],
 *                                            density_penalty_exponent) *
 *                                   (old_displacement_divs[q_point] *
 *                                      displacement_multiplier_phi_i_div *
 *                                      lambda_values[q_point] +
 *                                    2 * mu_values[q_point] *
 *                                      (displacement_multiplier_phi_i_symmgrad *
 *                                       old_displacement_symmgrads[q_point])));
 *   
 *                   /* Equation 5: This equation sets the lower slack
 *                    * variable equal to the unfiltered density,
 *                    * giving a minimum density of 0. */
 *                   cell_rhs(i) += fe_values.JxW(q_point) *
 *                                  (lower_slack_multiplier_phi_i *
 *                                   (old_unfiltered_density_values[q_point] -
 *                                    old_lower_slack_values[q_point]));
 *   
 *                   /* Equation 6: This equation sets the upper slack
 *                    * variable equal to one minus the unfiltered
 *                    * density. */
 *                   cell_rhs(i) += fe_values.JxW(q_point) *
 *                                  (upper_slack_multiplier_phi_i *
 *                                   (1 - old_unfiltered_density_values[q_point] -
 *                                    old_upper_slack_values[q_point]));
 *   
 *                   /* Equation 7: This is the difference between the
 *                    * density and the filter applied to the
 *                    * unfiltered density. This being driven to 0 by
 *                    * the Newton steps ensures that the filter is
 *                    * applied correctly. */
 *                   cell_rhs(i) += fe_values.JxW(q_point) *
 *                                  (unfiltered_density_multiplier_phi_i *
 *                                   (old_density_values[q_point] -
 *                                    filtered_unfiltered_density_values[q_point]));
 *   
 *                   /* Equation 8: This along with equation 9 give the
 *                    * requirement that s*z = \alpha for the barrier
 *                    * size alpha, and gives complementary slackness
 *                    * from KKT conditions when \alpha goes to 0. */
 *                   cell_rhs(i) +=
 *                     -1 * fe_values.JxW(q_point) *
 *                     (lower_slack_phi_i *
 *                      (old_lower_slack_multiplier_values[q_point] -
 *                       barrier_size / old_lower_slack_values[q_point]));
 *   
 *                   /* Equation 9 */
 *                   cell_rhs(i) +=
 *                     -1 * fe_values.JxW(q_point) *
 *                     (upper_slack_phi_i *
 *                      (old_upper_slack_multiplier_values[q_point] -
 *                       barrier_size / old_upper_slack_values[q_point]));
 *                 }
 *             }
 *   
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               if (face->at_boundary() &&
 *                   face->boundary_id() == BoundaryIds::down_force)
 *                 {
 *                   fe_face_values.reinit(cell, face);
 *   
 *                   for (const auto face_q_point :
 *                        fe_face_values.quadrature_point_indices())
 *                     {
 *                       for (const auto i : fe_face_values.dof_indices())
 *                         {
 *                           Tensor<1, dim> traction;
 *                           traction[1] = -1.;
 *   
 *                           cell_rhs(i) +=
 *                             -1 *
 *                             (traction * fe_face_values[displacements<dim>].value(
 *                                           i, face_q_point)) *
 *                             fe_face_values.JxW(face_q_point);
 *   
 *                           cell_rhs(i) +=
 *                             (traction *
 *                              fe_face_values[displacement_multipliers<dim>].value(
 *                                i, face_q_point)) *
 *                             fe_face_values.JxW(face_q_point);
 *                         }
 *                     }
 *                 }
 *             }
 *   
 *           MatrixTools::local_apply_boundary_values(boundary_values,
 *                                                    local_dof_indices,
 *                                                    dummy_cell_matrix,
 *                                                    cell_rhs,
 *                                                    true);
 *   
 *           constraints.distribute_local_to_global(cell_rhs,
 *                                                  local_dof_indices,
 *                                                  test_rhs);
 *         }
 *   
 *       return test_rhs;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Computingthemeritfunction"></a> 
 * <h4>Computing the merit function</h4>
 * 
 
 * 
 * The algorithm we use herein uses a "watchdog" strategy to
 * determine where and how far to go from the current iterate.  We
 * base the watchdog strategy on an exact @f$l_1@f$ merit function. This
 * function calculates the exact @f$l_1@f$ merit of a given, putative,
 * next iterate.
 *   
 
 * 
 * The merit function consists of the sum of the objective function
 * (which is simply an integral of external forces (on the boundary
 * of the domain) times the displacement values of a test solution
 * (typically, the current solution plus some multiple of the Newton
 * update), and the @f$l_1@f$ norms of the Lagrange multiplier
 * components of residual vectors. The following code computes these
 * parts in turn:
 * 
 * @code
 *     template <int dim>
 *     double SANDTopOpt<dim>::calculate_exact_merit(
 *       const BlockVector<double> &test_solution)
 *     {
 *       TimerOutput::Scope t(timer, "merit function");
 *   
 * @endcode
 * 
 * Start with computing the objective function:
 * 
 * @code
 *       double objective_function_merit = 0;
 *       {
 *         MappingQ<dim>         mapping(1);
 *         const QGauss<dim>     quadrature_formula(fe.degree + 1);
 *         const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 *         FEValues<dim>         fe_values(mapping,
 *                                 fe,
 *                                 quadrature_formula,
 *                                 update_values | update_gradients |
 *                                   update_quadrature_points | update_JxW_values);
 *         FEFaceValues<dim>     fe_face_values(mapping,
 *                                          fe,
 *                                          face_quadrature_formula,
 *                                          update_values |
 *                                            update_quadrature_points |
 *                                            update_normal_vectors |
 *                                            update_JxW_values);
 *   
 *         const unsigned int n_face_q_points = face_quadrature_formula.size();
 *   
 *         std::vector<Tensor<1, dim>> displacement_face_values(n_face_q_points);
 *   
 *         for (const auto &cell : dof_handler.active_cell_iterators())
 *           {
 *             for (const auto &face : cell->face_iterators())
 *               {
 *                 if (face->at_boundary() &&
 *                     face->boundary_id() == BoundaryIds::down_force)
 *                   {
 *                     fe_face_values.reinit(cell, face);
 *                     fe_face_values[ValueExtractors::displacements<dim>]
 *                       .get_function_values(test_solution,
 *                                            displacement_face_values);
 *                     for (unsigned int face_q_point = 0;
 *                          face_q_point < n_face_q_points;
 *                          ++face_q_point)
 *                       {
 *                         Tensor<1, dim> traction;
 *                         traction[1] = -1.;
 *   
 *                         objective_function_merit +=
 *                           (traction * displacement_face_values[face_q_point]) *
 *                           fe_face_values.JxW(face_q_point);
 *                       }
 *                   }
 *               }
 *           }
 *       }
 *   
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         {
 *           objective_function_merit =
 *             objective_function_merit -
 *             barrier_size * cell->measure() *
 *               std::log(test_solution.block(
 *                 SolutionBlocks::density_lower_slack)[cell->active_cell_index()]);
 *           objective_function_merit =
 *             objective_function_merit -
 *             barrier_size * cell->measure() *
 *               std::log(test_solution.block(
 *                 SolutionBlocks::density_upper_slack)[cell->active_cell_index()]);
 *         }
 *   
 * @endcode
 * 
 * Then compute the residual and take the @f$l_1@f$ norms of the
 * components that correspond to Lagrange mulipliers. We add
 * those to the objective function computed above, and return
 * the sum at the bottom:
 * 
 * @code
 *       const BlockVector<double> test_rhs = calculate_test_rhs(test_solution);
 *   
 *       const double elasticity_constraint_merit =
 *         penalty_multiplier *
 *         test_rhs.block(SolutionBlocks::displacement_multiplier).l1_norm();
 *       const double filter_constraint_merit =
 *         penalty_multiplier *
 *         test_rhs.block(SolutionBlocks::unfiltered_density_multiplier).l1_norm();
 *       const double lower_slack_merit =
 *         penalty_multiplier *
 *         test_rhs.block(SolutionBlocks::density_lower_slack_multiplier).l1_norm();
 *       const double upper_slack_merit =
 *         penalty_multiplier *
 *         test_rhs.block(SolutionBlocks::density_upper_slack_multiplier).l1_norm();
 *   
 *       const double total_merit =
 *         objective_function_merit + elasticity_constraint_merit +
 *         filter_constraint_merit + lower_slack_merit + upper_slack_merit;
 *       return total_merit;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Findingasearchdirection"></a> 
 * <h4>Finding a search direction</h4>
 * 
 
 * 
 * Next up is the function that actually computes a search direction
 * starting at the current state (passed as the first argument) and
 * returns the resulting vector. To this end, the function first
 * calls the functions that assemble the linear system that
 * corresponds to the Newton system, and that solve it.
 * 
 
 * 
 * This function also updates the penalty multiplier in the merit
 * function, and then returns the largest scaled feasible step.
 * It uses the `calculate_max_step_sizes()` function to find the
 * largest feasible step that satisfies @f$s>0@f$ and @f$z>0@f$.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     BlockVector<double> SANDTopOpt<dim>::find_max_step()
 *     {
 *       assemble_system();
 *       BlockVector<double> step = solve();
 *   
 * @endcode
 * 
 * Next we are going to update penalty_multiplier.  In
 * essence, a larger penalty multiplier makes us consider the
 * constraints more.  Looking at the Hessian and gradient with
 * respect to the step we want to take with our decision
 * variables, and comparing that to the norm of our constraint
 * error gives us a way to ensure that our merit function is
 * "exact" - that is, it has a minimum in the same location
 * that the objective function does.  As our merit function is
 * exact for any penalty multiplier over some minimum value,
 * we only keep the computed value if it increases the penalty
 * multiplier.
 * 
 
 * 
 * 
 * @code
 *       const std::vector<unsigned int> decision_variables = {
 *         SolutionBlocks::density,
 *         SolutionBlocks::displacement,
 *         SolutionBlocks::unfiltered_density,
 *         SolutionBlocks::density_upper_slack,
 *         SolutionBlocks::density_lower_slack};
 *       double hess_part = 0;
 *       double grad_part = 0;
 *       for (const unsigned int decision_variable_i : decision_variables)
 *         {
 *           for (const unsigned int decision_variable_j : decision_variables)
 *             {
 *               Vector<double> temp_vector(step.block(decision_variable_i).size());
 *               system_matrix.block(decision_variable_i, decision_variable_j)
 *                 .vmult(temp_vector, step.block(decision_variable_j));
 *               hess_part += step.block(decision_variable_i) * temp_vector;
 *             }
 *           grad_part -= system_rhs.block(decision_variable_i) *
 *                        step.block(decision_variable_i);
 *         }
 *   
 *       const std::vector<unsigned int> equality_constraint_multipliers = {
 *         SolutionBlocks::displacement_multiplier,
 *         SolutionBlocks::unfiltered_density_multiplier,
 *         SolutionBlocks::density_lower_slack_multiplier,
 *         SolutionBlocks::density_upper_slack_multiplier};
 *       double constraint_norm = 0;
 *       for (const unsigned int multiplier_i : equality_constraint_multipliers)
 *         constraint_norm += system_rhs.block(multiplier_i).linfty_norm();
 *   
 *   
 *       double test_penalty_multiplier;
 *       if (hess_part > 0)
 *         test_penalty_multiplier =
 *           (grad_part + .5 * hess_part) / (.05 * constraint_norm);
 *       else
 *         test_penalty_multiplier = (grad_part) / (.05 * constraint_norm);
 *   
 *       penalty_multiplier = std::max(penalty_multiplier, test_penalty_multiplier);
 *   
 * @endcode
 * 
 * Based on all of this, we can now compute step sizes for the
 * primal and dual (Lagrange multiplier) variables. Once we
 * have these, we scale the components of the solution vector,
 * and that is what this function returns.
 * 
 * @code
 *       const std::pair<double, double> max_step_sizes =
 *         calculate_max_step_size(nonlinear_solution, step);
 *       const double step_size_s = max_step_sizes.first;
 *       const double step_size_z = max_step_sizes.second;
 *   
 *       step.block(SolutionBlocks::density) *= step_size_s;
 *       step.block(SolutionBlocks::displacement) *= step_size_s;
 *       step.block(SolutionBlocks::unfiltered_density) *= step_size_s;
 *       step.block(SolutionBlocks::displacement_multiplier) *= step_size_z;
 *       step.block(SolutionBlocks::unfiltered_density_multiplier) *= step_size_z;
 *       step.block(SolutionBlocks::density_lower_slack) *= step_size_s;
 *       step.block(SolutionBlocks::density_lower_slack_multiplier) *= step_size_z;
 *       step.block(SolutionBlocks::density_upper_slack) *= step_size_s;
 *       step.block(SolutionBlocks::density_upper_slack_multiplier) *= step_size_z;
 *   
 *       return step;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Computingascaledstep"></a> 
 * <h4>Computing a scaled step</h4>
 * 
 
 * 
 * The next function then implements a back-tracking algorithm for a
 * line search. It keeps shrinking step size until it finds a step
 * where the merit is decreased, and then returns the new location
 * based on the current state vector, and the direction to go into,
 * times the step length.
 * 
 * @code
 *     template <int dim>
 *     BlockVector<double>
 *     SANDTopOpt<dim>::compute_scaled_step(const BlockVector<double> &state,
 *                                          const BlockVector<double> &max_step,
 *                                          const double descent_requirement)
 *     {
 *       const double merit_derivative =
 *         (calculate_exact_merit(state + 1e-4 * max_step) -
 *          calculate_exact_merit(state)) /
 *         1e-4;
 *       double       step_size                 = 1;
 *       unsigned int max_linesearch_iterations = 10;
 *       for (unsigned int k = 0; k < max_linesearch_iterations; ++k)
 *         {
 *           if (calculate_exact_merit(state + step_size * max_step) <
 *               calculate_exact_merit(state) +
 *                 step_size * descent_requirement * merit_derivative)
 *             break;
 *           else
 *             step_size = step_size / 2;
 *         }
 *       return state + (step_size * max_step);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Checkingforconvergence"></a> 
 * <h4>Checking for convergence</h4>
 * 
 
 * 
 * The final auxiliary function in this block is the one that checks
 * to see if the KKT conditions are sufficiently met so that the
 * overall algorithm can lower the barrier size. It does so by
 * computing the @f$l_1@f$ norm of the residual, which is what
 * `calculate_test_rhs()` computes.
 * 
 * @code
 *     template <int dim>
 *     bool SANDTopOpt<dim>::check_convergence(const BlockVector<double> &state)
 *     {
 *       const BlockVector<double> test_rhs      = calculate_test_rhs(state);
 *       const double              test_rhs_norm = test_rhs.l1_norm();
 *   
 *       const double convergence_condition = 1e-2;
 *       const double target_norm           = convergence_condition * barrier_size;
 *   
 *       std::cout << "    Checking convergence. Current rhs norm is "
 *                 << test_rhs_norm << ", target is " << target_norm << std::endl;
 *   
 *       return (test_rhs_norm < target_norm);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Postprocessingthesolution"></a> 
 * <h3>Postprocessing the solution</h3>
 * 
 
 * 
 * The first of the postprocessing functions outputs information
 * in a VTU file for visualization. It looks long, but it's really
 * just the same as what was done in @ref step_22 "step-22", for example, just
 * with (a lot) more solution variables:
 * 
 * @code
 *     template <int dim>
 *     void SANDTopOpt<dim>::output_results(const unsigned int iteration) const
 *     {
 *       std::vector<std::string> solution_names(1, "density");
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         data_component_interpretation(
 *           1, DataComponentInterpretation::component_is_scalar);
 *       for (unsigned int i = 0; i < dim; ++i)
 *         {
 *           solution_names.emplace_back("displacement");
 *           data_component_interpretation.push_back(
 *             DataComponentInterpretation::component_is_part_of_vector);
 *         }
 *       solution_names.emplace_back("unfiltered_density");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *       for (unsigned int i = 0; i < dim; ++i)
 *         {
 *           solution_names.emplace_back("displacement_multiplier");
 *           data_component_interpretation.push_back(
 *             DataComponentInterpretation::component_is_part_of_vector);
 *         }
 *       solution_names.emplace_back("unfiltered_density_multiplier");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *       solution_names.emplace_back("low_slack");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *       solution_names.emplace_back("low_slack_multiplier");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *       solution_names.emplace_back("high_slack");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *       solution_names.emplace_back("high_slack_multiplier");
 *       data_component_interpretation.push_back(
 *         DataComponentInterpretation::component_is_scalar);
 *   
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(nonlinear_solution,
 *                                solution_names,
 *                                DataOut<dim>::type_dof_data,
 *                                data_component_interpretation);
 *       data_out.build_patches();
 *   
 *       std::ofstream output("solution" + std::to_string(iteration) + ".vtu");
 *       data_out.write_vtu(output);
 *     }
 *   
 *   
 * @endcode
 * 
 * The second of these functions outputs the solution as an `.stl`
 * file for 3d
 * printing. [STL](https://en.wikipedia.org/wiki/STL_(file_format))
 * files are made up of triangles and normal vectors, and we will
 * use it to show all of those cells with a density value larger
 * than zero by first extruding the mesh from a @f$z@f$ value of zero
 * to @f$z=0.25@f$, and then generating two triangles for each face of
 * the cells with a sufficiently large density value. The triangle
 * nodes must go counter-clockwise when looking from the outside,
 * and the normal vectors must be unit vectors pointing outwards,
 * which requires a few checks.
 * 
 * @code
 *     template <int dim>
 *     void SANDTopOpt<dim>::write_as_stl()
 *     {
 *       static_assert(dim == 2,
 *                     "This function is not implemented for anything "
 *                     "other than the 2d case.");
 *   
 *       std::ofstream stlfile;
 *       stlfile.open("bridge.stl");
 *   
 *       stlfile << "solid bridge\n" << std::scientific;
 *       double height = .25;
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           if (nonlinear_solution.block(
 *                 SolutionBlocks::density)[cell->active_cell_index()] > 0.5)
 *             {
 * @endcode
 * 
 * We have now found a cell with a density value larger
 * than zero. Let us start by writing out the bottom
 * and top faces. Owing to the ordering issue mentioned
 * above, we have to make sure that we understand
 * whether a cell has a right- or left-handed
 * coordinate system. We do this by interrogating the
 * directions of the two edges starting at vertex 0 and
 * whether they form a right-handed coordinate system.
 * 
 * @code
 *               const Tensor<1, dim> edge_directions[2] = {cell->vertex(1) -
 *                                                            cell->vertex(0),
 *                                                          cell->vertex(2) -
 *                                                            cell->vertex(0)};
 *               const Tensor<2, dim> edge_tensor(
 *                 {{edge_directions[0][0], edge_directions[0][1]},
 *                  {edge_directions[1][0], edge_directions[1][1]}});
 *               const bool is_right_handed_cell = (determinant(edge_tensor) > 0);
 *   
 *               if (is_right_handed_cell)
 *                 {
 *                   /* Write one side at z = 0. */
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << -1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(0)[0] << ' '
 *                           << cell->vertex(0)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << -1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(3)[0] << ' '
 *                           << cell->vertex(3)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *   
 *                   /* Write one side at z = height. */
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << 1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(0)[0] << ' '
 *                           << cell->vertex(0)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << height << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << 1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(3)[0] << ' '
 *                           << cell->vertex(3)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << height << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                 }
 *               else /* The cell has a left-handed set up */
 *                 {
 *                   /* Write one side at z = 0. */
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << -1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(0)[0] << ' '
 *                           << cell->vertex(0)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << -1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(3)[0] << ' '
 *                           << cell->vertex(3)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << 0.000000e+00 << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *   
 *                   /* Write one side at z = height. */
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << 1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(0)[0] << ' '
 *                           << cell->vertex(0)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << height << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                   stlfile << "   facet normal " << 0.000000e+00 << ' '
 *                           << 0.000000e+00 << ' ' << 1.000000e+00 << '\n';
 *                   stlfile << "      outer loop\n";
 *                   stlfile << "         vertex " << cell->vertex(1)[0] << ' '
 *                           << cell->vertex(1)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(2)[0] << ' '
 *                           << cell->vertex(2)[1] << ' ' << height << '\n';
 *                   stlfile << "         vertex " << cell->vertex(3)[0] << ' '
 *                           << cell->vertex(3)[1] << ' ' << height << '\n';
 *                   stlfile << "      endloop\n";
 *                   stlfile << "   endfacet\n";
 *                 }
 *   
 * @endcode
 * 
 * Next we need to deal with the four faces of the
 * cell, extended into the @f$z@f$ direction. However, we
 * only need to write these faces if either the face
 * is on the domain boundary, or if it is the
 * interface between a cell with density greater than
 * 0.5, and a cell with a density less than 0.5.
 * 
 * @code
 *               for (unsigned int face_number = 0;
 *                    face_number < GeometryInfo<dim>::faces_per_cell;
 *                    ++face_number)
 *                 {
 *                   const typename DoFHandler<dim>::face_iterator face =
 *                     cell->face(face_number);
 *   
 *                   if ((face->at_boundary()) ||
 *                       (!face->at_boundary() &&
 *                        (nonlinear_solution.block(
 *                           0)[cell->neighbor(face_number)->active_cell_index()] <
 *                         0.5)))
 *                     {
 *                       const Tensor<1, dim> normal_vector =
 *                         (face->center() - cell->center());
 *                       const double normal_norm = normal_vector.norm();
 *                       if ((face->vertex(0)[0] - face->vertex(0)[0]) *
 *                               (face->vertex(1)[1] - face->vertex(0)[1]) *
 *                               0.000000e+00 +
 *                             (face->vertex(0)[1] - face->vertex(0)[1]) * (0 - 0) *
 *                               normal_vector[0] +
 *                             (height - 0) *
 *                               (face->vertex(1)[0] - face->vertex(0)[0]) *
 *                               normal_vector[1] -
 *                             (face->vertex(0)[0] - face->vertex(0)[0]) * (0 - 0) *
 *                               normal_vector[1] -
 *                             (face->vertex(0)[1] - face->vertex(0)[1]) *
 *                               (face->vertex(1)[0] - face->vertex(0)[0]) *
 *                               normal_vector[0] -
 *                             (height - 0) *
 *                               (face->vertex(1)[1] - face->vertex(0)[1]) * 0 >
 *                           0)
 *                         {
 *                           stlfile << "   facet normal "
 *                                   << normal_vector[0] / normal_norm << ' '
 *                                   << normal_vector[1] / normal_norm << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      outer loop\n";
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      endloop\n";
 *                           stlfile << "   endfacet\n";
 *                           stlfile << "   facet normal "
 *                                   << normal_vector[0] / normal_norm << ' '
 *                                   << normal_vector[1] / normal_norm << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      outer loop\n";
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      endloop\n";
 *                           stlfile << "   endfacet\n";
 *                         }
 *                       else
 *                         {
 *                           stlfile << "   facet normal "
 *                                   << normal_vector[0] / normal_norm << ' '
 *                                   << normal_vector[1] / normal_norm << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      outer loop\n";
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "      endloop\n";
 *                           stlfile << "   endfacet\n";
 *                           stlfile << "   facet normal "
 *                                   << normal_vector[0] / normal_norm << ' '
 *                                   << normal_vector[1] / normal_norm << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "      outer loop\n";
 *                           stlfile << "         vertex " << face->vertex(0)[0]
 *                                   << ' ' << face->vertex(0)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' '
 *                                   << 0.000000e+00 << '\n';
 *                           stlfile << "         vertex " << face->vertex(1)[0]
 *                                   << ' ' << face->vertex(1)[1] << ' ' << height
 *                                   << '\n';
 *                           stlfile << "      endloop\n";
 *                           stlfile << "   endfacet\n";
 *                         }
 *                     }
 *                 }
 *             }
 *         }
 *       stlfile << "endsolid bridge";
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Therunfunctiondrivingtheoverallalgorithm"></a> 
 * <h3>The run() function driving the overall algorithm</h3>
 * 
 
 * 
 * This function finally provides the overall driver logic. It is,
 * in the grand scheme of things, a rather complicated function
 * primarily because the optimization algorithm is difficult: It
 * isn't just about finding a Newton direction like in @ref step_15 "step-15" and
 * then going a fixed distance in this direction any more, but
 * instead about (i) determining what the optimal log-barrier
 * penalty parameter should be in the current step, (ii) a
 * complicated algorithm to determine how far we want to go, and
 * other ingredients. Let us see how we can break this down into
 * smaller chunks in the following documentation.
 *   
 
 * 
 * The function starts out simple enough with first setting up the
 * mesh, the DoFHandler, and then the various linear algebra objects
 * necessary for the following:
 * 
 * @code
 *     template <int dim>
 *     void SANDTopOpt<dim>::run()
 *     {
 *       std::cout << "filter r is: " << filter_r << std::endl;
 *   
 *       {
 *         TimerOutput::Scope t(timer, "setup");
 *   
 *         create_triangulation();
 *   
 *         dof_handler.distribute_dofs(fe);
 *         DoFRenumbering::component_wise(dof_handler);
 *   
 *         setup_boundary_values();
 *         setup_block_system();
 *         setup_filter_matrix();
 *       }
 *   
 * @endcode
 * 
 * We then set a number of parameters that affect the
 * log-barrier and line search components of the optimization
 * algorithm:
 * 
 * @code
 *       barrier_size                  = 25;
 *       const double min_barrier_size = .0005;
 *   
 *       const unsigned int max_uphill_steps    = 8;
 *       const double       descent_requirement = .0001;
 *   
 *   
 * @endcode
 * 
 * Now start the principal iteration. The overall algorithm
 * works by using an outer loop in which we loop until either
 * (i) the log-barrier parameter has become small enough, or (ii)
 * we have reached convergence. In any case, we terminate if
 * end up with too large a number of iterations. This overall
 * structure is encoded as a `do { ... } while (...)` loop
 * where the convergence condition is at the bottom.
 * 
 * @code
 *       unsigned int       iteration_number = 0;
 *       const unsigned int max_iterations   = 10000;
 *   
 *       do
 *         {
 *           std::cout << "Starting outer step in iteration " << iteration_number
 *                     << " with barrier parameter " << barrier_size << std::endl;
 *   
 * @endcode
 * 
 * Within this outer loop, we have an inner loop in which we
 * try to find an update direction using the watchdog
 * algorithm described in the introduction.
 *         
 
 * 
 * The general idea of the watchdog algorithm itself is
 * this: For a maximum of `max_uphill_steps` (i.e., a loop
 * within the "inner loop" mentioned above) attempts, we use
 * `find_max_step()` to compute a Newton update step, and
 * add these up in the `nonlinear_solution` vector.  In each of
 * these attempts (starting from the place reached at the
 * end of the previous attempt), we check whether we have
 * reached a target value of the merit function described
 * above. The target value is computed based on where this
 * algorithm starts (the `nonlinear_solution` at the beginning of
 * the watchdog loop, saves as `watchdog_state`) and the
 * first proposed direction provided by `find_max_step()` in
 * the first go-around of this loop (the `k==0` case).
 * 
 * @code
 *           do
 *             {
 *               std::cout << "  Starting inner step in iteration "
 *                         << iteration_number
 *                         << " with merit function penalty multiplier "
 *                         << penalty_multiplier << std::endl;
 *   
 *               bool watchdog_step_found = false;
 *   
 *               const BlockVector<double> watchdog_state = nonlinear_solution;
 *               BlockVector<double>       first_step;
 *               double target_merit     = numbers::signaling_nan<double>();
 *               double merit_derivative = numbers::signaling_nan<double>();
 *   
 *               for (unsigned int k = 0; k < max_uphill_steps; ++k)
 *                 {
 *                   ++iteration_number;
 *                   const BlockVector<double> update_step = find_max_step();
 *   
 *                   if (k == 0)
 *                     {
 *                       first_step = update_step;
 *                       merit_derivative =
 *                         ((calculate_exact_merit(watchdog_state +
 *                                                 .0001 * first_step) -
 *                           calculate_exact_merit(watchdog_state)) /
 *                          .0001);
 *                       target_merit = calculate_exact_merit(watchdog_state) +
 *                                      descent_requirement * merit_derivative;
 *                     }
 *   
 *                   nonlinear_solution += update_step;
 *                   const double current_merit =
 *                     calculate_exact_merit(nonlinear_solution);
 *   
 *                   std::cout << "    current watchdog state merit is: "
 *                             << current_merit << "; target merit is "
 *                             << target_merit << std::endl;
 *   
 *                   if (current_merit < target_merit)
 *                     {
 *                       watchdog_step_found = true;
 *                       std::cout << "    found workable step after " << k + 1
 *                                 << " iterations" << std::endl;
 *                       break;
 *                     }
 *                 }
 *   
 *   
 * @endcode
 * 
 * The next part of the algorithm then depends on
 * whether the watchdog loop above succeeded. If it
 * did, then we are satisfied and no further action is
 * necessary: We just stay where we are. If, however,
 * we took the maximal number of unsuccessful steps in
 * the loop above, then we need to do something else,
 * and this is what the following code block does.
 *             
 
 * 
 * Specifically, from the final (unsuccessful) state
 * of the loop above, we seek one more update
 * direction and take what is called a "stretch
 * step". If that stretch state satisfies a condition
 * involving the merit function, then we go there. On
 * the other hand, if the stretch state is also
 * unacceptable (as all of the watchdog steps above
 * were), then we discard all of the watchdog steps
 * taken above and start over again where we had
 * started the watchdog iterations -- that place was
 * stored in the `watchdog_state` variable above. More
 * specifically, the conditions below first test
 * whether we take a step from `watchdog_state` in
 * direction `first_step`, or whether we can do one
 * more update from the stretch state to find a new
 * place. It is possible that neither of these is
 * actually better than the state we started from at
 * the beginning of the watchdog algorithm, but even
 * if that is so, that place clearly was a difficult
 * place to be in, and getting away to start the next
 * iteration from another place might be a useful
 * strategy to eventually converge.
 *             
 
 * 
 * We keep repeating the watchdog steps above along
 * with the logic below until this inner iteration is
 * finally converged (or if we run up against the
 * maximal number of iterations -- where we count the
 * number of linear solves as iterations and increment
 * the counter every time we call `find_max_step()`
 * since that is where the linear solve actually
 * happens). In any case, at the end of each of these
 * inner iterations we also output the solution in a
 * form suitable for visualization.
 * 
 
 * 
 * 
 * @code
 *               if (watchdog_step_found == false)
 *                 {
 *                   ++iteration_number;
 *                   const BlockVector<double> update_step = find_max_step();
 *                   const BlockVector<double> stretch_state =
 *                     compute_scaled_step(nonlinear_solution,
 *                                         update_step,
 *                                         descent_requirement);
 *   
 * @endcode
 * 
 * If we did not get a successful watchdog step,
 * we now need to decide between going back to
 * where we started, or using the final state.  We
 * compare the merits of both of these locations,
 * and then take a scaled step from whichever
 * location is better.  As the scaled step is
 * guaranteed to lower the merit, we will end up
 * keeping one of the two.
 * 
 * @code
 *                   if ((calculate_exact_merit(nonlinear_solution) <
 *                        calculate_exact_merit(watchdog_state)) ||
 *                       (calculate_exact_merit(stretch_state) < target_merit))
 *                     {
 *                       std::cout << "    Taking scaled step from end of watchdog"
 *                                 << std::endl;
 *                       nonlinear_solution = stretch_state;
 *                     }
 *                   else
 *                     {
 *                       std::cout
 *                         << "    Taking scaled step from beginning of watchdog"
 *                         << std::endl;
 *                       if (calculate_exact_merit(stretch_state) >
 *                           calculate_exact_merit(watchdog_state))
 *                         {
 *                           nonlinear_solution =
 *                             compute_scaled_step(watchdog_state,
 *                                                 first_step,
 *                                                 descent_requirement);
 *                         }
 *                       else
 *                         {
 *                           ++iteration_number;
 *                           nonlinear_solution = stretch_state;
 *                           const BlockVector<double> stretch_step =
 *                             find_max_step();
 *                           nonlinear_solution =
 *                             compute_scaled_step(nonlinear_solution,
 *                                                 stretch_step,
 *                                                 descent_requirement);
 *                         }
 *                     }
 *                 }
 *   
 *               output_results(iteration_number);
 *             }
 *           while ((iteration_number < max_iterations) &&
 *                  (check_convergence(nonlinear_solution) == false));
 *   
 *   
 * @endcode
 * 
 * At the end of the outer loop, we have to update the
 * barrier parameter, for which we use the following
 * formula. The rest of the function is then simply about
 * checking the outer loop convergence condition, and if
 * we decide to terminate computations, about writing the
 * final "design" as an STL file for use in 3d printing,
 * and to output some timing information.
 * 
 * @code
 *           const double barrier_size_multiplier = .8;
 *           const double barrier_size_exponent   = 1.2;
 *   
 *           barrier_size =
 *             std::max(std::min(barrier_size * barrier_size_multiplier,
 *                               std::pow(barrier_size, barrier_size_exponent)),
 *                      min_barrier_size);
 *   
 *           std::cout << std::endl;
 *         }
 *       while (((barrier_size > min_barrier_size) ||
 *               (check_convergence(nonlinear_solution) == false)) &&
 *              (iteration_number < max_iterations));
 *   
 *       write_as_stl();
 *     }
 *   } // namespace SAND
 *   
 * @endcode
 * 
 * 
 * <a name="Themainfunction"></a> 
 * <h3>The main function</h3>
 * 
 
 * 
 * The remainder of the code, the `main()` function, is as usual:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         SAND::SANDTopOpt<2> elastic_problem_2d;
 *         elastic_problem_2d.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;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
<a name="TestProblem"></a><h3>Test Problem</h3>
 
The algorithms used above are tested against a traditional topology optimization
 problem called the Messerschmitt-Bolkow-Blohm Beam (MBB Beam).
 
This problem considers the optimal 2-d structure that can be built on a
rectangle 6 units wide, and 1 unit tall. The bottom corners are fixed in place
in the @f$y@f$ direction using a zero Dirichlet boundary condition, and a downward
force is applied in the center of the top of the beam by enforcing a Neumann
boundary condition. The rest of the boundary is allowed to move, and has no
external force applied, which takes the form of a zero Neumann boundary
condition. In essence, we are asking the following question: How should we
design a bridge in a way so that if the bottom left and bottom right point of
the bridge are on rollers that allow these points to move horizontally but not
vertically, and so that the displacement in response to the vertical force in
the center is minimal.
 
While the total volume of the domain is 6 units, 3 units of material are allowed for
the structure. Because of the symmetry of the problem, it could be posed on a
rectangle of width 3 and height 1 by cutting the original domain in half, and
using zero Dirichlet boundary conditions in the @f$x@f$ direction along the cut
edge. That said, symmetry of the solution is a good indicator that the program
is working as expected, so we solved the problem on the whole domain,
as shown below. @cite Bendse2004
 
<div style="text-align:center;">
  <img src="https://www.dealii.org/images/steps/developer/step-79.mbbgeometry.png"
       alt="The MBB problem domain and boundary conditions">
</div>
 
 
Using the program discussed above, we find the minimum volume of the MBB Beam and the
individual components of the solution look as follows:
 
<div class="onecolumn" style="width: 80%; text-align: center;">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-79.filtereddensity.png"
         alt="Filtered Density Solution">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-79.unfiltereddensity.png"
         alt="Unfiltered Density Solution">
  </div>
</div>
 
 
These pictures show that what we find here is in accordance with what one
typically sees in other publications on the topic @cite Bendse2004.
Maybe more interestingly, the
result looks like a truss bridge (except that we apply the load at the top of
the trusses, rather than the bottom as in real truss bridges, akin to a "deck
truss" bridge), suggesting that the designs that have been used in bridge-building
for centuries are indeed based on ideas we can now show to be optimal
in some sense.
<a name="Possibilitiesforextensions"></a><h4>Possibilities for extensions</h4>
 
 
The results shown above took around 75 iterations to find, which is quite
concerning given the expense in solving the large linear systems in each
iteration. Looking at the evolution, it does look as though the convergence has
moments of happening quickly and moments of happening slowly. We believe this to
be due to both a lack of precision on when and how to decrease the boundary
values, as well as our choice of merit function being sub-optimal. In the future,
a LOQO barrier update replacing the monotone reduction, as well as a Markov
Filter in place of a merit function will decrease the number of necessary
iterations significantly.
 
The barrier decrease is most sensitive in the middle of the convergence, which
is problematic, as it seems like we need it to decrease quickly, then slowly,
then quickly again.
 
Secondly, the linear solver used here is just the sparse direct solver based on
the SparseDirectUMFPACK class. This works reasonably well on small problems,
but the formulation of the optimization problem detailed above has quite a large
number of variables and so the linear problem is not only large but also has a
lot of nonzero entries in many rows, even on meshes that overall are still
relatively coarse. As a consequence, the solver time dominates the
computations, and more sophisticated approaches at solving the linear system
are necessary.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-79.cc"
*/