* 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 <typename number2>
+ void shift (const number2 offset);
/**
* Exception
protected:
+ /**
+ * This function performs the actual scaling.
+ */
+ static void scale (typename std::vector<number>& coefficients,
+ const number factor);
+
+ /**
+ * This function performs the actual shift
+ */
+ template <typename number2>
+ static void shift (typename std::vector<number>& coefficients,
+ const number2 shift);
+
+ /**
+ * Multiply polynomial by a factor.
+ */
+ static void multiply (typename std::vector<number>& coefficients,
+ const number factor);
+
/**
* Coefficients of the polynomial
* $\sum_i a_i x^i$. This vector
/**
- * 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
*/
generate_complete_basis (const unsigned int degree);
private:
+ /**
+ * Coefficients for the interval $[0,1]$.
+ */
+ static typename std::vector<const typename std::vector<number> *> shifted_coefficients;
+
/**
* Vector with already computed
* coefficients. For each degree
* vectors in order to simplify
* programming multithread-safe.
*/
- static typename std::vector<const typename std::vector<number> *> coefficients;
+ static typename std::vector<const typename std::vector<number> *> recursive_coefficients;
/**
* Compute coefficients recursively.
// for the start.
template <typename number>
typename std::vector<const typename std::vector<number> *>
-Legendre<number>::coefficients(20,
- static_cast<const typename std::vector<number>*>(0));
+Legendre<number>::recursive_coefficients(
+ 20, static_cast<const typename std::vector<number>*>(0));
+template <typename number>
+typename std::vector<const typename std::vector<number> *>
+Legendre<number>::shifted_coefficients(
+ 20, static_cast<const typename std::vector<number>*>(0));
// have a lock that guarantees that at most one thread is changing and
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)
{
// 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<number>(*c0);
+ c1 = new std::vector<number>(*c1);
+
+ Polynomial<number>::shift(*c0, (long double) -1.);
+ Polynomial<number>::scale(*c0, 2.);
+ Polynomial<number>::shift(*c1, (long double) -1.);
+ Polynomial<number>::scale(*c1, 2.);
+ Polynomial<number>::multiply(*c1, sqrt(3.));
+ shifted_coefficients[0]=c0;
+ shifted_coefficients[1]=c1;
}
else
{
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<number>(*ck);
+ shift(*ck,(long double) -1.);
+ Polynomial<number>::scale(*ck, 2.);
+ Polynomial<number>::multiply(*ck, sqrt(2.*k+1.));
+ shifted_coefficients[k] = ck;
};
};
// of coefficients. do that in a MT
// safe way
coefficients_lock.acquire ();
- const std::vector<number> *p = coefficients[k];
+ const std::vector<number> *p = shifted_coefficients[k];
coefficients_lock.release ();
// return the object pointed
#include <base/polynomial.h>
+#include <base/exceptions.h>
template <typename number>
Polynomial<number>::Polynomial (const typename std::vector<number> &a):
}
+template <typename number>
+void
+Polynomial<number>::scale(typename std::vector<number>& coefficients,
+ const number factor)
+{
+ double f = 1.;
+ for (typename std::vector<number>::iterator c = coefficients.begin();
+ c != coefficients.end(); ++c)
+ {
+ *c *= f;
+ f *= factor;
+ }
+}
+
+
+
+template <typename number>
+void
+Polynomial<number>::scale(const number factor)
+{
+ scale (coefficients, factor);
+}
+
+
+
+template <typename number>
+void
+Polynomial<number>::multiply(typename std::vector<number>& coefficients,
+ const number factor)
+{
+ for (typename std::vector<number>::iterator c = coefficients.begin();
+ c != coefficients.end(); ++c)
+ *c *= factor;
+}
+
+
+
+template <typename number>
+template <typename number2>
+void
+Polynomial<number>::shift(typename std::vector<number>& coefficients,
+ const number2 offset)
+{
+ // Copy coefficients to a vector of
+ // accuracy given by the argument
+ std::vector<number2> 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<new_coefficients.size(); ++d)
+ {
+ const unsigned int n = d;
+ // Binomial coefficients are
+ // needed for the
+ // computation. The rightmost
+ // value is unity.
+ unsigned int binomial_coefficient = 1;
+
+ // Powers of the offset will be
+ // needed and computed
+ // successively.
+ number2 offset_power = offset;
+
+ // Compute (x+offset)^d
+ // and modify all values c_k
+ // with k<d.
+ // The coefficient in front of
+ // x^d is not modified in this step.
+ for (unsigned int k=0;k<d;++k)
+ {
+ // Recursion from Bronstein
+ // Make sure no remainders
+ // occur in integer
+ // division.
+ binomial_coefficient = (binomial_coefficient*(n-k))/(k+1);
+
+ new_coefficients[d-k-1] += new_coefficients[d]
+ * binomial_coefficient
+ * offset_power;
+ offset_power *= offset;
+ }
+ // The binomial coefficient
+ // should have gone through a
+ // whole row of Pascal's
+ // triangle.
+ Assert (binomial_coefficient == 1, ExcInternalError());
+ }
+ coefficients.assign(new_coefficients.begin(), new_coefficients.end());
+}
+
+
+template <typename number>
+template <typename number2>
+void
+Polynomial<number>::shift(const number2 offset)
+{
+ shift(coefficients, offset);
+}
+
+
// ------------------------------------------------------------ //
template class Polynomial<float>;
template class Polynomial<double>;
template class Polynomial<long double>;
+
+template void Polynomial<float>::shift(const float offset);
+template void Polynomial<float>::shift(const double offset);
+template void Polynomial<float>::shift(const long double offset);
+template void Polynomial<double>::shift(const double offset);
+template void Polynomial<double>::shift(const long double offset);
+template void Polynomial<long double>::shift(const long double offset);
*/
unsigned int memory_consumption () const;
- /**
- * Exception
- */
- DeclException0 (ExcMatrixNotInitialized);
/**
* Exception
*/
int, int,
<< "The entry with index <" << arg1 << ',' << arg2
<< "> does not exist.");
+
/**
* Exception
*/
DeclException0 (ExcMatrixNotSquare);
+
/**
* Exception
*/
const Polynomial<double>& 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<gauss.n_quadrature_points;++i)
{
- double x = 2.*gauss.point(i)(0)-1.;
+ double x = gauss.point(i)(0);
double P1 = p1.value(x);
double P2 = p2.value(x);
- sum += 2.*gauss.weight(i) * P1 * P2;
+ sum += gauss.weight(i) * P1 * P2;
}
return sum;
}
deallog.depth_console(0);
std::vector<Polynomial<double> > p;
+ std::vector<Polynomial<double> > 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<double>(i));
for (unsigned int i=0;i<p.size();++i)
p.clear();
for (unsigned int i=0;i<6;++i)
- p.push_back(LagrangeEquidistant(6, i));
+ {
+ p.push_back(LagrangeEquidistant(6, i));
+ q.push_back(LagrangeEquidistant(6, i));
+ }
// We add 1.0001 bacuse of bugs in
// the ostream classes
for (unsigned int j=0;j<p.size();++j)
deallog << 'P' << i << "(x" << j
<< ") =" << p[i].value((double) j/p.size())+1.0001 << std::endl;
+
+ for (unsigned int i=0;i<p.size();++i)
+ {
+ q[i].scale(.5);
+ for (unsigned int j=0;j<p.size();++j)
+ {
+ double x = (double) j/p.size();
+ if (fabs(q[i].value(2.*x)-p[i].value(x)) > 1.e-15)
+ deallog << "Polynomial " << i
+ << ": Values q(" << 2.*x
+ << ") and p(" << x
+ << ") differ after scale: " << q[i].value(2.*x)
+ << " != " << p[i].value(x)
+ << std::endl;
+ }
+ q[i].shift((double) 1.);
+ for (unsigned int j=0;j<p.size();++j)
+ {
+ double x = (double) j/p.size();
+ double diff = fabs(q[i].value(2.*x-1.)-p[i].value(x));
+ if (diff > 1.e-13)
+ deallog << "Polynomial " << i
+ << ": Values q(" << 2.*x-1.
+ << ") and p(" << x
+ << ") differ by 10^" << log(diff)/log(10)
+ << " after shift: " << q[i].value(2.*x-1.)
+ << " != " << p[i].value(x)
+ << std::endl;
+ }
+ }
}
DEAL::Legendre
-DEAL::P0 * P0 =2.00
+DEAL::P0 * P0 =1.00
DEAL::P1 * P0 =0.00
-DEAL::P1 * P1 =0.667
+DEAL::P1 * P1 =1.00
DEAL::P2 * P0 =0.00
DEAL::P2 * P1 =0.00
-DEAL::P2 * P2 =0.400
+DEAL::P2 * P2 =1.00
DEAL::P3 * P0 =0.00
DEAL::P3 * P1 =0.00
DEAL::P3 * P2 =0.00
-DEAL::P3 * P3 =0.286
+DEAL::P3 * P3 =1.00
DEAL::P4 * P0 =0.00
DEAL::P4 * P1 =0.00
DEAL::P4 * P2 =0.00
DEAL::P4 * P3 =0.00
-DEAL::P4 * P4 =0.222
+DEAL::P4 * P4 =1.00
DEAL::P5 * P0 =0.00
DEAL::P5 * P1 =0.00
DEAL::P5 * P2 =0.00
DEAL::P5 * P3 =0.00
DEAL::P5 * P4 =0.00
-DEAL::P5 * P5 =0.182
+DEAL::P5 * P5 =1.00
DEAL::P6 * P0 =0.00
DEAL::P6 * P1 =0.00
DEAL::P6 * P2 =0.00
DEAL::P6 * P3 =0.00
DEAL::P6 * P4 =0.00
DEAL::P6 * P5 =0.00
-DEAL::P6 * P6 =0.154
+DEAL::P6 * P6 =1.00
DEAL::P7 * P0 =0.00
DEAL::P7 * P1 =0.00
DEAL::P7 * P2 =0.00
DEAL::P7 * P4 =0.00
DEAL::P7 * P5 =0.00
DEAL::P7 * P6 =0.00
-DEAL::P7 * P7 =0.133
+DEAL::P7 * P7 =1.00
DEAL::P8 * P0 =0.00
DEAL::P8 * P1 =0.00
DEAL::P8 * P2 =0.00
DEAL::P8 * P5 =0.00
DEAL::P8 * P6 =0.00
DEAL::P8 * P7 =0.00
-DEAL::P8 * P8 =0.118
+DEAL::P8 * P8 =1.00
DEAL::P9 * P0 =0.00
DEAL::P9 * P1 =0.00
DEAL::P9 * P2 =0.00
DEAL::P9 * P6 =0.00
DEAL::P9 * P7 =0.00
DEAL::P9 * P8 =0.00
-DEAL::P9 * P9 =0.105
+DEAL::P9 * P9 =1.00
DEAL::P10 * P0 =0.00
DEAL::P10 * P1 =0.00
DEAL::P10 * P2 =0.00
DEAL::P10 * P7 =0.00
DEAL::P10 * P8 =0.00
DEAL::P10 * P9 =0.00
-DEAL::P10 * P10 =0.0952
+DEAL::P10 * P10 =1.00
DEAL::P11 * P0 =0.00
DEAL::P11 * P1 =0.00
DEAL::P11 * P2 =0.00
DEAL::P11 * P8 =0.00
DEAL::P11 * P9 =0.00
DEAL::P11 * P10 =0.00
-DEAL::P11 * P11 =0.0870
-DEAL::P12 * P0 =0.00
-DEAL::P12 * P1 =0.00
-DEAL::P12 * P2 =0.00
-DEAL::P12 * P3 =0.00
-DEAL::P12 * P4 =0.00
-DEAL::P12 * P5 =0.00
-DEAL::P12 * P6 =0.00
-DEAL::P12 * P7 =0.00
-DEAL::P12 * P8 =0.00
-DEAL::P12 * P9 =0.00
-DEAL::P12 * P10 =0.00
-DEAL::P12 * P11 =0.00
-DEAL::P12 * P12 =0.0800
-DEAL::P13 * P0 =0.00
-DEAL::P13 * P1 =0.00
-DEAL::P13 * P2 =0.00
-DEAL::P13 * P3 =0.00
-DEAL::P13 * P4 =0.00
-DEAL::P13 * P5 =0.00
-DEAL::P13 * P6 =0.00
-DEAL::P13 * P7 =0.00
-DEAL::P13 * P8 =0.00
-DEAL::P13 * P9 =0.00
-DEAL::P13 * P10 =0.00
-DEAL::P13 * P11 =0.00
-DEAL::P13 * P12 =0.00
-DEAL::P13 * P13 =0.0741
-DEAL::P14 * P0 =0.00
-DEAL::P14 * P1 =0.00
-DEAL::P14 * P2 =0.00
-DEAL::P14 * P3 =0.00
-DEAL::P14 * P4 =0.00
-DEAL::P14 * P5 =0.00
-DEAL::P14 * P6 =0.00
-DEAL::P14 * P7 =0.00
-DEAL::P14 * P8 =0.00
-DEAL::P14 * P9 =0.00
-DEAL::P14 * P10 =0.00
-DEAL::P14 * P11 =0.00
-DEAL::P14 * P12 =0.00
-DEAL::P14 * P13 =0.00
-DEAL::P14 * P14 =0.0690
+DEAL::P11 * P11 =1.00
DEAL::LagrangeEquidistant
DEAL::P0(x0) =2.00
DEAL::P0(x1) =1.00
DEAL:Lagrange:3d:Polyno::P9 = 1.3 gradient -2.6 -14. -37. 2nd -0.095 0.029 0.074 0.029 0.074 0.40 0.074 0.40 -0.16
DEAL:Lagrange:3d:Polyno::
DEAL:Legendre:1d:Tensor::P0 = 10. gradient 0.0 2nd 0.0
-DEAL:Legendre:1d:Tensor::P1 = 5.0 gradient 10. 2nd 0.0
-DEAL:Legendre:1d:Tensor::P2 = -1.2 gradient 15. 2nd 3.0
+DEAL:Legendre:1d:Tensor::P1 = 0.0 gradient 35. 2nd 0.0
+DEAL:Legendre:1d:Tensor::P2 = -11. gradient 0.0 2nd 27.
DEAL:Legendre:1d:Tensor::
DEAL:Legendre:1d:Polyno::P0 = 10. gradient 0.0 2nd 0.0
-DEAL:Legendre:1d:Polyno::P1 = 5.0 gradient 10. 2nd 0.0
-DEAL:Legendre:1d:Polyno::P2 = -1.2 gradient 15. 2nd 3.0
+DEAL:Legendre:1d:Polyno::P1 = 0.0 gradient 35. 2nd 0.0
+DEAL:Legendre:1d:Polyno::P2 = -11. gradient 0.0 2nd 27.
DEAL:Legendre:1d:Polyno::
DEAL:Legendre:2d:Tensor::P0 = 1.0e+02 gradient 0.0 0.0 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Tensor::P1 = 50. gradient 1.0e+02 0.0 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Tensor::P2 = -12. gradient 1.5e+02 -0.0 2nd 3.0 0.0 0.0 -0.0
-DEAL:Legendre:2d:Tensor::P3 = 20. gradient 0.0 1.0e+02 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Tensor::P4 = 10. gradient 20. 50. 2nd 0.0 1.0 1.0 0.0
-DEAL:Legendre:2d:Tensor::P5 = -2.5 gradient 30. -12. 2nd 0.60 1.5 1.5 -0.0
-DEAL:Legendre:2d:Tensor::P6 = -44. gradient -0.0 60. 2nd -0.0 0.0 0.0 3.0
-DEAL:Legendre:2d:Tensor::P7 = -22. gradient -44. 30. 2nd -0.0 0.60 0.60 1.5
-DEAL:Legendre:2d:Tensor::P8 = 5.5 gradient -66. -7.5 2nd -1.3 0.90 0.90 -0.38
+DEAL:Legendre:2d:Tensor::P1 = 0.0 gradient 3.5e+02 0.0 2nd 0.0 0.0 0.0 0.0
+DEAL:Legendre:2d:Tensor::P2 = -1.1e+02 gradient 0.0 -0.0 2nd 27. 0.0 0.0 -0.0
+DEAL:Legendre:2d:Tensor::P3 = -1.0e+02 gradient -0.0 3.5e+02 2nd -0.0 0.0 0.0 0.0
+DEAL:Legendre:2d:Tensor::P4 = -0.0 gradient -3.6e+02 0.0 2nd -0.0 12. 12. 0.0
+DEAL:Legendre:2d:Tensor::P5 = 1.2e+02 gradient -0.0 -3.9e+02 2nd -28. 0.0 0.0 -0.0
+DEAL:Legendre:2d:Tensor::P6 = 8.9 gradient 0.0 -8.0e+02 2nd 0.0 -0.0 -0.0 27.
+DEAL:Legendre:2d:Tensor::P7 = 0.0 gradient 31. -0.0 2nd 0.0 -28. -28. 0.0
+DEAL:Legendre:2d:Tensor::P8 = -10. gradient 0.0 9.0e+02 2nd 2.4 -0.0 -0.0 -30.
DEAL:Legendre:2d:Tensor::
DEAL:Legendre:2d:Polyno::P0 = 1.0e+02 gradient 0.0 0.0 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Polyno::P1 = 50. gradient 1.0e+02 0.0 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Polyno::P2 = -12. gradient 1.5e+02 -0.0 2nd 3.0 0.0 0.0 -0.0
-DEAL:Legendre:2d:Polyno::P3 = 20. gradient 0.0 1.0e+02 2nd 0.0 0.0 0.0 0.0
-DEAL:Legendre:2d:Polyno::P4 = 10. gradient 20. 50. 2nd 0.0 1.0 1.0 0.0
-DEAL:Legendre:2d:Polyno::P5 = -44. gradient -0.0 60. 2nd -0.0 0.0 0.0 3.0
+DEAL:Legendre:2d:Polyno::P1 = 0.0 gradient 3.5e+02 0.0 2nd 0.0 0.0 0.0 0.0
+DEAL:Legendre:2d:Polyno::P2 = -1.1e+02 gradient 0.0 -0.0 2nd 27. 0.0 0.0 -0.0
+DEAL:Legendre:2d:Polyno::P3 = -1.0e+02 gradient -0.0 3.5e+02 2nd -0.0 0.0 0.0 0.0
+DEAL:Legendre:2d:Polyno::P4 = -0.0 gradient -3.6e+02 0.0 2nd -0.0 12. 12. 0.0
+DEAL:Legendre:2d:Polyno::P5 = 8.9 gradient 0.0 -8.0e+02 2nd 0.0 -0.0 -0.0 27.
DEAL:Legendre:2d:Polyno::
DEAL:Legendre:3d:Tensor::P0 = 1.0e+03 gradient 0.0 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P1 = 5.0e+02 gradient 1.0e+03 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P2 = -1.2e+02 gradient 1.5e+03 -0.0 -0.0 2nd 3.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0
-DEAL:Legendre:3d:Tensor::P3 = 2.0e+02 gradient 0.0 1.0e+03 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P4 = 1.0e+02 gradient 2.0e+02 5.0e+02 0.0 2nd 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P5 = -25. gradient 3.0e+02 -1.2e+02 -0.0 2nd 0.60 1.5 0.0 1.5 -0.0 -0.0 0.0 -0.0 -0.0
-DEAL:Legendre:3d:Tensor::P6 = -4.4e+02 gradient -0.0 6.0e+02 -0.0 2nd -0.0 0.0 -0.0 0.0 3.0 0.0 -0.0 0.0 -0.0
-DEAL:Legendre:3d:Tensor::P7 = -2.2e+02 gradient -4.4e+02 3.0e+02 -0.0 2nd -0.0 0.60 -0.0 0.60 1.5 0.0 -0.0 0.0 -0.0
-DEAL:Legendre:3d:Tensor::P8 = 55. gradient -6.6e+02 -75. 0.0 2nd -1.3 0.90 -0.0 0.90 -0.38 -0.0 -0.0 -0.0 0.0
-DEAL:Legendre:3d:Tensor::P9 = 3.0e+02 gradient 0.0 0.0 1.0e+03 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P10 = 1.5e+02 gradient 3.0e+02 0.0 5.0e+02 2nd 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0
-DEAL:Legendre:3d:Tensor::P11 = -38. gradient 4.5e+02 -0.0 -1.2e+02 2nd 0.90 0.0 1.5 0.0 -0.0 -0.0 1.5 -0.0 -0.0
-DEAL:Legendre:3d:Tensor::P12 = 60. gradient 0.0 3.0e+02 2.0e+02 2nd 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0
-DEAL:Legendre:3d:Tensor::P13 = 30. gradient 60. 1.5e+02 1.0e+02 2nd 0.0 0.30 0.20 0.30 0.0 0.50 0.20 0.50 0.0
-DEAL:Legendre:3d:Tensor::P14 = -7.5 gradient 90. -38. -25. 2nd 0.18 0.45 0.30 0.45 -0.0 -0.12 0.30 -0.12 -0.0
-DEAL:Legendre:3d:Tensor::P15 = -1.3e+02 gradient -0.0 1.8e+02 -4.4e+02 2nd -0.0 0.0 -0.0 0.0 0.90 0.60 -0.0 0.60 -0.0
-DEAL:Legendre:3d:Tensor::P16 = -66. gradient -1.3e+02 90. -2.2e+02 2nd -0.0 0.18 -0.44 0.18 0.45 0.30 -0.44 0.30 -0.0
-DEAL:Legendre:3d:Tensor::P17 = 16. gradient -2.0e+02 -23. 55. 2nd -0.40 0.27 -0.66 0.27 -0.11 -0.075 -0.66 -0.075 0.0
-DEAL:Legendre:3d:Tensor::P18 = -3.6e+02 gradient -0.0 -0.0 9.0e+02 2nd -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 3.0
-DEAL:Legendre:3d:Tensor::P19 = -1.8e+02 gradient -3.6e+02 -0.0 4.5e+02 2nd -0.0 -0.0 0.90 -0.0 -0.0 0.0 0.90 0.0 1.5
-DEAL:Legendre:3d:Tensor::P20 = 46. gradient -5.5e+02 0.0 -1.1e+02 2nd -1.1 -0.0 1.3 -0.0 0.0 -0.0 1.3 -0.0 -0.38
-DEAL:Legendre:3d:Tensor::P21 = -73. gradient -0.0 -3.6e+02 1.8e+02 2nd -0.0 -0.0 0.0 -0.0 -0.0 0.90 0.0 0.90 0.60
-DEAL:Legendre:3d:Tensor::P22 = -36. gradient -73. -1.8e+02 90. 2nd -0.0 -0.36 0.18 -0.36 -0.0 0.45 0.18 0.45 0.30
-DEAL:Legendre:3d:Tensor::P23 = 9.1 gradient -1.1e+02 46. -22. 2nd -0.22 -0.55 0.27 -0.55 0.0 -0.11 0.27 -0.11 -0.075
-DEAL:Legendre:3d:Tensor::P24 = 1.6e+02 gradient 0.0 -2.2e+02 -4.0e+02 2nd 0.0 -0.0 -0.0 -0.0 -1.1 0.54 -0.0 0.54 -1.3
-DEAL:Legendre:3d:Tensor::P25 = 80. gradient 1.6e+02 -1.1e+02 -2.0e+02 2nd 0.0 -0.22 -0.40 -0.22 -0.55 0.27 -0.40 0.27 -0.66
-DEAL:Legendre:3d:Tensor::P26 = -20. gradient 2.4e+02 27. 49. 2nd 0.48 -0.33 -0.59 -0.33 0.14 -0.068 -0.59 -0.068 0.17
+DEAL:Legendre:3d:Tensor::P1 = 0.0 gradient 3.5e+03 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+DEAL:Legendre:3d:Tensor::P2 = -1.1e+03 gradient 0.0 -0.0 -0.0 2nd 27. 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0
+DEAL:Legendre:3d:Tensor::P3 = -1.0e+03 gradient -0.0 3.5e+03 -0.0 2nd -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 -0.0
+DEAL:Legendre:3d:Tensor::P4 = -0.0 gradient -3.6e+03 0.0 -0.0 2nd -0.0 12. -0.0 12. 0.0 0.0 -0.0 0.0 -0.0
+DEAL:Legendre:3d:Tensor::P5 = 1.2e+03 gradient -0.0 -3.9e+03 0.0 2nd -28. 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0
+DEAL:Legendre:3d:Tensor::P6 = 89. gradient 0.0 -8.0e+03 0.0 2nd 0.0 -0.0 0.0 -0.0 27. -0.0 0.0 -0.0 0.0
+DEAL:Legendre:3d:Tensor::P7 = 0.0 gradient 3.1e+02 -0.0 0.0 2nd 0.0 -28. 0.0 -28. 0.0 -0.0 0.0 -0.0 0.0
+DEAL:Legendre:3d:Tensor::P8 = -1.0e+02 gradient 0.0 9.0e+03 -0.0 2nd 2.4 -0.0 0.0 -0.0 -30. 0.0 0.0 0.0 -0.0
+DEAL:Legendre:3d:Tensor::P9 = -6.9e+02 gradient -0.0 -0.0 3.5e+03 2nd -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0
+DEAL:Legendre:3d:Tensor::P10 = -0.0 gradient -2.4e+03 -0.0 0.0 2nd -0.0 -0.0 12. -0.0 -0.0 0.0 12. 0.0 0.0
+DEAL:Legendre:3d:Tensor::P11 = 7.7e+02 gradient -0.0 0.0 -3.9e+03 2nd -19. -0.0 0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0
+DEAL:Legendre:3d:Tensor::P12 = 7.2e+02 gradient 0.0 -2.4e+03 -3.6e+03 2nd 0.0 -0.0 -0.0 -0.0 -0.0 12. -0.0 12. -0.0
+DEAL:Legendre:3d:Tensor::P13 = 0.0 gradient 2.5e+03 -0.0 -0.0 2nd 0.0 -8.3 -12. -8.3 -0.0 0.0 -12. 0.0 -0.0
+DEAL:Legendre:3d:Tensor::P14 = -8.0e+02 gradient 0.0 2.7e+03 4.0e+03 2nd 19. -0.0 -0.0 -0.0 0.0 -13. -0.0 -13. 0.0
+DEAL:Legendre:3d:Tensor::P15 = -62. gradient -0.0 5.6e+03 3.1e+02 2nd -0.0 0.0 0.0 0.0 -19. -28. 0.0 -28. 0.0
+DEAL:Legendre:3d:Tensor::P16 = -0.0 gradient -2.1e+02 0.0 0.0 2nd -0.0 19. 1.1 19. -0.0 -0.0 1.1 -0.0 0.0
+DEAL:Legendre:3d:Tensor::P17 = 69. gradient -0.0 -6.2e+03 -3.5e+02 2nd -1.7 0.0 0.0 0.0 21. 31. 0.0 31. -0.0
+DEAL:Legendre:3d:Tensor::P18 = -5.8e+02 gradient -0.0 -0.0 -5.4e+03 2nd -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 27.
+DEAL:Legendre:3d:Tensor::P19 = -0.0 gradient -2.0e+03 -0.0 -0.0 2nd -0.0 -0.0 -19. -0.0 -0.0 -0.0 -19. -0.0 0.0
+DEAL:Legendre:3d:Tensor::P20 = 6.5e+02 gradient -0.0 0.0 6.0e+03 2nd -16. -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 -30.
+DEAL:Legendre:3d:Tensor::P21 = 6.0e+02 gradient 0.0 -2.0e+03 5.6e+03 2nd 0.0 -0.0 0.0 -0.0 -0.0 -19. 0.0 -19. -28.
+DEAL:Legendre:3d:Tensor::P22 = 0.0 gradient 2.1e+03 -0.0 0.0 2nd 0.0 -7.0 19. -7.0 -0.0 -0.0 19. -0.0 -0.0
+DEAL:Legendre:3d:Tensor::P23 = -6.8e+02 gradient 0.0 2.3e+03 -6.2e+03 2nd 16. -0.0 0.0 -0.0 0.0 21. 0.0 21. 31.
+DEAL:Legendre:3d:Tensor::P24 = -52. gradient -0.0 4.7e+03 -4.8e+02 2nd -0.0 0.0 -0.0 0.0 -16. 43. -0.0 43. 2.4
+DEAL:Legendre:3d:Tensor::P25 = -0.0 gradient -1.8e+02 0.0 -0.0 2nd -0.0 16. -1.7 16. -0.0 0.0 -1.7 0.0 0.0
+DEAL:Legendre:3d:Tensor::P26 = 58. gradient -0.0 -5.2e+03 5.4e+02 2nd -1.4 0.0 -0.0 0.0 17. -48. -0.0 -48. -2.7
DEAL:Legendre:3d:Tensor::
DEAL:Legendre:3d:Polyno::P0 = 1.0e+03 gradient 0.0 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P1 = 5.0e+02 gradient 1.0e+03 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P2 = -1.2e+02 gradient 1.5e+03 -0.0 -0.0 2nd 3.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0
-DEAL:Legendre:3d:Polyno::P3 = 2.0e+02 gradient 0.0 1.0e+03 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P4 = 1.0e+02 gradient 2.0e+02 5.0e+02 0.0 2nd 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P5 = -4.4e+02 gradient -0.0 6.0e+02 -0.0 2nd -0.0 0.0 -0.0 0.0 3.0 0.0 -0.0 0.0 -0.0
-DEAL:Legendre:3d:Polyno::P6 = 3.0e+02 gradient 0.0 0.0 1.0e+03 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P7 = 1.5e+02 gradient 3.0e+02 0.0 5.0e+02 2nd 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0
-DEAL:Legendre:3d:Polyno::P8 = 60. gradient 0.0 3.0e+02 2.0e+02 2nd 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0
-DEAL:Legendre:3d:Polyno::P9 = -3.6e+02 gradient -0.0 -0.0 9.0e+02 2nd -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 3.0
+DEAL:Legendre:3d:Polyno::P1 = 0.0 gradient 3.5e+03 0.0 0.0 2nd 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+DEAL:Legendre:3d:Polyno::P2 = -1.1e+03 gradient 0.0 -0.0 -0.0 2nd 27. 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0
+DEAL:Legendre:3d:Polyno::P3 = -1.0e+03 gradient -0.0 3.5e+03 -0.0 2nd -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 -0.0
+DEAL:Legendre:3d:Polyno::P4 = -0.0 gradient -3.6e+03 0.0 -0.0 2nd -0.0 12. -0.0 12. 0.0 0.0 -0.0 0.0 -0.0
+DEAL:Legendre:3d:Polyno::P5 = 89. gradient 0.0 -8.0e+03 0.0 2nd 0.0 -0.0 0.0 -0.0 27. -0.0 0.0 -0.0 0.0
+DEAL:Legendre:3d:Polyno::P6 = -6.9e+02 gradient -0.0 -0.0 3.5e+03 2nd -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0
+DEAL:Legendre:3d:Polyno::P7 = -0.0 gradient -2.4e+03 -0.0 0.0 2nd -0.0 -0.0 12. -0.0 -0.0 0.0 12. 0.0 0.0
+DEAL:Legendre:3d:Polyno::P8 = 7.2e+02 gradient 0.0 -2.4e+03 -3.6e+03 2nd 0.0 -0.0 -0.0 -0.0 -0.0 12. -0.0 12. -0.0
+DEAL:Legendre:3d:Polyno::P9 = -5.8e+02 gradient -0.0 -0.0 -5.4e+03 2nd -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 27.
DEAL:Legendre:3d:Polyno::
#include <dofs/dof_accessor.h>
#include <grid/grid_generator.h>
#include <fe/fe_q.h>
+#include <fe/fe_dgp.h>
#include <fe/fe_dgq.h>
#include <fe/fe_system.h>
#include <fe/mapping_q1.h>
}
+template<int dim>
+void plot_FE_DGP_shape_functions()
+{
+ MappingQ1<dim> m;
+ FE_DGP<dim> p1(1);
+ plot_shape_functions(m, p1, "DGP1");
+ plot_face_shape_functions(m, p1, "DGP1");
+ test_compute_functions(m, p1, "DGP1");
+ FE_DGP<dim> p2(2);
+ plot_shape_functions(m, p2, "DGP2");
+ plot_face_shape_functions(m, p2, "DGP2");
+ test_compute_functions(m, p2, "DGP2");
+ FE_DGP<dim> p3(3);
+ plot_shape_functions(m, p3, "DGP3");
+ plot_face_shape_functions(m, p3, "DGP3");
+ test_compute_functions(m, p3, "DGP3");
+ FE_DGP<dim> p4(4);
+ plot_shape_functions(m, p4, "DGP4");
+ plot_face_shape_functions(m, p4, "DGP4");
+ test_compute_functions(m, p4, "DGP4");
+// FE_DGP<dim> p5(5);
+// plot_shape_functions(m, p5, "DGP5");
+// FE_DGP<dim> p6(6);
+// plot_shape_functions(m, p6, "DGP6");
+// FE_DGP<dim> p7(7);
+// plot_shape_functions(m, p7, "DGP7");
+// FE_DGP<dim> p8(8);
+// plot_shape_functions(m, p8, "DGP8");
+// FE_DGP<dim> p9(9);
+// plot_shape_functions(m, p9, "DGP9");
+// FE_DGP<dim> p10(10);
+// plot_shape_functions(m, p10, "DGP10");
+}
+
+
int
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.