grad(u) + grad(u)^T

include/sctl/sph_harm.hpp

@@ -280,7 +280,7 @@ template <class Real> class SphericalHarmonics{
             du[7] = ((dr[2]*f[1]-dr[1]*f[2])*invr3 - 3*dr[2]*dr[1]*fdotr*invr5) * scal;
             du[8] = (                  fdotr*invr3 - 3*dr[2]*dr[2]*fdotr*invr5) * scal;
 
-            Real p = (2*fdotr*invr3) * scal;
+            Real p = 0;//(2*fdotr*invr3) * scal;
 
             Real K[9];
             K[0] = du[0] + du[0] - p; K[1] = du[1] + du[3] - 0; K[2] = du[2] + du[6] - 0;
@@ -319,30 +319,6 @@ template <class Real> class SphericalHarmonics{
         stokes_evalDL(x, Df_);
         stokes_evalKL(x, n, Kf_);
 
-        {
-          Vector<Real> dx(COORD_DIM);
-          if (0) {
-            dx = x;
-          } else {
-            auto cross_prod = [](Vector<Real>& x, Vector<Real>& y) {
-              Vector<Real> z(COORD_DIM);
-              z[0] = x[1] * y[2] - x[2] * y[1];
-              z[1] = x[2] * y[0] - x[0] * y[2];
-              z[2] = x[0] * y[1] - x[1] * y[0];
-              return z;
-            };
-            Vector<Real> z(COORD_DIM); z = 0; z[2] = 1;
-            dx = cross_prod(x,z);
-            dx = cross_prod(dx,x) * -1;
-          }
-
-          Vector<Real> v0, v1;
-          Real dr = sqrt<Real>(dx[0]*dx[0]+dx[1]*dx[1]+dx[2]*dx[2]);
-          stokes_evalSL(x - dx * 1e-5 / dr, v0);
-          stokes_evalSL(x + dx * 1e-5 / dr, v1);
-          Kf_ = (v1 - v0) / 2e-5;
-        }
-
         auto errSL = (Sf-Sf_)/(Sf+0.01);
         auto errDL = (Df-Df_)/(Df+0.01);
         auto errKL = (Kf-Kf_)/(Kf+0.01);

include/sctl/sph_harm.txx

@@ -1152,22 +1152,22 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalKL(const Vector<R
             Xt_t += ( Xr) / R[i];
           }
           { // Set Vr_p, Vt_p, Vp_p, Wr_p, Wt_p, Wp_p, Xr_p, Xt_p, Xp_p
-            auto C0 = Ynm_1p;
-            auto C1 = Ynm_1p * csc_theta;
-            auto C2 = (Anm * Ynm_2p - Bnm * Ynm_0p) * csc_theta;
+            auto C0 = -Ynm_1p;
+            auto C1 = -Ynm_1p * csc_theta;
+            auto C2 = -(Anm * Ynm_2p - Bnm * Ynm_0p) * csc_theta;
             SetVecSH(Vr_p, Vt_p, Vp_p, Wr_p, Wt_p, Wp_p, Xr_p, Xt_p, Xp_p, C0, C1, C2);
 
-            Vr_p -= (-sin_theta * Vp                 ) * csc_theta / R[i];
-            Vt_p -= (-cos_theta * Vp                 ) * csc_theta / R[i];
-            Vp_p -= ( sin_theta * Vr + cos_theta * Vt) * csc_theta / R[i];
+            Vr_p += (-sin_theta * Vp                 ) * csc_theta / R[i];
+            Vt_p += (-cos_theta * Vp                 ) * csc_theta / R[i];
+            Vp_p += ( sin_theta * Vr + cos_theta * Vt) * csc_theta / R[i];
 
-            Wr_p -= (-sin_theta * Wp                 ) * csc_theta / R[i];
-            Wt_p -= (-cos_theta * Wp                 ) * csc_theta / R[i];
-            Wp_p -= ( sin_theta * Wr + cos_theta * Wt) * csc_theta / R[i];
+            Wr_p += (-sin_theta * Wp                 ) * csc_theta / R[i];
+            Wt_p += (-cos_theta * Wp                 ) * csc_theta / R[i];
+            Wp_p += ( sin_theta * Wr + cos_theta * Wt) * csc_theta / R[i];
 
-            Xr_p -= (-sin_theta * Xp                 ) * csc_theta / R[i];
-            Xt_p -= (-cos_theta * Xp                 ) * csc_theta / R[i];
-            Xp_p -= ( sin_theta * Xr + cos_theta * Xt) * csc_theta / R[i];
+            Xr_p += (-sin_theta * Xp                 ) * csc_theta / R[i];
+            Xt_p += (-cos_theta * Xp                 ) * csc_theta / R[i];
+            Xp_p += ( sin_theta * Xr + cos_theta * Xt) * csc_theta / R[i];
           }
         }
 
@@ -1256,24 +1256,18 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalKL(const Vector<R
           SX[2][2] = a * Xp_p;
         }
 
-        { //////////////////////
-          norm0[0] = 0;
-          norm0[1] = 1;
-          norm0[2] = 0;
-        }
-
         Complex<Real> KV[COORD_DIM][COORD_DIM], KW[COORD_DIM][COORD_DIM], KX[COORD_DIM][COORD_DIM];
-        KV[0][0] = SV[0][0] + SV[0][0]*0;   KV[0][1] = SV[0][1] + SV[0][1]*0;   KV[0][2] = SV[0][2] + SV[0][2]*0;
-        KV[1][0] = SV[1][0] + SV[1][0]*0;   KV[1][1] = SV[1][1] + SV[1][1]*0;   KV[1][2] = SV[1][2] + SV[1][2]*0;
-        KV[2][0] = SV[2][0] + SV[2][0]*0;   KV[2][1] = SV[2][1] + SV[2][1]*0;   KV[2][2] = SV[2][2] + SV[2][2]*0;
+        KV[0][0] = SV[0][0] + SV[0][0];   KV[0][1] = SV[0][1] + SV[1][0];   KV[0][2] = SV[0][2] + SV[2][0];
+        KV[1][0] = SV[1][0] + SV[0][1];   KV[1][1] = SV[1][1] + SV[1][1];   KV[1][2] = SV[1][2] + SV[2][1];
+        KV[2][0] = SV[2][0] + SV[0][2];   KV[2][1] = SV[2][1] + SV[1][2];   KV[2][2] = SV[2][2] + SV[2][2];
 
-        KW[0][0] = SW[0][0] + SW[0][0]*0;   KW[0][1] = SW[0][1] + SW[0][1]*0;   KW[0][2] = SW[0][2] + SW[0][2]*0;
-        KW[1][0] = SW[1][0] + SW[1][0]*0;   KW[1][1] = SW[1][1] + SW[1][1]*0;   KW[1][2] = SW[1][2] + SW[1][2]*0;
-        KW[2][0] = SW[2][0] + SW[2][0]*0;   KW[2][1] = SW[2][1] + SW[2][1]*0;   KW[2][2] = SW[2][2] + SW[2][2]*0;
+        KW[0][0] = SW[0][0] + SW[0][0];   KW[0][1] = SW[0][1] + SW[1][0];   KW[0][2] = SW[0][2] + SW[2][0];
+        KW[1][0] = SW[1][0] + SW[0][1];   KW[1][1] = SW[1][1] + SW[1][1];   KW[1][2] = SW[1][2] + SW[2][1];
+        KW[2][0] = SW[2][0] + SW[0][2];   KW[2][1] = SW[2][1] + SW[1][2];   KW[2][2] = SW[2][2] + SW[2][2];
 
-        KX[0][0] = SX[0][0] + SX[0][0]*0;   KX[0][1] = SX[0][1] + SX[0][1]*0;   KX[0][2] = SX[0][2] + SX[0][2]*0;
-        KX[1][0] = SX[1][0] + SX[1][0]*0;   KX[1][1] = SX[1][1] + SX[1][1]*0;   KX[1][2] = SX[1][2] + SX[1][2]*0;
-        KX[2][0] = SX[2][0] + SX[2][0]*0;   KX[2][1] = SX[2][1] + SX[2][1]*0;   KX[2][2] = SX[2][2] + SX[2][2]*0;
+        KX[0][0] = SX[0][0] + SX[0][0];   KX[0][1] = SX[0][1] + SX[1][0];   KX[0][2] = SX[0][2] + SX[2][0];
+        KX[1][0] = SX[1][0] + SX[0][1];   KX[1][1] = SX[1][1] + SX[1][1];   KX[1][2] = SX[1][2] + SX[2][1];
+        KX[2][0] = SX[2][0] + SX[0][2];   KX[2][1] = SX[2][1] + SX[1][2];   KX[2][2] = SX[2][2] + SX[2][2];
 
 
         write_coeff(KV[0][0]*norm0[0] + KV[0][1]*norm0[1] + KV[0][2]*norm0[2], n, m, 0, 0);