]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Legendre polynomials are now orthoNORMAL on [0,1]
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Jul 2002 18:06:28 +0000 (18:06 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Jul 2002 18:06:28 +0000 (18:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@6225 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h
deal.II/base/source/legendre.cc
deal.II/base/source/polynomial.cc
deal.II/lac/include/lac/sparse_matrix_ez.h
tests/base/polynomial1d.cc
tests/base/polynomial1d.checked
tests/base/polynomial_test.checked
tests/fe/shapes.cc

index f15e4a2c26ca1e3560dee41d757aafc450585728..7b3c298cb5088e5a5d77558b70f08c0753b5953c 100644 (file)
@@ -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 <typename number2>
+    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<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
@@ -192,12 +255,14 @@ class LagrangeEquidistant: public Polynomial<double>
 
 
 /**
- * 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<number>
     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
@@ -238,7 +308,7 @@ class Legendre : public Polynomial<number>
                                      * 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.
index 2344a45f0fd4a0d0c92e0b37df681ec68487fc69..120e85376d9d63412ff7040bb7f33232fa7b6541 100644 (file)
 // 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
@@ -60,13 +64,13 @@ Legendre<number>::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<number>::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<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
        {
@@ -110,19 +126,26 @@ Legendre<number>::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<number>(*ck);
+         shift(*ck,(long double) -1.);
+         Polynomial<number>::scale(*ck, 2.);
+         Polynomial<number>::multiply(*ck, sqrt(2.*k+1.));
+         shifted_coefficients[k] = ck;
        };
     };
 
@@ -145,7 +168,7 @@ Legendre<number>::get_coefficients (const unsigned int k)
                                   // 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
index 71b42690303fb48cfdfc73b11129ccb4eb5e33c2..cb93b56208d09e7adc89ddc584d6cb0947d9c02b 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <base/polynomial.h>
+#include <base/exceptions.h>
 
 template <typename number>
 Polynomial<number>::Polynomial (const typename std::vector<number> &a):
@@ -91,6 +92,108 @@ Polynomial<number>::value (const number         x,
 }
 
 
+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);
+}
+
+
 // ------------------------------------------------------------ //
 
 
@@ -348,3 +451,10 @@ generate_complete_basis (const unsigned int degree)
 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);
index 2bbac3b8b6f33dad4bc7f0b78342cd4e49d36234..f0b8a36fec74dab70e03660ea16e9d65efbc327e 100644 (file)
@@ -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
                                      */
index 89d864b0b1fa32539fb2b644ec4e41d0faf276f4..efc0ae8c31b435dc7e6e8d2f1c48a407571e8ba6 100644 (file)
@@ -24,15 +24,15 @@ double scalar_product (const Polynomial<double>& p1,
                       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;
 }
@@ -45,10 +45,11 @@ int main ()
   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)
@@ -61,7 +62,10 @@ int main ()
   
   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
@@ -69,5 +73,35 @@ int main ()
     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;
+       }
+    }
 }
 
index 31c5c15030532a749bf1ec48bc5b0647458ae46b..2f5d26cd6efd1a1a83b829a1eec7c00487b6a27f 100644 (file)
@@ -1,33 +1,33 @@
 
 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
@@ -35,7 +35,7 @@ DEAL::P7 * P3 =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
@@ -44,7 +44,7 @@ DEAL::P8 * P4 =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
@@ -54,7 +54,7 @@ DEAL::P9 * P5 =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
@@ -65,7 +65,7 @@ DEAL::P10 * P6 =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
@@ -77,49 +77,7 @@ DEAL::P11 * P7 =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
index a98dcef4337a70ecf8a6681dc2c07b52ff1e95f2..0bf5f06631fa872a0a0c400432bcc9ba49579d34 100644 (file)
@@ -64,66 +64,66 @@ DEAL:Lagrange:3d:Polyno::P8 = -65.  gradient        1.3e+02 -1.1e+02        54.             2nd     4.7     0.21    -
 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::
index 04091e7b3bd71d12f01fdb210efeab4082b0ff81..44560287d4a6d5322bce9cbc5bb98942d93588a6 100644 (file)
@@ -11,6 +11,7 @@
 #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>
@@ -387,6 +388,41 @@ void plot_FE_DGQ_shape_functions()
 }
 
 
+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()
 {
@@ -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.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.