]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Generalize jacobian_determinants_at_vertices to alternating_form_a_vertices.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 2009 04:39:12 +0000 (04:39 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 2009 04:39:12 +0000 (04:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@19092 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/geometry_info.h
deal.II/base/source/geometry_info.cc
deal.II/deal.II/source/grid/grid_tools.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/doc/news/changes.h
tests/base/geometry_info_6.cc

index 596c4fcc8faeae1428ce8816f35192a8607de996..4b3745b53982ec126e19a06c85c4300b53dbc415 100644 (file)
@@ -2041,15 +2041,68 @@ struct GeometryInfo
                                      const unsigned int i);
     
                                     /**
-                                     * For a (bi-, tri-)linear mapping from
-                                     * the reference cell, face, or edge to
-                                     * the object specified by the given
-                                     * vertices, compute the determinant of
-                                     * the Jacobian at the vertices. Note
-                                     * that it is the actual determinant, not
-                                     * its absolute value as often used in
-                                     * transforming integrals from one
-                                     * coordinate system to another.
+                                     * For a (bi-, tri-)linear
+                                     * mapping from the reference
+                                     * cell, face, or edge to the
+                                     * object specified by the given
+                                     * vertices, compute the
+                                     * alternating form of the
+                                     * transformed unit vectors
+                                     * vertices. For an object of
+                                     * dimensionality @p dim, there
+                                     * are @p dim vectors with @p
+                                     * spacedim components each, and
+                                     * the alternating form is a
+                                     * tensor of rank spacedim-dim
+                                     * that corresponds to the wedge
+                                     * product of the @p dim unit
+                                     * vectors, and it corresponds to
+                                     * the volume and normal vectors
+                                     * of the mapping from reference
+                                     * element to the element
+                                     * described by the vertices.
+                                     *
+                                     * For example, if
+                                     * dim==spacedim==2, then the
+                                     * alternating form is a scalar
+                                     * (because spacedim-dim=0) and
+                                     * its value equals $\mathbf
+                                     * v_1\wedge \mathbf v_2=\mathbf
+                                     * v_1\cdot\mathbf v_2$. If
+                                     * dim==spacedim==3, then the
+                                     * result is again a scalar with
+                                     * value $\mathbf v_1\wedge
+                                     * \mathbf v_2 \wedge \mathbf v_3
+                                     * = (\mathbf v_1\times \mathbf
+                                     * v_2)\cdot \mathbf v_3$, where
+                                     * $\mathbf v_1, \mathbf v_2,
+                                     * \mathbf v_3$ are the images of
+                                     * the unit vectors at a vertex
+                                     * of the unit dim-dimensional
+                                     * cell under transformation to
+                                     * the dim-dimensional cell in
+                                     * spacedim-dimensional space. In
+                                     * both cases, i.e. for dim==2 or
+                                     * 3, the result happens to equal
+                                     * the determinant of the
+                                     * Jacobian of the mapping from
+                                     * reference cell to cell in real
+                                     * space. Note that it is the
+                                     * actual determinant, not its
+                                     * absolute value as often used
+                                     * in transforming integrals from
+                                     * one coordinate system to
+                                     * another. In particular, if the
+                                     * object specified by the
+                                     * vertices is a parallelogram
+                                     * (i.e. a linear transformation
+                                     * of the reference cell) then
+                                     * the computed values are the
+                                     * same at all vertices and equal
+                                     * the (signed) area of the cell;
+                                     * similarly, for
+                                     * parallel-epipeds, it is the
+                                     * volume of the cell.
                                      *
                                      * This function is used in order to
                                      * determine how distorted a cell is (see
@@ -2057,10 +2110,11 @@ struct GeometryInfo
                                      * @ref GlossDistorted "distorted cells"
                                      * in the glossary).
                                      */
+    template <int spacedim>
     static
     void
-    jacobian_determinants_at_vertices (const Point<dim> (&vertices)[vertices_per_cell],
-                                      double (&determinants)[vertices_per_cell]);
+    alternating_form_at_vertices (const Point<spacedim> (&vertices)[vertices_per_cell],
+                                 Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell]);
     
                                      /**
                                      * For each face of the reference
index 2388cf460cc72495363eb97c40d34eca977c92a7..81e07c7b3c376f00db5667c0b96c48f6dc69bd80 100644 (file)
@@ -1625,35 +1625,96 @@ d_linear_shape_function_gradient (const Point<dim> &,
 
 
 
-template <int dim>
-void
-GeometryInfo<dim>::
-jacobian_determinants_at_vertices (const Point<dim> (&vertices)[vertices_per_cell],
-                                  double (&determinants)[vertices_per_cell])
+namespace internal
 {
-                                  // for each of the vertices, form the
-                                  // Jacobian matrix and compute its
-                                  // determinant
-                                  //
-                                  // note that the transformation is
-                                  //    \vec x = sum_i \vec v_i phi_i(\vec xi)
-                                  // and we have to take the gradient
-                                  // with respect to \vec \xi.
-  for (unsigned int i=0; i<vertices_per_cell; ++i)
+  namespace GeometryInfo
+  {
+    template <int dim>
+    inline
+    void
+    alternating_form_at_vertices
+    (const Point<dim> (&vertices)[dealii::GeometryInfo<dim>::vertices_per_cell],
+     Tensor<0,dim> (&determinants)[dealii::GeometryInfo<dim>::vertices_per_cell])
     {
-      Tensor<2,dim> jacobian;
-      for (unsigned int j=0; j<vertices_per_cell; ++j)
+                                      // for each of the vertices, form the
+                                      // Jacobian matrix and compute its
+                                      // determinant
+                                      //
+                                      // note that the transformation is
+                                      //    \vec x = sum_i \vec v_i phi_i(\vec xi)
+                                      // and we have to take the gradient
+                                      // with respect to \vec \xi.
+      for (unsigned int i=0; i<dealii::GeometryInfo<dim>::vertices_per_cell; ++i)
        {
-         Tensor<2,dim> x;
-         outer_product (x,
-                        vertices[j],
-                        d_linear_shape_function_gradient (unit_cell_vertex(i),
-                                                          j));
-         jacobian += x;
-       }
+         Tensor<2,dim> jacobian;
+         for (unsigned int j=0; j<dealii::GeometryInfo<dim>::vertices_per_cell; ++j)
+           {
+             Tensor<2,dim> x;
+             outer_product (x,
+                            vertices[j],
+                            dealii::GeometryInfo<dim>::
+                            d_linear_shape_function_gradient (dealii::GeometryInfo<dim>::
+                                                              unit_cell_vertex(i),
+                                                              j));
+             jacobian += x;
+           }
       
-      determinants[i] = determinant (jacobian);
+         determinants[i] = determinant (jacobian);
+       }
+    }
+
+                                    /**
+                                     * Alternating form for lines in 2d.
+                                     */
+    inline
+    void
+    alternating_form_at_vertices
+    (const Point<2> (&vertices)[dealii::GeometryInfo<1>::vertices_per_cell],
+     Tensor<1,2> (&forms)[dealii::GeometryInfo<1>::vertices_per_cell])
+    {
+                                      // for each of the vertices,
+                                      // form the scaled normal
+                                      // vector. since the mapping is
+                                      // linear, all normals are the
+                                      // same
+      const Point<2> d = vertices[1]-vertices[0];
+
+                                      // choose the right normal
+      forms[0] = forms[1] = Point<2> (d[1], -d[0]);
+    }
+
+
+                                    /**
+                                     * Alternating form for quads in 3d.
+                                     */
+    inline
+    void
+    alternating_form_at_vertices
+    (const Point<3> (&vertices)[dealii::GeometryInfo<2>::vertices_per_cell],
+     Tensor<1,3> (&forms)[dealii::GeometryInfo<2>::vertices_per_cell])
+    {
+                                      // for each of the vertices,
+                                      // form the scaled normal
+                                      // vector. since the mapping is
+                                      // linear, all normals are the
+                                      // same
+      (void)vertices;
+      (void)forms;
     }
