Minimal example to understand for mass and laplace matrices

main.cpp

/*
 * Minimal example to understand for mass and laplace matrices
 */

// The program starts with the usual include files, all of which you should
// have seen before by now:
#include <deal.II/base/utilities.h>
#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/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>

#include <fstream>
#include <iostream>

int main()
{
  using namespace dealii;

  FE_Q<2> fe(1);

  Triangulation<2> mesh;
  GridGenerator::subdivided_hyper_cube(mesh, 1);

  DoFHandler<2> dof_handler(mesh);
  dof_handler.distribute_dofs(fe);

  SparsityPattern sparsity_pattern;
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);

  std::ofstream gnu("sparsity_pattern.gnuplot");
  sparsity_pattern.print_gnuplot(gnu);

  std::ofstream svg("sparsity_pattern.svg");
  sparsity_pattern.print_svg(svg);

  SparseMatrix<double> mass_matrix;
  SparseMatrix<double> laplace_matrix;
  mass_matrix.reinit(sparsity_pattern);
  laplace_matrix.reinit(sparsity_pattern);

  MatrixCreator::create_mass_matrix(dof_handler,
                                    QGauss<2>(fe.degree + 1),
                                    mass_matrix);

  MatrixCreator::create_laplace_matrix(dof_handler,
                                       QGauss<2>(fe.degree + 1),
                                       laplace_matrix);

  std::cout << "mass matrix:" << std::endl;
  mass_matrix.print_formatted(std::cout, 2, false, 5, "0.0");

  std::cout << "laplace matrix:" << std::endl;
  laplace_matrix.print_formatted(std::cout, 2, false, 5, "0.0");
  return 0;
}

mass matrix:
 0.11  0.06  0.06  0.03 
 0.06  0.11  0.03  0.06 
 0.06  0.03  0.11  0.06 
 0.03  0.06  0.06  0.11 

laplace matrix:
 0.67 -0.17 -0.17 -0.33 
-0.17  0.67 -0.33 -0.17 
-0.17 -0.33  0.67 -0.17 
-0.33 -0.17 -0.17  0.67 

• Mass matrix line sum: 0.26
• Laplace matrix line sum: 0.0