]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
step 39 running
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 2 Dec 2009 04:12:15 +0000 (04:12 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 2 Dec 2009 04:12:15 +0000 (04:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@20195 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/include/numerics/mesh_worker_info.h
deal.II/deal.II/include/numerics/mesh_worker_info.templates.h
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.templates.h
deal.II/deal.II/source/numerics/mesh_worker_assembler.all_dimensions.cc
deal.II/examples/step-39/step-39.cc

index 33af1b0de02829ca27dbf684dbe2a3e6fd9010f0..b0077fd09ab7785e2c05d860910796fc37be6642 100644 (file)
@@ -160,7 +160,7 @@ namespace MeshWorker
                                          * communicating the cell and
                                          * face vectors.
                                          */
-       typedef NamedData<SmartPointer<BlockVector<number>,CellsAndFaces<number> > > DataVectors;
+       typedef NamedData<BlockVector<number>*> DataVectors;
        
                                         /**
                                          * The initialization
@@ -219,7 +219,7 @@ namespace MeshWorker
                                          */
        template <int dim>
        void initialize_info(DoFInfo<dim>& info,
-                            bool interior_face);
+                            bool interior_face) const;
        
                                         /**
                                          * Assemble the local values
@@ -962,6 +962,116 @@ namespace MeshWorker
     };
 
     
+//----------------------------------------------------------------------//
+    
+    template <typename number>
+    inline void
+    Functional<number>::initialize(unsigned int n)
+    {
+      results.resize(n);
+    }
+    
+    
+    template <typename number>
+    template<int dim>
+    inline void
+    Functional<number>::assemble(const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       results[i] += info.J[i];
+    }
+    
+
+    template <typename number>
+    template<int dim>
+    inline void
+    Functional<number>::assemble(const DoFInfo<dim>& info1,
+                                const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<results.size();++i)
+       {
+         results[i] += info1.J[i];
+         results[i] += info2.J[i];
+       }
+    }
+
+    
+    template <typename number>
+    inline number
+    Functional<number>::operator() (unsigned int i) const
+    {
+      AssertIndexRange(i, results.size());
+      return results[i];
+    }
+    
+//----------------------------------------------------------------------//
+
+    template <typename number>
+    inline void
+    CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
+    {
+      Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
+      if (sep)
+       {
+         Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
+         AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
+       }
+      
+      results = r;
+      separate_faces = sep;
+    }
+
+    
+    template <typename number>
+    template <int dim>
+    inline void
+    CellsAndFaces<number>::initialize_info(DoFInfo<dim>& info, bool) const
+    {
+      info.initialize_numbers(results(0)->n_blocks());
+    }
+
+
+    template <typename number>
+    template<int dim>
+    inline void
+    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
+    {
+      for (unsigned int i=0;i<info.J.size();++i)
+       {
+         if (separate_faces &&
+             info.face_number != deal_II_numbers::invalid_unsigned_int)
+           results(1)->block(i)(info.face->user_index()) += info.J[i];
+         else
+           results(0)->block(i)(info.cell->user_index()) += info.J[i];
+       }
+    }
+    
+
+    template <typename number>
+    template<int dim>
+    inline void
+    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
+                                   const DoFInfo<dim>& info2)
+    {
+      for (unsigned int i=0;i<info1.J.size();++i)
+       {
+         if (separate_faces)
+           {
+             const double J = info1.J[i] + info2.J[i];
+             results(1)->block(i)(info1.face->user_index()) += J;
+             if (info2.face != info1.face)
+               results(1)->block(i)(info2.face->user_index()) += J;
+           }
+         else
+           {
+             results(0)->block(i)(info1.cell->user_index()) += .5*info1.J[i];
+             results(0)->block(i)(info2.cell->user_index()) += .5*info2.J[i];
+           }
+       }
+    }
+    
+
+
 //----------------------------------------------------------------------//
 
     template <class VECTOR>
index 947a66b23a7c7aa53962a3308735a65ca91a5c4e..cef55748d31e4117085e7fcde0b2d64f7b3efa05 100644 (file)
@@ -47,7 +47,8 @@ namespace MeshWorker
        MatrixBlock<FullMatrix<number> >& M,
        unsigned int row, unsigned int col);
     public:
-      void initialize_vectors(const unsigned int n_vectors);
+      void initialize_numbers(const unsigned int n);
+      void initialize_vectors(const unsigned int n);
                                       /**
                                        * Allocate @p n local
                                        * matrices. Additionally,
@@ -542,19 +543,20 @@ namespace MeshWorker
                                        */
       template <typename T>
       IntegrationInfoBox(const T&);
-       
+
+      template <class VECTOR>
+      void initialize_data(const VECTOR*, const std::string& name,
+                          bool values, bool gradients, bool hessians);
+      
       template <class WORKER>
       void initialize(const WORKER&,
                      const FiniteElement<dim, spacedim>& el,
                      const Mapping<dim, spacedim>& mapping);
        
-      template <class WORKER, class VECTOR,class P>
-      void initialize(const WORKER&,
-                     const FiniteElement<dim, spacedim>& el,
-                     const Mapping<dim, spacedim>& mapping,
-                     const NamedData<SmartPointer<VECTOR,P> >& data);
 //      private:
-       
+
+      boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > data;
+      
       CellInfo cell_info;
       FaceInfo bdry_info;
       FaceInfo face_info;
@@ -565,6 +567,14 @@ namespace MeshWorker
 
 //----------------------------------------------------------------------//
 
+  template <typename number>
+  inline void
+  LocalResults<number>::initialize_numbers(unsigned int n)
+  {
+    J.resize(n);
+  }
+
+    
   template <typename number>
   inline void
   LocalResults<number>::initialize_vectors(unsigned int n)
@@ -623,6 +633,8 @@ namespace MeshWorker
   inline void
   LocalResults<number>::reinit(const BlockIndices& bi)
   {
+    for (unsigned int i=0;i<J.size();++i)
+      J[i] = 0.;
     for (unsigned int i=0;i<R.size();++i)
       R[i].reinit(bi);
     for (unsigned int i=0;i<M1.size();++i)
@@ -818,7 +830,29 @@ namespace MeshWorker
                  subface_info(t),
                  neighbor_info(t)
   {}
+  
+  
+  template <int dim, int sdim>
+  template <typename VECTOR>
+  void
+  IntegrationInfoBox<dim,sdim>::initialize_data(
+    const VECTOR* v, const std::string& name,
+    bool values, bool gradients, bool hessians)
+  {
+    boost::shared_ptr<VectorData<VECTOR, dim, sdim> >
+      p = boost::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> ());
+    p->add(name, values, gradients, hessians);
+    p->initialize(v, name);
+    data = p;
     
+    cell_info.initialize_data(data);
+    bdry_info.initialize_data(data);
+    face_info.initialize_data(data);
+    subface_info.initialize_data(data);
+    neighbor_info.initialize_data(data);
+  }
+
+  
   template <int dim, int sdim>
   template <class WORKER>
   void
@@ -836,33 +870,14 @@ namespace MeshWorker
     cell_info.initialize<FEValues<dim,sdim> >(el, mapping, integrator.cell_quadrature,
                         integrator.cell_flags);
     bdry_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.bdry_quadrature,
-                        integrator.face_flags);
+                        integrator.bdry_flags);
     face_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
                         integrator.face_flags);
     subface_info.initialize<FESubfaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
                            integrator.face_flags);
     neighbor_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
                             integrator.ngbr_flags);
