--- Begin Message ---
------------------------------------------------------------
revno: 5923
committer: Garth N. Wells <gnw20@xxxxxxxxx>
branch nick: dolfin-all
timestamp: Wed 2011-06-01 23:19:07 +0100
message:
Replace Legendre code with Boost Legendre code.
modified:
dolfin/math/Legendre.cpp
dolfin/math/Legendre.h
dolfin/math/basic.cpp
dolfin/math/basic.h
dolfin/quadrature/GaussQuadrature.cpp
dolfin/quadrature/GaussianQuadrature.cpp
dolfin/quadrature/LobattoQuadrature.cpp
dolfin/quadrature/RadauQuadrature.cpp
dolfin/swig/docstrings.i
site-packages/dolfin/cppimports.py
--
lp:~dolfin-core/dolfin/wells
https://code.launchpad.net/~dolfin-core/dolfin/wells
You are subscribed to branch lp:~dolfin-core/dolfin/wells.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/wells/+edit-subscription
=== modified file 'dolfin/math/Legendre.cpp'
--- dolfin/math/Legendre.cpp 2011-05-24 09:50:28 +0000
+++ dolfin/math/Legendre.cpp 2011-06-01 22:19:07 +0000
@@ -21,101 +21,30 @@
// Last changed: 2009-02-17
#include <cmath>
+#include <boost/math/special_functions/legendre.hpp>
+
#include <dolfin/common/constants.h>
+#include <dolfin/common/real.h>
#include <dolfin/log/dolfin_log.h>
-#include <dolfin/common/real.h>
#include "Legendre.h"
using namespace dolfin;
//-----------------------------------------------------------------------------
-Legendre::Legendre(uint n) : n(n), cache_x(0.0), cache(n + 1)
-{
- cache[0] = 1.0; //constant value
-
- // eval to initialize cache
- eval(n, -1.0);
-}
-//-----------------------------------------------------------------------------
-real Legendre::operator() (real x)
-{
- return eval(n, x);
-}
-//-----------------------------------------------------------------------------
-real Legendre::ddx(real x)
-{
- return ddx(n, x);
-}
-//-----------------------------------------------------------------------------
-real Legendre::d2dx(real x)
-{
- return d2dx(n, x);
-}
-//-----------------------------------------------------------------------------
-real Legendre::eval(uint nn, real x)
-{
- //recursive formula, BETA page 254
- //return ( (2.0*nn-1.0)*x*eval(nn-1, x) - (nn-1.0)*eval(nn-2, x) ) / nn;
-
- //The special cases
- if (n == 0)
- return 1.0;
- else if (n == 1)
- return x;
-
- //check cache
- if (x != cache_x)
- {
- cache[1] = x;
- for (uint i = 2; i <= n; ++i)
- {
- real ii(i);
- cache[i] = ( (2.0*ii-1.0)*x*cache[i-1] - (ii-1.0)*cache[i-2] ) / ii;
- }
- cache_x = x;
- }
-
- return cache[nn];
+real Legendre::eval(uint n, real x)
+{
+ return boost::math::legendre_p(n, x);
}
//-----------------------------------------------------------------------------
real Legendre::ddx(uint n, real x)
{
- // Special cases
- if (n == 0)
- return 0.0;
- else if (n == 1)
- return 1.0;
-
- // Avoid division by zero
- if (real_abs(x - 1.0) < real_epsilon())
- x -= 2.0*real_epsilon();
-
- if (real_abs(x + 1.0) < real_epsilon())
- x += 2.0*real_epsilon();
-
- // Formula, BETA page 254
- const real nn = real(n);
- return nn * (x*eval(n, x) - eval(n-1, x)) / (x*x - 1.0);
+ assert((x*x - 1.0) != 0.0);
+ return -boost::math::legendre_p(n, 1, x)/(std::sqrt(1.0 - x*x));
}
//-----------------------------------------------------------------------------
-real Legendre::d2dx(uint, real x)
+real Legendre::d2dx(uint n, real x)
{
- // Special case n = 0
- if (n == 0)
- return 0.0;
-
- // Special case n = 1
- if (n == 1)
- return 0.0;
-
- // Avoid division by zero
- if (real_abs(x - 1.0) < real_epsilon())
- x -= 2.0*real_epsilon();
- if (real_abs(x + 1.0) < real_epsilon())
- x += 2.0*real_epsilon();
-
- // Formula, BETA page 254
- const real nn = real(n);
- return (2.0*x*ddx(n, x) - nn*(nn+1)*eval(n, x)) / (1.0-x*x);
+ assert((x*x - 1.0) != 0.0);
+ return boost::math::legendre_p(n, 2, x)/(1.0 - x*x);
}
//-----------------------------------------------------------------------------
=== modified file 'dolfin/math/Legendre.h'
--- dolfin/math/Legendre.h 2011-05-24 09:50:28 +0000
+++ dolfin/math/Legendre.h 2011-06-01 22:19:07 +0000
@@ -1,4 +1,4 @@
-// Copyright (C) 2003-2009 Anders Logg
+// Copyright (C) 2011 Garth N. Wells
//
// This file is part of DOLFIN.
//
@@ -15,56 +15,32 @@
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
-// First added: 2003-06-03
-// Last changed: 2009-02-17
+// First added: 2011-06-01
+// Last changed:
#ifndef __LEGENDRE_H
#define __LEGENDRE_H
-#include <vector>
+#include <dolfin/common/real.h>
#include <dolfin/common/types.h>
-#include <dolfin/common/real.h>
namespace dolfin
{
- /// Legendre polynomial of given degree n on the interval [-1,1].
- ///
- /// P0(x) = 1
- /// P1(x) = x
- /// P2(x) = (3x^2 - 1) / 2
- /// ...
- ///
- /// The function values and derivatives are computed using
- /// three-term recurrence formulas.
+ /// Interface for computing Legendre polynomials via Boost.
class Legendre
{
public:
- Legendre(uint n);
-
- /// Evaluation at given point
- real operator() (real x);
-
- /// Evaluation of derivative at given point
- real ddx(real x);
-
- /// Evaluation of second derivative at given point
- real d2dx(real x);
-
- /// Evaluation of arbitrary order, nn <= n (useful ie in RadauQuadrature)
- real eval(uint nn, real x);
-
- real ddx(uint n, real x);
- real d2dx(uint n, real x);
-
-
- private:
-
- const uint n;
- real cache_x;
- std::vector<real> cache;
+ /// Evaluate polynomial of order n at point x
+ static real eval(uint n, real x);
+
+ /// Evaluate first derivative of polynomial of order n at point x
+ static real ddx(uint n, real x);
+
+ /// Evaluate second derivative of polynomial of order n at point x
+ static real d2dx(uint n, real x);
};
=== modified file 'dolfin/math/basic.cpp'
--- dolfin/math/basic.cpp 2011-05-11 13:13:14 +0000
+++ dolfin/math/basic.cpp 2011-06-01 22:19:07 +0000
@@ -35,11 +35,6 @@
}
//-----------------------------------------------------------------------------
-double dolfin::sqr(double x)
-{
- return x*x;
-}
-//-----------------------------------------------------------------------------
dolfin::uint dolfin::ipow(uint a, uint n)
{
uint p = a;
@@ -50,7 +45,7 @@
//-----------------------------------------------------------------------------
double dolfin::rand()
{
- if ( !rand_seeded )
+ if (!rand_seeded)
{
unsigned int s = static_cast<long int>(time(0));
std::srand(s);
@@ -66,12 +61,3 @@
rand_seeded = true;
}
//-----------------------------------------------------------------------------
-bool dolfin::near(double x, double x0)
-{
- return dolfin::between(x0, x, x0);
-}
-//-----------------------------------------------------------------------------
-bool dolfin::between(double x0, double x, double x1)
-{
- return (x0-DOLFIN_EPS < x && x < x1+DOLFIN_EPS);
-}
=== modified file 'dolfin/math/basic.h'
--- dolfin/math/basic.h 2011-05-11 13:13:14 +0000
+++ dolfin/math/basic.h 2011-06-01 22:19:07 +0000
@@ -26,10 +26,6 @@
namespace dolfin
{
-
- /// Return the square of x
- double sqr(double x);
-
/// Return a to the power n
uint ipow(uint a, uint n);
@@ -39,11 +35,6 @@
/// Seed random number generator
void seed(unsigned int s);
- /// Return true if x is within DOLFIN_EPS of x0
- bool near(double x, double x0);
-
- /// Return true if x is between x0 and x1 (inclusive)
- bool between(double x0, double x, double x1);
}
#endif
=== modified file 'dolfin/quadrature/GaussQuadrature.cpp'
--- dolfin/quadrature/GaussQuadrature.cpp 2011-05-11 13:13:14 +0000
+++ dolfin/quadrature/GaussQuadrature.cpp 2011-06-01 22:19:07 +0000
@@ -33,14 +33,13 @@
{
init();
- if (!check(2*n-1))
+ if (!check(2*n - 1))
error("Gauss quadrature not ok, check failed.");
}
//-----------------------------------------------------------------------------
std::string GaussQuadrature::str(bool verbose) const
{
std::stringstream s;
-
if (verbose)
{
s << str(false) << std::endl << std::endl;
@@ -59,9 +58,7 @@
}
}
else
- {
s << "<GaussQuadrature with " << points.size() << " points on [-1, 1]>";
- }
return s.str();
}
@@ -80,26 +77,24 @@
return;
}
- Legendre p(n);
real x, dx;
// Compute the points by Newton's method
- for (unsigned int i = 0; i <= ((n-1)/2); i++)
+ for (unsigned int i = 0; i <= ((n - 1)/2); i++)
{
-
// Initial guess
- x = cos(DOLFIN_PI*(double(i+1)-0.25)/(double(n)+0.5));
+ x = cos(DOLFIN_PI*(double(i + 1) - 0.25)/(double(n) + 0.5));
// Newton's method
do
{
- dx = - p(x) / p.ddx(x);
+ dx = -Legendre::eval(n, x)/Legendre::ddx(n, x);
x = x + dx;
} while (real_abs(dx) > real_epsilon());
// Save the value using the symmetry of the points
- points[i] = - x;
- points[n-1-i] = x;
+ points[i] = -x;
+ points[n - 1 - i] = x;
}
// Set middle node
=== modified file 'dolfin/quadrature/GaussianQuadrature.cpp'
--- dolfin/quadrature/GaussianQuadrature.cpp 2011-05-11 13:13:14 +0000
+++ dolfin/quadrature/GaussianQuadrature.cpp 2011-06-01 22:19:07 +0000
@@ -72,10 +72,9 @@
// Compute the matrix coefficients
for (unsigned int i = 0; i < n; i++)
{
- Legendre p(i);
for (unsigned int j = 0; j < n; j++)
{
- A_real[i + n*j] = p(points[j]);
+ A_real[i + n*j] = Legendre::eval(i, points[j]);
_A(i, j) = to_double(A_real[i + n*j]);
_b[i] = 0.0;
b_real[i] = 0.0;
@@ -118,11 +117,12 @@
// value of the integral of the Legendre polynomial of degree q.
// This value should be zero for q > 0 and 2 for q = 0
- Legendre p(q);
+ //Legendre p(q);
+ Legendre p;
real sum = 0.0;
for (unsigned int i = 0; i < points.size(); i++)
- sum += weights[i]*p(points[i]);
+ sum += weights[i]*p.eval(q, points[i]);
//info("Checking quadrature weights: %.2e.", fabs(sum));
=== modified file 'dolfin/quadrature/LobattoQuadrature.cpp'
--- dolfin/quadrature/LobattoQuadrature.cpp 2011-05-11 13:13:14 +0000
+++ dolfin/quadrature/LobattoQuadrature.cpp 2011-06-01 22:19:07 +0000
@@ -89,7 +89,8 @@
return;
}
- Legendre p(n - 1);
+ //Legendre p(n - 1);
+ Legendre p;
real x, dx;
// Set the first and last nodal points which are 0 and 1
@@ -106,7 +107,7 @@
// Newton's method
do
{
- dx = -p.ddx(x)/p.d2dx(x);
+ dx = -p.ddx(n - 1, x)/p.d2dx(n - 1, x);
x = x + dx;
} while (real_abs(dx) > real_epsilon());
=== modified file 'dolfin/quadrature/RadauQuadrature.cpp'
--- dolfin/quadrature/RadauQuadrature.cpp 2011-05-11 13:13:14 +0000
+++ dolfin/quadrature/RadauQuadrature.cpp 2011-06-01 22:19:07 +0000
@@ -80,7 +80,6 @@
return;
}
- Legendre p(n);
real x, dx, step, sign;
// Set size of stepping for seeking starting points
@@ -93,20 +92,20 @@
x = -1.0 + step;
// Set the sign at -1 + epsilon
- sign = ((p.eval(n - 1, x) + p(x)) > 0 ? 1.0 : -1.0);
+ sign = (Legendre::eval(n - 1, x) + Legendre::eval(n, x) > 0 ? 1.0 : -1.0);
// Compute the rest of the nodes by Newton's method
for (unsigned int i = 1; i < n; i++)
{
// Step to a sign change
- while ((p.eval(n - 1, x) + p(x))*sign > 0.0)
+ while ((Legendre::eval(n - 1, x) + Legendre::eval(n, x))*sign > 0.0)
x += step;
// Newton's method
do
{
- dx = -(p.eval(n-1, x) + p(x))/(p.ddx(n - 1, x) + p.ddx(x));
+ dx = -(Legendre::eval(n-1, x) + Legendre::eval(n, x))/(Legendre::ddx(n - 1, x) + Legendre::ddx(n, x));
x = x + dx;
} while (real_abs(dx) > real_epsilon());
@@ -114,8 +113,8 @@
points[i] = x;
// Fix step so that it's not too large
- if (step > (points[i] - points[i-1])/10.0)
- step = (points[i] - points[i-1])/10.0;
+ if (step > (points[i] - points[i - 1])/10.0)
+ step = (points[i] - points[i - 1])/10.0;
// Step forward
sign = -sign;
=== modified file 'dolfin/swig/docstrings.i'
--- dolfin/swig/docstrings.i 2011-05-28 17:32:04 +0000
+++ dolfin/swig/docstrings.i 2011-06-01 22:19:07 +0000
@@ -703,11 +703,11 @@
%feature("docstring") dolfin::GenericDofMap::tabulate_coordinates "
**Overloaded versions**
-* tabulate_coordinates\ **(coordinates, ufc_cell)**
+* tabulate_coordinates\ **(boost::multi_array<double, coordinates, ufc_cell)**
Tabulate the coordinates of all dofs on a cell (UFC cell version)
-* tabulate_coordinates\ **(coordinates, cell)**
+* tabulate_coordinates\ **(boost::multi_array<double, coordinates, cell)**
Tabulate the coordinates of all dofs on a cell (DOLFIN cell version)
";
@@ -808,11 +808,11 @@
%feature("docstring") dolfin::DofMap::tabulate_coordinates "
**Overloaded versions**
-* tabulate_coordinates\ **(coordinates, ufc_cell)**
+* tabulate_coordinates\ **(boost::multi_array<double, coordinates, ufc_cell)**
Tabulate the coordinates of all dofs on a cell (UFC cell version)
-* tabulate_coordinates\ **(coordinates, cell)**
+* tabulate_coordinates\ **(boost::multi_array<double, coordinates, cell)**
Tabulate the coordinates of all dofs on a cell (DOLFIN cell version)
";
@@ -11161,10 +11161,6 @@
";
// Documentation extracted from: (module=math, header=basic.h)
-%feature("docstring") dolfin::sqr "
-Return the square of x
-";
-
%feature("docstring") dolfin::ipow "
Return a to the power n
";
@@ -11177,14 +11173,6 @@
Seed random number generator
";
-%feature("docstring") dolfin::near "
-Return true if x is within DOLFIN_EPS of x0
-";
-
-%feature("docstring") dolfin::between "
-Return true if x is between x0 and x1 (inclusive)
-";
-
// Documentation extracted from: (module=math, header=Lagrange.h)
%feature("docstring") dolfin::Lagrange "
Lagrange polynomial (basis) with given degree q determined by n = q + 1 nodal points.
=== modified file 'site-packages/dolfin/cppimports.py'
--- site-packages/dolfin/cppimports.py 2011-05-24 19:44:06 +0000
+++ site-packages/dolfin/cppimports.py 2011-06-01 22:19:07 +0000
@@ -102,7 +102,6 @@
# Classes from dolfin/common
# ------------------------
from cpp import MPI
-from cpp import near, between
# Classes from dolfin/parameter
# -----------------------------
--- End Message ---
