<commit_before>/* ---------------------------------------------------------------------
 *
 * Copyright (C) 2000 - 2016 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * ---------------------------------------------------------------------

 *
 * Author: <NAME>, University of Heidelberg, 2000
 */


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>

#include <fstream>
#include <iostream>

// From the following include file we will import the declaration of
// H1-conforming finite element shape functions. This family of finite
// elements is called <code>FE_Q</code>, and was used in all examples before
// already to define the usual bi- or tri-linear elements, but we will now use
// it for bi-quadratic elements:
#include <deal.II/fe/fe_q.h>
// We will not read the grid from a file as in the previous example, but
// generate it using a function of the library. However, we will want to write
// out the locally refined grids (just the grid, not the solution) in each
// step, so we need the following include file instead of
// <code>grid_in.h</code>:
#include <deal.II/grid/grid_out.h>


// When using locally refined grids, we will get so-called <code>hanging
// nodes</code>. However, the standard finite element methods assumes that the
// discrete solution spaces be continuous, so we need to make sure that the
// degrees of freedom on hanging nodes conform to some constraints such that
// the global solution is continuous. We are also going to store the boundary
// conditions in this object. The following file contains a class which is
// used to handle these constraints:
#include <deal.II/lac/constraint_matrix.h>

// In order to refine our grids locally, we need a function from the library
// that decides which cells to flag for refinement or coarsening based on the
// error indicators we have computed. This function is defined here:
#include <deal.II/grid/grid_refinement.h>

// Finally, we need a simple way to actually compute the refinement indicators
// based on some error estimate. While in general, adaptivity is very
// problem-specific, the error indicator in the following file often yields
// quite nicely adapted grids for a wide class of problems.
#include <deal.II/numerics/error_estimator.h>

// Finally, this is as in previous programs:
using namespace dealii;


// @sect3{The <code>Step6</code> class template}

// The main class is again almost unchanged. Two additions, however, are made:
// we have added the <code>refine_grid</code> function, which is used to
// adaptively refine the grid (instead of the global refinement in the
// previous examples), and a variable which will hold the constraints. In
// addition, we have added a destructor to the class for reasons that will
// become clear when we discuss its implementation.
template <int dim>
class Step6
{
public:
  Step6 ();
  ~Step6 ();

  void run ();

private:
  void setup_system ();
  void assemble_system ();
  void solve ();
  void refine_grid ();
  void output_results (const unsigned int cycle) const;

  Triangulation<dim>           triangulation;
  const SphericalManifold<dim> boundary;

  FE_Q<dim>            fe;
  DoFHandler<dim>      dof_handler;

  // This is the new variable in the main class. We need an object which holds
  // a list of constraints to hold the hanging nodes and the boundary
  // conditions.
  ConstraintMatrix     constraints;

  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;

  Vector<double>       solution;
  Vector<double>       system_rhs;
};


// @sect3{Nonconstant coefficients}

// The implementation of nonconstant coefficients is copied verbatim from
// step-5:
template <int dim>
double coefficient (const Point<dim> &p)
{
  if (p.square() < 0.5*0.5)
    return 20;
  else
    return 1;
}



// @sect3{The <code>Step6</code> class implementation}

// @sect4{Step6::Step6}

// The constructor of this class is mostly the same as before, but this time
// we want to use the quadratic element. To do so, we only have to replace the
// constructor argument (which was <code>1</code> in all previous examples) by
// the desired polynomial degree (here <code>2</code>):
template <int dim>
Step6<dim>::Step6 ()
  :
  fe (2),
  dof_handler (triangulation)
{}


// @sect4{Step6::~Step6}

