]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Introduce new flags face_flip and face_rotation for cells in 3D and line_orientation...
authorleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Feb 2007 06:59:30 +0000 (06:59 +0000)
committerleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Feb 2007 06:59:30 +0000 (06:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@14444 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria_accessor.h
deal.II/deal.II/include/grid/tria_accessor.templates.h
deal.II/deal.II/include/grid/tria_faces.h
deal.II/deal.II/include/grid/tria_objects.h
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/grid/tria_objects.cc

index da1a4559ae1e66da644275076dc62494059354fb..c40c304a32c62358784f6758997a1a98d04c97a3 100644 (file)
@@ -788,13 +788,13 @@ class TriaObjectAccessor :  public TriaAccessor<dim>
                                       * (@p false). Which is the
                                       * standard direction is
                                       * documented with the
-                                      * Triangulation class. In
+                                      * GeometryInfo class. In
                                       * 1d and 2d, this is always
                                       * @p true, but in 3d it may be
                                       * different, see the respective
                                       * discussion in the
                                       * documentation of the
-                                      * Triangulation class.
+                                      * GeometryInfo class.
                                       *
                                       * This function is really only
                                       * for internal use in the
@@ -802,7 +802,124 @@ class TriaObjectAccessor :  public TriaAccessor<dim>
                                       * know what this is all about.
                                       */
     bool face_orientation (const unsigned int face) const;
+
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 180 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_flip (const unsigned int face) const;
+
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 90 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_rotation (const unsigned int face) const;
+    
+                                     /**
+                                      * Return whether the line with index @p
+                                      * line is oriented in standard
+                                      * direction. @p true indicates, that the
+                                      * line is oriented from vertex 0 to vertex
+                                      * 1, whereas it is the other way around
+                                      * otherwise. In 1d and 2d, this is always
+                                      * @p true, but in 3d it may be different,
+                                      * see the respective discussion in the
+                                      * documentation of the 
+                                      * GeometryInfo classe.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool line_orientation (const unsigned int line) const;
     
+                                     /**
+                                      * Set whether the quad with
+                                      * index @p face has its normal
+                                      * pointing in the standard
+                                      * direction (@p true) or
+                                      * whether it is the opposite
+                                      * (@p false). Which is the
+                                      * standard direction is
+                                      * documented with the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is only for
+                                      * internal use in the
+                                      * library. Setting this flag to
+                                      * any other value than the one
+                                      * that the triangulation has
+                                      * already set is bound to bring
+                                      * you desaster.
+                                     *
+                                     * It is only possible to set the
+                                     * face_orientation of cells in 3d
+                                     * (i.e. <code>celldim==3 &&
+                                     * dim==3</code>).
+                                      */
+    void set_face_orientation (const unsigned int face,
+                               const bool         orientation) const;
+
+                                    /**
+                                      * Set the flag indicating, what
+                                      * <code>face_flip()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * face_flip of cells in 3d
+                                     * (i.e. <code>celldim==3 &&
+                                     * dim==3</code>).
+                                      */
+    void set_face_flip (const unsigned int face,
+                       const bool         flip) const;
+
+                                    /**
+                                      * Set the flag indicating, what
+                                      * <code>face_rotation()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * face_rotation of cells in 3d
+                                     * (i.e. <code>celldim==3 &&
+                                     * dim==3</code>).
+                                      */
+    void set_face_rotation (const unsigned int face,
+                           const bool         rotation) const;
+
+                                     /**
+                                      * Set the flag indicating, what
+                                      * <code>line_orientation()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * line_orientation of faces in 3d
+                                     * (i.e. <code>celldim==2 &&
+                                     * dim==3</code>).
+                                      */
+    void set_line_orientation (const unsigned int line,
+                               const bool         orientation) const;
+
   private:
                                     /**
                                      *  Copy operator. This is normally
@@ -1307,13 +1424,13 @@ class TriaObjectAccessor<1, dim> :  public TriaAccessor<dim>
                                       * (@p false). Which is the
                                       * standard direction is
                                       * documented with the
-                                      * Triangulation class. In
+                                      * GeometryInfo class. In
                                       * 1d and 2d, this is always
                                       * @p true, but in 3d it may be
                                       * different, see the respective
                                       * discussion in the
                                       * documentation of the
-                                      * Triangulation class.
+                                      * GeometryInfo class.
                                       *
                                       * This function is really only
                                       * for internal use in the
@@ -1322,6 +1439,70 @@ class TriaObjectAccessor<1, dim> :  public TriaAccessor<dim>
                                       */
     bool face_orientation (const unsigned int face) const;
 
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 180 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_flip (const unsigned int face) const;
+
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 90 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_rotation (const unsigned int face) const;
+    
+                                     /**
+                                      * Return whether the line with index @p
+                                      * line is oriented in standard
+                                      * direction. @p true indicates, that the
+                                      * line is oriented from vertex 0 to vertex
+                                      * 1, whereas it is the other way around
+                                      * otherwise. In 1d and 2d, this is always
+                                      * @p true, but in 3d it may be different,
+                                      * see the respective discussion in the
+                                      * documentation of the 
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool line_orientation (const unsigned int line) const;
+    
+                                     /**
+                                      * Set the flag indicating, what
+                                      * <code>line_orientation()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * line_orientation of faces in 3d
+                                     * (i.e. <code>celldim==2 &&
+                                     * dim==3</code>).
+                                      */
+    void set_line_orientation (const unsigned int line,
+                               const bool         orientation) const;
+    
   private:
                                     /**
                                      *  Copy operator. This is normally
@@ -1857,13 +2038,13 @@ class TriaObjectAccessor<2, dim> :  public TriaAccessor<dim>
                                       * (@p false). Which is the
                                       * standard direction is
                                       * documented with the
-                                      * Triangulation class. In
+                                      * GeometryInfo class. In
                                       * 1d and 2d, this is always
                                       * @p true, but in 3d it may be
                                       * different, see the respective
                                       * discussion in the
                                       * documentation of the
-                                      * Triangulation class.
+                                      * GeometryInfo class.
                                       *
                                       * This function is really only
                                       * for internal use in the
@@ -1872,6 +2053,70 @@ class TriaObjectAccessor<2, dim> :  public TriaAccessor<dim>
                                       */
     bool face_orientation (const unsigned int face) const;
 
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 180 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_flip (const unsigned int face) const;
+
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 90 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_rotation (const unsigned int face) const;
+    
+                                     /**
+                                      * Return whether the line with index @p
+                                      * line is oriented in standard
+                                      * direction. @p true indicates, that the
+                                      * line is oriented from vertex 0 to vertex
+                                      * 1, whereas it is the other way around
+                                      * otherwise. In 1d and 2d, this is always
+                                      * @p true, but in 3d it may be different,
+                                      * see the respective discussion in the
+                                      * documentation of the 
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool line_orientation (const unsigned int line) const;
+    
+                                     /**
+                                      * Set the flag indicating, what
+                                      * <code>line_orientation()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * line_orientation of faces in 3d
+                                     * (i.e. <code>celldim==2 &&
+                                     * dim==3</code>).
+                                      */
+    void set_line_orientation (const unsigned int line,
+                               const bool         orientation) const;
+
   private:
                                     /**
                                      * Copy operator. This is
@@ -2394,13 +2639,13 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       * (@p false). Which is the
                                       * standard direction is
                                       * documented with the
-                                      * Triangulation class.  In
+                                      * GeometryInfo class.  In
                                       * 1d and 2d, this is always
                                       * @p true, but in 3d it may be
                                       * different, see the respective
                                       * discussion in the
                                       * documentation of the
-                                      * Triangulation class.
+                                      * GeometryInfo class.
                                       *
                                       * This function is really only
                                       * for internal use in the
@@ -2408,7 +2653,57 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       * know what this is all about.
                                       */
     bool face_orientation (const unsigned int face) const;
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 180 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_flip (const unsigned int face) const;
 
+                                     /**
+                                      * Return whether the face with index @p
+                                      * face is rotated by 90 degrees (@p true)
+                                      * or or not (@p false). In 1d and 2d, this
+                                      * is always @p false, but in 3d it may be
+                                      * different, see the respective discussion
+                                      * in the documentation of the
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool face_rotation (const unsigned int face) const;
+    
+                                     /**
+                                      * Return whether the line with index @p
+                                      * line is oriented in standard
+                                      * direction. @p true indicates, that the
+                                      * line is oriented from vertex 0 to vertex
+                                      * 1, whereas it is the other way around
+                                      * otherwise. In 1d and 2d, this is always
+                                      * @p true, but in 3d it may be different,
+                                      * see the respective discussion in the
+                                      * documentation of the 
+                                      * GeometryInfo class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
+                                      */
+    bool line_orientation (const unsigned int line) const;
+    
                                      /**
                                       * Set whether the quad with
                                       * index @p face has its normal
@@ -2418,7 +2713,7 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       * (@p false). Which is the
                                       * standard direction is
                                       * documented with the
-                                      * Triangulation class.
+                                      * GeometryInfo class.
                                       *
                                       * This function is only for
                                       * internal use in the
@@ -2430,6 +2725,32 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       */
     void set_face_orientation (const unsigned int face,
                                const bool         orientation) const;
+
+                                    /**
+                                      * Set the flag indicating, what
+                                      * <code>face_flip()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * face_orientation of cells in 3d
+                                     * (i.e. <code>celldim==3 &&
+                                     * dim==3</code>).
+                                      */
+    void set_face_flip (const unsigned int face,
+                       const bool         flip) const;
+
+                                    /**
+                                      * Set the flag indicating, what
+                                      * <code>face_rotation()</code> will
+                                      * return.
+                                     *
+                                     * It is only possible to set the
+                                     * face_orientation of cells in 3d
+                                     * (i.e. <code>celldim==3 &&
+                                     * dim==3</code>).
+                                      */
+    void set_face_rotation (const unsigned int face,
+                           const bool         rotation) const;
     
   private:
                                     /**
@@ -2747,15 +3068,13 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                       * for the index of the
                                       * child.
                                       *
-                                      * However, the function
-                                      * is more complicated in 3d,
-                                      * since there faces may have
-                                      * more than one orientation, and
-                                      * we have to use
-                                      * @p face_orientation for both
-                                      * this and the neighbor cell to
-                                      * figure out which cell we want
-                                      * to have.
+                                      * However, the function is more
+                                      * complicated in 3d, since there faces may
+                                      * have more than one orientation, and we
+                                      * have to use @p face_orientation, @p
+                                      * face_flip and @p face_rotation for both
+                                      * this and the neighbor cell to figure out
+                                      * which cell we want to have.
                                       *
                                       * This can lead to surprising
                                       * results: if we are sitting on
@@ -2787,6 +3106,12 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                       * orderings are only identical
                                       * if the face orientation is @p
                                       * true, and reversed otherwise.
+                                     *
+                                     * Similarly, effects of
+                                     * <tt>face_flip()==true</tt> and
+                                     * <tt>face_rotation()==true()</tt>, both
+                                     * of which indicate a non-standard face
+                                     * have to be considered.
                                       *
                                       * Fortunately, this is only very rarely
                                       * of concern, since one is usually only
index a34928c5eb18dca47a4fec9fda7de71a7ca3afd1..ca3258bd5b62b642f405e75371ad948de37688df 100644 (file)
@@ -302,6 +302,33 @@ TriaObjectAccessor<1,dim>::face_orientation (const unsigned int) const
 }
 
 
+template <int dim>
+inline
+bool
+TriaObjectAccessor<1,dim>::face_flip (const unsigned int) const
+{
+  return false;
+}
+
+
+template <int dim>
+inline
+bool
+TriaObjectAccessor<1,dim>::face_rotation (const unsigned int) const
+{
+  return false;
+}
+
+
+template <int dim>
+inline
+bool
+TriaObjectAccessor<1,dim>::line_orientation (const unsigned int) const
+{
+  return true;
+}
+
+
 
 template <int dim>
 inline
@@ -594,6 +621,53 @@ TriaObjectAccessor<2,dim>::face_orientation (const unsigned int) const
 }
 
 
+template <int dim>
+inline
+bool
+TriaObjectAccessor<2,dim>::face_flip (const unsigned int) const
+{
+  return false;
+}
+
+
+template <int dim>
+inline
+bool
+TriaObjectAccessor<2,dim>::face_rotation (const unsigned int) const
+{
+  return false;
+}
+
+
+template <>
+inline
+bool
+TriaObjectAccessor<2,3>::line_orientation (const unsigned int line) const
+{
+                                  // we cannot usae the quads() function here,
+                                  // since it returns a reference to
+                                  // TriaObjects<Quad>, but we need a
+                                  // (reference to) TriaObjectsQuad3D
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (line<GeometryInfo<2>::lines_per_cell,
+          ExcIndexRange (line, 0, GeometryInfo<2>::lines_per_cell));
+      Assert (this->present_index * GeometryInfo<2>::lines_per_cell + line
+             < this->tria->faces->quads.line_orientations.size(),
+             ExcInternalError());
+  return this->tria->faces->quads.
+    line_orientations[this->present_index * GeometryInfo<2>::lines_per_cell + line];
+}
+
+
+template <int dim>
+inline
+bool
+TriaObjectAccessor<2,dim>::line_orientation (const unsigned int) const
+{
+  return true;
+}
+
 
 template <int dim>
 inline
@@ -785,35 +859,43 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const
   Assert (i<12, ExcIndexRange(i,0,12));
 
                                    // get the line index by asking the
-                                   // quads. make sure we handle
-                                   // reverted faces correctly
+                                   // quads. first assume standard orientation
                                    //
-                                   // so set up a table that for each
+                                   // set up a table that for each
                                    // line describes a) from which
                                    // quad to take it, b) which line
                                    // therein it is if the face is
-                                   // oriented correctly, and c) if in
-                                   // the opposite direction
-  static const unsigned int lookup_table[12][3] =
-    { { 4, 0, 2 }, // take first four lines from bottom face
-      { 4, 1, 3 },
-      { 4, 2, 0 },
-      { 4, 3, 1 },
-
-      { 5, 0, 2 }, // second four lines from top face
-      { 5, 1, 3 },
-      { 5, 2, 0 },
-      { 5, 3, 1 },
-
-      { 0, 0, 2 }, // the rest randomly
-      { 1, 0, 2 },
-      { 0, 1, 3 },
-      { 1, 1, 3 }};
+                                   // oriented correctly
+  static const unsigned int lookup_table[12][2] =
+    { { 4, 0 }, // take first four lines from bottom face
+      { 4, 1 },
+      { 4, 2 },
+      { 4, 3 },
+
+      { 5, 0 }, // second four lines from top face
+      { 5, 1 },
+      { 5, 2 },
+      { 5, 3 },
+
+      { 0, 0 }, // the rest randomly
+      { 1, 0 },
+      { 0, 1 },
+      { 1, 1 }};
+
+                                  // respect non-standard faces by calling the
+                                  // reordering function from GeometryInfo
+
+  const unsigned int quad_index=lookup_table[i][0];
+  const unsigned int std_line_index=lookup_table[i][1];
+  
+  const unsigned int line_index=GeometryInfo<dim>::standard_to_real_face_line(
+    std_line_index,
+    face_orientation(quad_index),
+    face_flip(quad_index),
+    face_rotation(quad_index));
 
-  return (this->quad(lookup_table[i][0])
-          ->line_index(face_orientation(lookup_table[i][0]) ?
-                       lookup_table[i][1] :
-                       lookup_table[i][2]));
+  return (this->quad(quad_index)
+          ->line_index(line_index));
 }
 
 
@@ -916,15 +998,147 @@ face_orientation (const unsigned int face) const
   Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
   Assert (face<GeometryInfo<3>::faces_per_cell,
           ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
-      Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
-             < this->tria->levels[this->present_level]
-             ->hexes.face_orientations.size(),
-             ExcInternalError());
+  Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+         < this->tria->levels[this->present_level]
+         ->hexes.face_orientations.size(),
+         ExcInternalError());
+  
+  return (this->tria->levels[this->present_level]
+         ->hexes.face_orientations[this->present_index *
+                                   GeometryInfo<3>::faces_per_cell
+                                   + face]);
+}
+
+
+template <>
+inline
+bool
+TriaObjectAccessor<3, 3>::
+face_flip (const unsigned int face) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (face<GeometryInfo<3>::faces_per_cell,
+          ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
+  Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+         < this->tria->levels[this->present_level]
+         ->hexes.face_flips.size(),
+         ExcInternalError());
+  
+  return (this->tria->levels[this->present_level]
+         ->hexes.face_flips[this->present_index *
+                            GeometryInfo<3>::faces_per_cell
+                            + face]);
+}
+
+
+template <>
+inline
+bool
+TriaObjectAccessor<3, 3>::
+face_rotation (const unsigned int face) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (face<GeometryInfo<3>::faces_per_cell,
+          ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
+  Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+         < this->tria->levels[this->present_level]
+         ->hexes.face_rotations.size(),
+         ExcInternalError());
       
-      return (this->tria->levels[this->present_level]
-             ->hexes.face_orientations[this->present_index *
-                                       GeometryInfo<3>::faces_per_cell
-                                       + face]);
+  return (this->tria->levels[this->present_level]
+         ->hexes.face_rotations[this->present_index *
+                                GeometryInfo<3>::faces_per_cell
+                                + face]);
+}
+
+
+template <>
+inline
+bool
+TriaObjectAccessor<3, 3>::
+line_orientation (const unsigned int line) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (line<GeometryInfo<3>::lines_per_cell,
+          ExcIndexRange (line, 0, GeometryInfo<3>::lines_per_cell));
+
+                                   // get the line index by asking the
+                                   // quads. first assume standard orientation
+                                   //
+                                   // set up a table that for each
+                                   // line describes a) from which
+                                   // quad to take it, b) which line
+                                   // therein it is if the face is
+                                   // oriented correctly
+  static const unsigned int lookup_table[12][2] =
+    { { 4, 0 }, // take first four lines from bottom face
+      { 4, 1 },
+      { 4, 2 },
+      { 4, 3 },
+
+      { 5, 0 }, // second four lines from top face
+      { 5, 1 },
+      { 5, 2 },
+      { 5, 3 },
+
+      { 0, 0 }, // the rest randomly
+      { 1, 0 },
+      { 0, 1 },
+      { 1, 1 }};
+
+  const unsigned int quad_index=lookup_table[line][0];
+  const unsigned int std_line_index=lookup_table[line][1];
+  
+  const unsigned int line_index=GeometryInfo<dim>::standard_to_real_face_line(
+    std_line_index,
+    face_orientation(quad_index),
+    face_flip(quad_index),
+    face_rotation(quad_index));
+  
+                                  // now we got to the correct line and ask
+                                  // the quad for its line_orientation. however, if
+                                  // the face is rotated, it might be possible,
+                                  // that a standard orientation of the line
+                                  // with respect to the face corrsponds to a
+                                  // non-standard orientation for the line with
+                                  // respect to the cell.
+                                  //
+                                  // set up a table indicating if the two
+                                  // standard orientations coincide
+                                  //
+                                  // first index: two pairs of lines 0(lines
+                                  // 0/1) and 1(lines 2/3)
+                                  //
+                                  // second index: face_orientation; 0:
+                                  // opposite normal, 1: standard
+                                  //
+                                  // third index: face_flip; 0: standard, 1:
+                                  // face rotated by 180 degrees
+                                  //
+                                  // forth index: face_rotation: 0: standard,
+                                  // 1: face rotated by 90 degrees
+
+  static const bool bool_table[2][2][2][2] =
+    { { { { true, false },   // lines 0/1, face_orientation=false, face_flip=false, face_rotation=false and true
+         { false, true }},  // lines 0/1, face_orientation=false, face_flip=true, face_rotation=false and true
+       { { true, true },    // lines 0/1, face_orientation=true, face_flip=false, face_rotation=false and true
+         { false, false }}},// linea 0/1, face_orientation=true, face_flip=true, face_rotation=false and true
+
+      { { { true, true },    // lines 2/3 ...
+         { false, false }},
+       { { true, false },
+         { false, true }}}};
+  
+  
+  return (this->quad(quad_index)
+          ->line_orientation(line_index)
+         == bool_table[std_line_index/2]
+         [face_orientation(quad_index)]
+         [face_flip(quad_index)]
+         [face_rotation(quad_index)]);
 }
 
 
index a095481006e3d7e26880043e0315a9d577494aa8..8981bdbb098301bf9fbce3e4c156cc0a8d8158bd 100644 (file)
@@ -110,7 +110,7 @@ namespace internal
                                          *  The TriaObject containing the data of quads.
                                          */
 
-       TriaObjects<Quad> quads;
+       TriaObjectsQuad3D quads;
 
                                         /**
                                          *  The TriaObject containing the data of lines.
index f3065e94120755d93c7ff3a467b760ffd25dbb19..ce6ca946b5f2ee8c59206f9fa915c7d803ac727d 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+// TODO: See if we can unify the class hierarchy a bit using partial
+// specialization of classes here, e.g. declare a general class
+// TriaObjects<dim,G>.
+// To consider for this: in that case we would have to duplicate quite a few
+// things, e.g. TriaObjectsHex is derived from TriaObjects<Hexahedron> and
+// declares mainly additional data. This would have to be changed in case of a
+// partial specialization.
+
 namespace internal
 {
   namespace Triangulation
@@ -236,6 +244,14 @@ namespace internal
                                          * has.
                                          */
        std::vector<bool> face_orientations;
+                                        /**
+                                         * flip = rotation by 180 degrees
+                                         */
+       std::vector<bool> face_flips;
+                                        /**
+                                         * rotation by 90 degrees
+                                         */
+       std::vector<bool> face_rotations;
 
                                          /**
                                           *  Assert that enough space is
@@ -270,6 +286,63 @@ namespace internal
                                           */
         unsigned int memory_consumption () const;          
     };
