Add the capability to initialise spherical surfaces to a certain

coordinate radius.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/SphericalSurface/trunk@20 40f6ab95-0e4f-0410-8daa-ee8d7420be1d

5 files changed, 213 insertions, 3 deletions

diff --git a/param.ccl b/param.ccl
index 35efc0d..4967f9c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -67,3 +67,56 @@ BOOLEAN symmetric_y[42] "Reflection symmetry in the y direction"
 BOOLEAN symmetric_z[42] "Reflection symmetry in the z direction"
 {
 } no
+
+
+
+PRIVATE:
+
+# Place a surfaces at a certain location
+
+CCTK_REAL origin_x[42] "Origin for surface"
+{
+  * :: "origin"
+} 0.0
+
+CCTK_REAL origin_y[42] "Origin for surface"
+{
+  * :: "origin"
+} 0.0
+
+CCTK_REAL origin_z[42] "Origin for surface"
+{
+  * :: "origin"
+} 0.0
+
+
+
+BOOLEAN set_spherical[42] "Place surface at a certain radius"
+{
+} no
+
+CCTK_REAL radius[42] "Radius for surface"
+{
+  * :: "radius"
+} 0.0
+
+
+
+BOOLEAN set_elliptic[42] "Place surface at a certain radius"
+{
+} no
+
+CCTK_REAL radius_x[42] "Radius for surface"
+{
+  * :: "radius"
+} 0.0
+
+CCTK_REAL radius_y[42] "Radius for surface"
+{
+  * :: "radius"
+} 0.0
+
+CCTK_REAL radius_z[42] "Radius for surface"
+{
+  * :: "radius"
+} 0.0
diff --git a/schedule.ccl b/schedule.ccl
index b7feabc..ad7de7b 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -9,7 +9,7 @@ STORAGE: sf_shape_descriptors
 
 
 
-SCHEDULE SphericalSurfaceInfo_Setup AT basegrid
+SCHEDULE SphericalSurface_Setup AT basegrid
 {
   LANG: C
 } "Calculate surface coordinate descriptors"
@@ -19,3 +19,11 @@ SCHEDULE SphericalSurfaceInfo_Setup AT basegrid
 SCHEDULE GROUP SphericalSurface_HasBeenSet AT poststep
 {
 } "Set the spherical surfaces before this group, and use it afterwards"
