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
* @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
-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);
}
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
// 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],
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];
= 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
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);
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];
= 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;
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(),
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(),
<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)
//---------------------------- geometry_info_6.cc ---------------------------
-// check GeometryInfo::jacobian_determinants_at_vertices
+// check GeometryInfo::alternating_form_at_vertices
#include "../tests.h"
#include <base/logstream.h>
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());
}
}
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());
}
}
}
}
- 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());
}
}
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]
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]