]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix various things for 3d with misoriented faces.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Oct 2003 14:27:33 +0000 (14:27 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Oct 2003 14:27:33 +0000 (14:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@8095 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/source/fe/mapping_q1.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/numerics/error_estimator.cc

index 9249aec2ac273beb0dfb40478781b40e3469409a..38fa7318729a4d60a00886f984f5f6ec84fc9375 100644 (file)
@@ -676,20 +676,28 @@ class TriaObjectAccessor :  public TriaAccessor<dim>
     unsigned int number_of_children () const;
 
                                      /**
-                                      * Enquire some information about
-                                      * whether the face with the
-                                      * given number is in standard
-                                      * orientation or not. This
-                                      * information is only useful
-                                      * (and available) in 3d, see the
+                                      * Return whether the face 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
+                                      * @ref{Triangulation} 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
-                                      * respective class. For all
-                                      * other dimensions, this
-                                      * function should not be called
-                                      * and triggers and assertion if
-                                      * done so.
+                                      * @ref{Triangulation} class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
                                       */
-    bool get_face_orientation (const unsigned int face) const;
+    bool face_orientation (const unsigned int face) const;
     
   private:
                                     /**
@@ -1135,20 +1143,28 @@ class TriaObjectAccessor<1, dim> :  public TriaAccessor<dim>
     unsigned int number_of_children () const;
     
                                      /**
-                                      * Enquire some information about
-                                      * whether the face with the
-                                      * given number is in standard
-                                      * orientation or not. This
-                                      * information is only useful
-                                      * (and available) in 3d, see the
+                                      * Return whether the face 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
+                                      * @ref{Triangulation} 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
-                                      * respective class. For all
-                                      * other dimensions, this
-                                      * function should not be called
-                                      * and triggers and assertion if
-                                      * done so.
+                                      * @ref{Triangulation} class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
                                       */
-    bool get_face_orientation (const unsigned int face) const;
+    bool face_orientation (const unsigned int face) const;
 
   private:
                                     /**
@@ -1616,20 +1632,28 @@ class TriaObjectAccessor<2, dim> :  public TriaAccessor<dim>
     unsigned int number_of_children () const;
 
                                      /**
-                                      * Enquire some information about
-                                      * whether the face with the
-                                      * given number is in standard
-                                      * orientation or not. This
-                                      * information is only useful
-                                      * (and available) in 3d, see the
+                                      * Return whether the face 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
+                                      * @ref{Triangulation} 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
-                                      * respective class. For all
-                                      * other dimensions, this
-                                      * function should not be called
-                                      * and triggers and assertion if
-                                      * done so.
+                                      * @ref{Triangulation} class.
+                                      *
+                                      * This function is really only
+                                      * for internal use in the
+                                      * library unless you absolutely
+                                      * know what this is all about.
                                       */
-    bool get_face_orientation (const unsigned int face) const;
+    bool face_orientation (const unsigned int face) const;
 
   private:
                                     /**
@@ -2098,7 +2122,7 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
     unsigned int number_of_children () const;
 
                                      /**
-                                      * Return whether the quad with
+                                      * Return whether the face with
                                       * index @p{face} has its normal
                                       * pointing in the standard
                                       * direction (@p{true}) or
@@ -2106,6 +2130,12 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       * (@p{false}). Which is the
                                       * standard direction is
                                       * documented with the
+                                      * @ref{Triangulation} 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
                                       * @ref{Triangulation} class.
                                       *
                                       * This function is really only
@@ -2113,7 +2143,7 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
                                       * library unless you absolutely
                                       * know what this is all about.
                                       */
-    bool get_face_orientation (const unsigned int face) const;
+    bool face_orientation (const unsigned int face) const;
 
                                      /**
                                       * Set whether the quad with
index 98e5833172a45b63258b18e5ee7cfea02d3cfb7f..7eb16fed5cb639e9559e3b263d9ac2bf8f1fa34b 100644 (file)
@@ -25,6 +25,8 @@
 #include <cmath>
 
 
+//TODO[WB]: TriaObjectAccessor<dim,N> should never be used since we have specializations. Can we remove the implementations of functions, or what are they there for?
+
 /*------------------------ Functions: TriaAccessor ---------------------------*/
 
 
@@ -213,6 +215,16 @@ TriaObjectAccessor<1,dim>::max_refinement_depth () const
 
 
 
+template <int dim>
+inline
+bool
+TriaObjectAccessor<1,dim>::face_orientation (const unsigned int) const
+{
+  return true;
+}
+
+
+
 template <int dim>
 inline
 void
@@ -400,6 +412,16 @@ TriaObjectAccessor<2,dim>::max_refinement_depth () const
 
 
 
+template <int dim>
+inline
+bool
+TriaObjectAccessor<2,dim>::face_orientation (const unsigned int) const
+{
+  return true;
+}
+
+
+
 template <int dim>
 inline
 void
@@ -534,7 +556,7 @@ TriaObjectAccessor<3,dim>::line (const unsigned int i) const
       { 4, 3, 0 }};
 
   return (this->quad(lookup_table[i][0])
-          ->line(get_face_orientation(lookup_table[i][0]) ?
+          ->line(face_orientation(lookup_table[i][0]) ?
                  lookup_table[i][1] :
                  lookup_table[i][2]));
 }
@@ -592,7 +614,7 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const
       { 4, 3, 0 }};
 
   return (this->quad(lookup_table[i][0])
-          ->line_index(get_face_orientation(lookup_table[i][0]) ?
+          ->line_index(face_orientation(lookup_table[i][0]) ?
                        lookup_table[i][1] :
                        lookup_table[i][2]));
 }