-  }
-
-    
-  template <int dim, int sdim>
-  template <class WORKER, class VECTOR, class P>
-  void
-  IntegrationInfoBox<dim,sdim>::initialize(
-    const WORKER& integrator,
-    const FiniteElement<dim,sdim>& el,
-    const Mapping<dim,sdim>& mapping,
-    const NamedData<SmartPointer<VECTOR,P> >& data)
-  {
-    cell_info.initialize_data(data);
-    bdry_info.initialize_data(data);
-    face_info.initialize_data(data);
-    subface_info.initialize_data(data);
-    neighbor_info.initialize_data(data);
-
-    initialize(integrator, el, mapping);
-  }
+  }    
 }
 
 
index c92dbc15e52ccb3eb8a6109ab5006ca687820679..95c5274699c24cc090934277026a460c42d7122f 100644 (file)
@@ -106,7 +106,7 @@ namespace MeshWorker
                                         // For all components
        for (unsigned int j=0;j<values[i].size();++j)
          {
-           values[i][j].resize(el.n_components());
+           values[i][j].resize(quadrature.size());
          }
       }
     
@@ -119,7 +119,7 @@ namespace MeshWorker
                                         // For all components
        for (unsigned int j=0;j<gradients[i].size();++j)
          {
-           gradients[i][j].resize(el.n_components());
+           gradients[i][j].resize(quadrature.size());
          }
       }
     
