add grid refinement with optimization

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

diff --git a/deal.II/deal.II/include/grid/grid_refinement.h b/deal.II/deal.II/include/grid/grid_refinement.h
index b93e3bfb7e1d4745178a2ff3cd692fe567d408b5..1815633c5aca4c6020d71ac4cf0ffa187befed95 100644 (file)
--- a/deal.II/deal.II/include/grid/grid_refinement.h
+++ b/deal.II/deal.II/include/grid/grid_refinement.h
@@ -14,6 +14,9 @@
 #define __deal2__grid_refinement_h
 
 
+#include<lac/vector.h>
+#include<vector>
+
 // forward declarations
 template <int dim> class Triangulation;
 template <class T> class Vector;
@@ -214,7 +217,7 @@ class GridRefinement
                                                 const Vector<number> &criteria,
                                                 const double         top_fraction_of_cells,
                                                 const double         bottom_fraction_of_cells);
-
+    
                                     /**
                                      * Refine the triangulation by
                                      * flagging those cells which
@@ -248,6 +251,41 @@ class GridRefinement
                                                   const double         top_fraction,
                                                   const double         bottom_fraction);
 
+
+
+                                    /**
+                                     * Refine the triangulation by
+                                     * flagging certain cells to reach
+                                     * an optimal grid:
+                                     * We try to minimize the error
+                                     * multiplied with the number of
+                                     * cells in the new grid. All cells
+                                     * with large error indicator are
+                                     * refined to generate an optimal
+                                     * grid in the above sense.
+                                     * We assume that the error in one
+                                     * cell is reduced to a quarter
+                                     * after refinement.
+                                     * The new triangulation has three
+                                     * new cells for every flagged cell.
+                                     *
+                                     * Refer to the general doc of
+                                     * this class for more
+                                     * information.
+                                     *
+                                     * Note that this function takes
+                                     * a vector of @p{float}s, rather
+                                     * than the usual @p{double}s,
+                                     * since accuracy is not so much
+                                     * needed here and saving memory
+                                     * may be a good goal when using
+                                     * many cells.
+                                     */
+    
+    template<int dim, typename number>
+    static void refine_and_coarsen_optimize (Triangulation<dim>   &tria,
+                                            const Vector<number> &criteria);
+
                                     /**
                                      * Exception
                                      */
@@ -259,8 +297,28 @@ class GridRefinement
                                      * Exception
                                      */
     DeclException0 (ExcInvalidParameterValue);
+
+
+  private:
+    
+                                    /**
+                                     * Sorts the vector ind as
+                                     * an index vector of a in
+                                     * increasing order.
+                                     * This implementation of
+                                     * quicksort seems to be faster
+                                     * than the STL version and is
+                                     * needed in
+                                     * @p{refine_and_coarsen_optimize}
+                                     */
+    template<typename number>
+    static void qsort_index(const Vector<number>  &a,
+                           vector<unsigned int>  &ind,
+                           int                   l,
+                           int                   r);
 };
 
 
 
 #endif //__deal2__grid_refinement_h
+

diff --git a/deal.II/deal.II/source/grid/grid_refinement.cc b/deal.II/deal.II/source/grid/grid_refinement.cc
index df1e59b8b4f2e141f42e59cbed9df6d791d4e5cc..a6e36d378e70d48a4d941a42c82b4b785843cbb7 100644 (file)
--- a/deal.II/deal.II/source/grid/grid_refinement.cc
+++ b/deal.II/deal.II/source/grid/grid_refinement.cc
@@ -11,8 +11,6 @@
 //
 //----------------------------  grid_refinement.cc  ---------------------------
 
-
-#include <lac/vector.h>
 #include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
@@ -22,6 +20,55 @@
 #include <algorithm>
 #include <cmath>
 
+#include <math.h>
+#include <fstream>
+
+template<typename number>
+void GridRefinement::qsort_index(const Vector<number>  &a,
+                                vector<unsigned int> &ind,
+                                int                  l,
+                                int                  r)
+{
+  int i,j,t;
+  number v;
+  
+  if (r<=l)
+    return;
+
+  v = a(ind[r]);
+  i = l-1;
+  j = r;
+  do
+    {
+      do 
+       {
+         ++i;
+       }
+      while ((a(ind[i])>v)&&(i<r));
+      do
+       {
+         --j;
+       }
+      while ((a(ind[j])<v)&&(j>0));
+      
+      if (i<j)
+       {
+         t=ind[i];
+         ind[i] = ind[j];
+         ind[j] = t;
+       }
+      else
+       {
+         t = ind[i];
+         ind[i] = ind[r];
+         ind[r] = t;
+       }
+    }
+  while (i<j);
+  qsort_index(a,ind,l,i-1);
+  qsort_index(a,ind,i+1,r);  
+}
+
 
 
 
@@ -196,6 +243,56 @@ GridRefinement::refine_and_coarsen_fixed_fraction (Triangulation<dim>   &tria,
 
 
 
+template <int dim, typename number>
+void
+GridRefinement::refine_and_coarsen_optimize (Triangulation<dim>   &tria,
+                                            const Vector<number> &criteria)
+{
+  Assert (criteria.size() == tria.n_active_cells(),
+         ExcInvalidVectorSize(criteria.size(), tria.n_active_cells()));
+  
+                                  // get an increasing order on
+                                  // the error indicator
+  vector<unsigned int> tmp(criteria.size());
+  for (unsigned int i=0;i<criteria.size();++i)
+    tmp[i] = i;
+  
+  qsort_index(criteria,tmp,0,criteria.size()-1);
+  
+  double s0 = 0.75 * criteria(tmp[0]);
+  double E  = criteria.l1_norm();
+  
+  unsigned int N = criteria.size();
+  unsigned int M = 0;
+  
+                                  // The first M cells are refined
+                                  // to minimize the expected error
+                                  // multiplied with the expected
+                                  // number of cells.
+                                  // We assume that the error is
+                                  // decreased by 3/4 a_K if the cell
+                                  // K with error indicator a_K is
+                                  // refined.
+                                  // The expected number of cells is
+                                  // N+3*M (N is the current number
+                                  // of cells)
+  double min = (3.*(1.+M)+N) * (E-s0);
+  
+  unsigned int minArg = N-1;
+  
+  for (M=1;M<criteria.size();++M)
+    {
+      s0+= 0.75 * criteria(tmp[M]);
+      
+      if ( (3.*(1+M)+N)*(E-s0) <= min)
+       {
+         min = (3.*(1+M)+N)*(E-s0);
+         minArg = M;
+       }
+    }
+  refine(tria,criteria,criteria(tmp[minArg]));
+};
+
 // explicit instantiations
 template void GridRefinement
 ::refine (Triangulation<deal_II_dimension> &, const Vector<float> &, const double);
@@ -234,3 +331,12 @@ template void GridRefinement
                                     const double         top_fraction,
                                     const double         bottom_fraction);
 
+template void GridRefinement
+::refine_and_coarsen_optimize (Triangulation<deal_II_dimension> &,
+                              const Vector<float> &criteria);
+
+template void GridRefinement
+::refine_and_coarsen_optimize (Triangulation<deal_II_dimension> &,
+                              const Vector<double> &criteria);
+
+