From 58523d21e5dbc4729d832b82734a2e7d9cf4024a Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Fri, 30 Oct 2009 03:58:31 +0000 Subject: [PATCH] implement quadrature of dimension zero git-svn-id: https://svn.dealii.org/trunk@20007 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature.h | 42 +++--- deal.II/base/source/quadrature.cc | 186 ++++++++++++++++++------- 2 files changed, 157 insertions(+), 71 deletions(-) diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 4db976bc85..d9e944abeb 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -33,10 +33,20 @@ DEAL_II_NAMESPACE_OPEN * integration formulæ. Their names names prefixed by * Q. Refer to the list of derived classes for more details. * - * The schemes for higher dimensions are tensor products of the - * one-dimansional formulæ. Therefore, a three-dimensional 5-point - * Gauss formula has 125 quadrature points. + * The schemes for higher dimensions are typically tensor products of the + * one-dimensional formulæ, but refer to the section on implementation + * detail below. * + * In order to allow for dimension independent programming, a + * quadrature formula of dimension zero exists. Since an integral over + * zero dimensions is the evaluation at a single point, any + * constructor of such a formula initializes to a single quadrature + * point with weight one. Access to the weight is possible, while + * access to the quadrature point is not permitted, since a Point of + * dimension zero contains no information. The main purpose of these + * formulæ is their use in QProjector, which will create a useful + * formula of dimension one out of them. + * *

Mathematical background