@@ -132,7 +132,7 @@ namespace MeshWorker
                                         // For all components
        for (unsigned int j=0;j<hessians[i].size();++j)
          {
-           hessians[i][j].resize(el.n_components());
+           hessians[i][j].resize(quadrature.size());
          }
       }
   }
index 4b5654ee45390b957175d4c76473058ead35e2fa..cab71251f3d83558a077e55d58c060dd50cc6ebe 100644 (file)
@@ -122,6 +122,12 @@ namespace MeshWorker
       template <class STREAM, typename DATA>
       void print (STREAM& s, const NamedData<DATA>& v) const;
       
+                                      /**
+                                       * Print the number of selections to the stream.
+                                       */
+      template <class STREAM>
+      void print (STREAM& s) const;
+      
     protected:
                                       /**
                                        * Selection of the vectors
@@ -238,9 +244,11 @@ namespace MeshWorker
        unsigned int size) const;
   };
 
+  
 /**
  * Based on VectorSelector, this is the class that implements the
- * function VectorDataBase::fill() for a certain type of vector.
+ * function VectorDataBase::fill() for a certain type of vector, using
+ * NamedData to identify vectors by name.
  *
  * @author Guido Kanschat, 2009
  */
@@ -250,6 +258,7 @@ namespace MeshWorker
   {
     public:
       void initialize(const NamedData<VECTOR*>&);
+      void initialize(const VECTOR*, const std::string& name);
       
       virtual void fill(
        std::vector<std::vector<std::vector<double> > >& values,
@@ -262,8 +271,9 @@ namespace MeshWorker
        unsigned int start,
        unsigned int size) const;
     private:
-      NamedData<SmartPointer<VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
+      NamedData<SmartPointer<const VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
   };
+
   
 //----------------------------------------------------------------------//
 
@@ -357,6 +367,17 @@ namespace MeshWorker
   }
   
 
+  template <class STREAM>
+  inline void
+  VectorSelector::print(STREAM& s) const
+  {
+    s << "values: " << n_values()
+      << " gradients: " << n_gradients()
+      << " hessians: " << n_hessians()
+      << std::endl;
+  }
+
+
   template <class STREAM, typename DATA>
   inline void
   VectorSelector::print(STREAM& s, const NamedData<DATA>& v) const
index d0d3e6a6df610af9ca7f1d86c3d8b336867fe825..e691467949996e712d72000756944903fe45f4a6 100644 (file)
@@ -47,6 +47,16 @@ namespace MeshWorker
   }
   
   
+  template <class VECTOR, int dim, int spacedim>
+  void
+  VectorData<VECTOR, dim, spacedim>::initialize(const VECTOR* v, const std::string& name)
+  {
+    SmartPointer<const VECTOR,VectorData<VECTOR, dim, spacedim> > p = v;
+    data.add(p, name);
+    VectorSelector::initialize(data);
+  }
+  
+  
   template <class VECTOR, int dim, int spacedim>
   void
   VectorData<VECTOR, dim, spacedim>::fill(
index c57593172c86480211e2be8beb6100f2032382c4..b6c27e8d8de5fc78f66d47096ad98ca4485ca460 100644 (file)
@@ -22,107 +22,9 @@ namespace MeshWorker
   namespace Assembler
   {
 //----------------------------------------------------------------------//
-    
-    template <typename number>
-    inline void
-    Functional<number>::initialize(unsigned int n)
-    {
-      results.resize(n);
-    }
-    
-    
-    template <typename number>
-    template<int dim>
-    inline void
-    Functional<number>::assemble(const DoFInfo<dim>& info)
-    {
-      for (unsigned int i=0;i<results.size();++i)
-       results[i] += info.J[i];
-    }
-    
 
-    template <typename number>
-    template<int dim>
-    inline void
-    Functional<number>::assemble(const DoFInfo<dim>& info1,
-                                const DoFInfo<dim>& info2)
-    {
-      for (unsigned int i=0;i<results.size();++i)
-       {
-         results[i] += info1.J[i];
-         results[i] += info2.J[i];
-       }
-    }
-
-    
-    template <typename number>
-    inline number
-    Functional<number>::operator() (unsigned int i) const
-    {
-      AssertIndexRange(i, results.size());
-      return results[i];
-    }
-    
 //----------------------------------------------------------------------//
 
-    template <typename number>
-    inline void
-    CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
-    {
-      Assert(r.name(0) == "cells", DataVectors::ExcNameMismatch(0, "cells"));
-      if (sep)
-       Assert(r.name(1) == "faces", DataVectors::ExcNameMismatch(1, "faces"));
-      AssertDimension(r(0).n_blocks(), r(1).n_blocks());
-      
-      results = r;
-      separate_faces = sep;
-    }
-    
-    
-    template <typename number>
-    template<int dim>
-    inline void
-    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
-    {
-      for (unsigned int i=0;i<info.J.size();++i)
-       {
-         if (separate_faces &&
-             info.face_number != deal_II_numbers::invalid_unsigned_int)
-           results.vector(1).block(i)(info.face->user_index()) += info.J[i];
-         else
-           results.vector(0).block(i)(info.cell->user_index()) += info.J[i];
-       }
-    }
-    
-
-    template <typename number>
-    template<int dim>
-    inline void
-    CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
-                                   const DoFInfo<dim>& info2)
-    {
-      for (unsigned int i=0;i<info1.J.size();++i)
-       {
-         if (separate_faces)
-           {
-             const double J = info1.J[i] + info2.J[i];
-             results.vector(1).block(i)(info1.face->user_index()) += J;
-             if (info2.face != info1.face)
-               results.vector(1).block(i)(info2.face->user_index()) += J;
-           }
-         else
-           {
-             results.vector(0).block(i)(info1.cell->user_index()) += info1.J[i];
-             results.vector(0).block(i)(info2.cell->user_index()) += info2.J[i];
-           }
-       }
-    }
-    
-
-//----------------------------------------------------------------------//
-
-
-
   }
 }
 