+
+
+/**
+ * For quadrilaterals in 3D the data of TriaObjects needs to be extended, as we
+ * can obtain faces (quads) with lines in non-standard-orientation, therefore we
+ * declare a class TriaObjectsQuad3D, which additionaly contains a bool-vector
+ * of the line-orientations.
+ * @ingroup grid
+ */
+    
+    class TriaObjectsQuad3D: public TriaObjects<Quad>
+    {
+      public:
+
+                                        /**
+                                         * In effect, this field has
+                                         * <tt>4*n_quads</tt> elements,
+                                         * being the number of quads
+                                         * times the four lines each
+                                         * has.
+                                         */
+       std::vector<bool> line_orientations;
+
+                                         /**
+                                          *  Assert that enough space is
+                                          *  allocated to accomodate
+                                          *  <tt>new_objs</tt> new objects.
+                                          *  This function does not only call
+                                          *  <tt>vector::reserve()</tt>, but
+                                          *  does really append the needed
+                                          *  elements.
+                                          */
+        void reserve_space (const unsigned int new_objs);
+
+                                        /**
+                                         *  Clear all the data contained in this object.
+                                         */
+       void clear();
+       
+                                         /**
+                                          *  Check the memory consistency of the
+                                          *  different containers. Should only be
+                                          *  called with the prepro flag @p DEBUG
+                                          *  set. The function should be called from
+                                          *  the functions of the higher
+                                          *  TriaLevel classes.
+                                          */
+        void monitor_memory (const unsigned int true_dimension) const;
+
+                                         /**
+                                          * Determine an estimate for the
+                                          * memory consumption (in bytes)
+                                          * of this object.
+                                          */
+        unsigned int memory_consumption () const;          
+    };
+    
   }
 }
 
