]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add several functions to GeometryInfo related to mappings from the reference cell...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 28 Jun 2009 19:13:07 +0000 (19:13 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 28 Jun 2009 19:13:07 +0000 (19:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@18980 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/geometry_info.h
deal.II/base/source/geometry_info.cc
deal.II/doc/news/changes.h
tests/base/geometry_info_4.cc [new file with mode: 0644]
tests/base/geometry_info_4/cmp/generic [new file with mode: 0644]
tests/base/geometry_info_5.cc [new file with mode: 0644]
tests/base/geometry_info_5/cmp/generic [new file with mode: 0644]
tests/base/geometry_info_6.cc [new file with mode: 0644]
tests/base/geometry_info_6/cmp/generic [new file with mode: 0644]

index 5e810ea517eed4858ef27176d556f1323be81592..21275a4136bc17e8c7794f3fac1d8115e176da29 100644 (file)
@@ -2020,6 +2020,57 @@ struct GeometryInfo
     double
     distance_to_unit_cell (const Point<dim> &p);
 
+                                    /**
+                                     * Compute the value of the $i$-th
+                                     * $d$-linear (i.e. (bi-,tri-)linear)
+                                     * shape function at location $\xi$.
+                                     */
+    static
+    double
+    d_linear_shape_function (const Point<dim> &xi,
+                            const unsigned int i);
+
+                                    /**
+                                     * Compute the gradient of the $i$-th
+                                     * $d$-linear (i.e. (bi-,tri-)linear)
+                                     * shape function at location $\xi$.
+                                     */
+    static
+    Tensor<1,dim>
+    d_linear_shape_function_gradient (const Point<dim> &xi,
+                                     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.
+                                     *
+                                     * This function is used in order to
+                                     * determine how distorted a cell is: if
+                                     * all of the determinants are of roughly
+                                     * equal value and on the order of
+                                     * $h^\text{dim}$ then the cell is
+                                     * well-shaped. For example, a square
+                                     * cell or face has determinants equal to
+                                     * $h^\text{dim}$ whereas a strongly
+                                     * sheared parallelogram has a
+                                     * determinant much smaller. Similarly, a
+                                     * cell with very unequal edge lengths
+                                     * will have widely varying
+                                     * determinants. Finally, an inverted
+                                     * cell will have negative determinants.
+                                     */
+    static
+    void
+    jacobian_determinants_at_vertices (const Point<dim> (&vertices)[vertices_per_cell],
+                                      double (&determinants)[vertices_per_cell]);
+    
                                      /**
                                      * For each face of the reference
                                      * cell, this field stores the
@@ -2136,6 +2187,26 @@ template <>
 const unsigned int GeometryInfo<4>::opposite_face[faces_per_cell];
 #endif
 
+
+template <>
+Tensor<1,1>
+GeometryInfo<1>::
+d_linear_shape_function_gradient (const Point<1> &xi,
+                                 const unsigned int i);
+template <>
+Tensor<1,2>
+GeometryInfo<2>::
+d_linear_shape_function_gradient (const Point<2> &xi,
+                                 const unsigned int i);
+template <>
+Tensor<1,3>
+GeometryInfo<3>::
+d_linear_shape_function_gradient (const Point<3> &xi,
+                                 const unsigned int i);
+
+
+
+
 /* -------------- inline functions ------------- */
 
 namespace internal
@@ -2475,9 +2546,11 @@ GeometryInfo<3>::cell_to_child_coordinates (const Point<3>         &p,
          ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
 
   Point<3> point=p;
-                                  // there might be a cleverer way to do this, but since this function is
-                                  // called in very few places for initialization purposes only, I don't
-                                  // care at the moment
+                                  // there might be a cleverer way to do
+                                  // this, but since this function is called
+                                  // in very few places for initialization
+                                  // purposes only, I don't care at the
+                                  // moment
   switch (refine_case)
     {
       case RefinementCase<3>::cut_x:
@@ -2580,9 +2653,11 @@ GeometryInfo<3>::child_to_cell_coordinates (const Point<3>         &p,
          ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
 
   Point<3> point=p;
-                                  // there might be a cleverer way to do this, but since this function is
-                                  // called in very few places for initialization purposes only, I don't
-                                  // care at the moment
+                                  // there might be a cleverer way to do
+                                  // this, but since this function is called
+                                  // in very few places for initialization
+                                  // purposes only, I don't care at the
+                                  // moment
   switch (refine_case)
     {
       case RefinementCase<3>::cut_x:
index 646b779dc73b0393b691d242c0a3b14bf549736f..770514fbaeb90a04842117cae40264bb1a53718b 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -13,6 +13,7 @@
 
 
 #include <base/geometry_info.h>
+#include <base/tensor.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1435,6 +1436,228 @@ GeometryInfo<dim>::distance_to_unit_cell (const Point<dim> &p)
 }
 
 
+
+template <int dim>
+double
+GeometryInfo<dim>::
+d_linear_shape_function (const Point<dim> &xi,
+                        const unsigned int i)
+{
+  Assert (i < GeometryInfo<dim>::vertices_per_cell,
+         ExcIndexRange (i, 0, GeometryInfo<dim>::vertices_per_cell));
+    
+  switch (dim)
+    {
+      case 1:
+      {
+       const double x = xi[0];
+       switch (i)
+         {
+           case 0:
+                 return 1-x;
+           case 1:
+                 return x;
+         }
+      }
+       
+      case 2:
+      {
+       const double x = xi[0];
+       const double y = xi[1];
+       switch (i)
+         {
+           case 0:
+                 return (1-x)*(1-y);
+           case 1:
+                 return x*(1-y);
+           case 2:
+                 return (1-x)*y;
+           case 3:
+                 return x*y;
+         }
+      }
+
+      case 3:
+      {
+       const double x = xi[0];
+       const double y = xi[1];
+       const double z = xi[2];
+       switch (i)
+         {
+           case 0:
+                 return (1-x)*(1-y)*(1-z);
+           case 1:
+                 return x*(1-y)*(1-z);
+           case 2:
+                 return (1-x)*y*(1-z);
+           case 3:
+                 return x*y*(1-z);
+           case 4:
+                 return (1-x)*(1-y)*z;
+           case 5:
+                 return x*(1-y)*z;
+           case 6:
+                 return (1-x)*y*z;
+           case 7:
+                 return x*y*z;
+         }
+      }
+
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+  return -1e9;
+}
+
+
+
+template <>
+Tensor<1,1>
+GeometryInfo<1>::
+d_linear_shape_function_gradient (const Point<1> &xi,
+                                 const unsigned int i)
+{
+  Assert (i < GeometryInfo<1>::vertices_per_cell,
+         ExcIndexRange (i, 0, GeometryInfo<1>::vertices_per_cell));
+    
+  const double x = xi[0];
+  switch (i)
+    {
+      case 0:
+           return Point<1>(-1.);
+      case 1:
+           return Point<1>(1.);
+    }
+
+  return Point<1>(-1e9);
+}
+
+
+
+template <>
+Tensor<1,2>
+GeometryInfo<2>::
+d_linear_shape_function_gradient (const Point<2> &xi,
+                                 const unsigned int i)
+{
+  Assert (i < GeometryInfo<2>::vertices_per_cell,
+         ExcIndexRange (i, 0, GeometryInfo<2>::vertices_per_cell));
+    
+  const double x = xi[0];
+  const double y = xi[1];
+  switch (i)
+    {
+      case 0:
+           return Point<2>(-(1-y),-(1-x));
+      case 1:
+           return Point<2>(1-y,-x);
+      case 2:
+           return Point<2>(-y, 1-x);
+      case 3:
+           return Point<2>(y,x);
+    }
+  return Point<2> (-1e9, -1e9);
+}
+
+
+
+template <>
+Tensor<1,3>
+GeometryInfo<3>::
+d_linear_shape_function_gradient (const Point<3> &xi,
+                                 const unsigned int i)
+{
+  Assert (i < GeometryInfo<3>::vertices_per_cell,
+         ExcIndexRange (i, 0, GeometryInfo<3>::vertices_per_cell));
+    
+  const double x = xi[0];
+  const double y = xi[1];
+  const double z = xi[2];
+  switch (i)
+    {
+      case 0:
+           return Point<3>(-(1-y)*(1-z),
+                           -(1-x)*(1-z),
+                           -(1-x)*(1-y));
+      case 1:
+           return Point<3>((1-y)*(1-z),
+                           -x*(1-z),
+                           -x*(1-y));
+      case 2:
+           return Point<3>(-y*(1-z),
+                           (1-x)*(1-z),
+                           -(1-x)*y);
+      case 3:
+           return Point<3>(y*(1-z),
+                           x*(1-z),
+                           -x*y);
+      case 4:
+           return Point<3>(-(1-y)*z,
+                           -(1-x)*z,
+                           (1-x)*(1-y));
+      case 5:
+           return Point<3>((1-y)*z,
+                           -x*z,
+                           x*(1-y));
+      case 6:
+           return Point<3>(-y*z,
+                           (1-x)*z,
+                           (1-x)*y);
+      case 7:
+           return Point<3>(y*z, x*z, x*y);
+    }
+
+  return Point<3> (-1e9, -1e9, -1e9);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+GeometryInfo<dim>::
+d_linear_shape_function_gradient (const Point<dim> &xi,
+                                 const unsigned int i)
+{
+  Assert (false, ExcNotImplemented());
+  return Tensor<1,dim>();
+}
+
+
+
+
+
+template <int dim>
+void
+GeometryInfo<dim>::
+jacobian_determinants_at_vertices (const Point<dim> (&vertices)[vertices_per_cell],
+                                  double (&determinants)[vertices_per_cell])
+{
+                                  // 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)
+    {
+      Tensor<2,dim> jacobian;
+      for (unsigned int j=0; j<vertices_per_cell; ++j)
+       {
+         Tensor<2,dim> x;
+         outer_product (x,
+                        vertices[j],
+                        d_linear_shape_function_gradient (unit_cell_vertex(i),
+                                                          j));
+         jacobian += x;
+       }
+      
+      determinants[i] = determinant (jacobian);
+    }
+}
+
+
 template class GeometryInfo<1>;
 template class GeometryInfo<2>;
 template class GeometryInfo<3>;
index 2ec5653fe0504a5c006656fbfcbd4546fd21944e..280e2b2d5a53cfa9e2ab2b0c7e3563404d05c0fc 100644 (file)
@@ -93,6 +93,27 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li>
+  <p>
+  New: The GeometryInfo::jacobian_determinants_at_vertices can be used
+  to investigate the degree of distortion of cells.
+  <br>
+  (WB 2009/06/28)
+  </p>
+  </li>
+
+  <li>
+  <p>
+  New: The GeometryInfo::d_linear_shape_function and
+  GeometryInfo::d_linear_shape_function_gradient functions can be used
+  to represent the $d$-linear shape functions that are frequently
+  used to map the reference cell to real cells (though the
+  Mapping class hierarchy also allows to use higher order mappings).
+  <br>
+  (WB 2009/06/28)
+  </p>
+  </li>
+
   <li>
   <p>
   New: The determinant() function is now implemented for rank-2 Tensor
diff --git a/tests/base/geometry_info_4.cc b/tests/base/geometry_info_4.cc
new file mode 100644 (file)
index 0000000..14f77eb
--- /dev/null
@@ -0,0 +1,82 @@
+//----------------------------  geometry_info_4.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2003, 2004, 2005, 2006, 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  geometry_info_4.cc  ---------------------------
+
+
+// check GeometryInfo::bilinear_shape_function
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/geometry_info.h>
+
+#include <fstream>
+#include <cstdlib>
+
+
+template <int dim>
+void test ()
+{
+  deallog << "Checking in " << dim << "d" << std::endl;
+
+                                  // check phi_i(v_j) = delta_{ij}
+  for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      {
+       const double
+         phi_i = GeometryInfo<dim>::bilinear_shape_function(GeometryInfo<dim>::unit_cell_vertex(v),i);
+       
+       deallog << phi_i << std::endl;
+       Assert (phi_i == (i==v ? 1 : 0),
+               ExcInternalError());
+      }
+
+                                  // check that
+                                  //    sum_i phi_i(x) == 1
+                                  // at all points. do so at every
+                                  // vertex, and then at the center
+  for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+    {
+      double s = 0;
+      for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+       s += GeometryInfo<dim>::bilinear_shape_function(GeometryInfo<dim>::unit_cell_vertex(v),i);
+      Assert (s == 1, ExcInternalError());
+
+      deallog << "Sum of shape functions: " << s << std::endl;
+    }
+  {
+    Point<dim> center;
+    for (unsigned int i=0; i<dim; ++i)
+      center[i] = 0.5;
+    
+    double s = 0;
+    for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+      s += GeometryInfo<dim>::bilinear_shape_function(center,i);
+    Assert (s == 1, ExcInternalError());
+
+    deallog << "Sum of shape functions: " << s << std::endl;
+  }
+}
+
+
+int main () 
+{
+  std::ofstream logfile("geometry_info_4/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+
+  return 0;
+}
diff --git a/tests/base/geometry_info_4/cmp/generic b/tests/base/geometry_info_4/cmp/generic
new file mode 100644 (file)
index 0000000..bc9aaa4
--- /dev/null
@@ -0,0 +1,105 @@
+
+DEAL::Checking in 1d
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Checking in 2d
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Checking in 3d
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
+DEAL::Sum of shape functions: 1.00000
diff --git a/tests/base/geometry_info_5.cc b/tests/base/geometry_info_5.cc
new file mode 100644 (file)
index 0000000..9f59a4a
--- /dev/null
@@ -0,0 +1,83 @@
+//----------------------------  geometry_info_5.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2003, 2004, 2005, 2006, 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  geometry_info_5.cc  ---------------------------
+
+
+// check GeometryInfo::bilinear_shape_function_gradient
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/geometry_info.h>
+
+#include <fstream>
+#include <cstdlib>
+
+
+template <int dim>
+void test ()
+{
+  deallog << "Checking in " << dim << "d" << std::endl;
+
+                                  // check phi_i(v_j) = delta_{ij}
+  for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      {
+       const Tensor<1,dim>
+         phi_i_grad
+         = GeometryInfo<dim>::bilinear_shape_function_gradient(GeometryInfo<dim>::unit_cell_vertex(v),i);
+       
+       deallog << phi_i_grad << std::endl;
+      }
+
+                                  // check that
+                                  //    sum_i phi_i(x) == const
+                                  // at all points by verifying that the
+                                  // gradient of the sum of shape functions
+                                  // is zero. do so at every vertex, and then
+                                  // at the center
+  for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+    {
+      Tensor<1,dim> s;
+      for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+       s += GeometryInfo<dim>::bilinear_shape_function_gradient(GeometryInfo<dim>::unit_cell_vertex(v),i);
+      Assert (s.norm() == 0, ExcInternalError());
+
+      deallog << "Sum of shape functions: " << s << std::endl;
+    }
+  {
+    Point<dim> center;
+    for (unsigned int i=0; i<dim; ++i)
+      center[i] = 0.5;
+    
+    Tensor<1,dim> s;
+    for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_cell;++i)
+      s += GeometryInfo<dim>::bilinear_shape_function_gradient(center,i);
+    Assert (s.norm() == 0, ExcInternalError());
+
+    deallog << "Sum of shape functions: " << s << std::endl;
+  }
+}
+
+
+int main () 
+{
+  std::ofstream logfile("geometry_info_5/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+
+  return 0;
+}
diff --git a/tests/base/geometry_info_5/cmp/generic b/tests/base/geometry_info_5/cmp/generic
new file mode 100644 (file)
index 0000000..f7dde08
--- /dev/null
@@ -0,0 +1,105 @@
+
+DEAL::Checking in 1d
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Checking in 2d
+DEAL::-1.00000 -1.00000
+DEAL::-1.00000 -0.00000
+DEAL::-0.00000 -1.00000
+DEAL::-0.00000 -0.00000
+DEAL::1.00000 -0.00000
+DEAL::1.00000 -1.00000
+DEAL::0.00000 -0.00000
+DEAL::0.00000 -1.00000
+DEAL::-0.00000 1.00000
+DEAL::-0.00000 0.00000
+DEAL::-1.00000 1.00000
+DEAL::-1.00000 0.00000
+DEAL::0.00000 0.00000
+DEAL::0.00000 1.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 1.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Checking in 3d
+DEAL::-1.00000 -1.00000 -1.00000
+DEAL::-1.00000 -0.00000 -0.00000
+DEAL::-0.00000 -1.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -1.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::1.00000 -0.00000 -0.00000
+DEAL::1.00000 -1.00000 -1.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -1.00000 -0.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -0.00000 -1.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::-0.00000 1.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-1.00000 1.00000 -1.00000
+DEAL::-1.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -1.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 1.00000 -0.00000
+DEAL::1.00000 0.00000 -0.00000
+DEAL::1.00000 1.00000 -1.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -1.00000
+DEAL::-0.00000 -0.00000 1.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-1.00000 -1.00000 1.00000
+DEAL::-1.00000 -0.00000 0.00000
+DEAL::-0.00000 -1.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 1.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::1.00000 -0.00000 0.00000
+DEAL::1.00000 -1.00000 1.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -1.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 0.00000 1.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 1.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-1.00000 1.00000 1.00000
+DEAL::-1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 1.00000 1.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
diff --git a/tests/base/geometry_info_6.cc b/tests/base/geometry_info_6.cc
new file mode 100644 (file)
index 0000000..bd65c97
--- /dev/null
@@ -0,0 +1,146 @@
+//----------------------------  geometry_info_6.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2003, 2004, 2005, 2006, 2009 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  geometry_info_6.cc  ---------------------------
+
+
+// check GeometryInfo::jacobian_determinants_at_vertices
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/geometry_info.h>
+
+#include <fstream>
+#include <cstdlib>
+
+
+template <int dim>
+void test ()
+{
+  deallog << "Checking in " << dim << "d" << std::endl;
+
+                                  // check the determinant of the
+                                  // transformation for the reference
+                                  // cell. the determinant should be one in
+                                  // that case
+  {
+    Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+    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,
+                                                         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());
+      }
+  }
+
+                                  // try the same, but move squash the cell
+                                  // in the x-direction by a factor of 10
+  {
+    Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      {
+       vertices[v] = GeometryInfo<dim>::unit_cell_vertex(v);
+       vertices[v][0] /= 10;
+      }
+
+    double determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::jacobian_determinants_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());
+      }
+  }
+
+
+                                  // try the same, but move squash the cell
+                                  // in the x-direction by a factor of 10 and
+                                  // rotate it around the z-axis (unless in
+                                  // 1d)
+  {
+    Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      {
+       vertices[v] = GeometryInfo<dim>::unit_cell_vertex(v);
+       vertices[v][0] /= 10;
+
+       if (dim > 1)
+         {
+           std::swap (vertices[v][0], vertices[v][1]);
+           vertices[v][1] *= -1;
+         }
+      }
+
+    double determinants[GeometryInfo<dim>::vertices_per_cell];
+    GeometryInfo<dim>::jacobian_determinants_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());
+      }
+  }
+
+                                  // pinched cell
+  {
+    Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      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,
+                                                         determinants);
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      deallog << "Pinched cell: " << determinants[v]
+             << std::endl;
+  }
+
+
+                                  // inverted cell
+  {
+    Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      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,
+                                                         determinants);
+    for (unsigned int v=0;v<GeometryInfo<dim>::vertices_per_cell;++v)
+      deallog << "Inverted cell: " << determinants[v]
+             << std::endl;
+  }
+}
+
+
+int main () 
+{
+  std::ofstream logfile("geometry_info_6/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+
+  return 0;
+}
diff --git a/tests/base/geometry_info_6/cmp/generic b/tests/base/geometry_info_6/cmp/generic
new file mode 100644 (file)
index 0000000..542baef
--- /dev/null
@@ -0,0 +1,105 @@
+JobId unknown Sun Jun 28 18:42:42 2009
+DEAL::Checking in 1d
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Sum of shape functions: 0.00000
+DEAL::Checking in 2d
+DEAL::-1.00000 -1.00000
+DEAL::-1.00000 -0.00000
+DEAL::-0.00000 -1.00000
+DEAL::-0.00000 -0.00000
+DEAL::1.00000 -0.00000
+DEAL::1.00000 -1.00000
+DEAL::0.00000 -0.00000
+DEAL::0.00000 -1.00000
+DEAL::-0.00000 1.00000
+DEAL::-0.00000 0.00000
+DEAL::-1.00000 1.00000
+DEAL::-1.00000 0.00000
+DEAL::0.00000 0.00000
+DEAL::0.00000 1.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 1.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000
+DEAL::Checking in 3d
+DEAL::-1.00000 -1.00000 -1.00000
+DEAL::-1.00000 -0.00000 -0.00000
+DEAL::-0.00000 -1.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -1.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::-0.00000 -0.00000 -0.00000
+DEAL::1.00000 -0.00000 -0.00000
+DEAL::1.00000 -1.00000 -1.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -1.00000 -0.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -0.00000 -1.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::0.00000 -0.00000 -0.00000
+DEAL::-0.00000 1.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-1.00000 1.00000 -1.00000
+DEAL::-1.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::-0.00000 0.00000 -1.00000
+DEAL::-0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 1.00000 -0.00000
+DEAL::1.00000 0.00000 -0.00000
+DEAL::1.00000 1.00000 -1.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -0.00000
+DEAL::0.00000 0.00000 -1.00000
+DEAL::-0.00000 -0.00000 1.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::-1.00000 -1.00000 1.00000
+DEAL::-1.00000 -0.00000 0.00000
+DEAL::-0.00000 -1.00000 0.00000
+DEAL::-0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 1.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::1.00000 -0.00000 0.00000
+DEAL::1.00000 -1.00000 1.00000
+DEAL::0.00000 -0.00000 0.00000
+DEAL::0.00000 -1.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 0.00000 1.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-0.00000 1.00000 0.00000
+DEAL::-0.00000 0.00000 0.00000
+DEAL::-1.00000 1.00000 1.00000
+DEAL::-1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 1.00000 1.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000
+DEAL::Sum of shape functions: 0.00000 0.00000 0.00000

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.