* 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 <int dim>
class PolynomialSpace
* @p{Polynomial<double>}.
*/
template <class Pol>
- PolynomialSpace(const std::vector<Pol> &pols);
+ PolynomialSpace
(const std::vector<Pol> &pols);
/**
* Computes the value and the
*
* 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
* functions, see below, in a
* loop over all polynomials.
*/
- void compute (const Point<dim> &unit_point,
- std::vector<double> &values,
+ void compute (const Point<dim> &unit_point,
+ std::vector<double> &values,
std::vector<Tensor<1,dim> > &grads,
std::vector<Tensor<2,dim> > &grad_grads) const;
*
* Consider using @p{compute} instead.
*/
- Tensor<2,dim> compute_grad_grad(const unsigned int i,
- const Point<dim> &p) const;
+ Tensor<2,dim> compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
/**
* Return the number of
* 2d, and @p{N(N+1)(N+2)/6 in
* 3d.
*/
- unsigned int n() const;
+ unsigned int n () const;
/**
* Exception.
* @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
};
+/* -------------- 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 <int dim>
template <class Pol>
#include <base/tensor.h>
#include <base/point.h>
#include <base/polynomial.h>
-#include <base/smartpointer.h>
#include <vector>
* $[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 <int dim>
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 <class Pol>
- TensorProductPolynomials(const std::vector<Pol> &pols);
+ TensorProductPolynomials
(const std::vector<Pol> &pols);
/**
* Computes the value and the
* loop over all tensor product
* polynomials.
*/
- void compute(const Point<dim> &unit_point,
- std::vector<double> &values,
- std::vector<Tensor<1,dim> > &grads,
- std::vector<Tensor<2,dim> > &grad_grads) const;
+ void compute (const Point<dim> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &grad_grads) const;
/**
* Computes the value of the
* and in a much more efficient
* way.
*/
- Tensor<2,dim> compute_grad_grad(const unsigned int i,
- const Point<dim> &p) const;
+ Tensor<2,dim> compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
/**
* Returns the number of tensor
*/
const unsigned int n_tensor_pols;
- /**
- * @p{n_pols_to[n]=polynomials.size()^n}
- * Filled by the constructor.
- *
- * For internal use only.
- */
- std::vector<unsigned int> 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 <int dim>
+inline
+unsigned int
+TensorProductPolynomials<dim>::
+x_to_the_dim (const unsigned int x)
+{
+ unsigned int y = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ y *= x;
+ return y;
+}
+
+
+
template <int dim>
template <class Pol>
TensorProductPolynomials<dim>::
TensorProductPolynomials(const std::vector<Pol> &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<dim; ++i)
- n_pols_to[i+1]=n_pols_to[i]*n_pols;
- Assert(n_pols_to[dim]==n_tensor_pols, ExcInternalError());
-}
+ n_tensor_pols(x_to_the_dim(pols.size()))
+{}
PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
{
unsigned int n_pols = n;
- for (unsigned int i=1;i<dim;++i)
+ for (unsigned int i=1; i<dim; ++i)
{
n_pols *= (n+i);
n_pols /= (i+1);
}
-template <int dim>
+template <>
+void
+PolynomialSpace<1>::
+compute_index(const unsigned int n,
+ unsigned int (&index)[1]) const
+{
+ index[0] = n;
+}
+
+
+
+template <>
void
-PolynomialSpace<dim>::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<n_1d-iy-iz; ++ix)
- if (k++ == n)
- {
- nz = iz;
- ny = iy;
- nx = ix;
- return;
- }
+ for (unsigned int iy=0; iy<n_1d; ++iy)
+ if (n < k+n_1d-iy)
+ {
+ index[0] = n-k;
+ index[1] = iy;
+ return;
+ }
+ else
+ k+=n_1d-iy;
}
+
+template <>
+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<n_1d; ++iz)
+ for (unsigned int iy=0; iy<n_1d-iz; ++iy)
+ if (n < k+n_1d-iy-iz)
+ {
+ index[0] = n-k;
+ index[1] = iy;
+ index[2] = iz;
+ return;
+ }
+ else
+ k += n_1d-iy-iz;
+}
+
+
+
template <int dim>
double
-PolynomialSpace<dim>::compute_value(const unsigned int i,
- const Point<dim> & p) const
+PolynomialSpace<dim>::compute_value (const unsigned int i,
+ const Point<dim> &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<dim; ++d)
+ result *= polynomials[ix[d]].value(p(d));
return result;
}
-
+
+
template <int dim>
Tensor<1,dim>
-PolynomialSpace<dim>::compute_grad(const unsigned int i,
- const Point<dim> &p) const
+PolynomialSpace<dim>::compute_grad (const unsigned int i,
+ const Point<dim> &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<dim;++d)
+ for (unsigned int d=0; d<dim; ++d)
result[d] = 1.;
-
+
+ // get value and first derivative
std::vector<double> v(2);
- for (unsigned int d=0;d<dim;++d)
+ for (unsigned int d=0; d<dim; ++d)
{
polynomials[ix[d]].value(p(d), v);
result[d] *= v[1];
- for (unsigned int d1=0;d1<dim;++d1)
+ for (unsigned int d1=0; d1<dim; ++d1)
if (d1 != d)
result[d1] *= v[0];
}
template <int dim>
Tensor<2,dim>
-PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
- const Point<dim> &p) const
+PolynomialSpace<dim>::compute_grad_grad (const unsigned int i,
+ const Point<dim> &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<dim;++d)
- for (unsigned int d1=0;d1<dim;++d1)
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int d1=0; d1<dim; ++d1)
result[d][d1] = 1.;
+ // get value, first and second
+ // derivatives
std::vector<double> v(3);
- for (unsigned int d=0;d<dim;++d)
+ for (unsigned int d=0; d<dim; ++d)
{
polynomials[ix[d]].value(p(d), v);
result[d][d] *= v[2];
- for (unsigned int d1=0;d1<dim;++d1)
+ for (unsigned int d1=0; d1<dim; ++d1)
{
if (d1 != d)
{
result[d][d1] *= v[1];
result[d1][d] *= v[1];
- for (unsigned int d2=0;d2<dim;++d2)
+ for (unsigned int d2=0; d2<dim; ++d2)
if (d2 != d)
result[d1][d2] *= v[0];
}
template <int dim>
-void PolynomialSpace<dim>::compute(
- const Point<dim> &p,
- std::vector<double> &values,
- std::vector<Tensor<1,dim> > &grads,
- std::vector<Tensor<2,dim> > &grad_grads) const
+void
+PolynomialSpace<dim>::compute (const Point<dim> &p,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &grad_grads) const
{
const unsigned int n_1d=polynomials.size();
+template <>
+void
+TensorProductPolynomials<1>::
+compute_index (const unsigned int i,
+ unsigned int (&indices)[1]) const
+{
+ Assert (i<polynomials.size(), ExcInternalError());
+ indices[0] = i;
+}
-template <int dim>
-unsigned int TensorProductPolynomials<dim>::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<y; ++i)
- value*=x;
- return value;
+ const unsigned int n_pols = polynomials.size();
+ Assert (i<n_pols*n_pols, ExcInternalError());
+
+ indices[0] = i % n_pols;
+ indices[1] = i / n_pols;
+}
+
+
+
+template <>
+void
+TensorProductPolynomials<3>::
+compute_index (const unsigned int i,
+ unsigned int (&indices)[3]) const
+{
+ const unsigned int n_pols = polynomials.size();
+ Assert (i<n_pols*n_pols*n_pols, ExcInternalError());
+
+ indices[0] = i % n_pols;
+ indices[1] = (i/n_pols) % n_pols;
+ indices[2] = i / (n_pols*n_pols);
}
template <int dim>
double
-TensorProductPolynomials<dim>::compute_value(const unsigned int i,
- const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_value (const unsigned int i,
+ const Point<dim> &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<dim; ++d)
- value *= polynomials[(i/n_pols_to[d])%n_pols].value(p(d));
+ value *= polynomials[indices[d]].value(p(d));
return value;
}
template <int dim>
Tensor<1,dim>
-TensorProductPolynomials<dim>::compute_grad(const unsigned int i,
- const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_grad (const unsigned int i,
+ const Point<dim> &p) const
{
- const unsigned int n_pols=polynomials.size();
+ unsigned int indices[dim];
+ compute_index (i, indices);
// compute values and
// uni-directional derivatives at
// co-ordinate direction
std::vector<std::vector<double> > v(dim, std::vector<double> (2));
for (unsigned int d=0; d<dim; ++d)
- polynomials[(i/n_pols_to[d])%n_pols].value(p(d), v[d]);
+ polynomials[indices[d]].value(p(d), v[d]);
Tensor<1,dim> grad;
for (unsigned int d=0; d<dim; ++d)
- grad[d]=1.;
-
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int x=0; x<dim; ++x)
- grad[d]*=v[x][d==x];
+ {
+ grad[d] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grad[d] *= v[x][d==x];
+ }
return grad;
}
template <int dim>
Tensor<2,dim>
-TensorProductPolynomials<dim>::compute_grad_grad(const unsigned int i,
- const Point<dim> &p) const
+TensorProductPolynomials<dim>::compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const
{
- const unsigned int n_pols=polynomials.size();
-
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
std::vector<std::vector<double> > v(dim, std::vector<double> (3));
for (unsigned int d=0; d<dim; ++d)
- polynomials[(i/n_pols_to[d])%n_pols].value(p(d), v[d]);
+ polynomials[indices[d]].value(p(d), v[d]);
Tensor<2,dim> grad_grad;
-
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int d2=0; d2<dim; ++d2)
- grad_grad[d1][d2]=1.;
-
- for (unsigned int x=0; x<dim; ++x)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- {
- unsigned int derivative=0;
- if (d1==x || d2==x)
- {
- if (d1==d2)
- derivative=2;
- else
- derivative=1;
- }
- grad_grad[d1][d2]*=v[x][derivative];
- }
+ {
+ grad_grad[d1][d2] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x || d2==x)
+ {
+ if (d1==d2)
+ derivative=2;
+ else
+ derivative=1;
+ }
+ grad_grad[d1][d2] *= v[x][derivative];
+ }
+ }
return grad_grad;
}
template <int dim>
-void TensorProductPolynomials<dim>::compute(
- const Point<dim> &p,
- std::vector<double> &values,
- std::vector<Tensor<1,dim> > &grads,
- std::vector<Tensor<2,dim> > &grad_grads) const
+void
+TensorProductPolynomials<dim>::
+compute (const Point<dim> &p,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &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<double> > v(dim, n_pols);
- for (unsigned int d=0; d<v.size()[0]; ++d)
- for (unsigned int i=0; i<v.size()[1]; ++i)
+ // compute the values (and
+ // derivatives, if necessary) of
+ // all polynomials at this
+ // evaluation point
+ Table<2,std::vector<double> > v(dim, polynomials.size());
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<polynomials.size(); ++i)
{
- v(d,i).resize (v_size, 0.);
+ v(d,i).resize (n_values_and_derivatives, 0.);
polynomials[i].value(p(d), v(d,i));
};
- if (update_values)
+ for (unsigned int i=0; i<n_tensor_pols; ++i)
{
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- values[i]=1;
+ // first get the
+ // one-dimensional indices of
+ // this particular tensor
+ // product polynomial
+ unsigned int indices[dim];
+ compute_index (i, indices);
- for (unsigned int x=0; x<dim; ++x)
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- values[i]*=v[x][(i/n_pols_to[x])%n_pols][0];
- }
+ if (update_values)
+ {
+ values[i] = 1;
+ for (unsigned int x=0; x<dim; ++x)
+ values[i] *= v(x,indices[x])[0];
+ }
- if (update_grads)
- {
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- for (unsigned int d=0; d<dim; ++d)
- grads[i][d]=1.;
-
- for (unsigned int x=0; x<dim; ++x)
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- for (unsigned int d=0; d<dim; ++d)
- grads[i][d]*=v[x][(i/n_pols_to[x])%n_pols][d==x];
- }
-
- if (update_grad_grads)
- {
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- grad_grads[i][d1][d2]=1.;
-
- for (unsigned int x=0; x<dim; ++x)
- for (unsigned int i=0; i<n_tensor_pols; ++i)
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int d2=0; d2<dim; ++d2)
- {
- unsigned int derivative=0;
- if (d1==x || d2==x)
- {
- if (d1==d2)
- derivative=2;
- else
- derivative=1;
- }
- grad_grads[i][d1][d2]*=
- v[x][(i/n_pols_to[x])%n_pols][derivative];
- }
+ if (update_grads)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ grads[i][d] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grads[i][d] *= v(x,indices[x])[d==x];
+ }
+
+ if (update_grad_grads)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ grad_grads[i][d1][d2] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x || d2==x)
+ {
+ if (d1==d2)
+ derivative=2;
+ else
+ derivative=1;
+ }
+ grad_grads[i][d1][d2]
+ *= v(x,indices[x])[derivative];
+ }
+ }
}
}
+
template<int dim>
unsigned int
TensorProductPolynomials<dim>::n() const
projects / dealii-svn.git / commitdiff