Surface Descriptions
--------------------

Author:     Erik Schnetter
Date:       September 2003
References: http://www.cactuscode.org/archives_html/developers/msg00928.html


Last Thursday the idea was brought up to have a common representation of 
closed 2-spheres in Cactus, so that thorns using such spheres can 
interoperate.  Possible applications include apparent horizons, event 
horizons, radiation extraction, inner an outer boundary conditions, and 
probably many more.

The requirements that I have towards such a thorn are:
1. must be able to represent several such surfaces
2. the surfaces may have symmetries, e.g. only one quadrant of such a 
surface might exist
3. it should be implemented "the Cactus way", so that input, output, and 
parallelisation come for free.

These items mean that certain representation are not possible.  I 
suggest settling for a r/theta/phi representation, which seems to be 
good enough.  Note that this is just a representation of the surface 
itself, which is nowhere singular, as r is well-defined everywhere.  
All thorns performing calculations can internally switch to other 
representations.  I believe r/theta/phi is therefore an acceptable 
choice.

I've written a suggestion of what such a thorn could look like.  Have a 
look at the enclosed files.  I'm quite sure that you will have many 
objections to my choice of variable names and indentation.  Once you've 
piped them to /dev/null, please respond with constructive criticism.  
My opinion is to start with something simple and add features once they 
are needed, other people's opinions might differ.

My implementation has two drawbacks, that I would be willing to accept:

1. Cactus requires that all surfaces have the same number of grid 
points.  That's life.  I'm sure Tom will be pleased to implement the 
multi-model stuff at some time in the future, which is when we will be 
able to activate the same thorn multiple times with different numbers 
of grid points.  Until then this seems to be the best way in the Cactus 
world.

2. The possible surface symmetries make it non-trivial to calculate the 
grid spacings dtheta and dphi.  I have therefore implemented aliased 
functions that perform these calculations.

3. IOASCII does not allow outputting just a single surface, it will 
output all of them.  We can live with that, or change IOASCII.  A good 
thing would be to output different surfaces to different files.

Enjoy,
- -erik

PS: Obvious suggestions would be to add grid functions for theta and 
phi, or for the grid point coordinates x, y, and z.

===============================================================================

CVS info   : $Header: /cactus/CactusWebSite/Development/Specs/SurfaceDescriptions.txt,v 1.1 2003/10/09 15:48:12 tradke Exp $

Cactus Code Thorn Surface
Thorn Author(s)     : Erik Schnetter <schnetter@aei.mpg.de>
Thorn Maintainer(s) : Erik Schnetter <schnetter@aei.mpg.de>
--------------------------------------------------------------------------

Purpose of the thorn:

Provide a common description for closed 2-spheres, so that thorns
using such spheres can interoperate.  Possible applications include
apparent horizons, event horizons, radiation extraction, inner and
outer boundary conditions, and many more.

# Interface definition for thorn Surface
# $Header: /cactus/CactusWebSite/Development/Specs/SurfaceDescriptions.txt,v 1.1 2003/10/09 15:48:12 tradke Exp $

IMPLEMENTS: Surface



# Commenting names
CCTK_STRING name[nsurfaces] TYPE=scalar

# Free or used
CCTK_INT used[nsurfaces] TYPE=scalar



# The coordinate radius as function of theta and phi.
# The grid points are staggered about the origin in the theta direction,
# and not staggered in the phi direction.
# (When there is a symmetry in the z direction, then the grid points are
#  also not staggered about the equator.)
CCTK_REAL radius[nsurfaces] TYPE=array DIM=2 SIZE=ntheta,nphi GHOSTSIZE=nghoststheta,nghostsphi TIMELEVELS=3



# The coordinate centre of the sphere
CCTK_REAL centre_x[nsurfaces] TYPE=scalar TIMELEVELS=3
CCTK_REAL centre_y[nsurfaces] TYPE=scalar TIMELEVELS=3
CCTK_REAL centre_z[nsurfaces] TYPE=scalar TIMELEVELS=3



# Allocate and free surfaces
CCTK_INT FUNCTION \
Surface_Allocate (CCTK_POINTER_TO_CONST IN cctkGH, \
                  CCTK_STRING IN thename)

FUNCTION \
Surface_Free (CCTK_POINTER_TO_CONST IN cctkGH, \
              CCTK_INT IN n)

# Grid spacing of a surface
FUNCTION \
Surface_spacing (CCTK_INT IN n, \
                 CCTK_REAL OUT theta0, \
                 CCTK_REAL OUT phi0, \
                 CCTK_REAL OUT dtheta, \
                 CCTK_REAL OUT dphi)

# Coordinates of surface grid points
FUNCTION \
Surface_rthetaphi (CCTK_POINTER_TO_CONST IN cctkGH, \
                   CCTK_INT IN n, \
                   CCTK_INT IN tl, \
                   CCTK_INT IN i, \
                   CCTK_INT IN j, \
                   CCTK_REAL OUT r, \
                   CCTK_REAL OUT theta, \
                   CCTK_REAL OUT phi)

