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

/*    $Id$       */
/*    Version: $Name$                                          */
/*                                                                */
/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

                                 // @sect3{Include files}

				 // The first few files have already
				 // been covered in previous examples
				 // and will thus not be further
				 // commented on.
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <base/logstream.h>
#include <base/utilities.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <grid/tria.h>
#include <dofs/hp_dof_handler.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/hp_fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>

#include <fstream>
#include <iostream>

#include <fe/fe_q.h>
#include <grid/grid_out.h>
#include <dofs/dof_constraints.h>
#include <grid/grid_refinement.h>
#include <numerics/error_estimator.h>

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

template <int dim>
class LaplaceProblem 
{
  public:
    LaplaceProblem ();
    ~LaplaceProblem ();

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

    Triangulation<dim>   triangulation;

    hp::DoFHandler<dim>      dof_handler;
    hp::FECollection<dim>    fe_collection;
    hp::QCollection<dim>     quadrature_collection;

    ConstraintMatrix     hanging_node_constraints;

    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;

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


template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
		dof_handler (triangulation)
{
  for (unsigned int degree=1; degree<5; ++degree)
    {
      fe_collection.push_back (FE_Q<dim>(degree));
      quadrature_collection.push_back (QGauss<dim>(degree+2));
    }
}


template <int dim>
LaplaceProblem<dim>::~LaplaceProblem () 
{
  dof_handler.clear ();
}

template <int dim>
void LaplaceProblem<dim>::setup_system ()
{
  dof_handler.distribute_dofs (fe_collection);

  sparsity_pattern.reinit (dof_handler.n_dofs(),
			   dof_handler.n_dofs(),
			   dof_handler.max_couplings_between_dofs());
  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);

  solution.reinit (dof_handler.n_dofs());
  system_rhs.reinit (dof_handler.n_dofs());

  hanging_node_constraints.clear ();
  DoFTools::make_hanging_node_constraints (dof_handler,
					   hanging_node_constraints);

  hanging_node_constraints.close ();

  hanging_node_constraints.condense (sparsity_pattern);
  sparsity_pattern.compress();

  system_matrix.reinit (sparsity_pattern);
}

template <int dim>
void LaplaceProblem<dim>::assemble_system () 
{
  hp::FEValues<dim> hp_fe_values (fe_collection,
				  quadrature_collection, 
				  update_values    |  update_gradients |
				  update_q_points  |  update_JxW_values);

  typename hp::DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
  for (; cell!=endc; ++cell)
    {
      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
      FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
      Vector<double>       cell_rhs (dofs_per_cell);

      std::vector<unsigned int> local_dof_indices (dofs_per_cell);

      cell_matrix = 0;
      cell_rhs = 0;

      hp_fe_values.reinit (cell);

      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
      
      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
	for (unsigned int i=0; i<dofs_per_cell; ++i)
	  {
	    for (unsigned int j=0; j<dofs_per_cell; ++j)
	      cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
				   fe_values.shape_grad(j,q_point) *
				   fe_values.JxW(q_point));

	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
			    1.0 *
			    fe_values.JxW(q_point));
	  }

      cell->get_dof_indices (local_dof_indices);
      for (unsigned int i=0; i<dofs_per_cell; ++i)
	{
	  for (unsigned int j=0; j<dofs_per_cell; ++j)
	    system_matrix.add (local_dof_indices[i],
			       local_dof_indices[j],
			       cell_matrix(i,j));
	  
	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
	}
    }

  hanging_node_constraints.condense (system_matrix);
  hanging_node_constraints.condense (system_rhs);
  std::map<unsigned int,double> boundary_values;
  VectorTools::interpolate_boundary_values (dof_handler,
					    0,
					    ZeroFunction<dim>(),
					    boundary_values);
  MatrixTools::apply_boundary_values (boundary_values,
				      system_matrix,
				      solution,
				      system_rhs);
}