index 92605d8832a56bad0e353cb433967e5013b421bc..a554bb5ce8a0fb55e397d16085d0829741167fcc 100644 (file)
@@ -282,7 +282,18 @@ int TriaObjectAccessor<2, dim>::vertex_index (const unsigned int corner) const
 {
   Assert (corner<4, ExcIndexRange(corner,0,4));
 
-  return line(corner%2)->vertex_index(corner/2);
+                                  // table used to switch the vertices, if the
+                                  // line orientation is wrong,
+                                  //
+                                  // first index: line orientation 0: false or
+                                  // 1: true=standard
+                                  //
+                                  // second index: vertex index to be switched
+                                  // (or not)
+  
+  static const unsigned int switch_table[2][2]={{1,0},{0,1}};
+  
+  return line(corner%2)->vertex_index(switch_table[line_orientation(corner%2)][corner/2]);
 }
 
 
@@ -677,6 +688,41 @@ unsigned int TriaObjectAccessor<2, dim>::number_of_children () const
 }
 
 
+template <int dim>
+void
+TriaObjectAccessor<2, dim>::
+set_line_orientation (const unsigned int ,
+                      const bool         ) const
+{
+  Assert(false,
+        ExcMessage("The line orientation of Quads can only be set in 3D. In 2D the lines are alway in standard orientation."));}
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void
+TriaObjectAccessor<2, 3>::
+set_line_orientation (const unsigned int line,
+                      const bool         orientation) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (line<GeometryInfo<dim>::lines_per_face,
+          ExcIndexRange (line, 0, GeometryInfo<3>::lines_per_face));
+  Assert (this->present_index * GeometryInfo<dim>::lines_per_face + line
+          < this->tria->faces->quads.line_orientations.size(),
+          ExcInternalError());
+          
+  this->tria->faces
+    ->quads.line_orientations[this->present_index *
+                              GeometryInfo<dim>::lines_per_face
+                              +
+                              line]  = orientation;
+}
+
+#endif
 
 
 /*------------------------ Functions: TriaObjectAccessor ---------------------------*/