FUNCTION \
Surface_xyz (CCTK_POINTER_TO_CONST IN cctkGH, \
             CCTK_INT IN n, \
             CCTK_INT IN tl, \
             CCTK_INT IN i, \
             CCTK_INT IN j, \
             CCTK_REAL OUT ARRAY xyz)

# Parameter definitions for thorn Surface
# $Header: /cactus/CactusWebSite/Development/Specs/SurfaceDescriptions.txt,v 1.1 2003/10/09 15:48:12 tradke Exp $



# The following parameters are needed to declare grid arrays
# in interface.ccl files, and hence must be public.
RESTRICTED:

# The constant "20" must be the same here and below in symmetric_*.
CCTK_INT nsurfaces "Number of surfaces"
{
  0:20 :: ""
} 0



CCTK_INT ntheta "Number of grid points in the theta direction"
{
  3:*:2 :: ""
} 19

CCTK_INT nphi "Number of grid points in the theta direction"
{
  3:* :: ""
} 38



CCTK_INT nghoststheta "Number of ghost zones in the theta direction"
{
  0:* :: ""
} 1

CCTK_INT nghostsphi "Number of ghost zones in the phi direction"
{
  0:* :: ""
} 1



PRIVATE:

# Indicators of sphere symmetries
# 0: no symmetry
# 1: there is a symmetriy, so that only that part of the sphere in the
#    nonnegative half-space is stored.
BOOLEAN symmetric_x[20] "Is there a symmetry in the x direction?"
{
} no

BOOLEAN symmetric_y[20] "Is there a symmetry in the y direction?"
{
} no

BOOLEAN symmetric_z[20] "Is there a symmetry in the z direction?"
{
} no

# Schedule definitions for thorn Surface
# $Header: /cactus/CactusWebSite/Development/Specs/SurfaceDescriptions.txt,v 1.1 2003/10/09 15:48:12 tradke Exp $

STORAGE: name used

SCHEDULE Surface_Init AT basegrid
{
  LANG: C
} "Initialise administrative stuff"

/* $Header: /cactus/CactusWebSite/Development/Specs/SurfaceDescriptions.txt,v 1.1 2003/10/09 15:48:12 tradke Exp $ */

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"


  
static int initialized = 0;

static CCTK_REAL pi;

static int name_index;
static int used_index;

static int radius_index;
static int centre_index[3];



void
Surface_Init (CCTK_ARGUMENTS)
{
  int n;
  
  pi = acos(-1.0);
  
  name_index = CCTK_VarIndex ("Surface::name");
  assert (name_index>=0);
  used_index = CCTK_VarIndex ("Surface::used");
  assert (used_index>=0);
  
  radius_index = CCTK_VarIndex ("Surface::radius");
  assert (radius_index>=0);
  
  centre_index[0] = CCTK_VarIndex ("Surface::centre_x");
  assert (centre_index[0]>=0);
  centre_index[1] = CCTK_VarIndex ("Surface::centre_y");
  assert (centre_index[1]>=0);
  centre_index[2] = CCTK_VarIndex ("Surface::centre_z");
  assert (centre_index[2]>=0);
  
  for (n=0; n<nsurfaces; ++n) {
    used[n] = 0;
  }
}



CCTK_INT
Surface_Allocate (cGH const * restrict const cctkGH,
                  char const * restrict const thename)
{
  CCTK_INT * restrict used;
  CCTK_INT * restrict name;
  int n;
  
  assert (thename);
  
  used = CCTK_VarDataPtrI (cctkGH, used_index);
  assert (used);
  name = CCTK_VarDataPtrI (cctkGH, name_index);
  assert (name);
  
  for (n=0; n<nsurfaces; ++n) {
    if (! used[n]) {
      
      used[n] = 1;
      name[n] = strdup (thename);
      return n;
      
    }
  }
  
  return -1;
}



void
Surface_Free (cGH const * restrict const cctkGH,
              CCTK_INT const n)
{
  CCTK_INT * restrict used;
  CCTK_INT * restrict name;
  
  assert (n>=0 && n<nsurfaces);
  
  used = CCTK_VarDataPtrI (cctkGH, used_index);
  assert (used);
  name = CCTK_VarDataPtrI (cctkGH, name_index);
  assert (name);
  
  assert (used[n]);
  
  free (name[n]);
  name[n] = 0;
  used[n] = 0;
}