@@ -672,6 +694,29 @@ TriaObjectAccessor<3,dim>::max_refinement_depth () const
 
 
 
+template <int dim>
+inline
+bool
+TriaObjectAccessor<3, dim>::
+face_orientation (const unsigned int face) const
+{
+  Assert (used(), typename 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());
+          
+  return (this->tria->levels[this->present_level]
+          ->hexes.face_orientations[this->present_index *
+                                    GeometryInfo<3>::faces_per_cell
+                                    +
+                                    face]);
+}
+
+
+
 template <int dim>
 inline
 void
index 92fb98e3b304b01bc02cc69681190863b839be18..aaae80db024fa0dc0467490bac7613f14565febd 100644 (file)
@@ -62,7 +62,7 @@ namespace internal
                                                // face or subface
         case 3:
               return ((face_no +
-                       (cell->get_face_orientation(face_no) == true ?
+                       (cell->face_orientation(face_no) == true ?
                         0 : GeometryInfo<dim>::faces_per_cell))
                       * n_quadrature_points);
 
@@ -90,6 +90,11 @@ namespace internal
     Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
             ExcInternalError());
     Assert (n_quadrature_points > 0, ExcInternalError());
+
+    Assert (cell->has_children() == false,
+            ExcMessage ("You can't use subface data for cells that are "
+                        "already refined. Iterate over their children "
+                        "instead in these cases."));
     
     switch (dim)
       {
@@ -103,10 +108,11 @@ namespace internal
         case 3:
               return (((face_no * GeometryInfo<dim>::subfaces_per_face +
                         subface_no)
-                       + (cell->get_face_orientation(face_no) == true ?
+                       + (cell->face_orientation(face_no) == true ?
                           0 :
                           GeometryInfo<dim>::faces_per_cell *
-                          GeometryInfo<dim>::subfaces_per_face))
+                          GeometryInfo<dim>::subfaces_per_face)
+                       )
                       * n_quadrature_points);
         default:
               Assert (false, ExcInternalError());
index d3663b6e9b17d4816e26dc5bfe82bc5496bf5c85..9afb26c52e9711389ce768a193f8c3fecf047961 100644 (file)
@@ -6124,12 +6124,12 @@ Triangulation<3>::execute_refinement ()
                                              // opposite direction, so
                                              // store that up front
             const bool face_orientation[6]