@@ -704,30 +750,21 @@ int TriaObjectAccessor<3, dim>::vertex_index (const unsigned int corner) const
 {
   Assert (corner<8, ExcIndexRange(corner,0,8));
 
-                                  // get the corner indices by asking
-                                  // either the bottom or the top
-                                  // face for its vertices. make sure
-                                  // we take into account that the
-                                  // face might have non-standard
-                                  // orientation; if this is not the
-                                  // case, transpose vertices 1 and 2
-                                  // of this face
-  static const unsigned int vertex_translation[4] = 
-    { 0, 2, 1, 3 };
-  if (corner<4)
-    {
-      if (face_orientation(4) == true)
-        return quad(4)->vertex_index(corner);
-      else
-        return quad(4)->vertex_index(vertex_translation[corner]);
-    }
-  else
-    {
-      if (face_orientation(5) == true)
-        return quad(5)->vertex_index(corner-4);
-      else
-        return quad(5)->vertex_index(vertex_translation[corner-4]);
-    }
+                                  // get the corner indices by asking either
+                                  // the bottom or the top face for its
+                                  // vertices. handle non-standard faces by
+                                  // calling the vertex reordering function
+                                  // from GeometryInfo
+
+                                  // bottom face (4) for first four vertices,
+                                  // top face (5) for the rest
+  const unsigned int face_index=4+corner/4;
+
+  return quad(face_index)->vertex_index(GeometryInfo<dim>
+                                       ::standard_to_real_face_vertex(corner%4,
+                                                                      face_orientation(face_index),
+                                                                      face_flip(face_index),
+                                                                      face_rotation(face_index)));
 }
 
 
