From 1b21bf26cd4bce7c1fe7dd5ed3e31ab2a13db6f3 Mon Sep 17 00:00:00 2001
From: richter <richter@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 10 Oct 2000 09:48:11 +0000
Subject: [PATCH] add grid refinement with optimization

git-svn-id: https://svn.dealii.org/trunk@3401 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/grid/grid_refinement.h    |  60 +++++++++-
 .../deal.II/source/grid/grid_refinement.cc    | 110 +++++++++++++++++-
 2 files changed, 167 insertions(+), 3 deletions(-)

diff --git a/deal.II/deal.II/include/grid/grid_refinement.h b/deal.II/deal.II/include/grid/grid_refinement.h
index b93e3bfb7e..1815633c5a 100644
--- 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 df1e59b8b4..a6e36d378e 100644
--- 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);
+
+
-- 
2.39.5