From: Guido Kanschat <dr.guido.kanschat@gmail.com> Date: Mon, 8 Jul 2002 13:39:43 +0000 (+0000) Subject: missing functions in PolynomialSpace coded X-Git-Tag: v8.0.0~17753 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2969c26fed6f451d3fad5603ac52ec2ce8377498;p=dealii.git missing functions in PolynomialSpace coded git-svn-id: https://svn.dealii.org/trunk@6223 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h index 2cd7a37a1f..586642a612 100644 --- a/deal.II/base/include/base/polynomial_space.h +++ b/deal.II/base/include/base/polynomial_space.h @@ -152,6 +152,20 @@ class PolynomialSpace */ const unsigned int n_pols; + /** + * Compute numbers in x, y and z + * direction. Given an index + * @p{n} in the d-dimensional + * polynomial space, compute the + * indices i,j,k such that + * @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; + /** * Static function used in the * constructor to compute the diff --git a/deal.II/base/source/polynomial_space.cc b/deal.II/base/source/polynomial_space.cc index 460a9f9d9d..1c01c060b1 100644 --- a/deal.II/base/source/polynomial_space.cc +++ b/deal.II/base/source/polynomial_space.cc @@ -30,34 +30,103 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n) }; +template <int dim> +void +PolynomialSpace<dim>::compute_index(unsigned int n, + unsigned int& nx, + unsigned int& ny, + unsigned int& nz) const +{ + 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; + } +} + 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 { - Assert(false, ExcNotImplemented()); - return 0.; + 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)); + 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 { - Assert(false, ExcNotImplemented()); - return Tensor<1,dim>(); + unsigned int ix[3]; + compute_index(i,ix[0],ix[1],ix[2]); + + Tensor<1,dim> result; + for (unsigned int d=0;d<dim;++d) + result[d] = 1.; + + std::vector<double> v(2); + 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) + if (d1 != d) + result[d1] *= v[0]; + } + return result; } 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 { - Assert(false, ExcNotImplemented()); - return Tensor<2,dim>(); + unsigned int ix[3]; + compute_index(i,ix[0],ix[1],ix[2]); + + Tensor<2,dim> result; + for (unsigned int d=0;d<dim;++d) + for (unsigned int d1=0;d1<dim;++d1) + result[d][d1] = 1.; + + std::vector<double> v(3); + 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) + { + if (d1 != d) + { + result[d][d1] *= v[1]; + result[d1][d] *= v[1]; + for (unsigned int d2=0;d2<dim;++d2) + if (d2 != d) + result[d1][d2] *= v[0]; + } + } + } + return result; }