// Here comes the added destructor of the class. The reason why we want to add
// it is to pick up an issue we already started to discuss in step-1:
// the <code>boundary</code> object was defined before and not after the
// <code>triangulation</code> object. Of course we could have left this order
// unchanged, but we would like to show what happens if the order is reversed
// since this produces a rather nasty side-effect and results in an error which
// is difficult to track down if one does not know what happens.
//
// Basically what happens is the following: when we set a manifold description
// to the triangulation using the function call
// <code>triangulation.set_manifold (0, boundary)</code>, the
// <code>Triangulation</code> object also stores a pointer to the
// <code>Manifold</code> object in use. Since this pointer is used until either
// another <code>Manifold</code> object is set as boundary description or until
// the <code>Triangulation</code> object is destroyed, it would be unwise if we
// would allow the <code>boundary</code> to be deleted before the
// <code>triangulation</code>. To disallow this, the <code>triangulation</code>
// increases a counter inside the <code>boundary</code> which counts how many
// objects use it (this is what the
// <code>Subscriptor</code>/<code>SmartPointer</code> class pair is used for,
// in case you want something like this for your own programs; see step-7 for
// a more complete discussion of this topic). The <code>boundary</code> will
// refuse its destruction if that counter is larger than zero, since then some
// other objects might rely on its persistence. An exception will then be
// thrown and the program will usually abort upon the attempt to destroy
// <code>boundary</code>.
//
// To be fair, such exceptions about still used objects are not particularly
// popular among programmers using deal.II, since they only tell us that
// something is wrong, namely that some other object is still using the object
// that is presently being destroyed, but most of the time not who this using
// object is. It is therefore often rather time-consuming to find out where the
// problem exactly is, although it is then usually straightforward to remedy
// the situation. However, we believe that the effort to find invalid
// references to objects that do no longer exist is less if the problem is
// detected once the reference becomes invalid, rather than when non-existent
// objects are actually accessed again, since then usually only invalid data
// is accessed, but no error is immediately raised.
//
// Coming back to the present situation, if we did not write this destructor,
// the compiler will generate code that triggers exactly the behavior sketched
// above. The reason is that member variables of the <code>Step6</code> class
// are destroyed bottom-up (i.e., in reverse order of their declaration in the
// class), as always in C++. Thus, the boundary object will be
// destroyed before the triangulation object, since its declaration is below
// the one of the triangulation. This triggers the situation above, and an<commit_msg>Improve some wording.<commit_after>/* ---------------------------------------------------------------------
 *
 * Copyright (C) 2000 - 2016 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * ---------------------------------------------------------------------

 *
 * Author: <NAME>, University of Heidelberg, 2000
 */


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>

#include <fstream>
#include <iostream>

// From the following include file we will import the declaration of
// H1-conforming finite element shape functions. This family of finite
// elements is called <code>FE_Q</code>, and was used in all examples before
// already to define the usual bi- or tri-linear elements, but we will now use
// it for bi-quadratic elements:
#include <deal.II/fe/fe_q.h>
// We will not read the grid from a file as in the previous example, but
// generate it using a function of the library. However, we will want to write
// out the locally refined grids (just the grid, not the solution) in each
// step, so we need the following include file instead of
// <code>grid_in.h</code>:
#include <deal.II/grid/grid_out.h>


// When using locally refined grids, we will get so-called <code>hanging
// nodes</code>. However, the standard finite element methods assumes that the
// discrete solution spaces be continuous, so we need to make sure that the
// degrees of freedom on hanging nodes conform to some constraints such that
// the global solution is continuous. We are also going to store the boundary
// conditions in this object. The following file contains a class which is
// used to handle these constraints:
#include <deal.II/lac/constraint_matrix.h>

// In order to refine our grids locally, we need a function from the library
// that decides which cells to flag for refinement or coarsening based on the
// error indicators we have computed. This function is defined here:
#include <deal.II/grid/grid_refinement.h>

// Finally, we need a simple way to actually compute the refinement indicators
// based on some error estimate. While in general, adaptivity is very
// problem-specific, the error indicator in the following file often yields
// quite nicely adapted grids for a wide class of problems.
#include <deal.II/numerics/error_estimator.h>

// Finally, this is as in previous programs:
using namespace dealii;


// @sect3{The <code>Step6</code> class template}

// The main class is again almost unchanged. Two additions, however, are made:
// we have added the <code>refine_grid</code> function, which is used to
// adaptively refine the grid (instead of the global refinement in the
// previous examples), and a variable which will hold the constraints. In
// addition, we have added a destructor to the class for reasons that will
// become clear when we discuss its implementation.
template <int dim>
class Step6
{
public:
  Step6 ();
  ~Step6 ();

