From: guido Date: Mon, 8 Jul 2002 18:06:28 +0000 (+0000) Subject: Legendre polynomials are now orthoNORMAL on [0,1] X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=18043b4519777b22291bd81ef1b60555f2c1306d;p=dealii-svn.git Legendre polynomials are now orthoNORMAL on [0,1] git-svn-id: https://svn.dealii.org/trunk@6225 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index f15e4a2c26..7b3c298cb5 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -91,6 +91,50 @@ class Polynomial : public Subscriptor * separately. */ unsigned int degree () const; + + /** + * Scale the abscissa of the + * polynomial. Given the + * polynomial $p(t)$ and the + * scaling $t = ax$, then the + * result of this operation is + * the polynomial $q$, such that + * $q(x) = p(t)$. + * + * The operation is performed in + * place. + */ + void scale (const number factor); + + /** + * Shift the abscissa oft the + * polynomial. Given the + * polynomial $p(t)$ and the + * shift $t = x + a$, then the + * result of this operation is + * the polynomial $q$, such that + * $q(x) = p(t)$. + * + * The template parameter allows + * to compute the new + * coefficients with higher + * accuracy, since all + * computations are performed + * with type @p{number2}. This + * may be necessary, since this + * operation involves a big + * number of additions. On a Sun + * Sparc Ultra with Solaris 2.8, + * the difference between + * @p{double} and @p{long double} + * was not significant, though. + * + * + * The operation is performed in + * place. + */ + template + void shift (const number2 offset); /** * Exception @@ -104,6 +148,25 @@ class Polynomial : public Subscriptor protected: + /** + * This function performs the actual scaling. + */ + static void scale (typename std::vector& coefficients, + const number factor); + + /** + * This function performs the actual shift + */ + template + static void shift (typename std::vector& coefficients, + const number2 shift); + + /** + * Multiply polynomial by a factor. + */ + static void multiply (typename std::vector& coefficients, + const number factor); + /** * Coefficients of the polynomial * $\sum_i a_i x^i$. This vector @@ -192,12 +255,14 @@ class LagrangeEquidistant: public Polynomial /** - * Legendre polynomials of arbitrary order on @p{[-1,1]}. + * Legendre polynomials of arbitrary order on @p{[0,1]}. * * Constructing a Legendre polynomial of order @p{k}, the coefficients * will be computed by the three-term recursion formula. The * coefficients are stored in a static data vector to be available - * when needed next time. + * when needed next time. Since the recursion is performed for the + * interval $[-1,1]$, the polynomials are shifted to $[0,1]$ by the + * @p{scale} and @p{shift} functions of @p{Polynomial}, afterwards. * * @author Guido Kanschat, 2000 */ @@ -228,6 +293,11 @@ class Legendre : public Polynomial generate_complete_basis (const unsigned int degree); private: + /** + * Coefficients for the interval $[0,1]$. + */ + static typename std::vector *> shifted_coefficients; + /** * Vector with already computed * coefficients. For each degree @@ -238,7 +308,7 @@ class Legendre : public Polynomial * vectors in order to simplify * programming multithread-safe. */ - static typename std::vector *> coefficients; + static typename std::vector *> recursive_coefficients; /** * Compute coefficients recursively. diff --git a/deal.II/base/source/legendre.cc b/deal.II/base/source/legendre.cc index 2344a45f0f..120e85376d 100644 --- a/deal.II/base/source/legendre.cc +++ b/deal.II/base/source/legendre.cc @@ -30,8 +30,12 @@ // for the start. template typename std::vector *> -Legendre::coefficients(20, - static_cast*>(0)); +Legendre::recursive_coefficients( + 20, static_cast*>(0)); +template +typename std::vector *> +Legendre::shifted_coefficients( + 20, static_cast*>(0)); // have a lock that guarantees that at most one thread is changing and @@ -60,13 +64,13 @@ Legendre::compute_coefficients (const unsigned int k_) k=1; // check: does the information // already exist? - if ((coefficients.size() < k+1) || - ((coefficients.size() >= k+1) && - (coefficients[k] == 0))) + if ((recursive_coefficients.size() < k+1) || + ((recursive_coefficients.size() >= k+1) && + (recursive_coefficients[k] == 0))) // no, then generate the // respective coefficients { - coefficients.resize (k+1, 0); + recursive_coefficients.resize (k+1, 0); if (k<=1) { @@ -87,8 +91,20 @@ Legendre::compute_coefficients (const unsigned int k_) // now make these arrays // const - coefficients[0] = c0; - coefficients[1] = c1; + recursive_coefficients[0] = c0; + recursive_coefficients[1] = c1; + // Compute polynomials + // orthogonal on [0,1] + c0 = new std::vector(*c0); + c1 = new std::vector(*c1); + + Polynomial::shift(*c0, (long double) -1.); + Polynomial::scale(*c0, 2.); + Polynomial::shift(*c1, (long double) -1.); + Polynomial::scale(*c1, 2.); + Polynomial::multiply(*c1, sqrt(3.)); + shifted_coefficients[0]=c0; + shifted_coefficients[1]=c1; } else { @@ -110,19 +126,26 @@ Legendre::compute_coefficients (const unsigned int k_) const number b = a*(2*k-1); const number c = a*(k-1); - (*ck)[k] = b*(*coefficients[k-1])[k-1]; - (*ck)[k-1] = b*(*coefficients[k-1])[k-2]; + (*ck)[k] = b*(*recursive_coefficients[k-1])[k-1]; + (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2]; for (unsigned int i=1 ; i<= k-2 ; ++i) - (*ck)[i] = b*(*coefficients[k-1])[i-1] - -c*(*coefficients[k-2])[i]; + (*ck)[i] = b*(*recursive_coefficients[k-1])[i-1] + -c*(*recursive_coefficients[k-2])[i]; - (*ck)[0] = -c*(*coefficients[k-2])[0]; + (*ck)[0] = -c*(*recursive_coefficients[k-2])[0]; // finally assign the newly // created vector to the // const pointer in the // coefficients array - coefficients[k] = ck; + recursive_coefficients[k] = ck; + // and compute the + // coefficients for [0,1] + ck = new std::vector(*ck); + shift(*ck,(long double) -1.); + Polynomial::scale(*ck, 2.); + Polynomial::multiply(*ck, sqrt(2.*k+1.)); + shifted_coefficients[k] = ck; }; }; @@ -145,7 +168,7 @@ Legendre::get_coefficients (const unsigned int k) // of coefficients. do that in a MT // safe way coefficients_lock.acquire (); - const std::vector *p = coefficients[k]; + const std::vector *p = shifted_coefficients[k]; coefficients_lock.release (); // return the object pointed diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index 71b4269030..cb93b56208 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -13,6 +13,7 @@ #include +#include template Polynomial::Polynomial (const typename std::vector &a): @@ -91,6 +92,108 @@ Polynomial::value (const number x, } +template +void +Polynomial::scale(typename std::vector& coefficients, + const number factor) +{ + double f = 1.; + for (typename std::vector::iterator c = coefficients.begin(); + c != coefficients.end(); ++c) + { + *c *= f; + f *= factor; + } +} + + + +template +void +Polynomial::scale(const number factor) +{ + scale (coefficients, factor); +} + + + +template +void +Polynomial::multiply(typename std::vector& coefficients, + const number factor) +{ + for (typename std::vector::iterator c = coefficients.begin(); + c != coefficients.end(); ++c) + *c *= factor; +} + + + +template +template +void +Polynomial::shift(typename std::vector& coefficients, + const number2 offset) +{ + // Copy coefficients to a vector of + // accuracy given by the argument + std::vector new_coefficients(coefficients.size()); + new_coefficients.assign(coefficients.begin(), coefficients.end()); + + // Traverse all coefficients from + // c_1. c_0 will be modified by + // higher degrees, only. + for (unsigned int d=1; d +template +void +Polynomial::shift(const number2 offset) +{ + shift(coefficients, offset); +} + + // ------------------------------------------------------------ // @@ -348,3 +451,10 @@ generate_complete_basis (const unsigned int degree) template class Polynomial; template class Polynomial; template class Polynomial; + +template void Polynomial::shift(const float offset); +template void Polynomial::shift(const double offset); +template void Polynomial::shift(const long double offset); +template void Polynomial::shift(const double offset); +template void Polynomial::shift(const long double offset); +template void Polynomial::shift(const long double offset); diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index 2bbac3b8b6..f0b8a36fec 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -627,10 +627,6 @@ class SparseMatrixEZ : public Subscriptor */ unsigned int memory_consumption () const; - /** - * Exception - */ - DeclException0 (ExcMatrixNotInitialized); /** * Exception */ @@ -638,10 +634,12 @@ class SparseMatrixEZ : public Subscriptor int, int, << "The entry with index <" << arg1 << ',' << arg2 << "> does not exist."); + /** * Exception */ DeclException0 (ExcMatrixNotSquare); + /** * Exception */ diff --git a/tests/base/polynomial1d.cc b/tests/base/polynomial1d.cc index 89d864b0b1..efc0ae8c31 100644 --- a/tests/base/polynomial1d.cc +++ b/tests/base/polynomial1d.cc @@ -24,15 +24,15 @@ double scalar_product (const Polynomial& p1, const Polynomial& p2) { unsigned int degree = (p1.degree() + p2.degree())/2 + 1; - QGauss<1> gauss(degree); + QGauss<1> gauss(degree+3); double sum = 0.; for (unsigned int i=0;i > p; + std::vector > q; deallog << "Legendre" << std::endl; - for (unsigned int i=0;i<15;++i) + for (unsigned int i=0;i<12;++i) p.push_back (Legendre(i)); for (unsigned int i=0;i p2(2); + plot_shape_functions(m, p2, "DGP2"); + plot_face_shape_functions(m, p2, "DGP2"); + test_compute_functions(m, p2, "DGP2"); + FE_DGP p3(3); + plot_shape_functions(m, p3, "DGP3"); + plot_face_shape_functions(m, p3, "DGP3"); + test_compute_functions(m, p3, "DGP3"); + FE_DGP p4(4); + plot_shape_functions(m, p4, "DGP4"); + plot_face_shape_functions(m, p4, "DGP4"); + test_compute_functions(m, p4, "DGP4"); +// FE_DGP p5(5); +// plot_shape_functions(m, p5, "DGP5"); +// FE_DGP p6(6); +// plot_shape_functions(m, p6, "DGP6"); +// FE_DGP p7(7); +// plot_shape_functions(m, p7, "DGP7"); +// FE_DGP p8(8); +// plot_shape_functions(m, p8, "DGP8"); +// FE_DGP p9(9); +// plot_shape_functions(m, p9, "DGP9"); +// FE_DGP p10(10); +// plot_shape_functions(m, p10, "DGP10"); +} + + int main() { @@ -398,7 +434,8 @@ main() plot_FE_Q_shape_functions<1>(); plot_FE_Q_shape_functions<2>(); - plot_FE_DGQ_shape_functions<2>(); +// plot_FE_DGP_shape_functions<1>(); + plot_FE_DGP_shape_functions<2>(); // plot_FE_Q_shape_functions<3>(); // FESystem test.