index f17955777cf517d7a9bb4c5782246968924823a7..2961469b4627c9c06f03833f6a9b18344923d67b 100644 (file)
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 #include <lac/precondition_block.h>
+#include <lac/block_vector.h>
 
                                 // Include files for setting up the
                                 // mesh
 #include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
 
                                 // Include files for FiniteElement
                                 // classes and DoFHandler.
@@ -102,6 +104,13 @@ class MatrixIntegrator : public Subscriptor
 };
 
 
+                                // On each cell, we integrate the
+                                // Dirichlet form. All local
+                                // integrations consist of nested
+                                // loops, first over all quadrature
+                                // points and the iner loops over the
+                                // degrees of freedom associated
+                                // with the shape functions.
 template <int dim>
 void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const
 {
@@ -115,7 +124,8 @@ void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::Ce
                             * fe.JxW(k);
 }
 
-
+                                // On boundary faces, we use the
+                                // Nitsche boundary condition
 template <int dim>
 void MatrixIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const
 {
@@ -153,7 +163,7 @@ void MatrixIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::Fa
     {
       Assert (info1.face == info2.face, ExcInternalError());
       Assert (info1.face->has_children(), ExcInternalError());
-      Assert (info1.cell->has_children(), ExcInternalError());
+//      Assert (info1.cell->has_children(), ExcInternalError());
       penalty1 *= 2;
     }
 const double penalty = penalty1 + penalty2;
@@ -224,6 +234,77 @@ void RHSIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceI
 {}
 
 
+template <int dim>
+class Estimator : public Subscriptor
+{
+  public:
+    void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const;
+    void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const;
+    void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+             typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2) const;
+};
+
+
+template <int dim>
+void Estimator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const
+{
+  const FEValuesBase<dim>& fe = info.fe();
+  
+  const std::vector<Tensor<2,dim> >& DDuh = info.hessians[0][0];
+  for (unsigned k=0;k<fe.n_quadrature_points;++k)
+    {
+      const double t = info.cell->diameter() * trace(DDuh[k]);
+      info.J[0] +=  t*t * fe.JxW(k);
+    }
+}
+
+
+template <int dim>
+void Estimator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const
+{
+  const FEFaceValuesBase<dim>& fe = info.fe();
+  
+  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) * info.face->measure() / info.cell->measure();
+  
+  for (unsigned k=0;k<fe.n_quadrature_points;++k)
+    info.J[0] += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
+                * fe.JxW(k);
+}
+
+
+template <int dim>
+void Estimator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+                         typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2) const
+{
+  const FEFaceValuesBase<dim>& fe = info1.fe();
+  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) * info1.face->measure() / info1.cell->measure();
+  const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
+  const double penalty = penalty1 + penalty2;
+  const double h = info1.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];
+      info1.J[0] += (penalty * diff1*diff1 + h * diff2*diff2)
+                   * fe.JxW(k);
+    }
+  info2.J[0] = info1.J[0];
+}
+
+
                                 // @sect3{The main class}
 template <int dim>
 class Step39