@@ -1649,6 +1686,52 @@ set_face_orientation (const unsigned int face,
                               face]  = orientation;
 }
 
+
+template <>
+void
+TriaObjectAccessor<3, 3>::
+set_face_flip (const unsigned int face,
+              const bool         flip) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (face<GeometryInfo<3>::faces_per_cell,
+          ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
+  Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+          < this->tria->levels[this->present_level]
+          ->hexes.face_flips.size(),
+          ExcInternalError());
+          
+  this->tria->levels[this->present_level]
+    ->hexes.face_flips[this->present_index *
+                      GeometryInfo<3>::faces_per_cell
+                      +
+                      face]  = flip;
+}
+
+
+template <>
+void
+TriaObjectAccessor<3, 3>::
+set_face_rotation (const unsigned int face,
+                  const bool         rotation) const
+{
+  const int dim=3;
+  Assert (used(), TriaAccessor<dim>::ExcCellNotUsed());
+  Assert (face<GeometryInfo<3>::faces_per_cell,
+          ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell));
+  Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face
+          < this->tria->levels[this->present_level]
+          ->hexes.face_rotations.size(),
+          ExcInternalError());
+          
+  this->tria->levels[this->present_level]
+    ->hexes.face_rotations[this->present_index *
+                          GeometryInfo<3>::faces_per_cell
+                          +
+                          face]  = rotation;
+}
+
 #endif
 // Remark: The following explicit instantiations needed to be moved to
 // this place here to work around a problem with gcc3.3 on Apple MacOSX.
