.

include/sctl/sph_harm.hpp

@@ -127,7 +127,7 @@ template <class Real> class SphericalHarmonics{
     static void StokesEvalDL(const Vector<Real>& S, SHCArrange arrange, Long p, const Vector<Real>& coord, bool interior, Vector<Real>& U);
 
 
-    static void test() {
+    static void test_stokes() {
       int p = 6;
       int dof = 3;
       int Nt = 100, Np = 200;
@@ -231,8 +231,8 @@ template <class Real> class SphericalHarmonics{
         }
       };
 
-      for (Long i = 0; i < 20; i++) { // Evaluate
-        Real R0 = (1.01 + i/20.0);
+      for (Long i = 0; i < 40; i++) { // Evaluate
+        Real R0 = (0.01 + i/20.0);
 
         Vector<Real> x(3);
         x[0] = drand48()-0.5;
@@ -262,7 +262,7 @@ template <class Real> class SphericalHarmonics{
       Clear();
     }
 
-    static void test_() {
+    static void test() {
       int p = 3;
       int dof = 1;
       int Nt = p+1, Np = 2*p+1;

include/sctl/sph_harm.txx

@@ -699,7 +699,7 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalSL(const Vector<R
         if (interior) {
           Real a,b;
           a = n / (Real)((2*n+1) * (2*n+3)) * pow<Real>(R[i], n+1);
-          b = (n+1) / (Real)(4*n+2) * (pow<Real>(R[i], n-1) - pow<Real>(R[i], n+1));
+          b = -(n+1) / (Real)(4*n+2) * (pow<Real>(R[i], n-1) - pow<Real>(R[i], n+1));
           SVr = a * Vr + b * Wr;
           SVt = a * Vt + b * Wt;
           SVp = a * Vp + b * Wp;
@@ -860,7 +860,7 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalDL(const Vector<R
         if (interior) {
           Real a,b;
           a = -2*n*(n+2) / (Real)((2*n+1) * (2*n+3)) * pow<Real>(R[i], n+1);
-          b = (n+1)*(n+2) / (Real)(2*n+1) * (pow<Real>(R[i], n+1) - pow<Real>(R[i], n-1));
+          b = -(n+1)*(n+2) / (Real)(2*n+1) * (pow<Real>(R[i], n+1) - pow<Real>(R[i], n-1));
           SVr = a * Vr + b * Wr;
           SVt = a * Vt + b * Wt;
           SVp = a * Vp + b * Wp;