From: Wolfgang Bangerth Date: Thu, 16 Jul 2009 04:39:12 +0000 (+0000) Subject: Generalize jacobian_determinants_at_vertices to alternating_form_a_vertices. X-Git-Tag: v8.0.0~7469 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0c87997a4beb909029b26a0ea0382efd4f436cfa;p=dealii.git Generalize jacobian_determinants_at_vertices to alternating_form_a_vertices. git-svn-id: https://svn.dealii.org/trunk@19092 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h index 596c4fcc8f..4b3745b539 100644 --- a/deal.II/base/include/base/geometry_info.h +++ b/deal.II/base/include/base/geometry_info.h @@ -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 static void - jacobian_determinants_at_vertices (const Point (&vertices)[vertices_per_cell], - double (&determinants)[vertices_per_cell]); + alternating_form_at_vertices (const Point (&vertices)[vertices_per_cell], + Tensor (&forms)[vertices_per_cell]); /** * For each face of the reference diff --git a/deal.II/base/source/geometry_info.cc b/deal.II/base/source/geometry_info.cc index 2388cf460c..81e07c7b3c 100644 --- a/deal.II/base/source/geometry_info.cc +++ b/deal.II/base/source/geometry_info.cc @@ -1625,35 +1625,96 @@ d_linear_shape_function_gradient (const Point &, -template -void -GeometryInfo:: -jacobian_determinants_at_vertices (const Point (&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 + inline + void + alternating_form_at_vertices + (const Point (&vertices)[dealii::GeometryInfo::vertices_per_cell], + Tensor<0,dim> (&determinants)[dealii::GeometryInfo::vertices_per_cell]) { - Tensor<2,dim> jacobian; - for (unsigned int j=0; j::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::vertices_per_cell; ++j) + { + Tensor<2,dim> x; + outer_product (x, + vertices[j], + dealii::GeometryInfo:: + d_linear_shape_function_gradient (dealii::GeometryInfo:: + 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 +template +void +GeometryInfo:: +alternating_form_at_vertices (const Point (&vertices)[vertices_per_cell], + Tensor (&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 diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index 3eb1644f95..b43ced022a 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -1259,14 +1259,14 @@ namespace internal // cell Point parent_vertices [GeometryInfo::vertices_per_cell]; - double parent_determinants + Tensor<0,dim> parent_determinants [GeometryInfo::vertices_per_cell]; for (unsigned int i=0; i::vertices_per_cell; ++i) parent_vertices[i] = cell->vertex(i); - GeometryInfo::jacobian_determinants_at_vertices (parent_vertices, - parent_determinants); + GeometryInfo::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 child_vertices [GeometryInfo::max_children_per_cell] [GeometryInfo::vertices_per_cell]; - double child_determinants + Tensor<0,dim> child_determinants [GeometryInfo::max_children_per_cell] [GeometryInfo::vertices_per_cell]; @@ -1302,8 +1302,8 @@ namespace internal = cell_mid_point; for (unsigned int c=0; cn_children(); ++c) - GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], - child_determinants[c]); + GeometryInfo::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; cn_children(); ++c) for (unsigned int i=0; i::vertices_per_cell; ++i) - objective += std::pow (child_determinants[c][i] - + objective += std::pow (static_cast(child_determinants[c][i]) - average_parent_jacobian/std::pow(2.,1.*dim), 2); @@ -1419,7 +1419,7 @@ namespace internal Point child_vertices [GeometryInfo::max_children_per_cell] [GeometryInfo::vertices_per_cell]; - double child_determinants + Tensor<0,dim> child_determinants [GeometryInfo::max_children_per_cell] [GeometryInfo::vertices_per_cell]; @@ -1439,14 +1439,14 @@ namespace internal = cell_mid_point; for (unsigned int c=0; cn_children(); ++c) - GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], - child_determinants[c]); + GeometryInfo::alternating_form_at_vertices (child_vertices[c], + child_determinants[c]); double min = child_determinants[0][0]; for (unsigned int c=0; cn_children(); ++c) for (unsigned int i=0; i::vertices_per_cell; ++i) - min = std::min (min, - child_determinants[c][i]); + min = std::min (min, + child_determinants[c][i]); if (test == 0) old_min_jacobian = min; diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 04514d8c9e..66f9338fcf 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -1154,9 +1154,9 @@ namespace for (unsigned int i=0; i::vertices_per_cell; ++i) vertices[i] = cell->vertex(i); - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, - determinants); + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, + determinants); for (unsigned int i=0; i::vertices_per_cell; ++i) if (determinants[i] <= 1e-9 * std::pow (cell->diameter(), @@ -1212,9 +1212,9 @@ namespace for (unsigned int i=0; i::vertices_per_cell; ++i) vertices[i] = cell->child(c)->vertex(i); - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, - determinants); + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, + determinants); for (unsigned int i=0; i::vertices_per_cell; ++i) if (determinants[i] <= 1e-9 * std::pow (cell->child(c)->diameter(), diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index cd6585b572..3a481d5e2c 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -136,7 +136,7 @@ inconvenience this causes.
  • - 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.
    (WB 2009/06/28) diff --git a/tests/base/geometry_info_6.cc b/tests/base/geometry_info_6.cc index bd65c974a9..03dafbf788 100644 --- a/tests/base/geometry_info_6.cc +++ b/tests/base/geometry_info_6.cc @@ -12,7 +12,7 @@ //---------------------------- geometry_info_6.cc --------------------------- -// check GeometryInfo::jacobian_determinants_at_vertices +// check GeometryInfo::alternating_form_at_vertices #include "../tests.h" #include @@ -36,14 +36,14 @@ void test () for (unsigned int v=0;v::vertices_per_cell;++v) vertices[v] = GeometryInfo::unit_cell_vertex(v); - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, determinants); for (unsigned int v=0;v::vertices_per_cell;++v) { deallog << "Reference cell: " << determinants[v] << std::endl; - Assert (determinants[v] == 1, ExcInternalError()); + Assert (static_cast(determinants[v]) == 1, ExcInternalError()); } } @@ -57,14 +57,14 @@ void test () vertices[v][0] /= 10; } - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, determinants); for (unsigned int v=0;v::vertices_per_cell;++v) { deallog << "Squashed cell: " << determinants[v] << std::endl; - Assert (determinants[v] == 0.1, ExcInternalError()); + Assert (static_cast(determinants[v]) == 0.1, ExcInternalError()); } } @@ -87,14 +87,14 @@ void test () } } - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, determinants); for (unsigned int v=0;v::vertices_per_cell;++v) { deallog << "Squashed+rotated cell: " << determinants[v] << std::endl; - Assert (determinants[v] == 0.1, ExcInternalError()); + Assert (static_cast(determinants[v]) == 0.1, ExcInternalError()); } } @@ -105,8 +105,8 @@ void test () vertices[v] = GeometryInfo::unit_cell_vertex(v); vertices[1] /= 10; - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, determinants); for (unsigned int v=0;v::vertices_per_cell;++v) deallog << "Pinched cell: " << determinants[v] @@ -121,8 +121,8 @@ void test () vertices[v] = GeometryInfo::unit_cell_vertex(v); std::swap (vertices[0], vertices[1]); - double determinants[GeometryInfo::vertices_per_cell]; - GeometryInfo::jacobian_determinants_at_vertices (vertices, + Tensor<0,dim> determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::alternating_form_at_vertices (vertices, determinants); for (unsigned int v=0;v::vertices_per_cell;++v) deallog << "Inverted cell: " << determinants[v]