]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Daniel Castanon Quiroz: Implement generating a torus and describing the...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Sep 2011 02:23:06 +0000 (02:23 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Sep 2011 02:23:06 +0000 (02:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@24288 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 6adb15686de02a8b583f40553a89bf43237b136a..2055fe1710718418a21b499e3796d3236690af71 100644 (file)
@@ -283,13 +283,17 @@ and DoF handlers embedded in higher dimensional space have been added.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: The GridGenerator::torus function and TorusBoundary class can
+create and describe the surface of a torus.
+<br>
+(Daniel Castanon Quiroz, 2011/09/08)
+
 <li> Fixed: FEFieldFunction class was not thread-safe because it keeps a
 cache on the side that was invalidated when different
 threads kept pouncing on it. This is now fixed.
 <br>
 (Patrick Sodr&eacute;, 2011/09/07)
 
-
 <li> New: There is now a function GridTools::volume computing the volume
 of a triangulation.
 <br>
index 0e5bc833575563ac31cb7303eadfec5738e9f176..9981b9563aea517ccbab524d9d6ac12f35803188 100644 (file)
@@ -47,9 +47,6 @@ template <typename number> class SparseMatrix;
  *
  *
  * @ingroup grid
- * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan
- * Nauber, Joerg Weimar, Yaqi Wang, Luca Heltai, 1998, 1999, 2000, 2001, 2002,
- * 2003, 2006, 2007, 2008, 2009, 2010, 2011.
  */
 class GridGenerator
 {
@@ -673,11 +670,11 @@ class GridGenerator
                                      * 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 
+                                     * and 3, respectively. Otherwise
                                      * all indicators are set to 0.
                                      *
                                      * @note The triangulation needs to be
@@ -692,7 +689,7 @@ class GridGenerator
                                  const unsigned int  n_cells = 0,
                                  const bool colorize = false);
 
-    
+
                                     /**
                                      * Produce a quarter hyper-shell,
                                      * i.e. the space between two
@@ -730,7 +727,7 @@ class GridGenerator
                                  const double        outer_radius,
                                  const unsigned int  n_cells = 0,
                                  const bool colorize = false);
-   
+
                                     /**
                                      * Produce a domain that is the space
                                      * between two cylinders in 3d, with
@@ -761,6 +758,30 @@ class GridGenerator
                                 const unsigned int  n_radial_cells = 0,
                                 const unsigned int  n_axial_cells = 0);
 
+
+
+                                     /**
+                                     * Produce the surface meshing of the
+                                     * torus. The axis of the torus is the
+                                     * $y$-axis while the plane of the torus
+                                     * is the $x$-$z$ plane. The boundary of
+                                     * this object can be described by the
+                                     * TorusBoundary class.
+                                     *
+                                     * @param R The radius of the circle,
+                                     * which forms the middle line of the
+                                     * torus containing the loop of
+                                     * cells. Must be greater than @p r.
+                                     *
+                                     * @param r The inner radius of the
+                                     * torus.
+                                     */
+
+    static void torus (Triangulation<2,3>&  tria,
+                      const double         R,
+                      const double         r);
+
+
                                     /**
                                      * This class produces a square
                                      * on the <i>xy</i>-plane with a
@@ -826,8 +847,8 @@ class GridGenerator
                                      * the third parameter. Previous
                                      * content of @p result will be
                                      * deleted.
-                                     * 
-                                     * This function is most often used 
+                                     *
+                                     * This function is most often used
                                      * to compose meshes for more
                                      * complicated geometries if the
                                      * geometry can be composed of
@@ -835,7 +856,7 @@ class GridGenerator
                                      * exist to generate coarse meshes.
                                      * For example, the channel mesh used
                                      * in step-35 could in principle be
-                                     * created using a mesh created by the 
+                                     * created using a mesh created by the
                                      * GridGenerator::hyper_cube_with_cylindrical_hole
                                      * function and several rectangles,
                                      * and merging them using the current
@@ -847,11 +868,11 @@ class GridGenerator
                                      * individual mesh building blocks are
                                      * GridTools::transform, GridTools::rotate,
                                      * and GridTools::scale).
-                                     * 
+                                     *
                                      * @note The two input triangulations
                                      * must be coarse meshes that have
-                                     * no refined cells. 
-                                     * 
+                                     * no refined cells.
+                                     *
                                      * @note The function copies the material ids
                                      * of the cells of the two input
                                      * triangulations into the output
@@ -863,7 +884,7 @@ class GridGenerator
                                      * then you will currently have to set
                                      * boundary indicators again by hand
                                      * in the output triangulation.
-                                     * 
+                                     *
                                      * @note For a related operation
                                      * on refined meshes when both
                                      * meshes are derived from the
index eb5f68f804063c5cfc06a84e0b10b1a2e0e87869..42ef1c8681d443d6817c3b6b6882452096fd8e1c 100644 (file)
@@ -388,7 +388,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * this class and the documentation of the
                                      * base class.
                                      */
-    virtual 
+    virtual
     Point<spacedim>
     get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
 
@@ -397,7 +397,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * this class and the documentation of the
                                      * base class.
                                      */
-    virtual 
+    virtual
     Point<spacedim>
     get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
 
@@ -410,7 +410,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * Calls
                                      * @p get_intermediate_points_between_points.
                                      */
-    virtual 
+    virtual
     void
     get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
                                     std::vector<Point<spacedim> > &points) const;
@@ -424,7 +424,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * Only implemented for <tt>dim=3</tt>
                                      * and for <tt>points.size()==1</tt>.
                                      */
-    virtual 
+    virtual
     void
     get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
                                     std::vector<Point<spacedim> > &points) const;
