<commit_before>// ---------------------------------------------------------------------
// $Id$
//
// Copyright (C) 2012 - 2013 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.
//
// ---------------------------------------------------------------------


/**
 * @file Test whether 
 */

#include "../tests.h"
#include <deal.II/base/function.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>

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

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

#include <deal.II/integrators/laplace.h>

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/local_integrator.h>
#include <deal.II/meshworker/simple.h>
#include <deal.II/meshworker/loop.h>

#include <deal.II/multigrid/sparse_matrix_collection.h>

using namespace dealii;
using namespace LocalIntegrators;


template <int dim>
class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
{
public:
  LaplaceMatrix();
  virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
		    MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const;
};


template <int dim>
LaplaceMatrix<dim>::LaplaceMatrix()
{}


template <int dim>
void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
{
  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
}


template <int dim>
void LaplaceMatrix<dim>::boundary(MeshWorker::DoFInfo<dim>& dinfo,
				       typename MeshWorker::IntegrationInfo<dim>& info) const
{
  const unsigned int deg = info.fe_values().get_fe().tensor_degree();
  Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(),
			      Laplace::compute_penalty(dinfo, dinfo, deg, deg));
}


template <int dim>
void LaplaceMatrix<dim>::face(
  MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
  MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const
{
  if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::H1))
    return;
  
  const unsigned int deg = info1.fe_values().get_fe().tensor_degree();
  
  if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hdiv) &&
	  !info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hcurl))
	Laplace::ip_tangential_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 
				      dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
				      info1.fe_values(), info2.fe_values(),
				      Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
      else
	Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 
			   dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
			   info1.fe_values(), info2.fe_values(),
			   Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
}


template <int dim>
void assemble_mg_matrix(DoFHandler<dim>& dof_handler,
   MeshWorker::LocalIntegrator<dim>& matrix_integrator, mg::SparseMatrixCollection<double>& mg)
{
  MGConstrainedDoFs mg_constraints;
  mg_constraints.clear();
  mg_constraints.initialize(dof_handler);
  
  mg.set_zero();

  MappingQ1<dim> mapping;
  
  MeshWorker::IntegrationInfoBox<dim> info_box;
  UpdateFlags update_flags = update_values | update_gradients | update_hessians;
  info_box.add_update_flags_all(update_flags);
  info_box.initialize(dof_handler.get_fe(), mapping);

  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
  
  MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
  assembler.initialize(mg_constraints);
  assembler.initialize(mg.matrix);
  assembler.initialize_interfaces(mg.matrix_in, mg.matrix_out);
  assembler.initialize_fluxes(mg.matrix_up, mg.matrix_down);
  
  MeshWorker::integration_loop<dim, dim> (
    dof_handler.begin_mg(), dof_handler.end_mg(),
    dof_info, info_box, matrix_integrator, assembler);

  const unsigned int nlevels = dof_handler.get_tria().n_levels();
  for (unsigned int level=0;level<nlevels;++level)
  {
    for(unsigned int i=0; i<dof_handler.n_dofs(level); ++i)
      if(mg.matrix[level].diag_element(i)==0)
        mg.matrix[level].set(i,i,1.);
  }
}


template <int dim>
void test(FiniteElement<dim> &fe)
{
  deallog << fe.get_name() << std::endl;

  Triangulation<dim> tr;
  std::vector<unsigned int> repititions (dim, 1);
  repititions[0] = 2;
  const Point<dim> p1 = (dim == 1 ? Point<dim> (-1.) : (dim == 2 ? Point<dim>(-1.,-1.) : Point<dim> (-1.,-1.,-1.)));
  const Point<dim> p2 = (dim == 1 ? Point<dim> (1.) : (dim == 2 ? Point<dim>(1.,1.) : Point<dim> (1.,1.,1.)));
  GridGenerator::subdivided_hyper_rectangle(tr, repititions, p1, p2);
  tr.begin_active()->set_refine_flag();
  tr.execute_coarsening_and_refinement();
  
  LaplaceMatrix<dim> matrix_integrator;
  
  MGDoFHandler<dim> dof(tr);
  dof.distribute_dofs(fe);
  mg::SparseMatrixCollection<double> mg;
  mg.resize(0,tr.n_levels()-1);
  mg.reinit(dof);
  assemble_mg_matrix(dof, matrix_integrator, mg);
}