@@ -239,7 +320,7 @@ class Step39
     void assemble_mg_matrix ();
     void assemble_right_hand_side ();
     void error ();
-    void estimate ();
+    double estimate ();
     void solve ();
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
@@ -254,7 +335,8 @@ class Step39
     SparseMatrix<double> matrix;
     Vector<double>       solution;
     Vector<double>       right_hand_side;
-
+    BlockVector<double>  estimates;
+    
     MGLevelObject<SparsityPattern> mg_sparsity;
     MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
     MGLevelObject<SparseMatrix<double> > mg_matrix;
@@ -268,11 +350,10 @@ Step39<dim>::Step39(const FiniteElement<dim>& fe)
                :
                fe(fe),
                mg_dof_handler(triangulation),
-               dof_handler(mg_dof_handler)
+               dof_handler(mg_dof_handler),
+               estimates(1)
 {
   GridGenerator::hyper_L(triangulation, -1, 1);
-  triangulation.begin()->set_refine_flag();
-  triangulation.execute_coarsening_and_refinement();
 }
 
 
@@ -474,7 +555,39 @@ Step39<dim>::error()
   QGauss<dim> quadrature(n_gauss_points);
   VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
                                    cell_errors, quadrature, VectorTools::L2_norm);
-  deallog << "Error " << cell_errors.l2_norm() << std::endl;
+  deallog << "Error    " << cell_errors.l2_norm() << std::endl;
+}
+
+
+template <int dim>
+double
+Step39<dim>::estimate()
+{
+  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);
+  
+  const Estimator<dim> local;
+  MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::CellsAndFaces<double>, Estimator<dim> >
+  integrator(local);
+  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+  integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
+  UpdateFlags update_flags = update_values | update_gradients;
+  integrator.add_update_flags(update_flags | update_hessians, true, true, true, true);
+//  integrator.add_update_flags(update_hessians, true, false, false, false);
+  integrator.add_update_flags(update_quadrature_points, false, true, false, false);
+
+  NamedData<BlockVector<double>* > out_data;
+  BlockVector<double>* est = &estimates;
+  out_data.add(est, "cells");
+  integrator.initialize(out_data, false);
+  MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
+  info_box.initialize_data(&solution, std::string("solution"), true, true, true);
+  info_box.initialize(integrator, fe, mapping);
+  MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator);
+  return estimates.block(0).l2_norm();
 }
 
 
@@ -483,10 +596,10 @@ void Step39<dim>::output_results (const unsigned int cycle) const
 {
                                   // Output of the solution in
                                   // gnuplot format.
-  std::string filename = "sol-";
-  filename += ('0' + cycle);
-  Assert (cycle < 10, ExcInternalError());
-  
+  char * fn = new char[100];
+  sprintf(fn, "sol-%02d", cycle);
+
+  std::string filename(fn);
   filename += ".gnuplot";
   std::cout << "Writing solution to <" << filename << ">..."
            << std::endl << std::endl;
@@ -495,6 +608,7 @@ void Step39<dim>::output_results (const unsigned int cycle) const
   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 ();
   
@@ -511,7 +625,17 @@ Step39<dim>::run(unsigned int n_steps)
     {
       deallog << "Step " << s << std::endl;
       if (s != 0)
-       triangulation.refine_global(1);
+       {
+         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 ();
+           }
+       }
       
       deallog << "Triangulation "
              << triangulation.n_active_cells() << " cells, "
@@ -531,8 +655,8 @@ Step39<dim>::run(unsigned int n_steps)
       assemble_right_hand_side();
       deallog << "Solve" << std::endl;
       solve();
-      deallog << "Error" << std::endl;
       error();
+      deallog << "Estimate " << estimate() << std::endl;
       output_results(s);
     }
 }
@@ -540,7 +664,7 @@ Step39<dim>::run(unsigned int n_steps)
 
 int main()
 {
-  FE_DGQ<2> fe1(1);
+  FE_DGQ<2> fe1(3);
   Step39<2> test1(fe1);
-  test1.run(7);
+  test1.run(13);
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.