-              = { hex->get_face_orientation (0),
-                  hex->get_face_orientation (1),
-                  hex->get_face_orientation (2),
-                  hex->get_face_orientation (3),
-                  hex->get_face_orientation (4),
-                  hex->get_face_orientation (5) };
+              = { hex->face_orientation (0),
+                  hex->face_orientation (1),
+                  hex->face_orientation (2),
+                  hex->face_orientation (3),
+                  hex->face_orientation (4),
+                  hex->face_orientation (5) };
                     
            const unsigned int line_indices[30]
              = {
@@ -6535,7 +6535,7 @@ Triangulation<3>::execute_refinement ()
             for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
               for (unsigned int s=0; s<GeometryInfo<dim>::subfaces_per_face; ++s)
                 new_hexes[GeometryInfo<dim>::child_cell_on_face(f,s)]
-                  ->set_face_orientation(f, hex->get_face_orientation(f));
+                  ->set_face_orientation(f, hex->face_orientation(f));
             
 
                                             /////////////////////////////////
@@ -6684,9 +6684,9 @@ Triangulation<3>::execute_refinement ()
                                                          // all
                                                          // wrong...)
                         const bool orient
-                          = (neighbor->get_face_orientation(nb_nb)
+                          = (neighbor->face_orientation(nb_nb)
                              &&
-                             hex->get_face_orientation(face));
+                             hex->face_orientation(face));
                         
                         static const unsigned int
                           child_switch_table[GeometryInfo<dim>::subfaces_per_face]
index fdd4a0e4360b9cf2c3c7ad8f0b7031634a33d433..0b69a294e8a0ef9a977b8cad93aa9dd5187f0c5d 100644 (file)
@@ -250,15 +250,6 @@ unsigned int TriaObjectAccessor<1, dim>::number_of_children () const
 
 
 
-template <int dim>
-bool
-TriaObjectAccessor<1,dim>::get_face_orientation (const unsigned int) const
-{
-  Assert (false, ExcInternalError());
-  return true;
-}
-
-
 /*------------------------ Functions: QuadAccessor ---------------------------*/
 
 template <int dim>
@@ -678,15 +669,6 @@ unsigned int TriaObjectAccessor<2, dim>::number_of_children () const
 
 
 
-template <int dim>
-bool
-TriaObjectAccessor<2,dim>::get_face_orientation (const unsigned int) const
-{
-  Assert (false, ExcInternalError());
-  return true;
-}
-
-
 
 /*------------------------ Functions: TriaObjectAccessor ---------------------------*/
 
@@ -715,14 +697,14 @@ int TriaObjectAccessor<3, dim>::vertex_index (const unsigned int corner) const
     { 0, 3, 2, 1 };
   if (corner<4)
     {
-      if (get_face_orientation(0) == true)
+      if (face_orientation(0) == true)
         return quad(0)->vertex_index(corner);
       else
         return quad(0)->vertex_index(vertex_translation[corner]);
     }
   else
     {
-      if (get_face_orientation(1) == true)
+      if (face_orientation(1) == true)
         return quad(1)->vertex_index(corner-4);
       else
         return quad(1)->vertex_index(vertex_translation[corner-4]);
@@ -1666,28 +1648,6 @@ unsigned int TriaObjectAccessor<3, dim>::number_of_children () const
 
 
 
-template <int dim>
-bool
-TriaObjectAccessor<3, dim>::
-get_face_orientation (const unsigned int face) const
-{
-  Assert (used(), typename 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());
-          
-  return (this->tria->levels[this->present_level]
-          ->hexes.face_orientations[this->present_index *
-                                    GeometryInfo<3>::faces_per_cell
-                                    +
-                                    face]);
-}
-
-
-
 template <int dim>
 void
 TriaObjectAccessor<3, dim>::
index a9c4bf9884997ecc2e6c3c681020895b03b59ac6..0e0d61db9a4f7a3f9efd4a3f5a3a28c03b14a69f 100644 (file)
@@ -373,25 +373,36 @@ estimate_some (const Mapping<dim>                  &mapping,
 {
   const unsigned int n_solution_vectors = solutions.size();
   
-                                  // make up a fe face values object for the
-                                  // restriction of the finite element function
-                                  // to a face, for the present cell and its
-                                  // neighbor. In principle we would only need
-                                  // one at a time, but this way we can
-                                  // have more fine grained access to what
-                                  // values really need to be computed (we
-                                  // need not compute all values on the
-                                  // neighbor cells, so using two objects
+                                  // make up a fe face values object
+                                  // for the restriction of the
+                                  // finite element function to a
+                                  // face, for the present cell and
+                                  // its neighbor. In principle we
+                                  // would only need one at a time,
+                                  // but this way we can have more
+                                  // fine grained access to what
+                                  // values really need to be
+                                  // computed (we need not compute
+                                  // all values on the neighbor
+                                  // cells, so using two objects
                                   // gives us a performance gain).
+                                   //
+                                   // in debug mode, make sure that
+                                   // some data matches, so compute
+                                   // quadrature points always
   FEFaceValues<dim> fe_face_values_cell (mapping,
                                         dof_handler.get_fe(),
                                         quadrature,
-                                        UpdateFlags(update_gradients      |
-                                                    update_JxW_values     |
-                                                    ((!neumann_bc.empty() ||
-                                                      (coefficients != 0))  ?
-                                                     update_q_points : 0) |
-                                                    update_normal_vectors)); 
+                                         UpdateFlags(
+                                           update_gradients      |
+                                           update_JxW_values     |
+                                           ((!neumann_bc.empty() ||
+                                             (coefficients != 0))  ?
+                                            update_q_points : 0) |
+#ifdef DEBUG
+                                           update_q_points |
+#endif
+                                           update_normal_vectors));
   FEFaceValues<dim> fe_face_values_neighbor (mapping,
                                             dof_handler.get_fe(),
                                             quadrature,
@@ -399,7 +410,12 @@ estimate_some (const Mapping<dim>                  &mapping,
   FESubfaceValues<dim> fe_subface_values (mapping,
                                          dof_handler.get_fe(),
                                          quadrature,
-                                         update_gradients);
+                                         UpdateFlags(
+                                            update_gradients
+#ifdef DEBUG
+                                            | update_q_points
+#endif
+                                            ));
 
 
   active_cell_iterator cell = dof_handler.begin_active();
@@ -987,43 +1003,57 @@ integrate_over_irregular_face (const DoFHandler<dim>               &dof_handler,
     {
                                       // get an iterator pointing to the
                                       // cell behind the present subface
-      const bool orientation_flag
-        = (dim != 3 ?
-           true :
-           cell->get_face_orientation(face_no) == true
-           &&
-           neighbor->get_face_orientation(neighbor_neighbor) == true);
       static const unsigned int subface_translation[4]
         = { 0, 3, 2, 1 };
+                                       // see whether face and
+                                       // the neighbor's
+                                       // counterface share the
+                                       // same indexing of
+                                       // children. if not so,
+                                       // translate child
+                                       // indices
+      const bool face_orientations_match
+        = (neighbor->face_orientation(neighbor_neighbor) ==
+           cell->face_orientation(face_no));
       const unsigned int neighbor_child_index
         = (GeometryInfo<dim>::
            child_cell_on_face(neighbor_neighbor,
-                              (orientation_flag ?
+                              (face_orientations_match ?
                                subface_no :
                                subface_translation[subface_no])));
       const active_cell_iterator neighbor_child
        = neighbor->child(neighbor_child_index);
-      Assert (neighbor_child->face(neighbor_neighbor) ==
-             cell->face(face_no)->child(subface_no),
-             ExcInternalError());
       Assert (!neighbor->child(neighbor_child_index)->has_children(),
              ExcInternalError());
       
-                                      // restrict the finite element on the
-                                      // present cell to the subface and
-                                      // store the gradient of the solution
-                                      // in psi
+                                      // restrict the finite element
+                                      // on the present cell to the
+                                      // subface
       fe_subface_values.reinit (cell, face_no, subface_no);
 
+                                      // restrict the finite element
+                                      // on the neighbor cell to the
+                                      // common @p{subface}.
+      fe_face_values.reinit (neighbor_child, neighbor_neighbor);
+
+                                       // make sure that quadrature
+                                       // points match
+      for (unsigned int q=0; q<n_q_points; ++q)
+        Assert ((fe_face_values.quadrature_point(q) -
+                 fe_subface_values.quadrature_point(q)).square()
+                <
+                1.e-15 * (fe_face_values.quadrature_point(q).square() + 
+                          fe_subface_values.quadrature_point(q).square()),
+                ExcInternalError());
+
+                                       // store the gradient of the
+                                      // solution in psi
       for (unsigned int n=0; n<n_solution_vectors; ++n)
        fe_subface_values.get_function_grads (*solutions[n], per_thread_data.psi[n]);
 
-                                      // restrict the finite element on the
-                                      // neighbor cell to the common @p{subface}.
-                                      // store the gradient in @p{neighbor_psi}
-      
-      fe_face_values.reinit (neighbor_child, neighbor_neighbor);
-
+                                       // store the gradient from the
+                                      // neighbor's side in
+                                      // @p{neighbor_psi}
       for (unsigned int n=0; n<n_solution_vectors; ++n)
        fe_face_values.get_function_grads (*solutions[n], per_thread_data.neighbor_psi[n]);
       
@@ -1070,7 +1100,7 @@ integrate_over_irregular_face (const DoFHandler<dim>               &dof_handler,
          if (coefficients->n_components == 1)
            {
              coefficients->value_list (fe_face_values.get_quadrature_points(),
-                                            per_thread_data.coefficient_values1);
+                                        per_thread_data.coefficient_values1);
              for (unsigned int n=0; n<n_solution_vectors; ++n)
                for (unsigned int component=0; component<n_components; ++component)
                  for (unsigned int point=0; point<n_q_points; ++point)

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.