From 39184813db91ac0e62ecaa9f2f852b2d93c222bd Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 21 Apr 2003 16:11:54 +0000 Subject: [PATCH] Add anisotropic polynomials. git-svn-id: https://svn.dealii.org/trunk@7416 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/base/tensor_product_polynomials.h | 232 +++++++++++++++ .../base/source/tensor_product_polynomials.cc | 268 ++++++++++++++++++ 2 files changed, 500 insertions(+) diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h index 7dbcd9d70a..066d66520e 100644 --- a/deal.II/base/include/base/tensor_product_polynomials.h +++ b/deal.II/base/include/base/tensor_product_polynomials.h @@ -222,6 +222,223 @@ class TensorProductPolynomials + +/** + * Anisotropic tensor product of given polynomials. + * + * Given one-dimensional polynomials @{Px1}, @{Px2}, ... in + * x-direction, @{Py1}, @{Py2}, ... in y-direction, and so on, this + * class generates polynomials of the form @p{ Qijk(x,y,z) = + * Pxi(x)Pyj(y)Pzk(z)}. If the base polynomials are mutually + * orthogonal on the interval $[-1,1]$ or $[0,d], then the tensor + * product polynomials are orthogonal on $[-1,1]^d$ or $[0,1]^d$, + * respectively. + * + * 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{Px1(x)Py1(y)}, @p{Px2(x)Py1(y)}, + * @p{Px3(x)Py1(y)}, ..., @p{Px1(x)Py2(y)}, @p{Px2(x)Py2(y)}, + * @p{Px3(x)Py2(y)}, ..., and likewise in 3d. + * + * @author Wolfgang Bangerth 2003 + */ +template +class AnisotropicPolynomials +{ + public: + /** + * Constructor. @p{pols} is a + * table of one-dimensional + * polynomials. The number of + * rows in this table should be + * equal to the space dimension, + * with the elements of each row + * giving the polynomials that + * shall be used in this + * particular coordinate + * direction. These polynomials + * may vary between coordinates, + * as well as their number. + */ + AnisotropicPolynomials (const std::vector > > &pols); + + /** + * Computes the value and the + * first and second derivatives + * of each tensor product + * polynomial at @p{unit_point}. + * + * The size of the vectors must + * either be equal @p{0} or equal + * @p{n_tensor_pols}. In the + * first case, the function will + * not compute these values. + * + * If you need values or + * derivatives of all tensor + * product 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 tensor product + * polynomials. + */ + void compute (const Point &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const; + + /** + * Computes the value of the + * @p{i}th tensor product + * polynomial at + * @p{unit_point}. Here @p{i} is + * given in tensor product + * numbering. + * + * Note, that using this function + * within a loop over all tensor + * product polynomials is not + * efficient, because then each + * point value of the underlying + * (one-dimensional) polynomials + * is (unnecessarily) computed + * several times. Instead use + * the @p{compute} function, see + * above, with + * @p{values.size()==n_tensor_pols} + * to get the point values of all + * tensor polynomials all at once + * and in a much more efficient + * way. + */ + double compute_value (const unsigned int i, + const Point &p) const; + + /** + * Computes the grad of the + * @p{i}th tensor product + * polynomial at + * @p{unit_point}. Here @p{i} is + * given in tensor product + * numbering. + * + * Note, that using this function + * within a loop over all tensor + * product polynomials is not + * efficient, because then each + * derivative value of the + * underlying (one-dimensional) + * polynomials is (unnecessarily) + * computed several times. + * Instead use the @p{compute} + * function, see above, with + * @p{grads.size()==n_tensor_pols} + * to get the point value of all + * tensor polynomials all at once + * and in a much more efficient + * way. + */ + Tensor<1,dim> compute_grad (const unsigned int i, + const Point &p) const; + + /** + * Computes the second + * derivative (grad_grad) of the + * @p{i}th tensor product + * polynomial at + * @p{unit_point}. Here @p{i} is + * given in tensor product + * numbering. + * + * Note, that using this function + * within a loop over all tensor + * product polynomials is not + * efficient, because then each + * derivative value of the + * underlying (one-dimensional) + * polynomials is (unnecessarily) + * computed several times. + * Instead use the @p{compute} + * function, see above, with + * @p{grad_grads.size()==n_tensor_pols} + * to get the point value of all + * tensor polynomials all at once + * and in a much more efficient + * way. + */ + Tensor<2,dim> compute_grad_grad (const unsigned int i, + const Point &p) const; + + /** + * Returns the number of tensor + * product polynomials. It is the + * product of the number of + * polynomials in each coordinate + * direction. + */ + unsigned int n() const; + + /** + * Exception. + */ + DeclException3 (ExcDimensionMismatch2, + int, int, int, + << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3); + /** + * Exception + */ + DeclException1 (ExcInvalidDim, + int, + << "The number of rows in this table must be equal to the " + << "space dimension, but is " << arg1); + + private: + /** + * Copy of the vector @p{pols} of + * polynomials given to the + * constructor. + */ + const std::vector > > polynomials; + + /** + * Number of tensor product + * polynomials. This is + * @p{Nx*Ny*Nz}, or with terms + * dropped if the number of space + * dimensions is less than 3. + */ + const unsigned int n_tensor_pols; + + /** + * 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; + + /** + * Given the input to the + * constructor, compute + * @p{n_tensor_pols}. + */ + static + unsigned int + get_n_tensor_pols (const std::vector > > &pols); +}; + + + + /* -------------- declaration of explicit specializations --- */ template <> @@ -265,4 +482,19 @@ TensorProductPolynomials(const std::vector &pols) +template <> +void +AnisotropicPolynomials<1>::compute_index(const unsigned int n, + unsigned int (&index)[1]) const; +template <> +void +AnisotropicPolynomials<2>::compute_index(const unsigned int n, + unsigned int (&index)[2]) const; +template <> +void +AnisotropicPolynomials<3>::compute_index(const unsigned int n, + unsigned int (&index)[3]) const; + + + #endif diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc index d60ecdea40..83e6dbc8c6 100644 --- a/deal.II/base/source/tensor_product_polynomials.cc +++ b/deal.II/base/source/tensor_product_polynomials.cc @@ -18,6 +18,8 @@ +/* ------------------- TensorProductPolynomials -------------- */ + template <> void TensorProductPolynomials<1>:: @@ -240,7 +242,273 @@ TensorProductPolynomials::n() const return n_tensor_pols; } + + + +/* ------------------- AnisotropicPolynomials -------------- */ + + +template +AnisotropicPolynomials:: +AnisotropicPolynomials(const std::vector > > &pols) + : + polynomials (pols), + n_tensor_pols(get_n_tensor_pols(pols)) +{ + Assert (pols.size() == dim, ExcInvalidDim(pols.size())); + for (unsigned int d=0; d 0, + ExcMessage ("The number of polynomials must be larger than zero " + "for all coordinate directions.")); +} + + + + +template <> +void +AnisotropicPolynomials<1>:: +compute_index (const unsigned int i, + unsigned int (&indices)[1]) const +{ + Assert (i +void +AnisotropicPolynomials<2>:: +compute_index (const unsigned int i, + unsigned int (&indices)[2]) const +{ + Assert (i < polynomials[0].size()*polynomials[1].size(), + ExcInternalError()); + + indices[0] = i % polynomials[0].size(); + indices[1] = i / polynomials[0].size(); +} + + + +template <> +void +AnisotropicPolynomials<3>:: +compute_index (const unsigned int i, + unsigned int (&indices)[3]) const +{ + Assert (i < polynomials[0].size()*polynomials[1].size()*polynomials[2].size(), + ExcInternalError()); + + indices[0] = i % polynomials[0].size(); + indices[1] = (i/polynomials[0].size()) % polynomials[1].size(); + indices[2] = i / (polynomials[0].size()*polynomials[1].size()); +} + + + +template +double +AnisotropicPolynomials::compute_value (const unsigned int i, + const Point &p) const +{ + unsigned int indices[dim]; + compute_index (i, indices); + + double value=1.; + for (unsigned int d=0; d +Tensor<1,dim> +AnisotropicPolynomials::compute_grad (const unsigned int i, + const Point &p) const +{ + unsigned int indices[dim]; + compute_index (i, indices); + + // compute values and + // uni-directional derivatives at + // the given point in each + // co-ordinate direction + std::vector > v(dim, std::vector (2)); + for (unsigned int d=0; d grad; + for (unsigned int d=0; d +Tensor<2,dim> +AnisotropicPolynomials::compute_grad_grad (const unsigned int i, + const Point &p) const +{ + 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 +void +AnisotropicPolynomials:: +compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const +{ + 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; + + + // compute the values (and + // derivatives, if necessary) of + // all polynomials at this + // evaluation point + std::vector > > v(dim); + for (unsigned int d=0; d +unsigned int +AnisotropicPolynomials::n() const +{ + return n_tensor_pols; +} + + +template +unsigned int +AnisotropicPolynomials:: +get_n_tensor_pols (const std::vector > > &pols) +{ + unsigned int y = 1; + for (unsigned int d=0; d; template class TensorProductPolynomials<2>; template class TensorProductPolynomials<3>; + +template class AnisotropicPolynomials<1>; +template class AnisotropicPolynomials<2>; +template class AnisotropicPolynomials<3>; -- 2.39.5