+
+
+
+SCHEDULE SphericalSurface_Set AT poststep BEFORE SphericalSurface_HasBeenSet
+{
+  LANG: C
+  OPTIONS: global
+} "Set surface radii"
diff --git a/src/make.code.defn b/src/make.code.defn
index af99905..d394b7c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = setup.c
+SRCS = radius.c setup.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/radius.c b/src/radius.c
new file mode 100644
index 0000000..a670a48
--- /dev/null
+++ b/src/radius.c
@@ -0,0 +1,149 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+
+
+static
+CCTK_REAL min (CCTK_REAL const x, CCTK_REAL const y)
+{
+  return x < y ? x : y;
+}
+
+static
+CCTK_REAL max (CCTK_REAL const x, CCTK_REAL const y)
+{
+  return x > y ? x : y;
+}
+
+
+
+void SphericalSurface_Set (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_REAL const pi = 3.1415926535897932384626433832795028841971693993751;
+  
+  int n;
+  int i, j;
+  
+  
+  
+  for (n=0; n<nsurfaces; ++n) {
+    
+    if (set_spherical[n]) {
+      
+      sf_valid[n] = +1;
+      
+      sf_area[n] = 4 * pi * pow (radius[n], 2);
+      
+      sf_mean_radius[n] = radius[n];
+      
+      sf_centroid_x[n] = origin_x[n];
+      sf_centroid_y[n] = origin_y[n];
+      sf_centroid_z[n] = origin_z[n];
+      
+      sf_quadrupole_xx[n] = 0.0;
+      sf_quadrupole_xy[n] = 0.0;
+      sf_quadrupole_xz[n] = 0.0;
+      sf_quadrupole_yy[n] = 0.0;
+      sf_quadrupole_yz[n] = 0.0;
+      sf_quadrupole_zz[n] = 0.0;
+      
+      sf_min_radius[n] = radius[n];
+      sf_max_radius[n] = radius[n];
+      
+      sf_min_x[n] = origin_x[n] - radius[n];
+      sf_min_y[n] = origin_y[n] - radius[n];
+      sf_min_z[n] = origin_z[n] - radius[n];
+      sf_max_x[n] = origin_x[n] + radius[n];
+      sf_max_y[n] = origin_y[n] + radius[n];
+      sf_max_z[n] = origin_z[n] + radius[n];
+      
+      for (j=0; j<nphi[n]; ++j) {
+        for (i=0; i<ntheta[n]; ++i) {
+          sf_radius[n] = radius[n];
+        }
+      }
+      
+      sf_origin_x[n] = origin_x[n];
+      sf_origin_y[n] = origin_y[n];
+      sf_origin_z[n] = origin_z[n];
+      
+    } else if (set_elliptic[n]) {
+      
+      /*
+       * Equation for an ellipsoid:
+       *    x^2 / rx^2 + y^2 / ry^2 + z^2/ rz^2 = 1
+       * where rx, ry, rz are the radii in the x, y, and z directions.
+       *
+       * With
+       *    x = r (sin theta) (cos phi)
+       *    y = r (sin theta) (sin phi)
+       *    z = r (cos theta)
+       *
+       * we get the solution
+       *    1/r^2 = x^2/rx^2 + y^2/ry^2 + z^2/rz^2
+       *
+       * With the form
+       *    x A x = 1
+       * we obtain
+       *    A = diag [1/rx^2, 1/ry^2, 1/rz^2]
+       */
+      
+      CCTK_REAL const rx2 = pow (radius_x[n], 2);
+      CCTK_REAL const ry2 = pow (radius_y[n], 2);
+      CCTK_REAL const rz2 = pow (radius_z[n], 2);
+      
+      sf_valid[n] = +1;
+      
+      sf_area[n] = 0 * 4 * pi * pow (radius[n], 2);
+      
+      sf_mean_radius[n] = 0 * radius[n];
+      
+      sf_centroid_x[n] = origin_x[n];
+      sf_centroid_y[n] = origin_y[n];
+      sf_centroid_z[n] = origin_z[n];
+      
+      sf_quadrupole_xx[n] = 1.0 / rz2;
+      sf_quadrupole_xy[n] = 0.0;
+      sf_quadrupole_xz[n] = 0.0;
+      sf_quadrupole_yy[n] = 1.0 / ry2;
+      sf_quadrupole_yz[n] = 0.0;
+      sf_quadrupole_zz[n] = 1.0 / rx2;
+      
+      sf_min_radius[n] = min (radius_x[n], min (radius_y[n], radius_z[n]));
+      sf_max_radius[n] = max (radius_x[n], max (radius_y[n], radius_z[n]));
+      
+      sf_min_x[n] = origin_x[n] - radius_x[n];
+      sf_min_y[n] = origin_y[n] - radius_y[n];
+      sf_min_z[n] = origin_z[n] - radius_z[n];
+      sf_max_x[n] = origin_x[n] + radius_x[n];
+      sf_max_y[n] = origin_y[n] + radius_y[n];
+      sf_max_z[n] = origin_z[n] + radius_z[n];
+      
+      for (j=0; j<nphi[n]; ++j) {
+        for (i=0; i<ntheta[n]; ++i) {
+          CCTK_REAL const theta = sf_origin_theta[n] + i * sf_delta_theta[n];
+          CCTK_REAL const phi = sf_origin_phi[n] + j * sf_delta_phi[n];
+          CCTK_REAL const x2 = pow (sin(theta) * cos(phi), 2);
+          CCTK_REAL const y2 = pow (sin(theta) * sin(phi), 2);
+          CCTK_REAL const z2 = pow (cos(theta)           , 2);
+          sf_radius[n] = 1.0 / sqrt (x2 / rx2 + y2 / ry2 + z2 / rz2);
+        }
+      }
+    }
+    
+    sf_origin_x[n] = origin_x[n];
+    sf_origin_y[n] = origin_y[n];
+    sf_origin_z[n] = origin_z[n];
+    
+  } /* for n */
+}
diff --git a/src/setup.c b/src/setup.c
index f30d1ce..5c929ff 100644
--- a/src/setup.c
+++ b/src/setup.c
@@ -9,7 +9,7 @@
 
 
 
-void SphericalSurfaceInfo_Setup (CCTK_ARGUMENTS)
+void SphericalSurface_Setup (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;