@@ -2183,62 +2266,51 @@ neighbor_child_on_subface (const unsigned int face,
           ExcMessage ("The neighbor must be on the same level as this cell!"));
   Assert (this->neighbor(face)->has_children() == true,
           ExcMessage ("The neighbor must have children!"));
-  
-  static const unsigned int subface_translation[4]
-    = { 0, 2, 1, 3 };
-                                   // see whether face and
-                                   // the neighbor's
-                                   // counterface share the
-                                   // same indexing of
-                                   // children. if not so,
-                                   // translate child
-                                   // indices
+                                  // this function returns they neighbor's
+                                  // child on a given face and
+                                  // subface. However, irrespective of the
+                                  // actual face orientation, flip and rotation
+                                  // of the cells face it returns that child,
+                                  // which is neighboring the subface which
+                                  // would be subface number @p subface, if the
+                                  // face was in standard orientation. thus,
+                                  // our first task is to extract the
+                                  // information, which real_subface_no
+                                  // corresponds to that for the actual face,
+                                  // which might have non-standard orientation,
+                                  // flip and rotation.
+
+                                  // to this end: the subfaces are numbered
+                                  // according to the vertices of the
+                                  // face. thus we can use the
+                                  // standard_to_real_face_vertex function in
+                                  // GeometryInfo<dim> to translate the given
+                                  // subface to the real_subface_no
+  unsigned int real_subface_no=GeometryInfo<3>::
+                              standard_to_real_face_vertex(subface,
+                                                           this->face_orientation(face),
+                                                           this->face_flip(face),
+                                                           this->face_rotation(face));
+                                  // get the neighbor's number for the given
+                                  // face and the neighbor
   const unsigned int neighbor_neighbor
     = this->neighbor_of_neighbor (face);
