[cig-commits] r20621 - mc/3D/CitcomS/trunk/lib

[becker at geodynamics.org]

Author: becker
Date: 2012-08-22 14:21:52 -0700 (Wed, 22 Aug 2012)
New Revision: 20621

Modified:
   mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
Log:
Implemented improved tracer fix by Rene.

"a new method of finding the element of a tracer that is too close to
an element boundary. Up to now such a tracer was shifted by a constant
epsilon theta and phi. If the element boundary is parallel to this
theta/phi direction it is not guaranteed to work well (thus the
number_of_tries check), and additionally I got the problem that
sometimes all elements refuse this tracer. Eh checked in a workaround
for this (r15742), deleting orphan tracers in Tracer_setup.c. Because
I did not know this, I created my own bugfix moving the tracers an
epsilon amount orthogonal to all boundaries that are too close. In
order to save computing time I use the already computed vectors for
the element boundaries (this assumes that the element boundaries are
nearly orthogonal to each other, but unless somebody tries to change
CitcomS elements that should work fine). The shift happens now in
cartesian coordinates since the boundary-vectors are cartesian and the
radius-coordinate of the tracer is normalized prior to this check
anyway, so I just need to re-normalize it after the shift. For now I
did not touch all the now (hopefully) useless security checks but as
far as I can see they do no harm either, so we can remove them later."




Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c	2012-08-22 19:06:20 UTC (rev 20620)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c	2012-08-22 21:21:52 UTC (rev 20621)
@@ -2315,12 +2315,10 @@
     double cross3[4];
     double cross4[4];
     double rad1,rad2,rad3,rad4;
-    double theta, phi;
+    double theta, phi,rad;
     double tiny, eps;
     double x,y,z;
 
-    double myatan();
-
     /* make vectors from node to node */
 
     makevector(v12,rnode2[1],rnode2[2],rnode2[3],rnode1[1],rnode1[2],rnode1[3]);
@@ -2357,7 +2355,7 @@
     /*  Hopefully, this doesn't happen often, may be expensive                  */
 
     tiny=1e-15;
-    eps=1e-6;
+    eps=1e-3;
 
     if (number_of_tries>3)
         {
@@ -2373,34 +2371,37 @@
         }
 
     if (fabs(rad1)<=tiny||fabs(rad2)<=tiny||fabs(rad3)<=tiny||fabs(rad4)<=tiny)
-        {
-            x=test_point[1];
-            y=test_point[2];
-            z=test_point[3];
-            theta=myatan(sqrt(x*x+y*y),z);
-            phi=myatan(y,x);
+      {
+	if (fabs(rad1) <= tiny){
+	  test_point[1] += v12[1] * eps; 
+	  test_point[2] += v12[2] * eps; 
+	  test_point[3] += v12[3] * eps; 
+	}
+	if (fabs(rad2) <= tiny){
+	  test_point[1] += v23[1] * eps; 
+	  test_point[2] += v23[2] * eps; 
+	  test_point[3] += v23[3] * eps; 
+	}
+	if (fabs(rad3) <= tiny){
+	  test_point[1] += v34[1] * eps; 
+	  test_point[2] += v34[2] * eps; 
+	  test_point[3] += v34[3] * eps; 
+	}
+	if (fabs(rad4) <= tiny){
+	  test_point[1] += v41[1] * eps; 
+	  test_point[2] += v41[2] * eps; 
+	  test_point[3] += v41[3] * eps; 
+	}
+	rad = sqrt(test_point[1]*test_point[1]+test_point[2]*test_point[2]+test_point[3]*test_point[3]);
+	test_point[1] /= rad;
+	test_point[2] /= rad;
+	test_point[3] /= rad;
+	
+	number_of_tries++;
+	goto try_again;
+	
+      }
 
-            if (theta<=M_PI/2.0)
-                {
-                    theta=theta+eps;
-                }
-            else
-                {
-                    theta=theta-eps;
-                }
-            phi=phi+eps;
-            x=sin(theta)*cos(phi);
-            y=sin(theta)*sin(phi);
-            z=cos(theta);
-            test_point[1]=x;
-            test_point[2]=y;
-            test_point[3]=z;
-
-            number_of_tries++;
-            goto try_again;
-
-        }
-
     icheck=0;
     if (rad1>0.0&&rad2>0.0&&rad3>0.0&&rad4>0.0) icheck=1;