  void run ();

private:
  void setup_system ();
  void assemble_system ();
  void solve ();
  void refine_grid ();
  void output_results (const unsigned int cycle) const;

  Triangulation<dim>           triangulation;
  const SphericalManifold<dim> boundary;

  FE_Q<dim>            fe;
  DoFHandler<dim>      dof_handler;

  // This is the new variable in the main class. We need an object which holds
  // a list of constraints to hold the hanging nodes and the boundary
  // conditions.
  ConstraintMatrix     constraints;

  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;

  Vector<double>       solution;
  Vector<double>       system_rhs;
};


// @sect3{Nonconstant coefficients}

// The implementation of nonconstant coefficients is copied verbatim from
// step-5:
template <int dim>
double coefficient (const Point<dim> &p)
{
  if (p.square() < 0.5*0.5)
    return 20;
  else
    return 1;
}



// @sect3{The <code>Step6</code> class implementation}

// @sect4{Step6::Step6}

// The constructor of this class is mostly the same as before, but this time
// we want to use the quadratic element. To do so, we only have to replace the
// constructor argument (which was <code>1</code> in all previous examples) by
// the desired polynomial degree (here <code>2</code>):
template <int dim>
Step6<dim>::Step6 ()
  :
  fe (2),
  dof_handler (triangulation)
{}


// @sect4{Step6::~Step6}

// Here comes the added destructor of the class. The reason why we want to add
// it is to pick up an issue we already started to discuss in step-1:
// the <code>boundary</code> object was defined before and not after the
// <code>triangulation</code> object. Of course we could have left this order
// unchanged, but we would like to show what happens if the order is reversed
// since this produces a rather nasty side-effect and results in an error which
// is difficult to track down if one does not know what happens.
//
// Basically what happens is the following: when we set a manifold description
// to the triangulation using the function call
// <code>triangulation.set_manifold (0, boundary)</code>, the
// <code>Triangulation</code> object also stores a pointer to the
// <code>Manifold</code> object in use. Since this pointer is used until either
// another <code>Manifold</code> object is set as boundary description or until
// the <code>Triangulation</code> object is destroyed, it would be unwise if we
// would allow the <code>boundary</code> to be deleted before the
// <code>triangulation</code>. To disallow this, the <code>triangulation</code>
// increases a counter inside the <code>boundary</code> which counts how many
// objects use it (this is what the
// <code>Subscriptor</code>/<code>SmartPointer</code> class pair is used for,
// in case you want something like this for your own programs; see step-7 for
// a more complete discussion of this topic). The <code>boundary</code> will
// refuse its destruction if that counter is larger than zero, since then some
// other objects might rely on its persistence. An exception will then be
// thrown and the program will usually abort upon the attempt to destroy
// <code>boundary</code>.
//
// To be fair, such exceptions about still used objects are not particularly
// popular among programmers using deal.II, since they only tell us that
// something is wrong, namely that <i>some</i> other object is still using the object
// that is presently being destroyed, but most of the time not <i>which</i> object is
// still using it. It is therefore often rather time-consuming to find out where the
// problem exactly is, although it is then usually straightforward to remedy
// the situation. However, we believe that the effort to find invalid
// references to objects that do no longer exist is less if the problem is
// detected once the reference becomes invalid, rather than when non-existent
// objects are actually accessed again, since then usually only invalid data
// is accessed, but no error is immediately raised.
//
// Coming back to the present situation, if we did not write this destructor,
// the compiler will generate code that triggers exactly the behavior sketched
// above. The reason is that member variables of the <code>Step6</code> class
// are destroyed bottom-up (i.e., in reverse order of their declaration in the
// class), as always in C++. Thus, the boundary object will be
// destroyed before the triangulation object, since its declaration is below
// the one of the triangulation. This triggers the situation above, and an
{'fraction_non_alphanumeric': 0.10532646048109966, 'fraction_numerical': 0.005383734249713631, 'mean_word_length': 3.543585740307052, 'pattern_counts': {'":': 0, '<': 187, '<?xml version=': 0, '>': 185, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6115413', 'n_tokens_mistral': 4949, 'n_tokens_neox': 4720, 'n_words': 2244}