From 6c334f9bfd91eb7c02684fa4130dd19d91cbce6b Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Mon, 24 Jun 2013 14:14:08 +0000
Subject: [PATCH] step 39 test with local refinement on a single processor

git-svn-id: https://svn.dealii.org/trunk@29879 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/mpi/step-39.cc | 109 +++++++++++++++++++++++--------------------
 1 file changed, 59 insertions(+), 50 deletions(-)

diff --git a/tests/mpi/step-39.cc b/tests/mpi/step-39.cc
index baf365bb0b..a147f43539 100644
--- a/tests/mpi/step-39.cc
+++ b/tests/mpi/step-39.cc
@@ -20,7 +20,7 @@
 #include <deal.II/lac/trilinos_precondition.h>
 
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/distributed/grid_refinement.h>
 #include <deal.II/grid/filtered_iterator.h>
 
 #include <deal.II/fe/fe_q.h>
@@ -576,50 +576,59 @@ namespace Step39
   }
 
 
-  // template <int dim>
-  // double
-  // InteriorPenaltyProblem<dim>::estimate()
-  // {
-  //   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);
+  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);
 
-  //   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);
+    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);
 
-  //   NamedData<TrilinosWrappers::MPI::Vector* > solution_data;
-  //   solution_data.add(&solution, "solution");
+    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);
 
-  //   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);
+    NamedData<TrilinosWrappers::MPI::Vector* > solution_data;
+    solution_data.add(&ghost, "solution");
 
-  //   info_box.add_update_flags_boundary(update_quadrature_points);
-  //   info_box.initialize(fe, mapping, solution_data);
+    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);
 
-  //   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+    info_box.add_update_flags_boundary(update_quadrature_points);
+    info_box.initialize(fe, mapping, solution_data);
 
-  //   MeshWorker::Assembler::CellsAndFaces<double> assembler;
-  //   NamedData<BlockTrilinosWrappers::MPI::Vector* > out_data;
-  //   BlockTrilinosWrappers::MPI::Vector *est = &estimates;
-  //   out_data.add(est, "cells");
-  //   assembler.initialize(out_data, false);
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
 
-  //   Estimator<dim> integrator;
-  //   MeshWorker::integration_loop<dim, dim> (
-  //     dof_handler.begin_active(), dof_handler.end(),
-  //     dof_info, info_box,
-  //     integrator, assembler);
+    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();
-  // }
+    triangulation.load_user_indices(old_user_indices);
+    return estimates.block(0).l2_norm();
+  }
 
 
   // template <int dim>
@@ -673,20 +682,20 @@ namespace Step39
   // {
   //   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);
   // }
 
@@ -698,15 +707,15 @@ namespace Step39
     for (unsigned int s=0; s<n_steps; ++s)
       {
         deallog << "Step " << s << std::endl;
-        //if (estimates.block(0).size() == 0)
+        if (estimates.block(0).size() == 0)
           triangulation.refine_global(1);
-	  //else
-          //{
-	  // GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
-	  //                                                   estimates.block(0),
-	  //                                                   0.5, 0.0);
-	  //triangulation.execute_coarsening_and_refinement ();
-          //}
+	else
+        {
+	  parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
+										    estimates.block(0),
+										    0.5, 0.0);
+	  triangulation.execute_coarsening_and_refinement ();
+        }
 
         deallog << "Triangulation "
                 << triangulation.n_active_cells() << " cells, "
@@ -727,7 +736,7 @@ namespace Step39
         deallog << "Solve" << std::endl;
         solve();
 	//error();
-        //deallog << "Estimate " << estimate() << std::endl;
+        deallog << "Estimate " << estimate() << std::endl;
         //output_results(s);
       }
   }
-- 
2.39.5