void
Surface_spacing (CCTK_INT const n,
                 CCTK_REAL * restrict const theta0,
                 CCTK_REAL * restrict const phi0,
                 CCTK_REAL * restrict const dtheta,
                 CCTK_REAL * restrict const dphi)
{
  DECLARE_CCTK_PARAMETERS;
  
  /* Check arguments */
  assert (n>=0 && n<nsurfaces);
  assert (theta0);
  assert (phi0);
  assert (dtheta);
  assert (dphi);
  
  /* Calculate the grid spacing in the theta direction:
     depends on the symmetry in the z direction */
  if (! symmetric_z[n]) {
    
    /* no symmetry in the z direction
       extent: north pole to south pole
       north pole: theta = nghosts-1/2
       south pole: theta = (ntheta-1) - (nghosts-1/2) */
    *dtheta = pi / (ntheta - 2*nghoststheta + 1);
    *theta0 = (- nghoststheta + 0.5) * *dtheta;
    
  } else {
    
    /* symmetry in the z direction (bitant mode in z direction)
       extent: north pole to equator
       north pole: theta = nghosts-1/2
       equator:    theta = (ntheta-1) - nghosts */
    *dtheta = pi / (ntheta - 2*nghoststheta + 0.5);
    *theta0 = (- nghoststheta + 0.5) * *dtheta;
    
  }
  
  /* Calculate the grid spacing in the phi direction:
     depends on the symmetry in the x and y directions */
  if (! symmetric_x[n]) {
    if (! symmetric_y[n]) {
      
      /* no symmetries in the x and y directions (full mode)
         extent: 0 meridian to 2pi meridian
         0 meridian:   phi = nghosts
         2pi meridian: phi = nphi - nghosts
         (note that 2pi is a ghost point) */
      *dphi = 2*pi / (nphi - 2*nghosts);
      *phi0 = (- nghostsphi) * *dphi;
      
    } else {
      
      /* symmetry only in the y direction (bitant mode in y direction)
         extent: -pi/2 meridian to +pi/2 meridian
         -pi/2 meridian: phi = nghosts
         +pi/2 meridian:  phi = (nphi-1) - nghosts */
      *dphi = pi / (nphi+1 - 2*nghosts);
      *phi0 = -pi/2 + (- nghostsphi) * *dphi;
      
    }
  } else {
    if (! symmetric_y[n]) {
      
      /* symmetry only in the x direction (bitant mode in x direction)
         extent: pi meridian to 2pi meridian
         pi meridian:  phi = nghosts
         2pi meridian: phi = (nphi+1) - nghosts */
      *dphi = pi / (nphi+1 - 2*nghosts);
      *phi0 = pi + (- nghostsphi) * *dphi;
      
    } else {
      
      /* symmetry in the x and y directions (quadrant mode in xy direction)
         extent: 3pi/2 meridian to 2pi meridian
         3pi/2 meridian: phi = nghosts
         2pi meridian:   phi = (nphi+1) - nghosts */
      *dphi = pi/2 / (nphi+1 - 2*nghosts);
      *phi0 = 3*pi/2 + (- nghostsphi) * *dphi;
      
    }
  }
}



void
Surface_rthetaphi (cGH const * restrict const cctkGH,
                   CCTK_INT const n,
                   CCTK_INT const tl,
                   CCTK_INT const i,
                   CCTK_INT const j,
                   CCTK_REAL * restrict const r,
                   CCTK_REAL * restrict const theta,
                   CCTK_REAL * restrict const phi)
{
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_REAL const * restrict radius;
  CCTK_REAL theta0, phi0, dtheta, dphi;
  
  /* Check arguments */
  assert (cctkGH);
  assert (n>=0 && n<nsurfaces);
  assert (tl>=0 && tl<CCTK_ActiveTimeLevelsVI(radius_index));
  assert (i>=0 && i<ntheta);
  assert (j>=0 && j<nphi);
  assert (r);
  assert (theta);
  assert (phi);
  
  /* Calculate grid spacings */
  Surface_spacings (n, &theta0, &phi0, &dtheta, &dphi);
  
  /* Get grid variables */
  radius =  CCTK_VarDataPtrI (cctkGH, tl, radius_index);
  assert (radius);
  
  /* Calculate the radius: it is stored */
  *r = radius[i * ntheta + j];
  
  /* Calculate theta and phi: from spacings */
  *theta = theta0 + i * dtheta;
  *phi   = phi0   + j * dphi;
}



void
Surface_xyz (cGH const * restrict const cctkGH,
             CCTK_INT const n,
             CCTK_INT const tl,
             CCTK_INT const i,
             CCTK_INT const j,
             CCTK_REAL * restrict const xyz)
{
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_REAL const * restrict centre[3];
  CCTK_REAL r, theta, phi;
  int d;
  
  /* Check arguments */
  assert (cctkGH);
  assert (n>=0 && n<nsurfaces);
  assert (tl>=0 && tl<CCTK_ActiveTimeLevelsVI(radius_index));
  assert (i>=0 && i<ntheta);
  assert (j>=0 && j<nphi);
  assert (xyz);
  
  /* Get grid variables */
  for (d=0; d<3; ++d) {
    centre[d] =  CCTK_VarDataPtrI (cctkGH, tl, centre_index[d]);
    assert (centre[d]);
  }
  
  /* Calculate r, theta, and phi */
  Surface_rthetaphi (cctkGH, n, tl, i, j, &r, &theta, &phi);
  
  /* Calculate x, y, and z */
  *xyz[0] = centre[0] + r * sin(theta) * cos(phi);
  *xyz[1] = centre[1] + r * sin(theta) * sin(phi);
  *xyz[2] = centre[2] + r * cos(theta);
}