From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 24 Jun 2013 16:08:05 +0000 (+0000)
Subject: reduce example
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e051f651277cc8c077a07e56b0f69eaa1c39e13c;p=dealii-svn.git

reduce example

git-svn-id: https://svn.dealii.org/trunk@29883 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/tests/mpi/flux_edge_01.cc b/tests/mpi/flux_edge_01.cc
index 070faebecf..5b68d0e1f3 100644
--- a/tests/mpi/flux_edge_01.cc
+++ b/tests/mpi/flux_edge_01.cc
@@ -1,13 +1,16 @@
-/* Author: Guido Kanschat, Texas A&M University, 2009 */
-
-/*    $Id: step-39.cc 28552 2013-02-25 15:18:15Z kanschat $       */
-/*                                                                */
-/*    Copyright (C) 2010-2013 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.                        */
+//----------------------------------------------------------------------
+//    $Id: mesh_worker_01.cc 28614 2013-02-27 15:20:34Z kanschat $
+//
+//    Copyright (C) 2013 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.
+//
+//----------------------------------------------------------------------
+
+// Test problem in MGTools::make_flux_sparsity_pattern_edge
 
 #include "../tests.h"
 #include <deal.II/lac/sparse_matrix.h>
@@ -55,283 +58,6 @@ namespace Step39
 {
   using namespace dealii;
 
-  Functions::SlitSingularityFunction<2> exact_solution;
-
-
-
-
-  template <int dim>
-  class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim> &dinfo,
-	      typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void boundary(MeshWorker::DoFInfo<dim> &dinfo,
-		  typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void face(MeshWorker::DoFInfo<dim> &dinfo1,
-	      MeshWorker::DoFInfo<dim> &dinfo2,
-	      typename MeshWorker::IntegrationInfo<dim> &info1,
-	      typename MeshWorker::IntegrationInfo<dim> &info2) const;
-  };
-
-
-  template <int dim>
-  void MatrixIntegrator<dim>::cell(
-    MeshWorker::DoFInfo<dim> &dinfo,
-    typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
-  }
-
-
-  template <int dim>
-  void MatrixIntegrator<dim>::boundary(
-    MeshWorker::DoFInfo<dim> &dinfo,
-    typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
-    LocalIntegrators::Laplace::nitsche_matrix(
-      dinfo.matrix(0,false).matrix, info.fe_values(0),
-      LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg));
-  }
-
-  template <int dim>
-  void MatrixIntegrator<dim>::face(
-    MeshWorker::DoFInfo<dim> &dinfo1,
-    MeshWorker::DoFInfo<dim> &dinfo2,
-    typename MeshWorker::IntegrationInfo<dim> &info1,
-    typename MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
-    LocalIntegrators::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(0), info2.fe_values(0),
-      LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
-  }
-
-  template <int dim>
-  class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void face(MeshWorker::DoFInfo<dim> &dinfo1,
-	      MeshWorker::DoFInfo<dim> &dinfo2,
-	      typename MeshWorker::IntegrationInfo<dim> &info1,
-	      typename MeshWorker::IntegrationInfo<dim> &info2) const;
-  };
-
-
-  template <int dim>
-  void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &, typename MeshWorker::IntegrationInfo<dim> &) const
-  {}
-
-
-  template <int dim>
-  void RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
-    Vector<double> &local_vector = dinfo.vector(0).block(0);
-
-    std::vector<double> boundary_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
-
-    const unsigned int deg = fe.get_fe().tensor_degree();
-    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-        local_vector(i) += (- fe.shape_value(i,k) * penalty * boundary_values[k]
-                            + (fe.normal_vector(k) * fe.shape_grad(i,k)) * boundary_values[k])
-                           * fe.JxW(k);
-  }
-
-
-  template <int dim>
-  void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
-                                MeshWorker::DoFInfo<dim> &,
-                                typename MeshWorker::IntegrationInfo<dim> &,
-                                typename MeshWorker::IntegrationInfo<dim> &) const
-  {}
-
-
-  template <int dim>
-  class Estimator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void face(MeshWorker::DoFInfo<dim> &dinfo1,
-	      MeshWorker::DoFInfo<dim> &dinfo2,
-	      typename MeshWorker::IntegrationInfo<dim> &info1,
-	      typename MeshWorker::IntegrationInfo<dim> &info2) const;
-  };
-
-
-  template <int dim>
-  void Estimator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
-
-    const std::vector<Tensor<2,dim> > &DDuh = info.hessians[0][0];
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      {
-        const double t = dinfo.cell->diameter() * trace(DDuh[k]);
-        dinfo.value(0) +=  t*t * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
-
-  template <int dim>
-  void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
-
-    std::vector<double> boundary_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
-
-    const std::vector<double> &uh = info.values[0][0];
-
-    const unsigned int deg = fe.get_fe().tensor_degree();
-    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      dinfo.value(0) += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
-                        * fe.JxW(k);
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
-
-
-  template <int dim>
-  void Estimator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
-                            MeshWorker::DoFInfo<dim> &dinfo2,
-                            typename MeshWorker::IntegrationInfo<dim> &info1,
-                            typename MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const FEValuesBase<dim> &fe = info1.fe_values();
-    const std::vector<double> &uh1 = info1.values[0][0];
-    const std::vector<double> &uh2 = info2.values[0][0];
-    const std::vector<Tensor<1,dim> > &Duh1 = info1.gradients[0][0];
-    const std::vector<Tensor<1,dim> > &Duh2 = info2.gradients[0][0];
-
-    const unsigned int deg = fe.get_fe().tensor_degree();
-    const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
-    const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
-    const double penalty = penalty1 + penalty2;
-    const double h = dinfo1.face->measure();
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      {
-        double diff1 = uh1[k] - uh2[k];
-        double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
-        dinfo1.value(0) += (penalty * diff1*diff1 + h * diff2*diff2)
-                           * fe.JxW(k);
-      }
-    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
-    dinfo2.value(0) = dinfo1.value(0);
-  }
-
-
-
-  template <int dim>
-  class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
-  {
-  public:
-    void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
-    void face(MeshWorker::DoFInfo<dim> &dinfo1,
-	      MeshWorker::DoFInfo<dim> &dinfo2,
-	      typename MeshWorker::IntegrationInfo<dim> &info1,
-	      typename MeshWorker::IntegrationInfo<dim> &info2) const;
-  };
-
-
-  template <int dim>
-  void ErrorIntegrator<dim>::cell(
-    MeshWorker::DoFInfo<dim> &dinfo,
-    typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
-    std::vector<Tensor<1,dim> > exact_gradients(fe.n_quadrature_points);
-    std::vector<double> exact_values(fe.n_quadrature_points);
-
-    exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients);
-    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
-
-    const std::vector<Tensor<1,dim> > &Duh = info.gradients[0][0];
-    const std::vector<double> &uh = info.values[0][0];
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      {
-        double sum = 0;
-        for (unsigned int d=0; d<dim; ++d)
-          {
-            const double diff = exact_gradients[k][d] - Duh[k][d];
-            sum += diff*diff;
-          }
-        const double diff = exact_values[k] - uh[k];
-        dinfo.value(0) +=  sum * fe.JxW(k);
-        dinfo.value(1) +=  diff*diff * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-    dinfo.value(1) = std::sqrt(dinfo.value(1));
-  }
-
-
-  template <int dim>
-  void ErrorIntegrator<dim>::boundary(
-    MeshWorker::DoFInfo<dim> &dinfo,
-    typename MeshWorker::IntegrationInfo<dim> &info) const
-  {
-    const FEValuesBase<dim> &fe = info.fe_values();
-
-    std::vector<double> exact_values(fe.n_quadrature_points);
-    exact_solution.value_list(fe.get_quadrature_points(), exact_values);
-
-    const std::vector<double> &uh = info.values[0][0];
-
-    const unsigned int deg = fe.get_fe().tensor_degree();
-    const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      {
-        const double diff = exact_values[k] - uh[k];
-        dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
-      }
-    dinfo.value(0) = std::sqrt(dinfo.value(0));
-  }
-
-
-  template <int dim>
-  void ErrorIntegrator<dim>::face(
-    MeshWorker::DoFInfo<dim> &dinfo1,
-    MeshWorker::DoFInfo<dim> &dinfo2,
-    typename MeshWorker::IntegrationInfo<dim> &info1,
-    typename MeshWorker::IntegrationInfo<dim> &info2) const
-  {
-    const FEValuesBase<dim> &fe = info1.fe_values();
-    const std::vector<double> &uh1 = info1.values[0][0];
-    const std::vector<double> &uh2 = info2.values[0][0];
-
-    const unsigned int deg = fe.get_fe().tensor_degree();
-    const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
-    const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
-    const double penalty = penalty1 + penalty2;
-
-    for (unsigned k=0; k<fe.n_quadrature_points; ++k)
-      {
-        double diff = uh1[k] - uh2[k];
-        dinfo1.value(0) += (penalty * diff*diff)
-                           * fe.JxW(k);
-      }
-    dinfo1.value(0) = std::sqrt(dinfo1.value(0));
-    dinfo2.value(0) = dinfo1.value(0);
-  }
-
-
-
-
   template <int dim>
   class InteriorPenaltyProblem
   {
@@ -344,13 +70,6 @@ namespace Step39
 
   private:
     void setup_system ();
-    void assemble_matrix ();
-    void assemble_mg_matrix ();
-    void assemble_right_hand_side ();
-    void error ();
-    double estimate ();
-    void solve ();
-    void output_results (const unsigned int cycle) const;
 
     parallel::distributed::Triangulation<dim>        triangulation;
     const MappingQ1<dim>      mapping;
@@ -360,9 +79,6 @@ namespace Step39
     IndexSet locally_relevant_set;
     
     TrilinosWrappers::SparseMatrix matrix;
-    TrilinosWrappers::MPI::Vector       solution;
-    TrilinosWrappers::MPI::Vector       right_hand_side;
-    BlockVector<double>  estimates;
 
     MGLevelObject<TrilinosWrappers::SparseMatrix> mg_matrix;
 
@@ -394,9 +110,7 @@ namespace Step39
     dof_handler.distribute_mg_dofs (fe);
     
     DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_set);