@@ -453,7 +453,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * and the documentation of the
                                      * base class.
                                      */
-    virtual 
+    virtual
     void
     get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
                             typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
@@ -467,7 +467,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                     /**
                                      * Return the radius of the ball.
                                      */
-    double 
+    double
     get_radius () const;
 
                                     /**
@@ -734,6 +734,109 @@ class HalfHyperShellBoundary : public HyperShellBoundary<dim>
 };
 
 
+/**
+ * Class describing the boundary of the torus. The axis of the torus is the
+ * $y$-axis while the plane of the torus is the $x$-$z$ plane. A torus of this
+ * kind can be generated by GridGenerator::torus.
+ *
+ * This class is only implemented for
+ * <tt>dim=2</tt>,<tt>spacedim=3</tt>, that is, just the surface.
+ */
+template <int dim, int spacedim>
+class TorusBoundary : public Boundary<dim,spacedim>
+{
+  public:
+                                     /**
+                                     * Constructor.<tt>R</tt> has to be
+                                     * greater than <tt>r</tt>.
+                                     */
+    TorusBoundary (const double R, const double r);
+
+//Boundary Refinenment  Functions
+                                     /**
+                                     * Construct a new point on a line.
+                                     */
+    virtual Point<spacedim>
+    get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
+
+                                     /**
+                                     * Construct a new point on a quad.
+                                     */
+    virtual Point<spacedim>
+    get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
+
+                                     /**
+                                     * Construct a new points on a line.
+                                     */
+    virtual void   get_intermediate_points_on_line (
+      const typename Triangulation< dim, spacedim >::line_iterator &  line,
+      std::vector< Point< spacedim > > &       points) const ;
+
+                                     /**
+                                     * Construct a new points on a quad.
+                                     */
+    virtual void  get_intermediate_points_on_quad (
+      const typename Triangulation< dim, spacedim >::quad_iterator & quad,
+      std::vector< Point< spacedim > > &       points ) const ;
+
+                                     /**
+                                     * Get the normal from cartesian
+                                     * coordinates. This normal does not have
+                                     * unit length.
+                                     */
+    virtual void get_normals_at_vertices (
+      const typename Triangulation< dim, spacedim >::face_iterator &face,
+      typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const ;
+
+  private:
+                                    //Handy functions
+                                     /**
+                                     * Function that corrects the value and
+                                     * sign of angle, that is, given
+                                     * <tt>angle=tan(abs(y/x))<tt>; if <tt>
+                                     * (y > 0) && (x < 0) <tt> then
+                                     * <tt>correct_angle = Pi - angle<tt>,
+                                     * etc.
+                                     */
+
+    double           get_correct_angle(const double angle,const double x,const double y) const;
+
+                                     /**
+                                     * Get the cartesian coordinates of the
+                                     * Torus, i.e., from <tt>(theta,phi)<tt>
+                                     * to <tt>(x,y,z)<tt>.
+                                     */
+    Point<spacedim>  get_real_coord(const Point<dim>& surfP) const ;
+
+                                     /**
+                                     * Get the surface coordinates of the
+                                     * Torus, i.e., from <tt>(x,y,z)<tt> to
+                                     * <tt>(theta,phi)<tt>.
+                                     */
+    Point<dim>       get_surf_coord(const Point<spacedim>& p) const ;
+
+                                     /**
+                                     * Get the normal from surface
+                                     * coordinates. This normal does not have
+                                     * unit length.
+                                     */
+    Point<spacedim>  get_surf_norm_from_sp(const Point<dim>& surfP)      const ;
+
+                                     /**
+                                     * Get the normal from cartesian
+                                     * coordinates. This normal does not have
+                                     * unit length.
+                                     */
+    Point<spacedim>  get_surf_norm(const Point<spacedim>& p) const ;
+
+                                    /**
+                                     * Inner and outer radii of the shell.
+                                     */
+    const double R;
+    const double r;
+};
+
+
 
 /* -------------- declaration of explicit specializations ------------- */
 
@@ -787,6 +890,7 @@ HalfHyperShellBoundary<1>::
 get_normals_at_vertices (const Triangulation<1,1>::face_iterator &,
                         Boundary<1,1>::FaceVertexNormals &) const;
 
+
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
index 515f92b1f2d1c11a905de9980b7a0b2cdc93186e..1e4bd92f11bf658578d78b1e271607735c166c84 100644 (file)
@@ -232,6 +232,141 @@ GridGenerator::moebius (
 
 
 
+void
+GridGenerator::torus (Triangulation<2,3>&  tria,
+                     const double         R,
+                     const double         r)
+{
+  Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+
+  const unsigned int dim=2;
+  const unsigned int spacedim=3;
+  std::vector<Point<spacedim> > vertices (16);
+
+  vertices[0]=Point<spacedim>(R-r,0,0);
+  vertices[1]=Point<spacedim>(R,-r,0);
+  vertices[2]=Point<spacedim>(R+r,0,0);
+  vertices[3]=Point<spacedim>(R, r,0);
+  vertices[4]=Point<spacedim>(0,0,R-r);
+  vertices[5]=Point<spacedim>(0,-r,R);
+  vertices[6]=Point<spacedim>(0,0,R+r);
+  vertices[7]=Point<spacedim>(0,r,R);
+  vertices[8]=Point<spacedim>(-(R-r),0,0);
+  vertices[9]=Point<spacedim>(-R,-r,0);
+  vertices[10]=Point<spacedim>(-(R+r),0,0);
+  vertices[11]=Point<spacedim>(-R, r,0);
+  vertices[12]=Point<spacedim>(0,0,-(R-r));
+  vertices[13]=Point<spacedim>(0,-r,-R);
+  vertices[14]=Point<spacedim>(0,0,-(R+r));
+  vertices[15]=Point<spacedim>(0,r,-R);
+
+  std::vector<CellData<dim> > cells (16);
+                                  //Right Hand Orientation
+  cells[0].vertices[0] =  0;
+  cells[0].vertices[1] =  4;
+  cells[0].vertices[2] =  7;
+  cells[0].vertices[3] =  3;
+  cells[0].material_id = 0;
+
+  cells[1].vertices[0] =  1;
+  cells[1].vertices[1] =  5;
+  cells[1].vertices[2] =  4;
+  cells[1].vertices[3] =  0;
+  cells[1].material_id = 0;
+
+  cells[2].vertices[0] =  2;
+  cells[2].vertices[1] =  6;
+  cells[2].vertices[2] =  5;
+  cells[2].vertices[3] =  1;
+  cells[2].material_id = 0;
+
+  cells[3].vertices[0] =  3;
+  cells[3].vertices[1] =  7;
+  cells[3].vertices[2] =  6;
+  cells[3].vertices[3] =  2;
+  cells[3].material_id = 0;
+
+  cells[4].vertices[0] =  4;
+  cells[4].vertices[1] =  8;
+  cells[4].vertices[2] =  11;
+  cells[4].vertices[3] =  7;
+  cells[4].material_id = 0;
+
+  cells[5].vertices[0] =  5;
+  cells[5].vertices[1] =  9;
+  cells[5].vertices[2] =  8;
+  cells[5].vertices[3] =  4;
+  cells[5].material_id = 0;
+
+  cells[6].vertices[0] =  6;
+  cells[6].vertices[1] =  10;
+  cells[6].vertices[2] =  9;
+  cells[6].vertices[3] =  5;
+  cells[6].material_id = 0;
+
+  cells[7].vertices[0] =  7;
+  cells[7].vertices[1] =  11;
+  cells[7].vertices[2] =  10;
+  cells[7].vertices[3] =  6;
+  cells[7].material_id = 0;
+
+  cells[8].vertices[0] =  8;
+  cells[8].vertices[1] =  12;
+  cells[8].vertices[2] =  15;
+  cells[8].vertices[3] =  11;
+  cells[8].material_id = 0;
+
+  cells[9].vertices[0] =  9;
+  cells[9].vertices[1] =  13;
+  cells[9].vertices[2] =  12;
+  cells[9].vertices[3] =  8;
+  cells[9].material_id = 0;
+
+  cells[10].vertices[0] =  10;
+  cells[10].vertices[1] =  14;
+  cells[10].vertices[2] =  13;
+  cells[10].vertices[3] =  9;
+  cells[10].material_id = 0;
+
+  cells[11].vertices[0] =  11;
+  cells[11].vertices[1] =  15;
+  cells[11].vertices[2] =  14;
+  cells[11].vertices[3] =  10;
+  cells[11].material_id = 0;
+
+  cells[12].vertices[0] =  12;
+  cells[12].vertices[1] =  0;
+  cells[12].vertices[2] =  3;
+  cells[12].vertices[3] =  15;
+  cells[12].material_id = 0;
+
+  cells[13].vertices[0] =  13;
+  cells[13].vertices[1] =  1;
+  cells[13].vertices[2] =  0;
+  cells[13].vertices[3] =  12;
+  cells[13].material_id = 0;
+
+  cells[14].vertices[0] =  14;
+  cells[14].vertices[1] =  2;
+  cells[14].vertices[2] =  1;
+  cells[14].vertices[3] =  13;
+  cells[14].material_id = 0;
+
+  cells[15].vertices[0] =  15;
+  cells[15].vertices[1] =  3;
+  cells[15].vertices[2] =  2;
+  cells[15].vertices[3] =  14;
+  cells[15].material_id = 0;
+
+                                  // Must call this to be able to create a
+                                  // correct triangulation in dealii, read
+                                  // GridReordering<> doc
+  GridReordering<dim,spacedim>::reorder_cells (cells);
+  tria.create_triangulation_compatibility (vertices, cells, SubCellData());
+}
+
+
+
 
 // Implementation for 2D only
 template<>
@@ -1697,7 +1832,7 @@ GridGenerator::half_hyper_shell (Triangulation<2>   &tria,
     };
 
   tria.create_triangulation (vertices, cells, SubCellData());
-  
+
   if (colorize)
   {
       Triangulation<2>::cell_iterator cell = tria.begin();
@@ -1706,7 +1841,7 @@ GridGenerator::half_hyper_shell (Triangulation<2>   &tria,
        cell->face(2)->set_boundary_indicator(1);
       }
       tria.begin()->face(0)->set_boundary_indicator(3);
-      
+
       tria.last()->face(1)->set_boundary_indicator(2);
   }
 }
@@ -1776,8 +1911,8 @@ void GridGenerator::quarter_hyper_shell (Triangulation<2>   &tria,
       cells[i].material_id = 0;
     };
 
-  tria.create_triangulation (vertices, cells, SubCellData());   
-  
+  tria.create_triangulation (vertices, cells, SubCellData());
+
   if (colorize)
   {
       Triangulation<2>::cell_iterator cell = tria.begin();
@@ -1786,7 +1921,7 @@ void GridGenerator::quarter_hyper_shell (Triangulation<2>   &tria,
        cell->face(2)->set_boundary_indicator(1);
       }
       tria.begin()->face(0)->set_boundary_indicator(3);
-      
+
       tria.last()->face(1)->set_boundary_indicator(2);
   }
 }
