From: wolf Date: Mon, 21 Apr 2003 16:09:21 +0000 (+0000) Subject: Clean up some real old and unintelligible cruft. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c060d0364c910fb54f4ebfba64c8ca3cfbde9b9a;p=dealii-svn.git Clean up some real old and unintelligible cruft. git-svn-id: https://svn.dealii.org/trunk@7415 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h index ad06c30f7c..0151d48dae 100644 --- a/deal.II/base/include/base/polynomial_space.h +++ b/deal.II/base/include/base/polynomial_space.h @@ -33,7 +33,7 @@ * of 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 less than or equal @p{n}. * - * @author Guido Kanschat, 2002 + * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003 */ template class PolynomialSpace @@ -55,7 +55,7 @@ class PolynomialSpace * @p{Polynomial}. */ template - PolynomialSpace(const std::vector &pols); + PolynomialSpace (const std::vector &pols); /** * Computes the value and the @@ -65,9 +65,12 @@ class PolynomialSpace * * 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. + * @p{n()}. In the first case, + * the function will not compute + * these values, i.e. you + * indicate what you want to have + * computed by resizing those + * vectors which you want filled. * * If you need values or * derivatives of all polynomials @@ -79,8 +82,8 @@ class PolynomialSpace * functions, see below, in a * loop over all polynomials. */ - void compute (const Point &unit_point, - std::vector &values, + void compute (const Point &unit_point, + std::vector &values, std::vector > &grads, std::vector > &grad_grads) const; @@ -112,8 +115,8 @@ class PolynomialSpace * * Consider using @p{compute} instead. */ - Tensor<2,dim> compute_grad_grad(const unsigned int i, - const Point &p) const; + Tensor<2,dim> compute_grad_grad (const unsigned int i, + const Point &p) const; /** * Return the number of @@ -127,7 +130,7 @@ class PolynomialSpace * 2d, and @p{N(N+1)(N+2)/6 in * 3d. */ - unsigned int n() const; + unsigned int n () const; /** * Exception. @@ -161,10 +164,8 @@ class PolynomialSpace * @p{p_n(x,y,z) = * p_i(x)p_j(y)p_k(z)}. */ - void compute_index(unsigned int n, - unsigned int& nx, - unsigned int& ny, - unsigned int& nz) const; + void compute_index (const unsigned int n, + unsigned int (&index)[dim]) const; /** * Static function used in the @@ -175,6 +176,21 @@ class PolynomialSpace }; +/* -------------- declaration of explicit specializations --- */ + +template <> +void PolynomialSpace<1>::compute_index(const unsigned int n, + unsigned int (&index)[1]) const; +template <> +void PolynomialSpace<2>::compute_index(const unsigned int n, + unsigned int (&index)[2]) const; +template <> +void PolynomialSpace<3>::compute_index(const unsigned int n, + unsigned int (&index)[3]) const; + + + +/* -------------- inline and template functions ------------- */ template template diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h index 92d05e3caa..7dbcd9d70a 100644 --- a/deal.II/base/include/base/tensor_product_polynomials.h +++ b/deal.II/base/include/base/tensor_product_polynomials.h @@ -19,7 +19,6 @@ #include #include #include -#include #include @@ -34,7 +33,13 @@ * $[0,d], then the tensor product polynomials are orthogonal on * $[-1,1]^d$ or $[0,1]^d$, respectively. * - * @author Ralf Hartmann, 2000, documentation Guido Kanschat + * Indexing is as following: the order of dim-dimensional polynomials + * is x-coordinates running fastest, then y-coordinate, etc. The first + * few polynomials are thus @p{P1(x)P1(y)}, @p{P2(x)P1(y)}, + * @p{P3(x)P1(y)}, ..., @p{P1(x)P2(y)}, @p{P2(x)P2(y)}, + * @p{P3(x)P2(y)}, ..., and likewise in 3d. + * + * @author Ralf Hartmann, Guido Kanschat, 2000, Wolfgang Bangerth 2003 */ template class TensorProductPolynomials @@ -42,13 +47,15 @@ class TensorProductPolynomials public: /** * Constructor. @p{pols} is a - * vector of pointers to - * one-dimensional polynomials - * and will be copied into the - * member variable @p{polynomials}. + * vector of objects that should + * be derived or otherwise + * convertible to one-dimensional + * polynomial objects and will be + * copied into the member + * variable @p{polynomials}. */ template - TensorProductPolynomials(const std::vector &pols); + TensorProductPolynomials (const std::vector &pols); /** * Computes the value and the @@ -74,10 +81,10 @@ class TensorProductPolynomials * loop over all tensor product * polynomials. */ - void compute(const Point &unit_point, - std::vector &values, - std::vector > &grads, - std::vector > &grad_grads) const; + void compute (const Point &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const; /** * Computes the value of the @@ -158,8 +165,8 @@ class TensorProductPolynomials * and in a much more efficient * way. */ - Tensor<2,dim> compute_grad_grad(const unsigned int i, - const Point &p) const; + Tensor<2,dim> compute_grad_grad (const unsigned int i, + const Point &p) const; /** * Returns the number of tensor @@ -191,42 +198,70 @@ class TensorProductPolynomials */ const unsigned int n_tensor_pols; - /** - * @p{n_pols_to[n]=polynomials.size()^n} - * Filled by the constructor. - * - * For internal use only. - */ - std::vector n_pols_to; + /** + * Each tensor product polynomial + * @รพ{i} is a product of + * one-dimensional polynomials in + * each space direction. Compute + * the indices of these + * one-dimensional polynomials + * for each space direction, + * given the index @p{i}. + */ + void compute_index (const unsigned int i, + unsigned int (&indices)[dim]) const; /** * 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. + * @p{dim} for unsigned int @p{x}. + * Used in the constructor. */ - static unsigned int power(const unsigned int x, const unsigned int y); + static + unsigned int x_to_the_dim (const unsigned int x); }; +/* -------------- declaration of explicit specializations --- */ + +template <> +void +TensorProductPolynomials<1>::compute_index(const unsigned int n, + unsigned int (&index)[1]) const; +template <> +void +TensorProductPolynomials<2>::compute_index(const unsigned int n, + unsigned int (&index)[2]) const; +template <> +void +TensorProductPolynomials<3>::compute_index(const unsigned int n, + unsigned int (&index)[3]) const; + + +/* ---------------- template and inline functions ---------- */ + +template +inline +unsigned int +TensorProductPolynomials:: +x_to_the_dim (const unsigned int x) +{ + unsigned int y = 1; + for (unsigned int d=0; d template TensorProductPolynomials:: TensorProductPolynomials(const std::vector &pols) : polynomials (pols.begin(), pols.end()), - n_tensor_pols(power(pols.size(), dim)), - n_pols_to(dim+1) -{ - const unsigned int n_pols=polynomials.size(); - - n_pols_to[0]=1; - for (unsigned int i=0; i::compute_n_pols (const unsigned int n) { unsigned int n_pols = n; - for (unsigned int i=1;i::compute_n_pols (const unsigned int n) } -template +template <> +void +PolynomialSpace<1>:: +compute_index(const unsigned int n, + unsigned int (&index)[1]) const +{ + index[0] = n; +} + + + +template <> void -PolynomialSpace::compute_index(unsigned int n, - unsigned int& nx, - unsigned int& ny, - unsigned int& nz) const +PolynomialSpace<2>:: +compute_index (const unsigned int n, + unsigned int (&index)[2]) const { + // there should be a better way to + // write this function (not + // linear in n_1d), someone + // should think about this... const unsigned int n_1d=polynomials.size(); 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; ix +void +PolynomialSpace<3>:: +compute_index (const unsigned int n, + unsigned int (&index)[3]) const +{ + // there should be a better way to + // write this function (not + // quadratic in n_1d), someone + // should think about this... + // + // (ah, and yes: the original + // algorithm was even cubic!) + const unsigned int n_1d=polynomials.size(); + unsigned int k=0; + for (unsigned int iz=0; iz double -PolynomialSpace::compute_value(const unsigned int i, - const Point & p) const +PolynomialSpace::compute_value (const unsigned int i, + const Point &p) const { - unsigned int ix = 0; - unsigned int iy = 0; - unsigned int iz = 0; - compute_index(i,ix,iy,iz); - - double result = polynomials[ix].value(p(0)); - if (dim>1) - result *= polynomials[iy].value(p(1)); - if (dim>2) - result *= polynomials[iz].value(p(2)); + unsigned int ix[dim]; + compute_index(i,ix); + + // take the product of the + // polynomials in the various space + // directions + double result = 1.; + for (unsigned int d=0; d Tensor<1,dim> -PolynomialSpace::compute_grad(const unsigned int i, - const Point &p) const +PolynomialSpace::compute_grad (const unsigned int i, + const Point &p) const { - unsigned int ix[3]; - compute_index(i,ix[0],ix[1],ix[2]); + unsigned int ix[dim]; + compute_index(i,ix); Tensor<1,dim> result; - for (unsigned int d=0;d v(2); - for (unsigned int d=0;d::compute_grad(const unsigned int i, template Tensor<2,dim> -PolynomialSpace::compute_grad_grad(const unsigned int i, - const Point &p) const +PolynomialSpace::compute_grad_grad (const unsigned int i, + const Point &p) const { - unsigned int ix[3]; - compute_index(i,ix[0],ix[1],ix[2]); + unsigned int ix[dim]; + compute_index(i,ix); Tensor<2,dim> result; - for (unsigned int d=0;d v(3); - for (unsigned int d=0;d::compute_grad_grad(const unsigned int i, template -void PolynomialSpace::compute( - const Point &p, - std::vector &values, - std::vector > &grads, - std::vector > &grad_grads) const +void +PolynomialSpace::compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const { const unsigned int n_1d=polynomials.size(); diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc index 6c9a62c1d3..d60ecdea40 100644 --- a/deal.II/base/source/tensor_product_polynomials.cc +++ b/deal.II/base/source/tensor_product_polynomials.cc @@ -18,29 +18,60 @@ +template <> +void +TensorProductPolynomials<1>:: +compute_index (const unsigned int i, + unsigned int (&indices)[1]) const +{ + Assert (i -unsigned int TensorProductPolynomials::power(const unsigned int x, - const unsigned int y) + + +template <> +void +TensorProductPolynomials<2>:: +compute_index (const unsigned int i, + unsigned int (&indices)[2]) const { - unsigned int value=1; - for (unsigned int i=0; i +void +TensorProductPolynomials<3>:: +compute_index (const unsigned int i, + unsigned int (&indices)[3]) const +{ + const unsigned int n_pols = polynomials.size(); + Assert (i double -TensorProductPolynomials::compute_value(const unsigned int i, - const Point &p) const +TensorProductPolynomials::compute_value (const unsigned int i, + const Point &p) const { - const unsigned int n_pols=polynomials.size(); + unsigned int indices[dim]; + compute_index (i, indices); double value=1.; for (unsigned int d=0; d::compute_value(const unsigned int i, template Tensor<1,dim> -TensorProductPolynomials::compute_grad(const unsigned int i, - const Point &p) const +TensorProductPolynomials::compute_grad (const unsigned int i, + const Point &p) const { - const unsigned int n_pols=polynomials.size(); + unsigned int indices[dim]; + compute_index (i, indices); // compute values and // uni-directional derivatives at @@ -59,15 +91,15 @@ TensorProductPolynomials::compute_grad(const unsigned int i, // co-ordinate direction std::vector > v(dim, std::vector (2)); for (unsigned int d=0; d grad; for (unsigned int d=0; d::compute_grad(const unsigned int i, template Tensor<2,dim> -TensorProductPolynomials::compute_grad_grad(const unsigned int i, - const Point &p) const +TensorProductPolynomials::compute_grad_grad (const unsigned int i, + const Point &p) const { - const unsigned int n_pols=polynomials.size(); - + unsigned int indices[dim]; + compute_index (i, indices); + std::vector > v(dim, std::vector (3)); for (unsigned int d=0; d grad_grad; - for (unsigned int d1=0; d1::compute_grad_grad(const unsigned int i, template -void TensorProductPolynomials::compute( - const Point &p, - std::vector &values, - std::vector > &grads, - std::vector > &grad_grads) const +void +TensorProductPolynomials:: +compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const { - const unsigned int n_pols=polynomials.size(); - - Assert(values.size()==n_tensor_pols || values.size()==0, - ExcDimensionMismatch2(values.size(), n_tensor_pols, 0)); - Assert(grads.size()==n_tensor_pols|| grads.size()==0, - ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0)); - Assert(grad_grads.size()==n_tensor_pols|| grad_grads.size()==0, - ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0)); - - unsigned int v_size=0; - bool update_values=false, update_grads=false, update_grad_grads=false; - if (values.size()==n_tensor_pols) - { - update_values=true; - v_size=1; - } - if (grads.size()==n_tensor_pols) - { - update_grads=true; - v_size=2; - } - if (grad_grads.size()==n_tensor_pols) - { - update_grad_grads=true; - v_size=3; - } + Assert (values.size()==n_tensor_pols || values.size()==0, + ExcDimensionMismatch2(values.size(), n_tensor_pols, 0)); + Assert (grads.size()==n_tensor_pols|| grads.size()==0, + ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0)); + Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0, + ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0)); + + const bool update_values = (values.size() == n_tensor_pols), + update_grads = (grads.size()==n_tensor_pols), + update_grad_grads = (grad_grads.size()==n_tensor_pols); + + // check how many + // values/derivatives we have to + // compute + unsigned int n_values_and_derivatives = 0; + if (update_values) + n_values_and_derivatives = 1; + if (update_grads) + n_values_and_derivatives = 2; + if (update_grad_grads) + n_values_and_derivatives = 3; + - Table<2,std::vector > v(dim, n_pols); - for (unsigned int d=0; d > v(dim, polynomials.size()); + for (unsigned int d=0; d unsigned int TensorProductPolynomials::n() const