template <int dim>
void LaplaceProblem<dim>::solve () 
{
  SolverControl           solver_control (1000, 1e-12);
  SolverCG<>              cg (solver_control);

  PreconditionSSOR<> preconditioner;
  preconditioner.initialize(system_matrix, 1.2);

  cg.solve (system_matrix, solution, system_rhs,
	    preconditioner);

  hanging_node_constraints.distribute (solution);
}

template <int dim>
void LaplaceProblem<dim>::refine_grid ()
{
  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

  KellyErrorEstimator<dim>::estimate (dof_handler,
				      QGauss<dim-1>(3),
				      typename FunctionMap<dim>::type(),
				      solution,
				      estimated_error_per_cell);

  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
						   estimated_error_per_cell,
						   0.3, 0.03);

  triangulation.execute_coarsening_and_refinement ();
}

template <int dim>
void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
{
  Assert (cycle < 10, ExcNotImplemented());
  
  {
    const std::string filename = "grid-" +
				 Utilities::int_to_string (cycle, 2) +
				 ".eps";
    std::ofstream output (filename.c_str());
    
    GridOut grid_out;
    grid_out.write_eps (triangulation, output);
  }
  
  {
    const std::string filename = "solution-" +
				 Utilities::int_to_string (cycle, 2) +
				 ".gnuplot";
    DataOut<dim,hp::DoFHandler<dim> > data_out;

    data_out.attach_dof_handler (dof_handler);
    data_out.add_data_vector (solution, "solution");
    data_out.build_patches ();
  
    std::ofstream output (filename.c_str());
    data_out.write_gnuplot (output);
  }
}


template <>
void LaplaceProblem<2>::create_coarse_grid ()
{
  const unsigned int dim = 2;
  
  static const Point<2> vertices_1[]
    = {  Point<2> (-1.,   -1.),
         Point<2> (-1./2, -1.),
         Point<2> (0.,    -1.),
         Point<2> (+1./2, -1.),
         Point<2> (+1,    -1.),
	     
         Point<2> (-1.,   -1./2.),
         Point<2> (-1./2, -1./2.),
         Point<2> (0.,    -1./2.),
         Point<2> (+1./2, -1./2.),
         Point<2> (+1,    -1./2.),
	     
         Point<2> (-1.,   0.),
         Point<2> (-1./2, 0.),
         Point<2> (+1./2, 0.),
         Point<2> (+1,    0.),
	     
         Point<2> (-1.,   1./2.),
         Point<2> (-1./2, 1./2.),
         Point<2> (0.,    1./2.),
         Point<2> (+1./2, 1./2.),
         Point<2> (+1,    1./2.),
	     
         Point<2> (-1.,   1.),
         Point<2> (-1./2, 1.),
         Point<2> (0.,    1.),			  
         Point<2> (+1./2, 1.),
         Point<2> (+1,    1.)    };
  const unsigned int
    n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
  const std::vector<Point<dim> > vertices (&vertices_1[0],
                                           &vertices_1[n_vertices]);
  static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
    = {{0, 1, 5, 6},
       {1, 2, 6, 7},
       {2, 3, 7, 8},
       {3, 4, 8, 9},
       {5, 6, 10, 11},
       {8, 9, 12, 13},
       {10, 11, 14, 15},
       {12, 13, 17, 18},
       {14, 15, 19, 20},
       {15, 16, 20, 21},
       {16, 17, 21, 22},
       {17, 18, 22, 23}};
  const unsigned int
    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);

  std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
  for (unsigned int i=0; i<n_cells; ++i) 
    {
      for (unsigned int j=0;
           j<GeometryInfo<dim>::vertices_per_cell;
           ++j)
        cells[i].vertices[j] = cell_vertices[i][j];
      cells[i].material_id = 0;
    }

  triangulation.create_triangulation (vertices,
                                    cells,
                                    SubCellData());
  triangulation.refine_global (1);
}



