marcfehling / mg-ev-estimator Goto Github PK

#include <deal.II/base/convergence_table.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/parameter_handler.h>

#include <deal.II/distributed/cell_weights.h>
#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_q_iso_q1.h>
#include <deal.II/fe/fe_system.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/hp/fe_values.h>

#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_creator.h>

using namespace dealii;

const bool manual_constraint = false;
const unsigned int fe_degree = 1;

/**
  Test contraints between a locally refined cell and a cell with
  two subdivisions.

  The two DoFs in the center of the shared face are constrained
  by the other ones:
     2 11:  0.5
     2 15:  0.5
    12 11:  0.5
    12 15:  0.5
 */
int
main(int argc, char *argv[])
{
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

  const unsigned int dim = 2;

  // 1) create system
  hp::FECollection<dim> fe;
  fe.push_back(FE_Q_iso_Q1<dim>(fe_degree));
  fe.push_back(FE_Q_iso_Q1<dim>(fe_degree + 1));

  hp::QCollection<dim> quadrature_collection;
  quadrature_collection.push_back(QIterated<dim>(QGauss<1>(2), fe_degree));
  quadrature_collection.push_back(QIterated<dim>(QGauss<1>(2), fe_degree + 1));

  hp::QCollection<dim> quadrature;

  for (const auto &f : fe)
    quadrature.push_back(Quadrature<dim>(f.get_unit_support_points()));

  Triangulation<dim> tria;
  GridGenerator::subdivided_hyper_rectangle(tria,
                                            {2, 1},
                                            Point<dim>(0.0, 0.0),
                                            Point<dim>(2.0, 1.0));

  tria.begin()->set_refine_flag();
  tria.execute_coarsening_and_refinement();

  DoFHandler<dim> dof_handler(tria);

  for (const auto &cell : dof_handler.active_cell_iterators())
    if (cell->center()[0] < 1.0)
      cell->set_active_fe_index(0);
    else
      cell->set_active_fe_index(1);

  dof_handler.distribute_dofs(fe);

  MappingQ1<dim> mapping;

  // 2) match support points with DoF indices
  hp::FEValues<dim> fe_values_hp(fe, quadrature, update_quadrature_points);

  std::vector<types::global_dof_index> dof_indices;

  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      fe_values_hp.reinit(cell);
      const auto &fe_values = fe_values_hp.get_present_fe_values();

      dof_indices.resize(fe_values.dofs_per_cell);

      cell->get_dof_indices(dof_indices);

      for (const auto q : fe_values.quadrature_point_indices())
        std::cout << dof_indices[q] << ": " << fe_values.quadrature_point(q)
                  << std::endl;

      std::cout << std::endl;
    }

  std::cout << dof_handler.n_dofs() << std::endl;

  // set up constraints
  AffineConstraints<double> constraints;
  if(true)
    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  else
    {
      constraints.add_line(12);
      constraints.add_entry(12, 2, 1.0);
    }
  constraints.close();

  constraints.print(std::cout);
  std::cout << std::endl;

  DynamicSparsityPattern dsp(dof_handler.n_dofs());
      DoFTools::make_sparsity_pattern(dof_handler,
                                      dsp,
                                      constraints);

  SparsityPattern sparsity_pattern;
      sparsity_pattern.copy_from(dsp);

  SparseMatrix<double> matrix;
      matrix.reinit(sparsity_pattern);
      MatrixCreator::create_laplace_matrix<dim, dim>(
        dof_handler,
        quadrature_collection,
        matrix,
        nullptr,
        constraints);

  FullMatrix<double> full_matrix(dof_handler.n_dofs(),dof_handler.n_dofs());
  full_matrix.copy_from(matrix);

  full_matrix.print_formatted(std::cout, 2, false, 7, "0.00");
  std::cout << std::endl;

  LAPACKFullMatrix<double> lapack_full_matrix;
  lapack_full_matrix.copy_from(full_matrix);
  lapack_full_matrix.compute_eigenvalues();

  std::vector<double> eigenvalues;

  for(unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
    eigenvalues.push_back(lapack_full_matrix.eigenvalue(i).real());
  
  std::sort(eigenvalues.begin(), eigenvalues.end());

  for(const auto i : eigenvalues)
    std::cout << i << " ";
  std::cout <<std::endl;

}