-    solution.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-    right_hand_side.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
-
+    
     CompressedSimpleSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());
     DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
     matrix.reinit(dof_handler.locally_owned_dofs(), c_sparsity, MPI_COMM_WORLD, true);
@@ -438,268 +152,6 @@ namespace Step39
   }
 
 
-  template <int dim>
-  void
-  InteriorPenaltyProblem<dim>::assemble_matrix()
-  {
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    UpdateFlags update_flags = update_values | update_gradients;
-    info_box.add_update_flags_all(update_flags);
-    info_box.initialize(fe, mapping);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::MatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
-    assembler.initialize(matrix);
-    
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-    
-    MatrixIntegrator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim>(
-      begin, end,
-      dof_info, info_box,
-      integrator, assembler);
-    
-    matrix.compress(VectorOperation::add);
-  }
-
-
-  template <int dim>
-  void
-  InteriorPenaltyProblem<dim>::assemble_mg_matrix()
-  {
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    UpdateFlags update_flags = update_values | update_gradients;
-    info_box.add_update_flags_all(update_flags);
-    info_box.initialize(fe, mapping);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::MGMatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
-    assembler.initialize(mg_matrix);
-    assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
-    
-    FilteredIterator<typename DoFHandler<dim>::level_cell_iterator>
-      begin(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.begin());
-    FilteredIterator<typename DoFHandler<dim>::level_cell_iterator>
-      end(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.end());
-
-    MatrixIntegrator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim> (
-      begin, end,
-      dof_info, info_box,
-      integrator, assembler);
-    
-    for (unsigned int level=mg_matrix.min_level(); level <= mg_matrix.max_level();++level)
-      {
-	mg_matrix[level].compress(VectorOperation::add);
-	if (level > mg_matrix.min_level())
-	  {
-	    mg_matrix_dg_up[level].compress(VectorOperation::add);
-	    mg_matrix_dg_down[level].compress(VectorOperation::add);
-	  }
-      }
-  }
-
-
-  template <int dim>
-  void
-  InteriorPenaltyProblem<dim>::assemble_right_hand_side()
-  {
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
-    info_box.add_update_flags_all(update_flags);
-    info_box.initialize(fe, mapping);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::ResidualSimple<TrilinosWrappers::MPI::Vector> assembler;
-    NamedData<TrilinosWrappers::MPI::Vector*> data;
-    TrilinosWrappers::MPI::Vector *rhs = &right_hand_side;
-    data.add(rhs, "RHS");
-    assembler.initialize(data);
-
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-    
-
-    RHSIntegrator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim>(begin, end,
-      dof_info, info_box,
-      integrator, assembler);
-
-    right_hand_side.compress(VectorOperation::add);
-    right_hand_side *= -1.;
-  }
-
-
-  template <int dim>
-  void
-  InteriorPenaltyProblem<dim>::solve()
-  {
-    SolverControl control(1000, 1.e-12);
-    SolverCG<TrilinosWrappers::MPI::Vector > solver(control);
-
-    MGTransferPrebuilt<TrilinosWrappers::MPI::Vector > mg_transfer;
-    mg_transfer.build_matrices(dof_handler);
-
-    SolverControl coarse_solver_control (1000, 1e-10, false, false);
-    SolverCG<TrilinosWrappers::MPI::Vector> coarse_solver(coarse_solver_control);
-    PreconditionIdentity identity;
-    TrilinosWrappers::SparseMatrix coarse_matrix = mg_matrix[0];
-    MGCoarseGridLACIteration<SolverCG<TrilinosWrappers::MPI::Vector>,TrilinosWrappers::MPI::Vector>
-      coarse_grid_solver(coarse_solver, coarse_matrix, identity);
-
-    typedef TrilinosWrappers::PreconditionJacobi Smoother;
-    MGSmootherPrecondition<TrilinosWrappers::SparseMatrix, Smoother, TrilinosWrappers::MPI::Vector> mg_smoother;
-    mg_smoother.initialize(mg_matrix, .7);
-    mg_smoother.set_steps(2);
-
-    MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgmatrix(&mg_matrix);
-    MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgdown(&mg_matrix_dg_down);
-    MGMatrix<TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector> mgup(&mg_matrix_dg_up);
-
-    Multigrid<TrilinosWrappers::MPI::Vector > mg(dof_handler, mgmatrix,
-                                  coarse_grid_solver, mg_transfer,
-                                  mg_smoother, mg_smoother);
-    mg.set_edge_flux_matrices(mgdown, mgup);
-
-    PreconditionMG<dim, TrilinosWrappers::MPI::Vector,
-                   MGTransferPrebuilt<TrilinosWrappers::MPI::Vector > >
-                   preconditioner(dof_handler, mg, mg_transfer);
-    solver.solve(matrix, solution, right_hand_side, preconditioner);
-  }
-
-
-  template <int dim>
-  double
-  InteriorPenaltyProblem<dim>::estimate()
-  {
-    TrilinosWrappers::MPI::Vector ghost;
-    ghost.reinit(locally_relevant_set, MPI_COMM_WORLD);
-    ghost = solution;
-    
-    std::vector<unsigned int> old_user_indices;
-    triangulation.save_user_indices(old_user_indices);
-
-    estimates.block(0).reinit(triangulation.n_active_cells());
-    
-    unsigned int i=0;
-    for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
-         cell != triangulation.end(); ++cell,++i)
-      cell->set_user_index(i);
-
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
-    info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
-
-    NamedData<TrilinosWrappers::MPI::Vector* > solution_data;
-    solution_data.add(&ghost, "solution");
-
-    info_box.cell_selector.add("solution", false, false, true);
-    info_box.boundary_selector.add("solution", true, true, false);
-    info_box.face_selector.add("solution", true, true, false);
-
-    info_box.add_update_flags_boundary(update_quadrature_points);
-    info_box.initialize(fe, mapping, solution_data);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::CellsAndFaces<double> assembler;
-    NamedData<BlockVector<double>* > out_data;
-    BlockVector<double> *est = &estimates;
-    out_data.add(est, "cells");
-    assembler.initialize(out_data, false);
-    
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-    
-    Estimator<dim> integrator;
-    MeshWorker::integration_loop<dim, dim> (
-					    begin, end,
-					    dof_info, info_box,
-					    integrator, assembler);
-
-    triangulation.load_user_indices(old_user_indices);
-    return estimates.block(0).l2_norm();
-  }
-
-
-  // template <int dim>
-  // void
-  // InteriorPenaltyProblem<dim>::error()
-  // {
-  //   BlockVector<double> errors(2);
-  //   errors.block(0).reinit(triangulation.n_active_cells());
-  //   errors.block(1).reinit(triangulation.n_active_cells());
-  //   unsigned int i=0;
-  //   for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
-  //        cell != triangulation.end(); ++cell,++i)
-  //     cell->set_user_index(i);
-
-  //   MeshWorker::IntegrationInfoBox<dim> info_box;
-  //   const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
-  //   info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
-
-  //   NamedData<TrilinosWrappers::MPI::Vector* > solution_data;
-  //   solution_data.add(&solution, "solution");
-
-  //   info_box.cell_selector.add("solution", true, true, false);
-  //   info_box.boundary_selector.add("solution", true, false, false);
-  //   info_box.face_selector.add("solution", true, false, false);
-
-  //   info_box.add_update_flags_cell(update_quadrature_points);
-  //   info_box.add_update_flags_boundary(update_quadrature_points);
-  //   info_box.initialize(fe, mapping, solution_data);
-
-  //   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-  //   MeshWorker::Assembler::CellsAndFaces<double> assembler;
-  //   NamedData<BlockVector<double>* > out_data;
-  //   BlockVector<double> *est = &errors;
-  //   out_data.add(est, "cells");
-  //   assembler.initialize(out_data, false);
-
-  //   ErrorIntegrator<dim> integrator;
-  //   MeshWorker::integration_loop<dim, dim> (
-  //     dof_handler.begin_active(), dof_handler.end(),
-  //     dof_info, info_box,
-  //     integrator, assembler);
-
-  //   deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl;
-  //   deallog << "L2-error:     " << errors.block(1).l2_norm() << std::endl;
-  // }
-
-
-  // template <int dim>
-  // void InteriorPenaltyProblem<dim>::output_results (const unsigned int cycle) const
-  // {
-  //   char *fn = new char[100];
-  //   sprintf(fn, "step-39/sol-%02d", cycle);
-    
-  //   std::string filename(fn);
-  //   filename += ".gnuplot";
-  //   deallog << "Writing solution to <" << filename << ">..."
-  //           << std::endl << std::endl;
-  //   std::ofstream gnuplot_output (filename.c_str());
-    
-  //   DataOut<dim> data_out;
-  //   data_out.attach_dof_handler (dof_handler);
-  //   data_out.add_data_vector (solution, "u");
-  //   data_out.add_data_vector (estimates.block(0), "est");
-    
-  //   data_out.build_patches ();
-    
-  //   data_out.write_gnuplot(gnuplot_output);
-  // }
-
   template <int dim>
   void
   InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
