]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New: GridGenerator::quarter_hyper_shell() in 3d.
authorheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 5 Mar 2012 16:42:52 +0000 (16:42 +0000)
committerheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 5 Mar 2012 16:42:52 +0000 (16:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@25210 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 177f59bbe271fbf23d3128616b02caf40cf3bd12..92fc9e9ffa0c1f4f008f38a037485694127df0e6 100644 (file)
@@ -145,6 +145,10 @@ enabled due to a missing include file in file
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: GridGenerator::quarter_hyper_shell() in 3d.
+<br>
+(Thomas Geenen, 2012/03/05)
+
 <li> New: DataOut::write_vtu_in_parallel(). This routine uses MPI I/O to write
 a big vtu file in parallel.
 <br>
index 6f8d623cce42bf727fbe8430554e7a451cfbfe29..6406061e86bb817b4695305ee7ba15e67a6b3ea2 100644 (file)
@@ -691,33 +691,28 @@ class GridGenerator
 
                                     /**
                                      * Produce a quarter hyper-shell,
-                                     * i.e. the space between two
-                                     * circles in two space
-                                     * dimensions and the region
-                                     * between two spheres in 3d,
-                                     * with given inner and outer
-                                     * radius and a given number of
-                                     * elements for this initial
-                                     * triangulation. All components are
-                                     * restricted to be positive, so the
-                                     * opening angle is 90 degrees.
-                                     *
-                                     * If the number of
-                                     * initial cells is zero (as is
-                                     * the default), then it is
-                                     * computed adaptively such that
-                                     * the resulting elements have
-                                     * the least aspect ratio.
-                                     *
-                                     * If colorize is set to true, the
-                                     * inner, outer, left, and right
-                                     * boundary get indicator 0, 1, 2,
-                                     * and 3, respectively. Otherwise
-                                     * all indicators are set to 0.
+                                     * i.e. the space between two circles in
+                                     * two space dimensions and the region
+                                     * between two spheres in 3d, with given
+                                     * inner and outer radius and a given
+                                     * number of elements for this initial
+                                     * triangulation.
+                                     *
+                                     * If the number of initial cells is zero
+                                     * (as is the default), then it is
+                                     * computed adaptively such that the
+                                     * resulting elements have the least
+                                     * aspect ratio in 2d.
+                                     *
+                                     * If colorize is set to true, the inner,
+                                     * outer, left, and right boundary get
+                                     * indicator 0, 1, 2, and 3 in 2d,
+                                     * respectively. Otherwise all indicators
+                                     * are set to 0. In 3d indicator 2 is at
+                                     * the face x=0, 3 at y=0, 4 at z=0.
                                      *
                                      * @note The triangulation needs to be
-                                     * void upon calling this
-                                     * function. Only implemented in 2d so far.
+                                     * void upon calling this function.
                                      */
     template <int dim>
     static void quarter_hyper_shell (Triangulation<dim>   &tria,
@@ -988,6 +983,22 @@ class GridGenerator
                          const double inner_radius,
                          const double outer_radius);
 
+
+                                    /**
+                                     * Assign boundary number zero the inner
+                                     * shell boundary, one to the outer shell
+                                     * boundary, two to the face with x=0,
+                                     * three to the face with y=0, four to
+                                     * the face with z=0.
+                                     */
+    template<int dim>
+    static
+    void
+    colorize_quarter_hyper_shell(Triangulation<dim> & tria,
+                                const Point<dim>& center,
+                                const double inner_radius,
+                                const double outer_radius);
+    
                                     /**
                                      * Solve the Laplace equation for
                                      * @p laplace_transformation
index 73f6e366141b5323309d10813f0722923dfd5440..745ff90461be011435fff7b960bfce7165c0a0de 100644 (file)
@@ -2841,16 +2841,138 @@ GridGenerator::half_hyper_shell (Triangulation<3>& tria,
 
 }
 
+// Implementation for 3D only
+template <>
+void
+GridGenerator::colorize_quarter_hyper_shell(Triangulation<3> & tria,
+                                            const Point<3>& center,
+                                            const double inner_radius,
+                                            const double outer_radius)
+{
+
+  if (tria.n_cells() != 3)
+    AssertThrow (false, ExcNotImplemented());
 
+  double middle = (outer_radius-inner_radius)/2e0 + inner_radius;
+  double eps = 1e-3*middle;
+  Triangulation<3>::cell_iterator cell = tria.begin();
+
+  for (;cell!=tria.end();++cell)
+    for(unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+      {
+        if(!cell->face(f)->at_boundary())
+          continue;
+
+        double radius = cell->face(f)->center().norm() - center.norm();
+        if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2
+          {
+            cell->face(f)->set_boundary_indicator(2);
+            for (unsigned int j=0;j<GeometryInfo<3>::lines_per_face;++j)
+              if(cell->face(f)->line(j)->at_boundary())
+                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                  cell->face(f)->line(j)->set_boundary_indicator(2);
+          }
+        else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3
+          {
+            cell->face(f)->set_boundary_indicator(3);
+            for (unsigned int j=0;j<GeometryInfo<3>::lines_per_face;++j)
+              if(cell->face(f)->line(j)->at_boundary())
+                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                  cell->face(f)->line(j)->set_boundary_indicator(3);
+          }
+        else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4
+          {
+            cell->face(f)->set_boundary_indicator(4);
+            for (unsigned int j=0;j<GeometryInfo<3>::lines_per_face;++j)
+              if(cell->face(f)->line(j)->at_boundary())
+                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                  cell->face(f)->line(j)->set_boundary_indicator(4);
+          }
+        else if (radius < middle) // inner radius set boundary 0
+          {
+            cell->face(f)->set_boundary_indicator(0);
+            for (unsigned int j=0;j<GeometryInfo<3>::lines_per_face;++j)
+              if(cell->face(f)->line(j)->at_boundary())
+                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
+                  cell->face(f)->line(j)->set_boundary_indicator(0);
+          }
+        else if (radius > middle) // outer radius set boundary 1
+          {
+            cell->face(f)->set_boundary_indicator(1);
+            for (unsigned int j=0;j<GeometryInfo<3>::lines_per_face;++j)
+              if(cell->face(f)->line(j)->at_boundary())
+                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
+                  cell->face(f)->line(j)->set_boundary_indicator(1);
+          }
+        else
+          AssertThrow (false, ExcInternalError());
+      }
+}
+
+
+
+// Implementation for 3D only
 template <>
-void GridGenerator::quarter_hyper_shell (Triangulation<3> &,
-                                        const Point<3>&,
-                                        const double,
-                                        const double,
-                                        const unsigned int,
-                                        const bool)
+void GridGenerator::quarter_hyper_shell (Triangulation<3> & tria,
+                                         const Point<3>& center,
+                                         const double inner_radius,
+                                         const double outer_radius,
+                                         const unsigned int n,
+                                         const bool colorize)
 {
-  Assert (false, ExcNotImplemented());
+  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+          ExcInvalidRadii ());
+  if (n == 0 || n == 3)
+    {
+      const double a = inner_radius*std::sqrt(2.0)/2e0;
+      const double b = outer_radius*std::sqrt(2.0)/2e0;
+      const double c = a*std::sqrt(3.0)/2e0;
+      const double d = b*std::sqrt(3.0)/2e0;
+      const double e = outer_radius/2e0;
+      const double f = outer_radius*std::sqrt(3.0)/2e0;
+      const double g = f*std::sqrt(3.0)/2e0;
+      const double h = inner_radius/2e0;
+
+      std::vector<Point<3> > vertices;
+
+      vertices.push_back (center+Point<3>( 0,  inner_radius, 0)); //0
+      vertices.push_back (center+Point<3>( a,  a, 0));                 //1
+      vertices.push_back (center+Point<3>( b,  b, 0));                 //2
+      vertices.push_back (center+Point<3>( 0, outer_radius, 0));       //3
+      vertices.push_back (center+Point<3>( 0, a , a));                 //4
+      vertices.push_back (center+Point<3>( c, c , h));                 //5
+      vertices.push_back (center+Point<3>( d, d , e));                 //6
+      vertices.push_back (center+Point<3>( 0, b , b));                 //7
+      vertices.push_back (center+Point<3>( inner_radius, 0 , 0));      //8
+      vertices.push_back (center+Point<3>( outer_radius, 0 , 0));      //9
+      vertices.push_back (center+Point<3>( a, 0 , a)); //10
+      vertices.push_back (center+Point<3>( b, 0 , b)); //11
+      vertices.push_back (center+Point<3>( 0, 0 , inner_radius));      //12
+      vertices.push_back (center+Point<3>( 0, 0 , outer_radius));      //13
+
+      const int cell_vertices[3][8] = {
+        {0, 1, 3, 2, 4, 5, 7, 6},
+        {1, 8, 2, 9, 5, 10, 6, 11},
+        {4, 5, 7, 6, 12, 10, 13, 11},
+      };
+      std::vector<CellData<3> > cells(n);
+
+      for (unsigned int i=0; i<n; ++i)
+        {
+          for (unsigned int j=0; j<8; ++j)
+            cells[i].vertices[j] = cell_vertices[i][j];
+          cells[i].material_id = 0;
+        }
+
+      tria.create_triangulation ( vertices, cells, SubCellData());       // no boundary information
+    }
+  else
+    {
+      AssertThrow(false, ExcNotImplemented());
+    }
+
+  if (colorize)
+      colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
 }
 
 

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.