@@ -2813,13 +2948,13 @@ merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
          ExcMessage ("The input triangulations must be coarse meshes."));
   Assert (triangulation_2.n_levels() == 1,
          ExcMessage ("The input triangulations must be coarse meshes."));
-  
+
   // get the union of the set of vertices
   std::vector<Point<spacedim> > vertices = triangulation_1.get_vertices();
   vertices.insert (vertices.end(),
                   triangulation_2.get_vertices().begin(),
                   triangulation_2.get_vertices().end());
-  
+
   // now form the union of the set of cells
   std::vector<CellData<dim> > cells;
   cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells());
@@ -2832,8 +2967,8 @@ merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
     this_cell.material_id = cell->material_id();
     cells.push_back (this_cell);
   }
-  
-  // now do the same for the other other mesh. note that we have to 
+
+  // now do the same for the other other mesh. note that we have to
   // translate the vertex indices
   for (typename Triangulation<dim,spacedim>::cell_iterator
     cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell)
@@ -2844,7 +2979,7 @@ merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
     this_cell.material_id = cell->material_id();
     cells.push_back (this_cell);
   }
-  
+
   // throw out duplicated vertices from the two meshes
   // and create the triangulation
   SubCellData subcell_data;
index 2cf45eac6409940316b539c11dfa880a1b70780c..88074b4c1968929ad2ba3168dad01a39ca97d9e2 100644 (file)
@@ -1220,6 +1220,304 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
 
 
 
