* Given a vector of <i>n</i> one-dimensional polynomials
* <i>P<sub>0</sub></i> to <i>P<sub>n</sub></i>, where
* <i>P<sub>i</sub></i> has degree <i>i</i>, this class generates all
- * polynomials of the form <i> P<sub>ijk</sub>(x,y,z) =
+ * multi-dimensional polynomials of the form <i>
+ * P<sub>ijk</sub>(x,y,z) =
* P<sub>i</sub>(x)P<sub>j</sub>(y)P<sub>k</sub>(z)</i>, where the sum
* of <i>i</i>, <i>j</i> and <i>k</i> is less than or equal <i>n</i>.
*
- * @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 <int dim>
class PolynomialSpace
template <class Pol>
PolynomialSpace (const std::vector<Pol> &pols);
+ /**
+ * Prints the list of the indices
+ * to <tt>out</tt>.
+ */
+ void output_indices(std::ostream &out) const;
+
+ /**
+ * Sets the ordering of the
+ * polynomials. Requires
+ * <tt>index_map.size()==n()</tt>. Stores
+ * a copy of <tt>index_map</tt>.
+ */
+ void set_polynomial_ordering(const vector<unsigned int> &index_map);
+
/**
* Computes the value and the
* first and second derivatives
int, int, int,
<< "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3);
-
+ protected:
+
+ /**
+ * Compute numbers in x, y and z
+ * direction. Given an index
+ * <tt>n</tt> in the d-dimensional
+ * polynomial space, compute the
+ * indices i,j,k such that
+ * <i>p<sub>n</sub>(x,y,z) =
+ * p<sub>i</sub>(x)p<sub>j</sub>(y)p<sub>k</sub>(z)</i>.
+ */
+ void compute_index (const unsigned int n,
+ unsigned int (&index)[dim]) const;
+
private:
/**
* Copy of the vector <tt>pols</tt> of
const unsigned int n_pols;
/**
- * Compute numbers in x, y and z
- * direction. Given an index
- * <tt>n</tt> in the d-dimensional
- * polynomial space, compute the
- * indices i,j,k such that
- * <i>p<sub>n</sub>(x,y,z) =
- * p<sub>i</sub>(x)p<sub>j</sub>(y)p<sub>k</sub>(z)</i>.
+ * Index map for reordering the
+ * polynomials.
*/
- void compute_index (const unsigned int n,
- unsigned int (&index)[dim]) const;
+ std::vector<unsigned int> index_map;
+
+ /**
+ * Index map for reordering the
+ * polynomials.
+ */
+ std::vector<unsigned int> reverse_index_map;
/**
* Static function used in the
template <int dim>
template <class Pol>
-PolynomialSpace<dim>::
-PolynomialSpace (const std::vector<Pol> &pols)
+PolynomialSpace<dim>::PolynomialSpace (const std::vector<Pol> &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<n_pols; ++i)
+ {
+ index_map[i]=i;
+ reverse_index_map[i]=i;
+ }
+}
template<int dim>
// $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
#include <base/table.h>
+
template <int dim>
unsigned int
PolynomialSpace<dim>::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<index_map.size(),
+ ExcIndexRange(i,0,index_map.size()));
+ const unsigned int n=index_map[i];
index[0] = n;
}
template <>
void
PolynomialSpace<2>::
-compute_index (const unsigned int n,
+compute_index (const unsigned int i,
unsigned int (&index)[2]) const
{
+ Assert(i<index_map.size(),
+ ExcIndexRange(i,0,index_map.size()));
+ const unsigned int n=index_map[i];
// there should be a better way to
// write this function (not
// linear in n_1d), someone
template <>
void
PolynomialSpace<3>::
-compute_index (const unsigned int n,
+compute_index (const unsigned int i,
unsigned int (&index)[3]) const
{
+ Assert(i<index_map.size(),
+ ExcIndexRange(i,0,index_map.size()));
+ const unsigned int n=index_map[i];
// there should be a better way to
// write this function (not
// quadratic in n_1d), someone
}
+template <int dim>
+void
+PolynomialSpace<dim>::output_indices(std::ostream &out) const
+{
+ unsigned int ix[dim];
+ for (unsigned int i=0; i<n_pols; ++i)
+ {
+ compute_index(i,ix);
+ out << i << "\t";
+ for (unsigned int d=0; d<dim; ++d)
+ out << ix[d] << " ";
+ out << endl;
+ }
+}
+
+
+
+template <int dim>
+void
+PolynomialSpace<dim>::set_polynomial_ordering(
+ const vector<unsigned int> &imap)
+{
+ Assert(imap.size()==index_map.size(),
+ ExcDimensionMismatch(imap.size(), index_map.size()));
+
+ index_map=imap;
+ for (unsigned int i=0; i<index_map.size(); ++i)
+ reverse_index_map[index_map[i]]=i;
+}
+
+
template <int dim>
double
{
unsigned int ix[dim];
compute_index(i,ix);
-
// take the product of the
// polynomials in the various space
// directions
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<n_1d-iy-iz; ++ix)
- values[k++] = v[0][ix][0]
- * ((dim>1) ? 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)
for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
{
+ const unsigned int k2=reverse_index_map[k++];
for (unsigned int d=0;d<dim;++d)
- grads[k][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.);
- ++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.);
}
}
for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
{
+ const unsigned int k2=reverse_index_map[k++];
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int d2=0; d2<dim; ++d2)
{
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;
}
}
}