@@ -708,16 +160,13 @@ namespace Step39
     for (unsigned int s=0; s<n_steps; ++s)
       {
         deallog << "Step " << s << std::endl;
-        if (estimates.block(0).size() == 0)
+        if (s==0)
           triangulation.refine_global(1);
 	else
         {
 	  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 	    triangulation.begin_active(triangulation.n_levels()-1)->set_refine_flag();
 	      
-	  // parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
-	  // 									    estimates.block(0),
-	  // 									    0.5, 0.0);
 	  triangulation.execute_coarsening_and_refinement ();
         }
 
@@ -730,18 +179,6 @@ namespace Step39
         for (unsigned int l=0; l<triangulation.n_levels(); ++l)
           deallog << ' ' << dof_handler.n_dofs(l);
         deallog << std::endl;
-
-        deallog << "Assemble matrix" << std::endl;
-        assemble_matrix();
-        deallog << "Assemble multilevel matrix" << std::endl;
-        assemble_mg_matrix();
-        deallog << "Assemble right hand side" << std::endl;
-        assemble_right_hand_side();
-        deallog << "Solve" << std::endl;
-        solve();
-	//error();
-        deallog << "Estimate " << estimate() << std::endl;
-        //output_results(s);
       }
   }
 }
@@ -758,7 +195,7 @@ int main(int argc, char *argv[])
   
   try
     {
-      FE_DGQ<2> fe1(2);
+      FE_DGQ<2> fe1(0);
       InteriorPenaltyProblem<2> test1(fe1);
       test1.run(20);
     }