//---------------------------------------------------------------------------
-// $Id$
+// $Id$
// Version: $Name$
//
-// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
// to be a problem since we only need it on very rare occasions. if
// someone finds this is a bottleneck, feel free to replace it by a
// more fine-grained solution
-namespace
+namespace
{
Threads::ThreadMutex coefficients_lock;
}
Assert (coefficients.size() > 0, ExcEmptyObject());
Assert (values.size() > 0, ExcZero());
const unsigned int values_size=values.size();
-
-
+
+
// if we only need the value, then
// call the other function since
// that is significantly faster
{
*c *= f;
f *= factor;
- }
+ }
}
return *this;
}
-
+
template <typename number>
Polynomial<number>&
Polynomial<number>::operator *= (const Polynomial<number>& p)
unsigned int new_degree = this->degree() + p.degree();
std::vector<number> new_coefficients(new_degree+1, 0.);
-
+
for (unsigned int i=0; i<p.coefficients.size(); ++i)
for (unsigned int j=0; j<this->coefficients.size(); ++j)
new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i];
this->coefficients = new_coefficients;
-
+
return *this;
}
-
+
template <typename number>
Polynomial<number>&
Polynomial<number>::operator += (const Polynomial<number>& p)
return *this;
}
-
+
template <typename number>
Polynomial<number>&
Polynomial<number>::operator -= (const Polynomial<number>& p)
return *this;
}
-
+
template <typename number>
template <typename number2>
void
Polynomial<number>::shift(std::vector<number>& coefficients,
const number2 offset)
- {
+ {
#ifdef DEAL_II_LONG_DOUBLE_LOOP_BUG
AssertThrow (false,
ExcMessage("Sorry, but the compiler you are using has a bug that disallows "
// args. note that this code is
// actually unreachable
coefficients[0] = offset;
-#else
+#else
// Copy coefficients to a vector of
// accuracy given by the argument
std::vector<number2> new_coefficients(coefficients.begin(),
coefficients.end());
-
+
// Traverse all coefficients from
// c_1. c_0 will be modified by
// higher degrees, only.
// needed and computed
// successively.
number2 offset_power = offset;
-
+
// Compute (x+offset)^d
// and modify all values c_k
// with k<d.
}
-
+
template <typename number>
Polynomial<number>
Polynomial<number>::derivative () const
return Polynomial<number> (newcoefficients);
}
-
+
template <typename number>
Polynomial<number>
return Polynomial<number> (newcoefficients);
}
-
+
template <typename number>
void
result[n] = coefficient;
return result;
}
-
-
+
+
template <typename number>
Monomial<number>::Monomial (unsigned int n,
double coefficient)
: Polynomial<number>(make_vector(n, coefficient))
{}
-
-
+
+
template <typename number>
std::vector<Polynomial<number> >
Monomial<number>::generate_complete_basis (const unsigned int degree)
- std::vector<double>
+ std::vector<double>
LagrangeEquidistant::compute_coefficients (const unsigned int n,
const unsigned int support_point)
{
Assert(support_point<n_functions,
ExcIndexRange(support_point, 0, n_functions));
double const *x=0;
-
+
switch (n)
{
case 1:
0.0, 1.0
};
x=&x1[0];
- break;
+ break;
}
case 2:
{
0.0, -1.0, 22.0/3.0, -16.0, 32.0/3.0
};
x=&x4[0];
- break;
+ break;
}
case 5:
{
Assert(x!=0, ExcInternalError());
for (unsigned int i=0; i<n_functions; ++i)
a[i]=x[support_point*n_functions+i];
-
+
return a;
}
//----------------------------------------------------------------------//
-
+
std::vector<Polynomial<double> >
Lagrange::generate_complete_basis (const std::vector<Point<1> >& points)
{
std::vector<double> linear(2, 1.);
// We start with a constant polynomial
std::vector<double> one(1, 1.);
-
+
for (unsigned int i=0;i<p.size();++i)
{
// Construct interpolation formula
}
}
}
-
+
return p;
}
-
+
// ------------------ class Legendre --------------- //
{}
-
+
void
Legendre::compute_coefficients (const unsigned int k_)
{
#else
typedef long double SHIFT_TYPE;
#endif
-
+
unsigned int k = k_;
// first make sure that no other
// respective coefficients
{
recursive_coefficients.resize (k+1, 0);
-
+
if (k<=1)
{
// create coefficients
// orthogonal on [0,1]
c0 = new std::vector<double>(*c0);
c1 = new std::vector<double>(*c1);
-
+
Polynomial<double>::shift<SHIFT_TYPE> (*c0, -1.);
Polynomial<double>::scale(*c0, 2.);
Polynomial<double>::shift<SHIFT_TYPE> (*c1, -1.);
coefficients_lock.acquire ();
std::vector<double> *ck = new std::vector<double>(k+1);
-
+
const double a = 1./(k);
const double b = a*(2*k-1);
const double c = a*(k-1);
-
+
(*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)
std::vector<double> coefficients (3);
coefficients[0] = 0.0;
- coefficients[1] = -1.0 * std::sqrt (3);
- coefficients[2] = std::sqrt (3);
+ coefficients[1] = -1.0 * std::sqrt (3.);
+ coefficients[2] = std::sqrt (3.);
return coefficients;
}
std::vector<double> legendre_coefficients_tmp1 (p);
std::vector<double> legendre_coefficients_tmp2 (p - 1);
- coefficients[0] = -1.0 * std::sqrt (3);
- coefficients[1] = 2.0 * std::sqrt (3);
+ coefficients[0] = -1.0 * std::sqrt (3.);
+ coefficients[1] = 2.0 * std::sqrt (3.);
legendre_coefficients_tmp1[0] = 1.0;
for (unsigned int i = 2; i < p; ++i) {
for (unsigned int j = 0; j < i; ++j)
legendre_coefficients_tmp1[j] = coefficients[j];
- coefficients[0] = std::sqrt (2 * i + 1) * ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] / std::sqrt (2 * i - 1) + (1.0 - i) * legendre_coefficients_tmp2[0] / std::sqrt (2 * i - 3)) / i;
+ coefficients[0] = std::sqrt (2 * i + 1.) * ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] / std::sqrt (2 * i - 1.) + (1.0 - i) * legendre_coefficients_tmp2[0] / std::sqrt (2 * i - 3.)) / i;
for (unsigned int j = 1; j < i - 1; ++j)
- coefficients[j] = std::sqrt (2 * i + 1) * (std::sqrt (2 * i - 1) * (2.0 * legendre_coefficients_tmp1[j - 1] - legendre_coefficients_tmp1[j]) + (1.0 - i) * legendre_coefficients_tmp2[j] / std::sqrt (2 * i - 3)) / i;
+ coefficients[j] = std::sqrt (2 * i + 1.) * (std::sqrt (2 * i - 1.) * (2.0 * legendre_coefficients_tmp1[j - 1] - legendre_coefficients_tmp1[j]) + (1.0 - i) * legendre_coefficients_tmp2[j] / std::sqrt (2 * i - 3.)) / i;
- coefficients[i - 1] = std::sqrt (4 * i * i - 1) * (2.0 * legendre_coefficients_tmp1[i - 2] - legendre_coefficients_tmp1[i - 1]) / i;
- coefficients[i] = 2.0 * std::sqrt (4 * i * i - 1) * legendre_coefficients_tmp1[i - 1] / i;
+ coefficients[i - 1] = std::sqrt (4 * i * i - 1.) * (2.0 * legendre_coefficients_tmp1[i - 2] - legendre_coefficients_tmp1[i - 1]) / i;
+ coefficients[i] = 2.0 * std::sqrt (4 * i * i - 1.) * legendre_coefficients_tmp1[i - 1] / i;
}
for (int i = p; i > 0; --i)
// until we quit this function
Threads::ThreadMutex::ScopedLock lock(coefficients_lock);
- // The first 2 coefficients
+ // 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.size() >= k+1) &&
(recursive_coefficients[k] == 0)) )
// no, then generate the
// respective coefficients
{
recursive_coefficients.resize (k+1, 0);
-
+
if (k<=1)
{
// create coefficients
(*c2)[0] = 0.*a;
(*c2)[1] = -4.*a;
(*c2)[2] = 4.*a;
-
+
recursive_coefficients[2] = c2;
}
else
coefficients_lock.acquire ();
std::vector<double> *ck = new std::vector<double>(k+1);
-
+
const double a = 1.; //1./(2.*k);
(*ck)[0] = - a*(*recursive_coefficients[k-1])[0];
-
+
for (unsigned int i=1; i<=k-1; ++i)
(*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1]
- (*recursive_coefficients[k-1])[i] );
-
+
(*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1];
// for even degrees, we need
// to add a multiple of
(*ck)[1] += b*(*recursive_coefficients[2])[1];
(*ck)[2] += b*(*recursive_coefficients[2])[2];
- }
+ }
// finally assign the newly
// created vector to the
// const pointer in the