* * For each quadrature formula we denote by m, the maximal @@ -62,22 +72,10 @@ DEAL_II_NAMESPACE_OPEN * points in dim dimensions, where N is the constructor * parameter of QGauss. * - * For some programs it is necessary to have a quadrature object for - * faces. These programs fail to link if compiled for only one space - * dimension, since there quadrature rules for faces just don't make - * no sense. In order to allow these programs to be linked anyway, for - * class Quadrature@<0@> all functions are provided in the - * quadrature.cc file, but they will throw exceptions if actually - * called. The only function which is allowed to be called is the - * constructor taking one integer, which in this case ignores its - * parameter, and of course the destructor. Besides this, it is - * necessary to provide a class Point@<0@> to make the compiler - * happy. This class also does nothing. + * @note Instantiations for this template are provided for dimensions + * 0, 1, 2, and 3 (see the section on @ref Instantiations). * - * @note Instantiations for this template are provided for dimensions 1, 2, - * and 3 (see the section on @ref Instantiations in the manual). - * - * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2005 + * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2005, 2009 */ template class Quadrature : public Subscriptor @@ -413,11 +411,11 @@ Quadrature::weight (const unsigned int i) const template <> Quadrature<0>::Quadrature (const unsigned int); template <> -Quadrature<0>::Quadrature (const Quadrature<0> &); -template <> Quadrature<0>::Quadrature (const Quadrature<-1> &, const Quadrature<1> &); template <> +Quadrature<0>::Quadrature (const Quadrature<1> &); +template <> Quadrature<0>::~Quadrature (); template <> @@ -429,8 +427,6 @@ Quadrature<1>::Quadrature (const Quadrature<0> &); template <> const std::vector > & Quadrature<0>::get_points () const; -template <> -const std::vector & Quadrature<0>::get_weights () const; #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index ea04435538..80d46e8449 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 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 @@ -43,15 +43,8 @@ namespace template <> Quadrature<0>::Quadrature (const unsigned int) - : n_quadrature_points(0) -{} - - - -template <> -Quadrature<0>::Quadrature (const Quadrature<0> &) - : Subscriptor(), - n_quadrature_points(0) + : n_quadrature_points(1), + weights (1, 1.) {} @@ -111,21 +104,9 @@ template <> Quadrature<0>::Quadrature (const Quadrature<-1> &, const Quadrature<1> &) : - n_quadrature_points (0) -{ - Assert (false, ExcInternalError()); -} - - - -template <> -Quadrature<1>::Quadrature (const Quadrature<0> &, - const Quadrature<1> &) - : - n_quadrature_points (0) -{ - Assert (false, ExcInternalError()); -} + n_quadrature_points (1), + weights(1, 1.) +{} @@ -172,10 +153,55 @@ Quadrature::Quadrature (const SubQuadrature &q1, +template <> +Quadrature<1>::Quadrature (const SubQuadrature&, + const Quadrature<1>& q2) + : + n_quadrature_points (q2.size()), + quadrature_points (n_quadrature_points), + weights (n_quadrature_points, 0) +{ + unsigned int present_index = 0; + for (unsigned int i2=0; i2 0) + { + double sum = 0; + for (unsigned int i=0; i0.999999) && (sum<1.000001), ExcInternalError()); + } +#endif +} + + + +template <> +Quadrature<0>::Quadrature (const Quadrature<1> &) + : + n_quadrature_points (1), + weights (1, 1.) +{} + + template <> Quadrature<1>::Quadrature (const Quadrature<0> &) : - Subscriptor(), n_quadrature_points (numbers::invalid_unsigned_int), quadrature_points (), weights () @@ -279,16 +305,6 @@ Quadrature::get_weights () const -template <> -const std::vector & -Quadrature<0>::get_weights () const -{ - Assert (false, ExcInternalError()); - return weights; -} - - - template unsigned int Quadrature::memory_consumption () const @@ -424,10 +440,14 @@ QProjector::rotate (const Quadrature<2> &q, template <> void QProjector<1>::project_to_face (const Quadrature<0> &, - const unsigned int, - std::vector > &) + const unsigned int face_no, + std::vector > &q_points) { - Assert (false, ExcNotImplemented()); + const unsigned int dim=1; + Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim)); + AssertDimension (q_points.size(), 1); + + q_points[0] = Point((double) face_no); } @@ -520,12 +540,16 @@ QProjector<3>::project_to_face (const Quadrature<2> &quadrature, template <> void QProjector<1>::project_to_subface (const Quadrature<0> &, + const unsigned int face_no, const unsigned int, - const unsigned int, - std::vector > &, + std::vector >& q_points, const RefinementCase<0> &) { - Assert(false, ExcNotImplemented()); + const unsigned int dim=1; + Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim)); + AssertDimension (q_points.size(), 1); + + q_points[0] = Point((double) face_no); } @@ -697,10 +721,42 @@ QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, template <> Quadrature<1> -QProjector<1>::project_to_all_faces (const Quadrature<0> &) +QProjector<1>::project_to_all_faces (const Quadrature<0>& quadrature) { - Assert (false, ExcImpossibleInDim(1)); - return Quadrature<1>(0); + const unsigned int dim = 1; + + const unsigned int n_points = 1, + n_faces = GeometryInfo::faces_per_cell; + + // first fix quadrature points + std::vector > q_points; + q_points.reserve(n_points * n_faces); + std::vector > help(n_points); + + + // project to each face and append + // results + for (unsigned int face=0; face weights; + weights.reserve (n_points * n_faces); + for (unsigned int face=0; face(q_points, weights); } @@ -811,10 +867,44 @@ QProjector<3>::project_to_all_faces (const SubQuadrature &quadrature) template <> Quadrature<1> -QProjector<1>::project_to_all_subfaces (const Quadrature<0> &) +QProjector<1>::project_to_all_subfaces (const Quadrature<0>& quadrature) { - Assert (false, ExcImpossibleInDim(1)); - return Quadrature<1>(0); + const unsigned int dim = 1; + + const unsigned int n_points = 1, + n_faces = GeometryInfo::faces_per_cell, + subfaces_per_face = GeometryInfo::max_children_per_face; + + // first fix quadrature points + std::vector > q_points; + q_points.reserve (n_points * n_faces * subfaces_per_face); + std::vector > help(n_points); + + // project to each face and copy + // results + for (unsigned int face=0; face weights; + weights.reserve (n_points * n_faces * subfaces_per_face); + for (unsigned int face=0; face(q_points, weights); } -- 2.39.5