-  const bool face_orientations_match
-    = (this->neighbor(face)->face_orientation(neighbor_neighbor) ==
-       this->face_orientation(face));
+  const TriaIterator<3,CellAccessor<3> > neighbor=this->neighbor(face);
+  
+                                  // now use the info provided by GeometryInfo
+                                  // to extract the neighbors child number
   const unsigned int neighbor_child_index
-    = GeometryInfo<3>::child_cell_on_face(neighbor_neighbor, subface,
-                                         face_orientations_match);
+    = GeometryInfo<3>::child_cell_on_face(neighbor_neighbor, real_subface_no,
+                                         neighbor->face_orientation(neighbor_neighbor),
+                                         neighbor->face_flip(neighbor_neighbor),
+                                         neighbor->face_rotation(neighbor_neighbor));
   const TriaIterator<3,CellAccessor<3> > neighbor_child=
-    this->neighbor(face)->child(neighbor_child_index);
+    neighbor->child(neighbor_child_index);
 
-                                   // if the face on the side of the present
-                                   // cell is in the correct order, then make
-                                   // sure that the neighbor child cell we
+                                  // make sure that the neighbor child cell we
                                    // have found shares the desired subface.
-                                   //
-                                   // otherwise: if the face is turned
-                                   // inside out when viewed from the
-                                   // present cell, then the subface we are
-                                   // interested in is not
-                                   // this->face(face)->child(subface), but
-                                   // instead
-                                   // this->face(face)->child(subface'),
-                                   // where subface' is the translated
-                                   // subface number. this is so because we
-                                   // only store the face only once, so its
-                                   // children are defined in terms of its
-                                   // own (circular) orientation, not in
-                                   // terms of face_orientation as viewed
-                                   // from one of the adjacent cells. in
-                                   // that case, cell->face(f)->subface(sf)
-                                   // may yield unexpected results; in that
-                                   // case, a caller may need to adjust
-                                   // according to the face_orientation
-                                   // flag, though one in general only wants
-                                   // to loop over all subfaces, and not
-                                   // pick a particular one
-  Assert(((this->face_orientation(face) == true)
-          &&
-          (this->face(face)->child(subface) ==
-           neighbor_child->face(neighbor_neighbor)))
-         ||
-         ((this->face_orientation(face) == false)
-          &&
-          (this->face(face)->child(subface_translation[subface]) ==
-           neighbor_child->face(neighbor_neighbor))),
+  Assert((this->face(face)->child(real_subface_no) ==
+         neighbor_child->face(neighbor_neighbor)),
          ExcInternalError());
 
   return neighbor_child;
index bed78dfc3a75e5d2a9043cd2ac58c2646be0e4d9..07b1f2c1b797f9abaa5777d7128d31cef3954d01 100644 (file)
@@ -162,10 +162,41 @@ namespace internal
                                    new_size * GeometryInfo<3>::faces_per_cell
                                    - face_orientations.size(),
                                    true);