template <int dim>
void LaplaceProblem<dim>::run () 
{
  for (unsigned int cycle=0; cycle<5; ++cycle)
    {
      std::cout << "Cycle " << cycle << ':' << std::endl;

      if (cycle == 0)
	create_coarse_grid ();
      else
	refine_grid ();
      

      std::cout << "   Number of active cells:       "
		<< triangulation.n_active_cells()
		<< std::endl;

      setup_system ();

      std::cout << "   Number of degrees of freedom: "
		<< dof_handler.n_dofs()
		<< std::endl;
      
      assemble_system ();
      solve ();
      output_results (cycle);
    }
}

int main () 
{
  try
    {
      deallog.depth_console (0);

      LaplaceProblem<2> laplace_problem_2d;
      laplace_problem_2d.run ();
    }
  catch (std::exception &exc)<commit_msg>Make the program do something useful, and do it without bombing out...<commit_after>/* $Id$ */
/* Author: <NAME>, University of Heidelberg, 2000 */

/*    $Id$       */
/*    Version: $Name$                                          */
/*                                                                */
/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

                                 // @sect3{Include files}

				 // The first few files have already
				 // been covered in previous examples
				 // and will thus not be further
				 // commented on.
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <base/logstream.h>
#include <base/utilities.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <grid/tria.h>
#include <dofs/hp_dof_handler.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/hp_fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>

#include <fstream>
#include <iostream>

#include <fe/fe_q.h>
#include <grid/grid_out.h>
#include <dofs/dof_constraints.h>
#include <grid/grid_refinement.h>
#include <numerics/error_estimator.h>

#include <complex>

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

template <int dim>
class LaplaceProblem 
{
  public:
    LaplaceProblem ();
    ~LaplaceProblem ();

    void run ();
    
  private:
    void setup_system ();
    void assemble_system ();
    void solve ();
    void create_coarse_grid ();
    void refine_grid ();
    void estimate_smoothness (Vector<float> &smoothness_indicators) const;
    void output_results (const unsigned int cycle) const;

    Triangulation<dim>   triangulation;

    hp::DoFHandler<dim>      dof_handler;
    hp::FECollection<dim>    fe_collection;
    hp::QCollection<dim>     quadrature_collection;
    hp::QCollection<dim-1>   face_quadrature_collection;

    ConstraintMatrix     hanging_node_constraints;

    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;

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



template <int dim>
class RightHandSide : public Function<dim>
{
  public:
    RightHandSide () : Function<dim> () {};
    
    virtual double value (const Point<dim>   &p,
			  const unsigned int  component) const;
};


template <int dim>
double
RightHandSide<dim>::value (const Point<dim>   &p,
			   const unsigned int  /*component*/) const
{
  double product = 1;
  for (unsigned int d=0; d<dim; ++d)
    product *= (p[d]+1);
  return product;
}




template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
		dof_handler (triangulation)
{
  for (unsigned int degree=2; degree<7; ++degree)
    {
      fe_collection.push_back (FE_Q<dim>(degree));
      quadrature_collection.push_back (QGauss<dim>(degree+2));
      face_quadrature_collection.push_back (QGauss<dim-1>(degree+2));
    }
}


template <int dim>
LaplaceProblem<dim>::~LaplaceProblem () 
{
  dof_handler.clear ();
}

template <int dim>
void LaplaceProblem<dim>::setup_system ()
{
  dof_handler.distribute_dofs (fe_collection);

  sparsity_pattern.reinit (dof_handler.n_dofs(),
			   dof_handler.n_dofs(),
			   dof_handler.max_couplings_between_dofs());
  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);

  solution.reinit (dof_handler.n_dofs());
  system_rhs.reinit (dof_handler.n_dofs());

  hanging_node_constraints.clear ();
  DoFTools::make_hanging_node_constraints (dof_handler,
					   hanging_node_constraints);

  hanging_node_constraints.close ();

  hanging_node_constraints.condense (sparsity_pattern);
  sparsity_pattern.compress();

  system_matrix.reinit (sparsity_pattern);
}

template <int dim>
void LaplaceProblem<dim>::assemble_system () 
{
  hp::FEValues<dim> hp_fe_values (fe_collection,
				  quadrature_collection, 
				  update_values    |  update_gradients |
				  update_q_points  |  update_JxW_values);

  const RightHandSide<dim> rhs_function;
  
  typename hp::DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
  for (; cell!=endc; ++cell)
    {
      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
      FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
      Vector<double>       cell_rhs (dofs_per_cell);

      std::vector<unsigned int> local_dof_indices (dofs_per_cell);

      cell_matrix = 0;
      cell_rhs = 0;

      hp_fe_values.reinit (cell);
      
      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();

      std::vector<double>  rhs_values (fe_values.n_quadrature_points);
      rhs_function.value_list (fe_values.get_quadrature_points(),
			       rhs_values);
      
      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
	for (unsigned int i=0; i<dofs_per_cell; ++i)
	  {
	    for (unsigned int j=0; j<dofs_per_cell; ++j)
	      cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
				   fe_values.shape_grad(j,q_point) *
				   fe_values.JxW(q_point));

	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
			    rhs_values[q_point] *
			    fe_values.JxW(q_point));
	  }

      cell->get_dof_indices (local_dof_indices);
      for (unsigned int i=0; i<dofs_per_cell; ++i)
	{
	  for (unsigned int j=0; j<dofs_per_cell; ++j)
	    system_matrix.add (local_dof_indices[i],
			       local_dof_indices[j],
			       cell_matrix(i,j));
	  
	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
	}
    }

  hanging_node_constraints.condense (system_matrix);
  hanging_node_constraints.condense (system_rhs);
  std::map<unsigned int,double> boundary_values;
  VectorTools::interpolate_boundary_values (dof_handler,
					    0,
					    ZeroFunction<dim>(),
					    boundary_values);
  MatrixTools::apply_boundary_values (boundary_values,
				      system_matrix,
				      solution,
				      system_rhs);
}

