From: brian Date: Tue, 15 Oct 2002 23:27:02 +0000 (+0000) Subject: Added Hierarchical polynomial class. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3c1f50d12869e855370f1755c37ecbca1ddc2d9d;p=dealii-svn.git Added Hierarchical polynomial class. git-svn-id: https://svn.dealii.org/trunk@6659 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index 61b866cf3a..69b83ed706 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -678,10 +678,170 @@ namespace Polynomials v.push_back (Legendre(i)); return v; }; - + + + +// ------------------ class Hierarchical --------------- // + + + +// Reserve space for polynomials up to degree 19. Should be sufficient +// for the start. + template + std::vector *> + Hierarchical::recursive_coefficients( + 20, static_cast*>(0)); + + + + template + void + Hierarchical::compute_coefficients (const unsigned int k_) + { + unsigned int k = k_; + + // first make sure that no other + // thread intercepts the operation + // of this function + coefficients_lock.acquire (); + + // The first 2 coefficients + // are hard-coded + if (k==0) + k=1; + // check: does the information + // already exist? + if ( (recursive_coefficients.size() < k+1) || + ((recursive_coefficients.size() >= k+1) && (recursive_coefficients[k] == 0)) ) + // no, then generate the + // respective coefficients + { + recursive_coefficients.resize (k+1, 0); + + if (k<=1) + { + // create coefficients + // vectors for k=0 and k=1 + // + // allocate the respective + // amount of memory and + // later assign it to the + // coefficients array to + // make it const + std::vector *c0 = new std::vector(2); + (*c0)[0] = 1.0; + (*c0)[1] = -1.0; + + std::vector *c1 = new std::vector(2); + (*c1)[0] = 0.0; + (*c1)[1] = 1.0; + + // now make these arrays + // const + recursive_coefficients[0] = c0; + recursive_coefficients[1] = c1; + } + else if (k==2) + { + compute_coefficients(1); + + std::vector *c2 = new std::vector(3); + + (*c2)[0] = 0.; + (*c2)[1] = -4.; + (*c2)[2] = 4.; + + recursive_coefficients[2] = c2; + } + else + { + // for larger numbers, + // compute the coefficients + // recursively. to do so, + // we have to release the + // lock temporarily to + // allow the called + // function to acquire it + // itself + coefficients_lock.release (); + compute_coefficients(k-1); + coefficients_lock.acquire (); + + std::vector *ck = new std::vector(k+1); + + (*ck)[0] = -(*recursive_coefficients[k-1])[0]; + + for (unsigned int i=1; i<=k-1; ++i) + (*ck)[i] = ( 2.*(*recursive_coefficients[k-1])[i-1] + - (*recursive_coefficients[k-1])[i] ); + + (*ck)[k] = 2.*(*recursive_coefficients[k-1])[k-1]; + // for even degrees, we need + // to add a multiple of + // basis fcn phi_2 + if ( (k%2) == 0 ) + { + (*ck)[1] += (*recursive_coefficients[2])[1]; + (*ck)[2] += (*recursive_coefficients[2])[2]; + } + // finally assign the newly + // created vector to the + // const pointer in the + // coefficients array + recursive_coefficients[k] = ck; + }; + }; + + // now, everything is done, so + // release the lock again + coefficients_lock.release (); + }; + + + + template + const typename std::vector & + Hierarchical::get_coefficients (const unsigned int k) + { + // first make sure the coefficients + // get computed if so necessary + compute_coefficients (k); + + // then get a pointer to the array + // of coefficients. + const std::vector *p = recursive_coefficients[k]; + + // return the object pointed + // to. since this object does not + // change any more once computed, + // this is MT safe + return *p; + }; + + + + template + Hierarchical::Hierarchical (const unsigned int k) + : + Polynomial (get_coefficients(k)) + {} + + + + template + std::vector > + Hierarchical::generate_complete_basis (const unsigned int degree) + { + std::vector > v; + v.reserve(degree+1); + for (unsigned int i=0; i<=degree; ++i) + v.push_back (Hierarchical(i)); + return v; + }; } + // ------------------ explicit instantiations --------------- // namespace Polynomials @@ -698,4 +858,5 @@ namespace Polynomials template void Polynomial::shift(const long double offset); template class Legendre; + template class Hierarchical; } diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index 774be5c151..718fbf7f63 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -167,6 +167,14 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

base

    +
  1. + New: Added Hierarchical Polynomial (similar to Legendre class). Will + eventually be used in a hierarchical FiniteElement class similar to + FE_Q class. Included in Polynomials namespace. +
    + (Brian Carnes 2002/10/15) +

    +
  2. Changed: Because they became too many, the classes describing 1d polynomials are now in a namespace Polynomials. @@ -640,6 +648,7 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
    (WB 2002/06/04)

    +