/* ------------------------------------------------------------------------
 *
 * SPDX-License-Identifier: LGPL-2.1-or-later
 * Copyright (C) 1999 - 2024 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * Part of the source code is dual licensed under Apache-2.0 WITH
 * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
 * governing the source code and code contributions can be found in
 * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
 *
 * ------------------------------------------------------------------------
 */
 
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/mapping_q1.h>
 
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
 
#include <deal.II/dofs/dof_renumbering.h>
 
#include <fstream>
 
using namespace dealii;
 
 
 
void make_grid(Triangulation<2> &triangulation)
{
  const Point<2> center(1, 0);
  const double   inner_radius = 0.5, outer_radius = 1.0;
  GridGenerator::hyper_shell(
    triangulation, center, inner_radius, outer_radius, 5);
 
  for (unsigned int step = 0; step < 3; ++step)
    {
      for (const auto &cell : triangulation.active_cell_iterators())
        for (const auto v : cell->vertex_indices())
          {
            const double distance_from_center =
              center.distance(cell->vertex(v));
 
            if (std::fabs(distance_from_center - inner_radius) <=
                1e-6 * inner_radius)
              {
                cell->set_refine_flag();
                break;
              }
          }
 
      triangulation.execute_coarsening_and_refinement();
    }
 
  std::ofstream mesh_file("mesh.gnuplot");
  GridOut().write_gnuplot(triangulation, mesh_file);
}
 
 
 
void write_dof_locations(const DoFHandler<2> &dof_handler,
                         const std::string   &filename)
{
  const std::map<types::global_dof_index, Point<2>> dof_location_map =
    DoFTools::map_dofs_to_support_points(MappingQ1<2>(), dof_handler);
 
  std::ofstream dof_location_file(filename);
  DoFTools::write_gnuplot_dof_support_point_info(dof_location_file,
                                                 dof_location_map);
}
 
 
 
void distribute_dofs(DoFHandler<2> &dof_handler)
{
  const FE_Q<2> finite_element(1);
  dof_handler.distribute_dofs(finite_element);
 
  write_dof_locations(dof_handler, "dof-locations-1.gnuplot");
 
  DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
                                                  dof_handler.n_dofs());
 
  DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern);
 
  SparsityPattern sparsity_pattern;
  sparsity_pattern.copy_from(dynamic_sparsity_pattern);
 
  std::ofstream out("sparsity-pattern-1.svg");
  sparsity_pattern.print_svg(out);
}
 
 
 
void renumber_dofs(DoFHandler<2> &dof_handler)
{
  DoFRenumbering::Cuthill_McKee(dof_handler);
 
  write_dof_locations(dof_handler, "dof-locations-2.gnuplot");
 
 
  DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
                                                  dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern);
 
  SparsityPattern sparsity_pattern;
  sparsity_pattern.copy_from(dynamic_sparsity_pattern);
 
  std::ofstream out("sparsity-pattern-2.svg");
  sparsity_pattern.print_svg(out);
}
 
 
 
 
 
int main()
{
  cout << "\n";
  cout << "step-2():\n";
  cout << "  dealii example code.\n";
  cout << "\n";

  Triangulation<2> triangulation;
  make_grid(triangulation);
 
  DoFHandler<2> dof_handler(triangulation);
 
  distribute_dofs(dof_handler);
  renumber_dofs(dof_handler);

  cout << "\n";
  cout << "step-2():\n";
  cout << "  Normal end of execution.\n";

  return 0;
}