From 614752c04d832acdfe06e20cd32719d87bf3de1c Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Tue, 9 Mar 2004 14:39:31 +0000 Subject: [PATCH] New functions output_indices and set_polynomial_ordering. git-svn-id: https://svn.dealii.org/trunk@8698 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial_space.h | 78 ++++++++++++++++---- deal.II/base/source/polynomial_space.cc | 74 +++++++++++++++---- 2 files changed, 120 insertions(+), 32 deletions(-) diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h index b837533f08..c3b375d5fc 100644 --- a/deal.II/base/include/base/polynomial_space.h +++ b/deal.II/base/include/base/polynomial_space.h @@ -31,11 +31,19 @@ * Given a vector of n one-dimensional polynomials * P0 to Pn, where * Pi has degree i, this class generates all - * polynomials of the form Pijk(x,y,z) = + * multi-dimensional polynomials of the form + * Pijk(x,y,z) = * Pi(x)Pj(y)Pk(z), where the sum * of i, j and k is less than or equal n. * - * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003 + * The @ref{output_indices} function prints the ordering of the + * polynomials, i.e. for each multi-dimensional polynomial in the + * polynomial space it gives the indices i,j,k of the one-dimensional + * polynomials in x,y and z direction. The ordering of the + * multi-dimensional polynomials can be changed by using the + * @p{set_polynomial_ordering} function. + * + * @author Guido Kanschat, 2002, Wolfgang Bangerth, 2003, Ralf Hartmann 2004 */ template class PolynomialSpace @@ -58,6 +66,20 @@ class PolynomialSpace template PolynomialSpace (const std::vector &pols); + /** + * Prints the list of the indices + * to out. + */ + void output_indices(std::ostream &out) const; + + /** + * Sets the ordering of the + * polynomials. Requires + * index_map.size()==n(). Stores + * a copy of index_map. + */ + void set_polynomial_ordering(const vector &index_map); + /** * Computes the value and the * first and second derivatives @@ -152,7 +174,20 @@ class PolynomialSpace int, int, int, << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3); - + protected: + + /** + * Compute numbers in x, y and z + * direction. Given an index + * n in the d-dimensional + * polynomial space, compute the + * indices i,j,k such that + * pn(x,y,z) = + * pi(x)pj(y)pk(z). + */ + void compute_index (const unsigned int n, + unsigned int (&index)[dim]) const; + private: /** * Copy of the vector pols of @@ -169,16 +204,16 @@ class PolynomialSpace const unsigned int n_pols; /** - * Compute numbers in x, y and z - * direction. Given an index - * n in the d-dimensional - * polynomial space, compute the - * indices i,j,k such that - * pn(x,y,z) = - * pi(x)pj(y)pk(z). + * Index map for reordering the + * polynomials. */ - void compute_index (const unsigned int n, - unsigned int (&index)[dim]) const; + std::vector index_map; + + /** + * Index map for reordering the + * polynomials. + */ + std::vector reverse_index_map; /** * Static function used in the @@ -208,12 +243,23 @@ void PolynomialSpace<3>::compute_index(const unsigned int n, template template -PolynomialSpace:: -PolynomialSpace (const std::vector &pols) +PolynomialSpace::PolynomialSpace (const std::vector &pols) : polynomials (pols.begin(), pols.end()), - n_pols (compute_n_pols(polynomials.size())) -{} + n_pols (compute_n_pols(polynomials.size())), + index_map(n_pols), + reverse_index_map(n_pols) +{ + // per default set this index map + // to identity. This map can be + // changed by the user through the + // set_polynomial_ordering function + for (unsigned int i=0; i diff --git a/deal.II/base/source/polynomial_space.cc b/deal.II/base/source/polynomial_space.cc index aefd6e5ff3..eaa6ef02ff 100644 --- a/deal.II/base/source/polynomial_space.cc +++ b/deal.II/base/source/polynomial_space.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -17,6 +17,7 @@ #include + template unsigned int PolynomialSpace::compute_n_pols (const unsigned int n) @@ -34,9 +35,12 @@ PolynomialSpace::compute_n_pols (const unsigned int n) template <> void PolynomialSpace<1>:: -compute_index (const unsigned int n, +compute_index (const unsigned int i, unsigned int (&index)[1]) const { + Assert(i void PolynomialSpace<2>:: -compute_index (const unsigned int n, +compute_index (const unsigned int i, unsigned int (&index)[2]) const { + Assert(i void PolynomialSpace<3>:: -compute_index (const unsigned int n, +compute_index (const unsigned int i, unsigned int (&index)[3]) const { + Assert(i +void +PolynomialSpace::output_indices(std::ostream &out) const +{ + unsigned int ix[dim]; + for (unsigned int i=0; i +void +PolynomialSpace::set_polynomial_ordering( + const vector &imap) +{ + Assert(imap.size()==index_map.size(), + ExcDimensionMismatch(imap.size(), index_map.size())); + + index_map=imap; + for (unsigned int i=0; i double @@ -104,7 +145,6 @@ PolynomialSpace::compute_value (const unsigned int i, { unsigned int ix[dim]; compute_index(i,ix); - // take the product of the // polynomials in the various space // directions @@ -235,9 +275,10 @@ PolynomialSpace::compute (const Point &p, 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][0] : 1.) - * ((dim>2) ? v[2][iz][0] : 1.); + values[reverse_index_map[k++]] = + v[0][ix][0] + * ((dim>1) ? v[1][iy][0] : 1.) + * ((dim>2) ? v[2][iz][0] : 1.); } if (update_grads) @@ -248,11 +289,11 @@ PolynomialSpace::compute (const Point &p, 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; + grads[k2][d] = v[0][ix][(d==0) ? 1 : 0] + * ((dim>1) ? v[1][iy][(d==1) ? 1 : 0] : 1.) + * ((dim>2) ? v[2][iz][(d==2) ? 1 : 0] : 1.); } } @@ -264,6 +305,7 @@ PolynomialSpace::compute (const Point &p, for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy) for (unsigned int ix=0; ix::compute (const Point &p, const unsigned int j2 = ((d1==2) ? 1 : 0) + ((d2==2) ? 1 : 0); - grad_grads[k][d1][d2] = v[0][ix][j0] - * ((dim>1) ? v[1][iy][j1] : 1.) - * ((dim>2) ? v[2][iz][j2] : 1.); + grad_grads[k2][d1][d2] = + v[0][ix][j0] + * ((dim>1) ? v[1][iy][j1] : 1.) + * ((dim>2) ? v[2][iz][j2] : 1.); } - ++k; } } } -- 2.39.5