]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix computation of points in QGaussLobatto.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 12 Aug 2010 20:25:38 +0000 (20:25 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 12 Aug 2010 20:25:38 +0000 (20:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@21648 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/quadrature_lib.cc
deal.II/doc/news/changes.h

index 4370d4bf51dfc0c728b5e494cdd44268a7d9e6a4..0ee73d7fbf9ae74b3f622d0c3dfe0f35f5e40724 100644 (file)
@@ -230,6 +230,14 @@ compute_quadrature_points(const unsigned int q,
        double_eps * 5
        );
 
+                                  // The following implementation
+                                  // follows closely the one given in
+                                  // the appendix of the book by
+                                  // Karniadakis and Sherwin:
+                                  // Spectral/hp element methods for
+                                  // computational fluid dynamics
+                                  // (Oxford University Press, 2005)
+
                                   // we take the zeros of the Chebyshev
                                   // polynomial (alpha=beta=-0.5) as
                                   // initial values:
@@ -237,6 +245,7 @@ compute_quadrature_points(const unsigned int q,
     x[i] = - std::cos( (long double) (2*i+1)/(2*m) * numbers::PI );
 
   long double r, s, J_x, f, delta;
+
   for (unsigned int k=0; k<m; ++k)
     {
       r = x[k];
@@ -245,9 +254,9 @@ compute_quadrature_points(const unsigned int q,
 
       do
        {
-         s = 1.;
+         s = 0.;
          for (unsigned int i=0; i<k; ++i)
-           s /= r - x[i];
+           s += 1./(r - x[i]);
 
          J_x   =  0.5*(alpha + beta + m + 1)*JacobiP(r, alpha+1, beta+1, m-1);
          f     = JacobiP(r, alpha, beta, m);
index 257f1d3e4d4b7b16aeb62bb03b86eab5e9d8225f..e068cbfd194413b0428f1a538f4d1b2f9362d0fc 100644 (file)
@@ -164,13 +164,21 @@ inconvenience this causes.
 
 
 <ol>
+  <li><p> Fixed: The computation of quadrature points in the QGaussLobatto
+  class uses a Newton method that was wrongly implemented. While the
+  results were correct (at least for moderate orders), it required more
+  iterations than necessary. This is now fixed.
+  <br>
+  (Andrea Bonito 2010/08/12)
+  </p>
+
   <li><p> Changed: The DataOutBase::write_vtu function now writes data
   as 64-bit values, rather than 32-bit values.
   <br>
   (Scott Miller 2010/08/5)
   </p>
 
-  <li><p> New: MappingQ and MappingQEulerian now support order > 1 also in 
+  <li><p> New: MappingQ and MappingQEulerian now support order > 1 also in
   codimension one. Step-34 has been modified to show how this works.
   <br>
   (Luca Heltai 2010/07/23-27)
@@ -187,12 +195,12 @@ inconvenience this causes.
 <a name="lac"></a>
 <h3>lac</h3>
 
-<ol> 
+<ol>
   <li><p>Updated: Changes to the <code>SLEPcWrappers</code> to preserve
   compatibility with SLEPc version 3.1. Main new updated features are
   new solver classes: (i) Power; and (ii) Generalized Davidson.
-  <br> 
-  (Toby D. Young 2010/08/04) 
+  <br>
+  (Toby D. Young 2010/08/04)
   </p></li>
 
   <li><p>Fixed: SparseMatrix::precondition_SSOR and
@@ -216,7 +224,7 @@ inconvenience this causes.
   <br>
   (Andrew McBride 2010/07/29)
   </p>
-  
+
   <li>
   <p>
   Fixed: In some rather pathological cases, the function

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.