]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Yaqi's function for non-uniformly graded hyper rectangles.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Nov 2006 18:04:46 +0000 (18:04 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Nov 2006 18:04:46 +0000 (18:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@14196 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/source/grid/grid_generator.cc
deal.II/doc/news/changes.html

index e6697bf84305c0bfac4cccb7de3fcc14fa926bed..16190bbcb717626657562fd03e1a2b769c5c524a 100644 (file)
@@ -267,6 +267,44 @@ class GridGenerator
                                const Point<dim>                &p2,
                                const bool                       colorize=false);
 
+                                    /**
+                                     * Like the previous
+                                     * function. However, here the
+                                     * second argument does not
+                                     * denote the number of
+                                     * subdivisions in each
+                                     * coordinate direction, but a
+                                     * sequence of step sizes for
+                                     * each coordinate direction. The
+                                     * domain will therefore be
+                                     * subdivided into
+                                     * <code>step_sizes[i].size()</code>
+                                     * cells in coordinate direction
+                                     * <code>i</code>, with widths
+                                     * <code>step_sizes[i][j]</code>
+                                     * for the <code>j</code>th cell.
+                                     *
+                                     * This function is therefore the
+                                     * right one to generate graded
+                                     * meshes where cells are
+                                     * concentrated in certain areas,
+                                     * rather than a uniformly
+                                     * subdivided mesh as the
+                                     * previous function generates.
+                                     *
+                                     * The step sizes have to add up
+                                     * to the dimensions of the hyper
+                                     * rectangle specified by the
+                                     * points @p p1 and @p p2.
+                                     */
+    template <int dim>
+    void
+    subdivided_hyper_rectangle(Triangulation<dim>              &tria,
+                              const std::vector<std::vector<double> > &step_sizes,
+                              const Point<dim>                &p_1,
+                              const Point<dim>                &p_2,
+                              const bool                       colorize);
+    
                                     /**
                                      * A parallelogram. The first
                                      * corner point is the
index 9f34e64200f81faaac2f22a3a7e84448862a49ab..cfb251096d00331f3ca0d32d255269828293784b 100644 (file)
@@ -442,6 +442,202 @@ subdivided_hyper_rectangle (Triangulation<dim>              &tria,
 }
 
 
+
+template <int dim>
+void
+GridGenerator::
+subdivided_hyper_rectangle(Triangulation<dim>              &tria,
+                           const std::vector<std::vector<double> > &step_sz,
+                           const Point<dim>                &p_1,
+                           const Point<dim>                &p_2,
+                           const bool                       colorize)
+{
+                                  // contributed by Joerg R. Weimar
+                                  // (j.weimar@jweimar.de) 2003
+                                  // modified by Yaqi Wang 2006
+  Assert(step_sz.size() == dim, 
+        ExcInvalidRepetitionsDimension(dim));
+
+
+                                  // First, normalize input such that
+                                  // p1 is lower in all coordinate
+                                  // directions.
+
+                                   // and check the consistency of
+                                   // step sizes, i.e. that they all
+                                   // add up to the sizes specified by
+                                   // p_1 and p_2
+  Point<dim> p1(p_1);
+  Point<dim> p2(p_2);
+  std::vector< std::vector<double> > step_sizes(step_sz);
+
+  for (unsigned int i=0;i<dim;++i)
+    {
+      if (p1(i) > p2(i))
+       {
+         std::swap (p1(i), p2(i));
+         std::reverse (step_sizes[i].begin(), step_sizes[i].end());
+       }
+
+      double x = 0;
+      for (unsigned int j=0; j<step_sizes.at(i).size(); j++) 
+       x += step_sizes[i][j];
+      Assert(std::fabs(x - (p2(i)-p1(i))) <= 1e-12*std::fabs(x),
+            ExcInvalidRepetitions (i) );
+    }
+
+                                   // then generate the necessary
+                                   // points
+  std::vector<Point<dim> > points;
+  switch (dim)
+    {
+      case 1:
+      {
+       double x=0;
+       for (unsigned int i=0; i<=step_sizes[0].size(); ++i)
+         { 
+           points.push_back (Point<dim> (p1[0]+x));
+           x += step_sizes[0][i];
+         }
+       break;
+      }
+      
+      case 2:
+      {
+       double y=0;
+       for (unsigned int j=0; j<=step_sizes[1].size(); ++j)
+         {
+           double x=0;
+           for (unsigned int i=0; i<=step_sizes[0].size(); ++i)
+             {
+               points.push_back (Point<dim> (p1[0]+x,
+                                             p1[1]+y));
+               x += step_sizes[0][i];
+             }
+           y += step_sizes[1][j];
+         }
+       break;
+
+      }
+      case 3:
+      {
+       double z=0;
+       for (unsigned int k=0; k<=step_sizes[2].size(); ++k)
+         {
+           double y=0;
+           for (unsigned int j=0; j<=step_sizes[1].size(); ++j)
+             {
+               double x=0;
+               for (unsigned int i=0; i<=step_sizes[0].size(); ++i)
+                 {
+                   points.push_back (Point<dim> (p1[0]+x,
+                                                 p1[1]+y,
+                                                 p1[2]+z));
+                   x += step_sizes[0][i];
+                 }
+               y += step_sizes[1][j];
+             }
+           z += step_sizes[2][k];
+         }
+       break;
+      }
+      
+      default:
+            Assert (false, ExcNotImplemented());
+    }
+
+                                   // next create the cells
+                                  // Prepare cell data
+  std::vector<CellData<dim> > cells;
+  switch (dim)
+    {
+      case 1:
+      {
+        cells.resize (step_sizes[0].size());
+        for (unsigned int x=0; x<step_sizes[0].size(); ++x)
+          {
+            cells[x].vertices[0] = x;
+            cells[x].vertices[1] = x+1;
+            cells[x].material_id = 0;
+          }
+        break;
+      }
+      
+      case 2:
+      {
+        cells.resize (step_sizes[1].size()*step_sizes[0].size());
+        for (unsigned int y=0; y<step_sizes[1].size(); ++y)
+          for (unsigned int x=0; x<step_sizes[0].size(); ++x)
+            {
+              const unsigned int c = x+y*step_sizes[0].size();
+              cells[c].vertices[0] = y*(step_sizes[0].size()+1)+x;
+              cells[c].vertices[1] = y*(step_sizes[0].size()+1)+x+1;
+              cells[c].vertices[2] = (y+1)*(step_sizes[0].size()+1)+x;
+              cells[c].vertices[3] = (y+1)*(step_sizes[0].size()+1)+x+1;
+              cells[c].material_id = 0;
+            }
+        break;
+      }
+      
+      case 3:
+      {
+        const unsigned int n_x  = (step_sizes[0].size()+1);
+        const unsigned int n_xy = (step_sizes[0].size()+1)*(step_sizes[1].size()+1);
+        
+        cells.resize (step_sizes[2].size()*step_sizes[1].size()*step_sizes[0].size());
+        for (unsigned int z=0; z<step_sizes[2].size(); ++z)
+          for (unsigned int y=0; y<step_sizes[1].size(); ++y)
+            for (unsigned int x=0; x<step_sizes[0].size(); ++x)
+              {
+                const unsigned int c = x+y*step_sizes[0].size() +
+                                       z*step_sizes[0].size()*step_sizes[1].size();
+                cells[c].vertices[0] = z*n_xy + y*n_x + x;
+                cells[c].vertices[1] = z*n_xy + y*n_x + x+1;
+                cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
+                cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x;
+                cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
+                cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+                cells[c].material_id = 0;
+              }
+        break;
+        
+      }
+
+      default:
+            Assert (false, ExcNotImplemented());
+    }
+
+  tria.create_triangulation (points, cells, SubCellData());  
+
+  if (colorize)
+    {
+                                      // to colorize, run through all
+                                      // faces of all cells and set
+                                      // boundary indicator to the
+                                      // correct value if it was 0.
+     
+                                      // use a large epsilon to
+                                      // compare numbers to avoid
+                                      // roundoff problems.
+      double min_size = *std::min_element (step_sizes[0].begin(),
+                                          step_sizes[0].end());
+      for (unsigned int i=1; i<dim; ++i)
+       min_size = std::min (min_size,
+                            *std::min_element (step_sizes[i].begin(),
+                                               step_sizes[i].end()));
+      const double epsilon = 0.01 * min_size;
+    
+                                       // actual code is external since
+                                       // 1-D is different from 2/3D.
+      colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
+    }
+}
+
+
+
 #if deal_II_dimension == 1
 
 // Implementation for 1D only
@@ -1731,6 +1927,13 @@ GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
  const Point<deal_II_dimension>&,
  const Point<deal_II_dimension>&, bool);
 
+template void
+GridGenerator::subdivided_hyper_rectangle<deal_II_dimension> 
+(Triangulation<deal_II_dimension> &,
+ const std::vector<std::vector<double> >&,
+ const Point<deal_II_dimension>&,
+ const Point<deal_II_dimension>&, bool);
+
 template void
 GridGenerator::parallelogram<deal_II_dimension> (
   Triangulation<deal_II_dimension> &,
index 7130190b1ef83649fc92d810360ddcc1fc0e75e9..e0b9c9c9d22a2bd07e159a017a732bc4e1fdc03a 100644 (file)
@@ -418,7 +418,7 @@ inconvenience this causes.
   <li> <p> Extended: The <code class="class">QGaussLobatto</code> quadrature
        rule computes Legendre-Gauss-Lobatto nodes and quadrature weights.
        <br>
-       (F. Prill 2006/11/02)
+       (Florian Prill 2006/11/02)
        </p>
 
   <li> <p> Extended: The <code>contract</code> function family contracts two
@@ -717,7 +717,7 @@ inconvenience this causes.
   <li> <p> Improved: PETSc 2.3.1 (<i>build 16</i>) is now supported by the
        linear preconditioner classes. The new PETSc functions 
        <code class="member">PCFactorSetXXX</code> are used.<br>       
-       (F. Prill, 2006/08/04)
+       (Florian Prill, 2006/08/04)
        </p>
 
   <li> <p>New: The class <code class="class">TransposeMatrix</code>
@@ -922,11 +922,19 @@ inconvenience this causes.
 
 <ol>
 
+  <li> <p> New: There is now a function <code
+       class="member">GridGenerator::subdivided_hyper_rectangle</code>
+       that produces a non-uniformly subdivided rectangle, ideally
+       suited for graded meshes.
+       <br>
+       (Yaqi Wang 2006/11/15)
+       </p>
+
   <li> <p> Fixed: Corrected <code class="member">clone</code> method 
        of <code class="class">FE_DGQ</code> class for non-equidistant 
        support points.
        <br>
-       (F. Prill 2006/10/31)
+       (Florian Prill 2006/10/31)
        </p>
 
   <li> <p> Improved: The lookup mechanism in <code

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.