template <int dim>
void LaplaceProblem<dim>::solve () 
{
  SolverControl           solver_control (1000, 1e-12);
  SolverCG<>              cg (solver_control);

  PreconditionSSOR<> preconditioner;
  preconditioner.initialize(system_matrix, 1.2);

  cg.solve (system_matrix, solution, system_rhs,
	    preconditioner);

  hanging_node_constraints.distribute (solution);
}


unsigned int
int_pow (const unsigned int x,
	 const unsigned int n)
{
  unsigned int p=1;
  for (unsigned int i=0; i<n; ++i)
    p *= x;
  return p;
}


template <int dim>
void
LaplaceProblem<dim>::
estimate_smoothness (Vector<float> &smoothness_indicators) const
{
  const unsigned int N = 5;

  const unsigned n_fourier_modes = int_pow (N+1, dim);
  std::vector<Tensor<1,dim> > k_vectors (n_fourier_modes);
  for (unsigned int k=0; k<n_fourier_modes; ++k)
    switch (dim)
      {
	case 2:
	      k_vectors[k][0] = deal_II_numbers::PI * (k / (N+1));
	      k_vectors[k][1] = deal_II_numbers::PI * (k % (N+1));
	      break;
	default:
	      Assert (false, ExcNotImplemented());
      }
  
  QGauss<1>      base_quadrature (2);
  QIterated<dim> quadrature (base_quadrature, N);
  
  std::vector<Table<2,std::complex<double> > >
    fourier_transforms (fe_collection.size());
  for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
    {
      fourier_transforms[fe].reinit (n_fourier_modes,
				     fe_collection[fe].dofs_per_cell);

      for (unsigned int k=0; k<n_fourier_modes; ++k)
	for (unsigned int i=0; i<fe_collection[fe].dofs_per_cell; ++i)
	  {
	    std::complex<double> sum = 0;
	    for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
	      {
		const Point<dim> x_q = quadrature.point(q);
		sum += std::exp(std::complex<double>(0,1) *
				(k_vectors[k] * x_q)) *
		       fe_collection[fe].shape_value(i,x_q) *
		       quadrature.weight(q);
	      }
	    fourier_transforms[fe](k,i) = sum;
	  }
    }

  typename hp::DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
  std::vector<std::complex<double> > transformed_values (n_fourier_modes);
  for (unsigned int index=0; cell!=endc; ++cell, ++index)
    {
      Vector<double> dof_values (cell->get_fe().dofs_per_cell);
      cell->get_dof_values (solution, dof_values);

      for (unsigned int f=0; f<n_fourier_modes; ++f)
	{
	  transformed_values[f] = 0;
	  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
	    transformed_values[f] +=
	      fourier_transforms[cell->active_fe_index()](f,i)
	      *
	      dof_values(i);
	}

				       // now try to fit a curve C|k|^{-s} to
				       // the fourier transform of the
				       // solution on this cell
      double A[2][2] = { { 0,0 }, { 0,0 }};
      double F[2]    = { 0,0 };
      
      for (unsigned int f=0; f<n_fourier_modes; ++f)
	{
	  const double k_abs = std::sqrt(k_vectors[f] *
					 k_vectors[f]);

	  if (k_abs == 0)
	    continue;
	  
	  A[0][0] += 1;
	  A[1][0] += std::log (k_abs);
	  A[1][1] += std::pow (std::log (k_abs), 2.);

	  F[0] += std::log (std::abs (transformed_values[f]));
	  F[1] += std::log (std::abs (transformed_values[f])) *
		  std::log (k_abs);
	}
      A[0][1] = A[1][0];
      
      const double det = A[0][0] * A[1][1] - A[0][1] * A[0][1];
      const double s = (A[0][0]*F[1] - A[1][1]*F[0]) / det;
      
      smoothness_indicators(index) = s;
    }
}



  
template <int dim>
void LaplaceProblem<dim>::refine_grid ()
{
  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
  KellyErrorEstimator<dim>::estimate (dof_handler,
				      face_quadrature_collection,
				      typename FunctionMap<dim>::type(),
				      solution,
				      estimated_error_per_cell);

  Vector<float> smoothness_indicators (triangulation.n_active_cells());
  estimate_smoothness (smoothness_indicators);
  
  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
						   estimated_error_per_cell,
						   0.3, 0.03);

  float max_smoothness = 0,
	 min_smoothness = smoothness_indicators.linfty_norm();
  {
    typename hp::DoFHandler<dim>::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();
    for (unsigned int index=0; cell!=endc; ++cell, ++index)
      if (cell->refine_flag_set())
	{
	  max_smoothness = std::max (max_smoothness,
				     smoothness_indicators(index));
	  min_smoothness = std::min (min_smoothness,
				     smoothness_indicators(index));
	}
  }
  const float cutoff_smoothness = (max_smoothness + min_smoothness) / 2;
  {
    typename hp::DoFHandler<dim>::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();
    for (unsigned int index=0; cell!=endc; ++cell, ++index)
      if (cell->refine_flag_set()
	  &&
	  (smoothness_indicators(index) > cutoff_smoothness))
	{
	  cell->clear_refine_flag();
	  cell->set_active_fe_index (std::min (cell->active_fe_index() + 1,
					       fe_collection.size() - 1));
	}
  } 
  
  triangulation.execute_coarsening_and_refinement ();
}



