[cig-commits] commit: Fix rotations of initial data.

[Mercurial]

changeset:   121:fd0aa4755b4e
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Mar 27 16:06:14 2012 -0700
files:       initial.cxx
description:
Fix rotations of initial data.


diff -r af86ee9240e2 -r fd0aa4755b4e initial.cxx
--- a/initial.cxx	Mon Mar 26 10:55:06 2012 -0700
+++ b/initial.cxx	Tue Mar 27 16:06:14 2012 -0700
@@ -52,9 +52,11 @@ void initial(const Model &model, double 
   gsl_permutation_free(perm);
 
 
-  const double theta(0.1);
+  const double theta(pi);
+  // const double theta(0);
   FTensor::Tensor2<double,2,2>
-    rot(std::cos(theta),-std::sin(theta),std::sin(theta),std::cos(theta));
+    rot(std::cos(theta),std::sin(theta),-std::sin(theta),std::cos(theta));
+  FTensor::Tensor1<double,2> offset(0.5,0.5);
   const FTensor::Index<'a',2> a;
   const FTensor::Index<'b',2> b;
 
@@ -64,7 +66,7 @@ void initial(const Model &model, double 
         if(i<Nx && j<Ny)
           {
             FTensor::Tensor1<double,2> coord((i+0.5)*h,(j+0.5)*h), xy;
-            xy(a)=rot(a,b)*coord(b);
+            xy(a)=rot(a,b)*(coord(b)-offset(b)) + offset(a);
             double x(xy(0)), y(xy(1)), r2(x*x+y*y);
             switch(model)
               {
@@ -114,7 +116,7 @@ void initial(const Model &model, double 
         if(j<Ny)
           {
             FTensor::Tensor1<double,2> coord(i*h,(j+0.5)*h), xy;
-            xy(a)=rot(a,b)*coord(b);
+            xy(a)=rot(a,b)*(coord(b)-offset(b)) + offset(a);
             double x(xy(0)), y(xy(1)), r2(x*x+y*y);
             zx[i][j]=0;
             fx[i][j]=0;
@@ -124,6 +126,10 @@ void initial(const Model &model, double 
                 log_etax[i][j]=y*log_max_eta;
                 break;
               case Model::solCx:
+                {
+                  FTensor::Tensor1<double,2> f(0,pi*pi*sin(pi*y)*cos(pi*x));
+                  fx[i][j]=f(b)*rot(0,b);
+                }
                 distx[i][j]=x-middle;
                 if(x>middle)
                   {
@@ -152,10 +158,16 @@ void initial(const Model &model, double 
                       sign=-1;
                     }                      
                   const double csh(std::cosh(kx)), snh(std::sinh(kx)),
-                    sn(std::sin(kx));
-                  zx[i][j]=(tau_eta*snh - (v+tau_eta)*kx*csh
-                            - sign*rho*(kx*csh/2 - sn/2)/eta)
-                    *std::cos(pi*y)*eta;
+                    sn(std::sin(kx)), cs(std::cos(kx));
+
+                  FTensor::Tensor1<double,2>
+                    z((tau_eta*snh - (v+tau_eta)*kx*csh
+                       - sign*rho*(kx*csh/2 - sn/2)/eta)
+                      *std::cos(pi*y)*eta,
+                      (v*csh + (v+tau_eta)*kx*snh
+                       - sign*rho*(cs - csh - kx*snh)/(2*eta))
+                      *std::sin(pi*y)*eta);
+                  zx[i][j]=rot(0,b)*z(b);
                 }
                 break;
               case Model::sinker:
@@ -193,7 +205,7 @@ void initial(const Model &model, double 
         if(i<Nx)
           {
             FTensor::Tensor1<double,2> coord((i+0.5)*h, j*h), xy;
-            xy(a)=rot(a,b)*coord(b);
+            xy(a)=rot(a,b)*(coord(b)-offset(b)) + offset(a);
             double x(xy(0)), y(xy(1)), r2(x*x+y*y);
             zy[i][j]=0;
 
@@ -204,7 +216,10 @@ void initial(const Model &model, double 
                 log_etay[i][j]=y*log_max_eta;
                 break;
               case Model::solCx:
-                fy[i][j]=pi*pi*sin(pi*y)*cos(pi*x);
+                {
+                  FTensor::Tensor1<double,2> f(0,pi*pi*sin(pi*y)*cos(pi*x));
+                  fy[i][j]=f(b)*rot(1,b);
+                }
                 disty[i][j]=x-middle;
                 if(x>middle)
                   {
@@ -215,7 +230,6 @@ void initial(const Model &model, double 
                     log_etay[i][j]=log_min_eta;
                   }
                 {
-                  const double pi(atan(1.0)*4);
                   double kx, v, tau_eta, eta, sign;
                   if(x<middle)
                     {
@@ -233,11 +247,17 @@ void initial(const Model &model, double 
                       eta=max_eta;
                       sign=-1;
                     }                      
-                  const double csh(std::cosh(kx)), cs(std::cos(kx)),
-                    snh(std::sinh(kx));
-                  zy[i][j]=(v*csh + (v+tau_eta)*kx*snh
-                            - sign*rho*(cs - csh - kx*snh)/(2*eta))
-                    *std::sin(pi*y)*eta;
+                  const double snh(std::sinh(kx)), csh(std::cosh(kx)),
+                    sn(std::sin(kx)), cs(std::cos(kx));
+                    
+                  FTensor::Tensor1<double,2>
+                    z((tau_eta*snh - (v+tau_eta)*kx*csh
+                       - sign*rho*(kx*csh/2 - sn/2)/eta)
+                      *std::cos(pi*y)*eta,
+                      (v*csh + (v+tau_eta)*kx*snh
+                       - sign*rho*(cs - csh - kx*snh)/(2*eta))
+                      *std::sin(pi*y)*eta);
+                  zy[i][j]=rot(1,b)*z(b);
                 }
                 break;
               case Model::sinker:
@@ -266,7 +286,8 @@ void initial(const Model &model, double 
                     log_etay[i][j]=log_min_eta;
                     phi=-2*epsilon*A*r2_inclusion/z;
                     dphi=-phi/z;
-                    psi=-2*epsilon*(min_eta*z+A*r2_inclusion*r2_inclusion/(z*z*z));
+                    psi=-2*epsilon
+                      * (min_eta*z + A*r2_inclusion*r2_inclusion/(z*z*z));
                   }
                 v=(phi - z*conj(dphi) - conj(psi))/2.0;
                 zy[i][j]=v.imag();