int main()
{
  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
  std::ofstream logfile(logname.c_str());
  deallog.attach(logfile);
  deallog.depth_console (0);<commit_msg>print some output<commit_after>// ---------------------------------------------------------------------
// $Id$
//
// Copyright (C) 2012 - 2013 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.
//
// ---------------------------------------------------------------------


/**
 * @file Test whether 
 */

#include "../tests.h"
#include <deal.II/base/function.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>

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

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

#include <deal.II/integrators/laplace.h>

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/local_integrator.h>
#include <deal.II/meshworker/simple.h>
#include <deal.II/meshworker/loop.h>

#include <deal.II/multigrid/sparse_matrix_collection.h>

using namespace dealii;
using namespace LocalIntegrators;


template <int dim>
class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
{
public:
  LaplaceMatrix();
  virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
		    MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const;
};


template <int dim>
LaplaceMatrix<dim>::LaplaceMatrix()
{}


template <int dim>
void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
{
  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
}


template <int dim>
void LaplaceMatrix<dim>::boundary(MeshWorker::DoFInfo<dim>& dinfo,
				       typename MeshWorker::IntegrationInfo<dim>& info) const
{
  const unsigned int deg = info.fe_values().get_fe().tensor_degree();
  Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(),
			      Laplace::compute_penalty(dinfo, dinfo, deg, deg));
}


template <int dim>
void LaplaceMatrix<dim>::face(
  MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
  MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const
{
  if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::H1))
    return;
  
  const unsigned int deg = info1.fe_values().get_fe().tensor_degree();
  
  if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hdiv) &&
	  !info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hcurl))
	Laplace::ip_tangential_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 
				      dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
				      info1.fe_values(), info2.fe_values(),
				      Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
      else
	Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 
			   dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
			   info1.fe_values(), info2.fe_values(),
			   Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
}


template <int dim>
void assemble_mg_matrix(DoFHandler<dim>& dof_handler,
   MeshWorker::LocalIntegrator<dim>& matrix_integrator, mg::SparseMatrixCollection<double>& mg)
{
  MGConstrainedDoFs mg_constraints;
  mg_constraints.clear();
  mg_constraints.initialize(dof_handler);
  
  mg.set_zero();

  MappingQ1<dim> mapping;
  
  MeshWorker::IntegrationInfoBox<dim> info_box;
  UpdateFlags update_flags = update_values | update_gradients | update_hessians;
  info_box.add_update_flags_all(update_flags);
  info_box.initialize(dof_handler.get_fe(), mapping);

  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
  
  MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
  assembler.initialize(mg_constraints);
  assembler.initialize(mg.matrix);
  assembler.initialize_interfaces(mg.matrix_in, mg.matrix_out);
  assembler.initialize_fluxes(mg.matrix_up, mg.matrix_down);
  
  MeshWorker::integration_loop<dim, dim> (
    dof_handler.begin_mg(), dof_handler.end_mg(),
    dof_info, info_box, matrix_integrator, assembler);

  const unsigned int nlevels = dof_handler.get_tria().n_levels();
  for (unsigned int level=0;level<nlevels;++level)
  {
    for(unsigned int i=0; i<dof_handler.n_dofs(level); ++i)
      if(mg.matrix[level].diag_element(i)==0)
        mg.matrix[level].set(i,i,1.);
  }
}


template <int dim>
void test(FiniteElement<dim> &fe)
{
  deallog << fe.get_name() << std::endl;

  Triangulation<dim> tr;
  std::vector<unsigned int> repititions (dim, 1);
  repititions[0] = 2;
  const Point<dim> p1 = (dim == 1 ? Point<dim> (-1.) : (dim == 2 ? Point<dim>(-1.,-1.) : Point<dim> (-1.,-1.,-1.)));
  const Point<dim> p2 = (dim == 1 ? Point<dim> (1.) : (dim == 2 ? Point<dim>(1.,1.) : Point<dim> (1.,1.,1.)));
  GridGenerator::subdivided_hyper_rectangle(tr, repititions, p1, p2);
  tr.begin_active()->set_refine_flag();
  tr.execute_coarsening_and_refinement();
  
  LaplaceMatrix<dim> matrix_integrator;
  
  MGDoFHandler<dim> dof(tr);
  dof.distribute_dofs(fe);
  mg::SparseMatrixCollection<double> mg;
  mg.resize(0,tr.n_levels()-1);
  mg.reinit(dof);
  assemble_mg_matrix(dof, matrix_integrator, mg);

  for (unsigned int level=0;level<tr.nlevels();++level)
  {
    deallog << "mg" << std::endl;
    mg.matrix[level].print(deallog.get_file_stream());
    deallog << "in" << std::endl;
    mg.matrix_in[level].print(deallog.get_file_stream());
    deallog << "out" << std::endl;
    mg.matrix_out[level].print(deallog.get_file_stream());
    deallog << "up" << std::endl;
    mg.matrix_up[level].print(deallog.get_file_stream());
    deallog << "down" << std::endl;
    mg.matrix_down[level].print(deallog.get_file_stream());
  }
}

int main()
{
  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
  std::ofstream logfile(logname.c_str());
  deallog.attach(logfile);
  deallog.depth_console (0);
{'fraction_non_alphanumeric': 0.17505224240475808, 'fraction_numerical': 0.012940041793923806, 'mean_word_length': 4.5798206278026905, 'pattern_counts': {'":': 0, '<': 176, '<?xml version=': 0, '>': 145, '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': '2553060', 'n_tokens_mistral': 4639, 'n_tokens_neox': 4410, 'n_words': 777}