From 0f68474e4f753c78656e4d9d084a5ca6b27252b3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 28 Jun 2009 19:13:07 +0000 Subject: [PATCH] Add several functions to GeometryInfo related to mappings from the reference cell to real cells. git-svn-id: https://svn.dealii.org/trunk@18980 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/geometry_info.h | 87 ++++++++- deal.II/base/source/geometry_info.cc | 225 +++++++++++++++++++++- deal.II/doc/news/changes.h | 21 ++ tests/base/geometry_info_4.cc | 82 ++++++++ tests/base/geometry_info_4/cmp/generic | 105 ++++++++++ tests/base/geometry_info_5.cc | 83 ++++++++ tests/base/geometry_info_5/cmp/generic | 105 ++++++++++ tests/base/geometry_info_6.cc | 146 ++++++++++++++ tests/base/geometry_info_6/cmp/generic | 105 ++++++++++ 9 files changed, 952 insertions(+), 7 deletions(-) create mode 100644 tests/base/geometry_info_4.cc create mode 100644 tests/base/geometry_info_4/cmp/generic create mode 100644 tests/base/geometry_info_5.cc create mode 100644 tests/base/geometry_info_5/cmp/generic create mode 100644 tests/base/geometry_info_6.cc create mode 100644 tests/base/geometry_info_6/cmp/generic diff --git a/deal.II/base/include/base/geometry_info.h b/deal.II/base/include/base/geometry_info.h index 5e810ea517..21275a4136 100644 --- a/deal.II/base/include/base/geometry_info.h +++ b/deal.II/base/include/base/geometry_info.h @@ -2020,6 +2020,57 @@ struct GeometryInfo double distance_to_unit_cell (const Point &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 &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 &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 (&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: diff --git a/deal.II/base/source/geometry_info.cc b/deal.II/base/source/geometry_info.cc index 646b779dc7..770514fbae 100644 --- a/deal.II/base/source/geometry_info.cc +++ b/deal.II/base/source/geometry_info.cc @@ -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 +#include DEAL_II_NAMESPACE_OPEN @@ -1435,6 +1436,228 @@ GeometryInfo::distance_to_unit_cell (const Point &p) } + +template +double +GeometryInfo:: +d_linear_shape_function (const Point &xi, + const unsigned int i) +{ + Assert (i < GeometryInfo::vertices_per_cell, + ExcIndexRange (i, 0, GeometryInfo::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 +Tensor<1,dim> +GeometryInfo:: +d_linear_shape_function_gradient (const Point &xi, + const unsigned int i) +{ + Assert (false, ExcNotImplemented()); + return Tensor<1,dim>(); +} + + + + + +template +void +GeometryInfo:: +jacobian_determinants_at_vertices (const Point (&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 jacobian; + for (unsigned int j=0; j 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>; diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 2ec5653fe0..280e2b2d5a 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -93,6 +93,27 @@ inconvenience this causes.

base

    +
  1. +

    + New: The GeometryInfo::jacobian_determinants_at_vertices can be used + to investigate the degree of distortion of cells. +
    + (WB 2009/06/28) +

    +
  2. + +
  3. +

    + 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). +
    + (WB 2009/06/28) +

    +
  4. +
  5. 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 index 0000000000..14f77eb439 --- /dev/null +++ b/tests/base/geometry_info_4.cc @@ -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 +#include + +#include +#include + + +template +void test () +{ + deallog << "Checking in " << dim << "d" << std::endl; + + // check phi_i(v_j) = delta_{ij} + for (unsigned int i=0;i::vertices_per_cell;++i) + for (unsigned int v=0;v::vertices_per_cell;++v) + { + const double + phi_i = GeometryInfo::bilinear_shape_function(GeometryInfo::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::vertices_per_cell;++v) + { + double s = 0; + for (unsigned int i=0;i::vertices_per_cell;++i) + s += GeometryInfo::bilinear_shape_function(GeometryInfo::unit_cell_vertex(v),i); + Assert (s == 1, ExcInternalError()); + + deallog << "Sum of shape functions: " << s << std::endl; + } + { + Point center; + for (unsigned int i=0; i::vertices_per_cell;++i) + s += GeometryInfo::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 index 0000000000..bc9aaa403a --- /dev/null +++ b/tests/base/geometry_info_4/cmp/generic @@ -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 index 0000000000..9f59a4a67b --- /dev/null +++ b/tests/base/geometry_info_5.cc @@ -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 +#include + +#include +#include + + +template +void test () +{ + deallog << "Checking in " << dim << "d" << std::endl; + + // check phi_i(v_j) = delta_{ij} + for (unsigned int i=0;i::vertices_per_cell;++i) + for (unsigned int v=0;v::vertices_per_cell;++v) + { + const Tensor<1,dim> + phi_i_grad + = GeometryInfo::bilinear_shape_function_gradient(GeometryInfo::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::vertices_per_cell;++v) + { + Tensor<1,dim> s; + for (unsigned int i=0;i::vertices_per_cell;++i) + s += GeometryInfo::bilinear_shape_function_gradient(GeometryInfo::unit_cell_vertex(v),i); + Assert (s.norm() == 0, ExcInternalError()); + + deallog << "Sum of shape functions: " << s << std::endl; + } + { + Point center; + for (unsigned int i=0; i s; + for (unsigned int i=0;i::vertices_per_cell;++i) + s += GeometryInfo::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 index 0000000000..f7dde08ae2 --- /dev/null +++ b/tests/base/geometry_info_5/cmp/generic @@ -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 index 0000000000..bd65c974a9 --- /dev/null +++ b/tests/base/geometry_info_6.cc @@ -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 +#include + +#include +#include + + +template +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 vertices[GeometryInfo::vertices_per_cell]; + 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, + determinants); + for (unsigned int v=0;v::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 vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int v=0;v::vertices_per_cell;++v) + { + vertices[v] = GeometryInfo::unit_cell_vertex(v); + vertices[v][0] /= 10; + } + + double determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::jacobian_determinants_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()); + } + } + + + // 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 vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int v=0;v::vertices_per_cell;++v) + { + vertices[v] = GeometryInfo::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::vertices_per_cell]; + GeometryInfo::jacobian_determinants_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()); + } + } + + // pinched cell + { + Point vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int v=0;v::vertices_per_cell;++v) + vertices[v] = GeometryInfo::unit_cell_vertex(v); + vertices[1] /= 10; + + double determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::jacobian_determinants_at_vertices (vertices, + determinants); + for (unsigned int v=0;v::vertices_per_cell;++v) + deallog << "Pinched cell: " << determinants[v] + << std::endl; + } + + + // inverted cell + { + Point vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int v=0;v::vertices_per_cell;++v) + 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, + determinants); + for (unsigned int v=0;v::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 index 0000000000..542baef3cd --- /dev/null +++ b/tests/base/geometry_info_6/cmp/generic @@ -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 -- 2.39.5