template <int dim>
void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
{
  Assert (cycle < 10, ExcNotImplemented());
  
  {
    const std::string filename = "grid-" +
				 Utilities::int_to_string (cycle, 2) +
				 ".eps";
    std::ofstream output (filename.c_str());
    
    GridOut grid_out;
    grid_out.write_eps (triangulation, output);
  }
  
  {
    Vector<float> smoothness_indicators (triangulation.n_active_cells());
    estimate_smoothness (smoothness_indicators);


    Vector<float> fe_indices (triangulation.n_active_cells());
    {
      typename hp::DoFHandler<dim>::active_cell_iterator
	cell = dof_handler.begin_active(),
	endc = dof_handler.end();
      for (unsigned int index=0; cell!=endc; ++cell, ++index)
	{
	  
	  fe_indices(index) = cell->active_fe_index();
//	  smoothness_indicators(index) *= std::sqrt(cell->diameter());
	}
    }
    
    const std::string filename = "solution-" +
				 Utilities::int_to_string (cycle, 2) +
				 ".vtk";
    DataOut<dim,hp::DoFHandler<dim> > data_out;

    data_out.attach_dof_handler (dof_handler);
    data_out.add_data_vector (solution, "solution");
    data_out.add_data_vector (smoothness_indicators, "smoothness");
    data_out.add_data_vector (fe_indices, "fe_index");
    data_out.build_patches ();
  
    std::ofstream output (filename.c_str());
    data_out.write_vtk (output);
  }
}


