+
+/**
+ * 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 <int dim>
+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<std::vector<Polynomials::Polynomial<double> > > &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<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
+ * @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<dim> &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<dim> &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<dim> &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<std::vector<Polynomials::Polynomial<double> > > 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<std::vector<Polynomials::Polynomial<double> > > &pols);
+};
+
+
+
+
/* -------------- declaration of explicit specializations --- */
template <>
+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
+/* ------------------- TensorProductPolynomials -------------- */
+
template <>
void
TensorProductPolynomials<1>::
return n_tensor_pols;
}
+
+
+
+/* ------------------- AnisotropicPolynomials -------------- */
+
+
+template <int dim>
+AnisotropicPolynomials<dim>::
+AnisotropicPolynomials(const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
+ :
+ polynomials (pols),
+ n_tensor_pols(get_n_tensor_pols(pols))
+{
+ Assert (pols.size() == dim, ExcInvalidDim(pols.size()));
+ for (unsigned int d=0; d<dim; ++d)
+ Assert (pols[d].size() > 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<polynomials[0].size(), ExcInternalError());
+ indices[0] = i;
+}
+
+
+
+template <>
+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 <int dim>
+double
+AnisotropicPolynomials<dim>::compute_value (const unsigned int i,
+ const Point<dim> &p) const
+{
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ double value=1.;
+ for (unsigned int d=0; d<dim; ++d)
+ value *= polynomials[d][indices[d]].value(p(d));
+
+ return value;
+}
+
+
+template <int dim>
+Tensor<1,dim>
+AnisotropicPolynomials<dim>::compute_grad (const unsigned int i,
+ const Point<dim> &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<std::vector<double> > v(dim, std::vector<double> (2));
+ for (unsigned int d=0; d<dim; ++d)
+ polynomials[d][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 x=0; x<dim; ++x)
+ grad[d] *= v[x][d==x];
+ }
+
+ return grad;
+}
+
+
+template <int dim>
+Tensor<2,dim>
+AnisotropicPolynomials<dim>::compute_grad_grad (const unsigned int i,
+ const Point<dim> &p) const
+{
+ 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[d][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)
+ {
+ 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
+AnisotropicPolynomials<dim>::
+compute (const Point<dim> &p,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &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<std::vector<std::vector<double> > > v(dim);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ v[d].resize (polynomials[d].size());
+ for (unsigned int i=0; i<polynomials[d].size(); ++i)
+ {
+ v[d][i].resize (n_values_and_derivatives, 0.);
+ polynomials[d][i].value(p(d), v[d][i]);
+ };
+ }
+ for (unsigned int i=0; i<n_tensor_pols; ++i)
+ {
+ // first get the
+ // one-dimensional indices of
+ // this particular tensor
+ // product polynomial
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ 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 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 ? 1 : 0];
+ }
+
+ 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
+AnisotropicPolynomials<dim>::n() const
+{
+ return n_tensor_pols;
+}
+
+
+template <int dim>
+unsigned int
+AnisotropicPolynomials<dim>::
+get_n_tensor_pols (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
+{
+ unsigned int y = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ y *= pols[d].size();
+ return y;
+}
+
+
+
+/* ------------------- explicit instantiations -------------- */
template class TensorProductPolynomials<1>;
template class TensorProductPolynomials<2>;
template class TensorProductPolynomials<3>;
+
+template class AnisotropicPolynomials<1>;
+template class AnisotropicPolynomials<2>;
+template class AnisotropicPolynomials<3>;
projects / dealii-svn.git / commitdiff