+
+template <int dim, int spacedim>
+TorusBoundary<dim,spacedim>::TorusBoundary (const double R__,
+                                           const double r__)
+               :
+               R(R__),
+               r(r__)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+TorusBoundary<2,3>::TorusBoundary (const double R__,
+                                  const double r__)
+               :
+               R(R__),
+               r(r__)
+{
+  Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+}
+
+
+
+template <int dim, int spacedim>
+double
+TorusBoundary<dim,spacedim>::get_correct_angle(const double angle,
+                                              const double x,
+                                              const double y) const
+{
+  if(y>=0)
+    {
+      if (x >=0)
+       return angle;
+
+      return numbers::PI-angle;
+    }
+
+  if (x <=0)
+    return numbers::PI+angle;
+
+  return 2.0*numbers::PI-angle;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_real_coord (const Point<2>& surfP) const
+{
+  const double theta=surfP(0);
+  const double phi=surfP(1);
+
+  return Point<3> ((R+r*std::cos(phi))*std::cos(theta),
+                  r*std::sin(phi),
+                  (R+r*std::cos(phi))*std::sin(theta));
+}
+
+
+
+template <>
+Point<2>
+TorusBoundary<2,3>::get_surf_coord(const Point<3>& p) const
+{
+  const double phi=std::asin(std::abs(p(1))/r);
+  const double Rr_2=p(0)*p(0)+p(2)*p(2);
+
+  Point<2> surfP;
+  surfP(1)=get_correct_angle(phi,Rr_2-R*R,p(1));//phi
+
+  if(std::abs(p(0))<1.E-5)
+    {
+      if(p(2)>=0)
+       surfP(0) =  numbers::PI*0.5;
+      else
+       surfP(0) = -numbers::PI*0.5;
+    }
+  else
+    {
+      const double theta = std::atan(std::abs(p(2)/p(0)));
+      surfP(0)=  get_correct_angle(theta,p(0),p(2));
+    }
+
+  return surfP;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_new_point_on_line (const typename Triangulation<2,3>::line_iterator &line) const
+{
+                                  //Just get the average
+  Point<2>  p0=get_surf_coord(line->vertex(0));
+  Point<2>  p1=get_surf_coord(line->vertex(1));
+
+  Point<2>  middle(0,0);
+
+                                  //Take care for periodic conditions,
+                                  //For instance phi0= 0, phi1= 3/2*Pi  middle has to be 7/4*Pi not 3/4*Pi
+                                  //This also works for -Pi/2 + Pi, middle is 5/4*Pi
+  for(unsigned i=0; i<2;i++)
+    if(std::abs(p0(i)-p1(i))> numbers::PI)
+      middle(i)=2*numbers::PI;
+
+  middle+=  p0 + p1;
+  middle*=0.5;
+
+  Point<3> midReal=get_real_coord(middle);
+  return midReal;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_new_point_on_quad (const typename Triangulation<2,3>::quad_iterator &quad) const
+{
+                                  //Just get the average
+  Point<2> p[4];
+
+  for(unsigned i=0;i<4;i++)
+    p[i]=get_surf_coord(quad->vertex(i));
+
+  Point<2>  middle(0,0);
+
+                                  //Take care for periodic conditions, see get_new_point_on_line() above
+                                  //For instance phi0= 0, phi1= 3/2*Pi  middle has to be 7/4*Pi not 3/4*Pi
+                                  //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi
+  for(unsigned i=0;i<2;i++)
+    for(unsigned j=1;j<4;j++){
+      if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+       middle(i)+=2*numbers::PI;
+      }
+    }
+
+  for(unsigned i=0;i<4;i++)
+    middle+=p[i];
+
+  middle*= 0.25;
+
+  return get_real_coord(middle);
+}
+
+
+
+//Normal field without unit length
+template <>
+Point<3>
+TorusBoundary<2,3>:: get_surf_norm_from_sp(const Point<2>& surfP) const
+{
+
+  Point<3> n;
+  double theta=surfP[0];
+  double phi=surfP[1];
+
+  double f=R+r*std::cos(phi);
+
+  n[0]=r*std::cos(phi)*std::cos(theta)*f;
+  n[1]=r*std::sin(phi)*f;
+  n[2]=r*std::sin(theta)*std::cos(phi)*f;
+
+  return n;
+}
+
+
+
+//Normal field without unit length
+template <>
+Point<3>
+TorusBoundary<2,3>::get_surf_norm(const Point<3>& p) const{
+
+  Point<2> surfP=get_surf_coord(p);
+  return get_surf_norm_from_sp(surfP);
+
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_intermediate_points_on_line (const typename Triangulation<2, 3>::line_iterator &  line,
+                                std::vector< Point< 3 > > & points) const
+{
+                                  //Almost the same implementation as  StraightBoundary<2,3>
+  unsigned npoints=points.size();
+  if(npoints==0) return;
+
+  Point<2> p[2];
+
+  for(unsigned i=0;i<2;i++)
+    p[i]=get_surf_coord(line->vertex(i));
+
+  unsigned  offset[2];
+  offset[0]=0;
+  offset[1]=0;
+
+                                  //Take care for periodic conditions & negative angles,
+                                  //see get_new_point_on_line() above
+                                  //Because we dont have a symmetric interpolation (just the middle) we need to
+                                  //add 2*Pi to each almost zero and negative angles.
+  for(unsigned i=0;i<2;i++)
+    for(unsigned j=1;j<2;j++){
+      if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+       offset[i]++;
+       break;
+      }
+    }
+
+  for(unsigned  i=0;i<2;i++)
+    for(unsigned  j=0;j<2;j++)
+      if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+       p[j](i)+=2*numbers::PI*offset[i];
+
+
+  double dx=1.0/(npoints+1);
+  double x=dx;
+  Point<2>  target;
+  for(unsigned i=0; i<npoints;i++,x+=dx){
+    target=  (1-x)*p[0] + x*p[1];
+    points[i]=get_real_coord(target);
+  }
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_intermediate_points_on_quad (const typename Triangulation< 2, 3 >::quad_iterator &quad,
+                                std::vector< Point< 3 > > & points )const
+{
+                                  //Almost the same implementation as  StraightBoundary<2,3>
+  const unsigned int n=points.size(),
+                    m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+                                  // is n a square number
+  Assert(m*m==n, ExcInternalError());
+
+  const double ds=1./(m+1);
+  double y=ds;
+
+  Point<2>  p[4];
+
+  for(unsigned  i=0;i<4;i++)
+    p[i]=get_surf_coord(quad->vertex(i));
+
+  Point<2>  target;
+  unsigned  offset[2];
+  offset[0]=0;
+  offset[1]=0;
+
+                                  //Take care for periodic conditions & negative angles,
+                                  //see get_new_point_on_line() above
+                                  //Because we dont have a symmetric interpolation (just the middle) we need to
+                                  //add 2*Pi to each almost zero and negative angles.
+  for(unsigned  i=0;i<2;i++)
+    for(unsigned  j=1;j<4;j++){
+      if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+       offset[i]++;
+       break;
+      }
+    }
+
+  for(unsigned  i=0;i<2;i++)
+    for(unsigned  j=0;j<4;j++)
+      if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+       p[j](i)+=2*numbers::PI*offset[i];
+
+  for (unsigned  i=0; i<m; ++i, y+=ds)
+    {
+      double x=ds;
+      for (unsigned  j=0; j<m; ++j, x+=ds){
+       target=((1-x) * p[0] +
+               x     * p[1]) * (1-y) +
+              ((1-x) * p[2] +
+               x     * p[3]) * y;
+
+       points[i*m+j]=get_real_coord(target);
+      }
+    }
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_normals_at_vertices (const typename Triangulation<2,3 >::face_iterator &face,
+                        FaceVertexNormals &face_vertex_normals) const
+{
+  for(unsigned i=0; i<GeometryInfo<2>::vertices_per_face; i++)
+    face_vertex_normals[i]=get_surf_norm(face->vertex(i));
+}
+
+
+
 // explicit instantiations
 #include "tria_boundary_lib.inst"
 

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.