+  }
+}
+
+
+template <int dim>
+template <int spacedim>
+void
+GeometryInfo<dim>::
+alternating_form_at_vertices (const Point<spacedim> (&vertices)[vertices_per_cell],
+                             Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
+{
+                                  // forward to a template that we
+                                  // can specialize
+  internal::GeometryInfo::alternating_form_at_vertices (vertices, forms);
 }
 
 
@@ -1662,4 +1723,36 @@ template class GeometryInfo<2>;
 template class GeometryInfo<3>;
 template class GeometryInfo<4>;
 
+template
+void
+GeometryInfo<1>::
+alternating_form_at_vertices (const Point<1> (&vertices)[vertices_per_cell],
+                             Tensor<1-1,1> (&forms)[vertices_per_cell]);
+
+template
+void
+GeometryInfo<1>::
+alternating_form_at_vertices (const Point<2> (&vertices)[vertices_per_cell],
+                             Tensor<2-1,2> (&forms)[vertices_per_cell]);
+
+template
+void
+GeometryInfo<2>::
+alternating_form_at_vertices (const Point<2> (&vertices)[vertices_per_cell],
+                             Tensor<2-2,2> (&forms)[vertices_per_cell]);
+
+template
+void
+GeometryInfo<2>::
+alternating_form_at_vertices (const Point<3> (&vertices)[vertices_per_cell],
+                             Tensor<3-2,3> (&forms)[vertices_per_cell]);
+
+template
+void
+GeometryInfo<3>::
+alternating_form_at_vertices (const Point<3> (&vertices)[vertices_per_cell],
+                             Tensor<3-3,3> (&forms)[vertices_per_cell]);
+
+
+
 DEAL_II_NAMESPACE_CLOSE
index 3eb1644f959528ae891fa8c0dcf80be9e352def9..b43ced022a017c48776fa7b246bed04d476bc468 100644 (file)
@@ -1259,14 +1259,14 @@ namespace internal
                                         // cell
        Point<spacedim> parent_vertices
          [GeometryInfo<dim>::vertices_per_cell];
-       double parent_determinants
+       Tensor<0,dim> parent_determinants
          [GeometryInfo<dim>::vertices_per_cell];
 
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
          parent_vertices[i] = cell->vertex(i);
 
-       GeometryInfo<dim>::jacobian_determinants_at_vertices (parent_vertices,
-                                                             parent_determinants);
+       GeometryInfo<dim>::alternating_form_at_vertices (parent_vertices,
+                                                        parent_determinants);
 
        const double average_parent_jacobian
          = std::accumulate (&parent_determinants[0],
@@ -1283,7 +1283,7 @@ namespace internal
        Point<spacedim> child_vertices
          [GeometryInfo<dim>::max_children_per_cell]
          [GeometryInfo<dim>::vertices_per_cell];
-       double child_determinants
+       Tensor<0,dim> child_determinants
          [GeometryInfo<dim>::max_children_per_cell]
          [GeometryInfo<dim>::vertices_per_cell];
        
@@ -1302,8 +1302,8 @@ namespace internal
            = cell_mid_point;
        
        for (unsigned int c=0; c<cell->n_children(); ++c)
-         GeometryInfo<dim>::jacobian_determinants_at_vertices (child_vertices[c],
-                                                               child_determinants[c]);
+         GeometryInfo<dim>::alternating_form_at_vertices (child_vertices[c],
+                                                          child_determinants[c]);
 
                                         // on a uniformly refined
                                         // hypercube cell, the child
@@ -1318,7 +1318,7 @@ namespace internal
        double objective = 0;
        for (unsigned int c=0; c<cell->n_children(); ++c)
          for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-           objective += std::pow (child_determinants[c][i] -
+           objective += std::pow (static_cast<double>(child_determinants[c][i]) -
                                   average_parent_jacobian/std::pow(2.,1.*dim),
                                   2);
        
@@ -1419,7 +1419,7 @@ namespace internal
            Point<spacedim> child_vertices
              [GeometryInfo<dim>::max_children_per_cell]
              [GeometryInfo<dim>::vertices_per_cell];
-           double child_determinants
+           Tensor<0,dim> child_determinants
              [GeometryInfo<dim>::max_children_per_cell]
              [GeometryInfo<dim>::vertices_per_cell];
 
@@ -1439,14 +1439,14 @@ namespace internal
                = cell_mid_point;
        
            for (unsigned int c=0; c<cell->n_children(); ++c)
-             GeometryInfo<dim>::jacobian_determinants_at_vertices (child_vertices[c],
-                                                                   child_determinants[c]);
+             GeometryInfo<dim>::alternating_form_at_vertices (child_vertices[c],
+                                                              child_determinants[c]);
        
            double min = child_determinants[0][0];
            for (unsigned int c=0; c<cell->n_children(); ++c)
              for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-               min = std::min (min,
-                               child_determinants[c][i]);
+               min = std::min<double> (min,
+                                       child_determinants[c][i]);
 
            if (test == 0)
              old_min_jacobian = min;
index 04514d8c9e2e81d5b378ae8ea3221d0c83dde896..66f9338fcf0897ad64521597d29fa005d37cfeb1 100644 (file)
@@ -1154,9 +1154,9 @@ namespace
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
          vertices[i] = cell->vertex(i);
 
-       double determinants[GeometryInfo<dim>::vertices_per_cell];
-       GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
-                                                             determinants);
+       Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+       GeometryInfo<dim>::alternating_form_at_vertices (vertices,
+                                                        determinants);
 
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
          if (determinants[i] <= 1e-9 * std::pow (cell->diameter(),
@@ -1212,9 +1212,9 @@ namespace
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
          vertices[i] = cell->child(c)->vertex(i);
 
-       double determinants[GeometryInfo<dim>::vertices_per_cell];
-       GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
-                                                             determinants);
+       Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+       GeometryInfo<dim>::alternating_form_at_vertices (vertices,
+                                                        determinants);
 
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
          if (determinants[i] <= 1e-9 * std::pow (cell->child(c)->diameter(),
index cd6585b57232f5d95fb255d9a116e10d4ebe0b4b..3a481d5e2cc861ab75aeb0c1ed177d13a3a33241 100644 (file)
@@ -136,7 +136,7 @@ inconvenience this causes.
 
   <li>
   <p>
-  New: The GeometryInfo::jacobian_determinants_at_vertices can be used
+  New: The GeometryInfo::alternating_form_at_vertices can be used
   to investigate the degree of distortion of cells.
   <br>
   (WB 2009/06/28)
index bd65c974a9146bbe3f3e3cae5b5f1676d41acf5b..03dafbf788099053ebcb86744a64dfb6eef6a63e 100644 (file)
@@ -12,7 +12,7 @@
 //----------------------------  geometry_info_6.cc  ---------------------------
 
 
-// check GeometryInfo::jacobian_determinants_at_vertices
+// check GeometryInfo::alternating_form_at_vertices
 
 #include "../tests.h"
 #include <base/logstream.h>
@@ -36,14 +36,14 @@ void test ()
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       vertices[v] = GeometryInfo<dim>::unit_cell_vertex(v);
 
-    double determinants[GeometryInfo<dim>::vertices_per_cell];
-    GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+    Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::alternating_form_at_vertices (vertices,
                                                          determinants);
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       {
        deallog << "Reference cell: " << determinants[v]
                << std::endl;
-       Assert (determinants[v] == 1, ExcInternalError());
+       Assert (static_cast<double>(determinants[v]) == 1, ExcInternalError());
       }
   }
 
@@ -57,14 +57,14 @@ void test ()
        vertices[v][0] /= 10;
       }
 
-    double determinants[GeometryInfo<dim>::vertices_per_cell];
-    GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+    Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::alternating_form_at_vertices (vertices,
                                                          determinants);
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       {
        deallog << "Squashed cell: " << determinants[v]
                << std::endl;
-       Assert (determinants[v] == 0.1, ExcInternalError());
+       Assert (static_cast<double>(determinants[v]) == 0.1, ExcInternalError());
       }
   }
 
@@ -87,14 +87,14 @@ void test ()
          }
       }
 
-    double determinants[GeometryInfo<dim>::vertices_per_cell];
-    GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+    Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::alternating_form_at_vertices (vertices,
                                                          determinants);
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       {
        deallog << "Squashed+rotated cell: " << determinants[v]
                << std::endl;
-       Assert (determinants[v] == 0.1, ExcInternalError());
+       Assert (static_cast<double>(determinants[v]) == 0.1, ExcInternalError());
       }
   }
 
@@ -105,8 +105,8 @@ void test ()
       vertices[v] = GeometryInfo<dim>::unit_cell_vertex(v);
     vertices[1] /= 10;
     
-    double determinants[GeometryInfo<dim>::vertices_per_cell];
-    GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+    Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::alternating_form_at_vertices (vertices,
                                                          determinants);
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       deallog << "Pinched cell: " << determinants[v]
@@ -121,8 +121,8 @@ void test ()
       vertices[v] = GeometryInfo<dim>::unit_cell_vertex(v);
     std::swap (vertices[0], vertices[1]);
     
-    double determinants[GeometryInfo<dim>::vertices_per_cell];
-    GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+    Tensor<0,dim> determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::alternating_form_at_vertices (vertices,
                                                          determinants);
     for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
       deallog << "Inverted cell: " << determinants[v]

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.