template <>
void LaplaceProblem<2>::create_coarse_grid ()
{
  const unsigned int dim = 2;
  
  static const Point<2> vertices_1[]
    = {  Point<2> (-1.,   -1.),
         Point<2> (-1./2, -1.),
         Point<2> (0.,    -1.),
         Point<2> (+1./2, -1.),
         Point<2> (+1,    -1.),
	     
         Point<2> (-1.,   -1./2.),
         Point<2> (-1./2, -1./2.),
         Point<2> (0.,    -1./2.),
         Point<2> (+1./2, -1./2.),
         Point<2> (+1,    -1./2.),
	     
         Point<2> (-1.,   0.),
         Point<2> (-1./2, 0.),
         Point<2> (+1./2, 0.),
         Point<2> (+1,    0.),
	     
         Point<2> (-1.,   1./2.),
         Point<2> (-1./2, 1./2.),
         Point<2> (0.,    1./2.),
         Point<2> (+1./2, 1./2.),
         Point<2> (+1,    1./2.),
	     
         Point<2> (-1.,   1.),
         Point<2> (-1./2, 1.),
         Point<2> (0.,    1.),			  
         Point<2> (+1./2, 1.),
         Point<2> (+1,    1.)    };
  const unsigned int
    n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
  const std::vector<Point<dim> > vertices (&vertices_1[0],
                                           &vertices_1[n_vertices]);
  static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
    = {{0, 1, 5, 6},
       {1, 2, 6, 7},
       {2, 3, 7, 8},
       {3, 4, 8, 9},
       {5, 6, 10, 11},
       {8, 9, 12, 13},
       {10, 11, 14, 15},
       {12, 13, 17, 18},
       {14, 15, 19, 20},
       {15, 16, 20, 21},
       {16, 17, 21, 22},
       {17, 18, 22, 23}};
  const unsigned int
    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);

  std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
  for (unsigned int i=0; i<n_cells; ++i) 
    {
      for (unsigned int j=0;
           j<GeometryInfo<dim>::vertices_per_cell;
           ++j)
        cells[i].vertices[j] = cell_vertices[i][j];
      cells[i].material_id = 0;
    }

  triangulation.create_triangulation (vertices,
                                    cells,
                                    SubCellData());
  triangulation.refine_global (3);
}



template <int dim>
void LaplaceProblem<dim>::run () 
{
  for (unsigned int cycle=0; cycle<8; ++cycle)
    {
      std::cout << "Cycle " << cycle << ':' << std::endl;

      if (cycle == 0)
	create_coarse_grid ();
      else
	refine_grid ();
      

      std::cout << "   Number of active cells:       "
		<< triangulation.n_active_cells()
		<< std::endl;

      setup_system ();

      std::cout << "   Number of degrees of freedom: "
		<< dof_handler.n_dofs()
		<< std::endl;
      std::cout << "   Number of constraints       : "
		<< hanging_node_constraints.n_constraints()
		<< std::endl;
      
      assemble_system ();
      solve ();
      output_results (cycle);
    }
}

int main () 
{
  try
    {
      deallog.depth_console (0);

      LaplaceProblem<2> laplace_problem_2d;
      laplace_problem_2d.run ();
    }
  catch (std::exception &exc)
{'fraction_non_alphanumeric': 0.14392501646074596, 'fraction_numerical': 0.022347883341725085, 'mean_word_length': 2.7398609501738123, 'pattern_counts': {'":': 0, '<': 326, '<?xml version=': 0, '>': 266, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5903773', 'n_tokens_mistral': 9979, 'n_tokens_neox': 9061, 'n_words': 1760}