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
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
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:
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:
// $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
#include <base/geometry_info.h>
+#include <base/tensor.h>
DEAL_II_NAMESPACE_OPEN
}
+
+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>;
<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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+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