From 4a3639ddeb5da4f442c45c39fb8112c73158348b Mon Sep 17 00:00:00 2001 From: guido Date: Mon, 11 Feb 2002 12:58:09 +0000 Subject: [PATCH] new class PolynomialSpace git-svn-id: https://svn.dealii.org/trunk@5500 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial_space.h | 169 ++++++++++++ .../include/base/tensor_product_polynomials.h | 2 +- deal.II/base/source/legendre.cc | 3 +- deal.II/base/source/polynomial_space.cc | 186 +++++++++++++ .../base/source/tensor_product_polynomials.cc | 2 +- deal.II/deal.II/source/fe/mapping_q.cc | 2 +- tests/base/polynomial1d.cc | 18 ++ tests/base/polynomial1d.checked | 38 +++ tests/base/polynomial_test.cc | 258 +++++------------- tests/base/polynomial_test.checked | 224 ++++++++------- 10 files changed, 616 insertions(+), 286 deletions(-) create mode 100644 deal.II/base/include/base/polynomial_space.h create mode 100644 deal.II/base/source/polynomial_space.cc diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h new file mode 100644 index 0000000000..6fd8423abe --- /dev/null +++ b/deal.II/base/include/base/polynomial_space.h @@ -0,0 +1,169 @@ +//---------------------- polynomials.h ------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2002 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. +// +//---------------------- polynomials.h ------------- +#ifndef __deal2__polynomial_space_h +#define __deal2__polynomial_space_h + + +#include +#include +#include +#include +#include +#include + +#include + + +/** + * Polynomial space of degree at most n in higher dimensions. + * + * Given a vector of @{n} one-dimensional polynomials @{P0} to @{Pn}, + * where @{Pi} has degree @p{i}, this class generates all polynomials + * the form @p{ Pijk(x,y,z) = Pi(x)Pj(y)Pk(z)}, where the sum of + * @p{i}, @p{j} and @p{k} is below/equal @p{n}. + * + * @author Guido Kanschat, 2002 + */ +template +class PolynomialSpace +{ + public: + /** + * Constructor. @p{pols} is a + * vector of pointers to + * one-dimensional polynomials + * and will be copied into the + * member variable @p{polynomials}. + */ + template + PolynomialSpace(const typename std::vector &pols); + + /** + * Computes the value and the + * first and second derivatives + * of each polynomial at + * @p{unit_point}. + * + * The size of the vectors must + * either be equal @p{0} or equal + * @p{n()}. In the + * first case, the function will + * not compute these values. + * + * If you need values or + * derivatives of all polynomials + * then use this function, rather + * than using any of the + * @p{compute_value}, + * @p{compute_grad} or + * @p{compute_grad_grad} + * functions, see below, in a + * loop over all polynomials. + */ + void compute(const Point &unit_point, + std::vector &values, + typename std::vector > &grads, + typename std::vector > &grad_grads) const; + + /** + * Computes the value of the + * @p{i}th polynomial at + * @p{unit_point}. + * + * Consider using @p{compute} instead. + */ + double compute_value (const unsigned int i, + const Point &p) const; + + /** + * Computes the gradient of the + * @p{i}th polynomial at + * @p{unit_point}. + * + * Consider using @p{compute} instead. + */ + Tensor<1,dim> compute_grad (const unsigned int i, + const Point &p) const; + + /** + * Computes the second derivative + * (grad_grad) of the @p{i}th + * polynomial at + * @p{unit_point}. + * + * Consider using @p{compute} instead. + */ + Tensor<2,dim> compute_grad_grad(const unsigned int i, + const Point &p) const; + + /** + * Returns the number of tensor + * product polynomials. For $n$ + * 1d polynomials this is $n^dim$. + */ + unsigned int n() const; + + /** + * Exception. + */ + DeclException3 (ExcDimensionMismatch2, + int, int, int, + << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3); + + + private: + /** + * Copy of the vector @p{pols} of + * polynomials given to the + * constructor. + */ + std::vector > polynomials; + + /** + * Number of tensor product + * polynomials. For $n$ 1d + * polynomials this is $n^dim$. + */ + unsigned int n_pols; + + /** + * Computes @p{x} to the power of + * @p{y} for unsigned int @p{x} + * and @p{y}. It is a private + * function as it is only used in + * this class. + */ + static unsigned int power(const unsigned int x, const unsigned int y); +}; + + + +template +template +PolynomialSpace::PolynomialSpace( + const typename std::vector &pols): + polynomials (pols.begin(), pols.end()) +{ + const unsigned int n=polynomials.size(); + + n_pols = n; + for (unsigned int i=1;i +#include + + + +template +unsigned int PolynomialSpace::power(const unsigned int x, + const unsigned int y) +{ + unsigned int value=1; + for (unsigned int i=0; i +double +PolynomialSpace::compute_value(const unsigned int i, + const Point &p) const +{ + Assert(false, ExcNotImplemented()); + return 0.; +} + + +template +Tensor<1,dim> +PolynomialSpace::compute_grad(const unsigned int i, + const Point &p) const +{ + Assert(false, ExcNotImplemented()); + return Tensor<1,dim>(); +} + + +template +Tensor<2,dim> +PolynomialSpace::compute_grad_grad(const unsigned int i, + const Point &p) const +{ + Assert(false, ExcNotImplemented()); + return Tensor<2,dim>(); +} + + + + +template +void PolynomialSpace::compute( + const Point &p, + std::vector &values, + typename std::vector > &grads, + typename std::vector > &grad_grads) const +{ + const unsigned int n_1d=polynomials.size(); + + Assert(values.size()==n_pols || values.size()==0, + ExcDimensionMismatch2(values.size(), n_pols, 0)); + Assert(grads.size()==n_pols|| grads.size()==0, + ExcDimensionMismatch2(grads.size(), n_pols, 0)); + Assert(grad_grads.size()==n_pols|| grad_grads.size()==0, + ExcDimensionMismatch2(grad_grads.size(), n_pols, 0)); + + unsigned int v_size=0; + bool update_values=false, update_grads=false, update_grad_grads=false; + if (values.size()==n_pols) + { + update_values=true; + v_size=1; + } + if (grads.size()==n_pols) + { + update_grads=true; + v_size=2; + } + if (grad_grads.size()==n_pols) + { + update_grad_grads=true; + v_size=3; + } + + // Store data in a single + // vector. Access is by + // v[d][n][o] + // d: coordinate direction + // n: number of 1d polynomial + // o: order of derivative + std::vector > > + v(dim, + std::vector > (n_1d, + std::vector (v_size, 0.))); + + for (unsigned int d=0; d >& v_d=v[d]; + Assert(v_d.size()==n_1d, ExcInternalError()); + for (unsigned int i=0; i2) ? n_1d : 1);++iz) + for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy) + for (unsigned int ix=0; ix1) ? v[1][iy][0] : 1.) + * ((dim>2) ? v[2][iz][0] : 1.); + } + + if (update_grads) + { + unsigned int k = 0; + + for (unsigned int iz=0;iz<((dim>2) ? n_1d : 1);++iz) + for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy) + for (unsigned int ix=0; ix1) ? v[1][iy][(d==1) ? 1 : 0] : 1.) + * ((dim>2) ? v[2][iz][(d==2) ? 1 : 0] : 1.); + ++k; + } + } + + if (update_grad_grads) + { + unsigned int k = 0; + + for (unsigned int iz=0;iz<((dim>2) ? n_1d : 1);++iz) + for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy) + for (unsigned int ix=0; ix1) ? v[1][iy][j1] : 1.) + * ((dim>2) ? v[2][iz][j2] : 1.); + } + ++k; + } + } +} + + +template +unsigned int +PolynomialSpace::n() const +{ + return n_pols; +} + + +template class PolynomialSpace<1>; +template class PolynomialSpace<2>; +template class PolynomialSpace<3>; diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc index 56c71b7a8d..e6522ac3fd 100644 --- a/deal.II/base/source/tensor_product_polynomials.cc +++ b/deal.II/base/source/tensor_product_polynomials.cc @@ -202,7 +202,7 @@ void TensorProductPolynomials::compute( template unsigned int -TensorProductPolynomials::n_tensor_product_polynomials() const +TensorProductPolynomials::n() const { return n_tensor_pols; } diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index d062606494..bdad90bbfe 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -120,7 +120,7 @@ MappingQ::MappingQ (const unsigned int p): v.push_back(LagrangeEquidistant(degree,i)); tensor_pols = new TensorProductPolynomials (v); - Assert (n_shape_functions==tensor_pols->n_tensor_product_polynomials(), + Assert (n_shape_functions==tensor_pols->n(), ExcInternalError()); Assert(n_inner+n_outer==n_shape_functions, ExcInternalError()); diff --git a/tests/base/polynomial1d.cc b/tests/base/polynomial1d.cc index b40f43618a..1e9636d5f3 100644 --- a/tests/base/polynomial1d.cc +++ b/tests/base/polynomial1d.cc @@ -45,6 +45,9 @@ int main () deallog.depth_console(0); std::vector > p (15); + + deallog << "Legendre" << endl; + for (unsigned int i=0;i(i); @@ -52,4 +55,19 @@ int main () for (unsigned int j=0;j<=i;++j) deallog << 'P' << i << " * P" << j << " =" << scalar_product(p[i], p[j]) << std::endl; + + + deallog << "LagrangeEquidistant" << endl; + + p.resize(6); + for (unsigned int i=0;i +void check_poly(const Point& x, + const POLY& p) +{ + const unsigned int n = p.n(); + vector values (n); + vector > gradients(n); + vector > second(n); + + p.compute (x, values, gradients, second); + + for (unsigned int k=0;k +void +check_tensor (const vector >& v, + const Point& x) { - double eps=1e-14; - if ((i==j && std::fabs(value-1) p(v); + check_poly (x, p); + deallog.pop(); } -void Q3_4th_shape_function_values_and_grads_dim2( - const Point<2> &point, double &v_exact, - Tensor<1,2> &grad_exact, Tensor<2,2> &grad2); +template +void +check_poly (const vector >& v, + const Point& x) +{ + deallog.push("Polyno"); + PolynomialSpace p(v); + check_poly (x, p); + deallog.pop(); +} +void +check_dimensions (const vector >& p) +{ + deallog.push("1d"); + check_tensor(p, Point<1>(.5)); + check_poly(p, Point<1>(.5)); + deallog.pop(); + deallog.push("2d"); + check_tensor(p, Point<2>(.5, .2)); + check_poly(p, Point<2>(.5, .2)); + deallog.pop(); + deallog.push("3d"); + check_tensor(p, Point<3>(.5, .2, .3)); + check_poly(p, Point<3>(.5, .2, .3)); + deallog.pop(); +} int main() { @@ -47,185 +103,15 @@ int main() logfile.precision(4); deallog.attach(logfile); deallog.depth_console(0); - - std::vector values(1); - deallog << "LagrangeEquidistant polynoms:" << std::endl; - for (unsigned int order=1; order<=4; ++order) - { - deallog << "Polynomial p of order " << order << std::endl; - for (unsigned int s_point=0; s_point<=order; ++s_point) - { - LagrangeEquidistant polynom(order, s_point); - - // support points in vertices - for (unsigned int i=0; i<=order; ++i) - { - double x=static_cast(i)/order; - polynom.value(x, values); - deallog << " p_" << s_point << "(" << x << ")"; -// deallog << "=" << values[0]; - if (equals_delta_ij(values[0], s_point, i)) - deallog << " ok"; - else - deallog << " false"; - deallog << std::endl; - - // now also check - // whether the other - // @p{value} function - // returns the same - // result - if (polynom.value(x) != values[0]) - { - deallog << "The two `value' functions return different results!" - << std::endl; - abort (); - }; - } - } - } - - deallog << std::endl << "Test derivatives computed by the Horner scheme:" << std::endl; - LagrangeEquidistant pol(4, 2); - std::vector v_horner(6); - for (unsigned int i=0; i<=10; ++i) - { - double xi=i*0.1; - deallog << "x=" << xi << ", all derivatives: "; - std::vector v_exact(6); - - v_exact[0]=64.0*xi*xi*xi*xi-128.0*xi*xi*xi+76.0*xi*xi-12.0*xi; - v_exact[1]=256.0*xi*xi*xi-384.0*xi*xi+152.0*xi-12.0; - v_exact[2]=768.0*xi*xi-768.0*xi+152.0; - v_exact[3]=1536*xi-768; - v_exact[4]=1536; - v_exact[5]=0; - - pol.value(xi, v_horner); - - bool ok=true; - for (unsigned int i=0; i a_const(1,1.); - const Polynomial pol_const(a_const); - std::vector exact_values(5,0.); - exact_values[0]=1.; - std::vector computed_values(5); - - pol_const.value(0.24, computed_values); - bool ok=true; - for (unsigned int i=0; i1e-15) - ok=false; - } - - if (ok) - deallog << "ok"; - else - deallog << "false"; - - deallog << std::endl; - } - - - - - deallog << std::endl << "Test of TensorProductPolynomials:" << std::endl; - deallog << "2D Example:" << std::endl; - unsigned int p=3, - n_tensor_pols=(p+1)*(p+1); - std::vector > pols; - - for (unsigned int i=0; i<=p; ++i) - pols.push_back(LagrangeEquidistant(p, i)); - - TensorProductPolynomials<2> tp_pol(pols); - double v_exact; - Tensor<1,2> grad_exact; - Tensor<2,2> grad_grad_exact; - - Point<2> point(0.35,0.62); - // 4th shape function of Q3<2> is - // equivalent to its 1st shape - // function in lexicographical - // order. - Q3_4th_shape_function_values_and_grads_dim2(point, v_exact, grad_exact, grad_grad_exact); - - unsigned int i=1; - double v=tp_pol.compute_value(i, point); - Tensor<1,2> grad=tp_pol.compute_grad(i, point); - Tensor<2,2> grad_grad=tp_pol.compute_grad_grad(i, point); - - std::vector vs(n_tensor_pols); - std::vector > grads(n_tensor_pols); - std::vector > grad_grads(n_tensor_pols); - tp_pol.compute(point, vs, grads, grad_grads); - - - deallog << "v=" << v << std::endl; - deallog << "vs[" << i << "]=" << vs[i] << std::endl; - deallog << "v_exact=" << v_exact << std::endl; - deallog << "grad=" << grad << std::endl; - deallog << "grads[" << i << "]=" << grads[i] << std::endl; - deallog << "grad_exact=" << grad_exact << std::endl; - for (unsigned int j=0; j