+          face_flips.reserve (new_size * GeometryInfo<3>::faces_per_cell);
+          face_flips.insert (face_flips.end(),
+                                   new_size * GeometryInfo<3>::faces_per_cell
+                                   - face_flips.size(),
+                                   false);
+          face_rotations.reserve (new_size * GeometryInfo<3>::faces_per_cell);
+          face_rotations.insert (face_rotations.end(),
+                                   new_size * GeometryInfo<3>::faces_per_cell
+                                   - face_rotations.size(),
+                                   false);
         };
     }
 
 
+    void
+    TriaObjectsQuad3D::reserve_space (const unsigned int new_quads)
+    {
+      const unsigned int new_size = new_quads +
+                                    std::count_if (used.begin(),
+                                                   used.end(),
+                                                   std::bind2nd (std::equal_to<bool>(), true));
+
+                                       // see above...
+      if (new_size>cells.size())
+        {
+          TriaObjects<Quad>::reserve_space(new_quads);
+          line_orientations.reserve (new_size * GeometryInfo<2>::lines_per_cell);
+          line_orientations.insert (line_orientations.end(),
+                                   new_size * GeometryInfo<2>::lines_per_cell
+                                   - line_orientations.size(),
+                                   true);
+        };
+    }
+    
+
     template<>
     void
     TriaObjects<Line>::monitor_memory (const unsigned int) const
@@ -286,6 +317,33 @@ namespace internal
               == face_orientations.size(),
               ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell,
                                 face_orientations.size()));
+      Assert (cells.size() * GeometryInfo<3>::faces_per_cell
+              == face_flips.size(),
+              ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell,
+                                face_flips.size()));
+      Assert (cells.size() * GeometryInfo<3>::faces_per_cell
+              == face_rotations.size(),
+              ExcMemoryInexact (cells.size() * GeometryInfo<3>::faces_per_cell,
+                                face_rotations.size()));
+    }
+
+
+    void
+    TriaObjectsQuad3D::monitor_memory (const unsigned int) const
+    {
+                                       // check that we have not allocated
+                                       // too much memory. note that bool
+                                       // vectors allocate their memory in
+                                       // chunks of whole integers, so
+                                       // they may over-allocate by up to
+                                       // as many elements as an integer
+                                       // has bits
+      Assert (cells.size() * GeometryInfo<2>::lines_per_cell
+              == line_orientations.size(),
+              ExcMemoryInexact (cells.size() * GeometryInfo<2>::lines_per_cell,
+                                line_orientations.size()));
+      TriaObjects<Quad>::monitor_memory (3);
+      
     }
 
 
@@ -307,6 +365,16 @@ namespace internal
     {
       TriaObjects<Hexahedron>::clear();
       face_orientations.clear();
+      face_flips.clear();
+      face_rotations.clear();
+    }
+
+
+    void
+    TriaObjectsQuad3D::clear()
+    {
+      TriaObjects<Quad>::clear();
+      line_orientations.clear();
     }
     
     
@@ -327,7 +395,17 @@ namespace internal
     TriaObjectsHex::memory_consumption () const
     {
       return (MemoryConsumption::memory_consumption (face_orientations) +
-             this->TriaObjects<Hexahedron>::memory_consumption() );
+             MemoryConsumption::memory_consumption (face_flips) +
+             MemoryConsumption::memory_consumption (face_rotations) +
+             TriaObjects<Hexahedron>::memory_consumption() );
+    }
+
+
+    unsigned int
+    TriaObjectsQuad3D::memory_consumption () const
+    {
+      return (MemoryConsumption::memory_consumption (line_orientations) +
+             this->TriaObjects<Quad>::memory_consumption() );
     }
 
                                     // explicit instantiations

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.