Initial add to the repository.

These are the basic files needed to re-implement isosurfacer in C.
This has mostly been cleaned of C++ code, but hasn't been tested out
yet.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IsoSurfacer/trunk@2 bfcf8e34-485d-4d46-a995-1fd6fa6fb178

15 files changed, 2526 insertions, 0 deletions

diff --git a/COPYRIGHT b/COPYRIGHT
new file mode 100644
index 0000000..57b2896
--- /dev/null
+++ b/COPYRIGHT
@@ -0,0 +1,352 @@
+
+
+This thorn is (c) Copyright the authors listed in the sources files. This
+thorn is distributed under the GNU GPL with the specific exception that
+the GNU GPL does not migrate to code linking to or using infrastructure
+from this thorn.
+
+- The Cactus Team <cactus@cactuscode.org>
+
+**************************************************************************
+
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                          675 Mass Ave, Cambridge, MA 02139, USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	Appendix: How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) 19yy  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) 19yy name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
+
+
diff --git a/README b/README
new file mode 100644
index 0000000..fadcb89
--- /dev/null
+++ b/README
@@ -0,0 +1,20 @@
+Cactus Code Thorn IsoSurfacer
+Authors    : 
+Managed by : 
+Version    : 2.0
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn extracts isosurfaces of grid functions.
+
+2. Dependencies of the thorn
+
+This thorn additionally requires thorns PUGH.
+
+3. Thorn distribution
+
+This thorn is available to public.
+
+4. Additional information
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..0601d30
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,10 @@
+# Interface definition for thorn isosurfacer
+# $Header$
+
+implements: isosurfacer
+inherits: Grid
+
+# USES INCLUDE HEADER: sstream.hpp
+# USES INCLUDE HEADER: Server.hpp
+# USES INCLUDE HEADER: CactusCommands.hpp
+# USES INCLUDE HEADER: http_GHExtensions.hpp # no connection to HTTP yet
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..a4140dd
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,40 @@
+# Parameter definitions for thorn IsoSurfacer
+# $Header$
+
+private:
+
+STRING isosurfacer "What to isosurface and how. Format: {(functionName) (isolevel1, isolevel2, ...) (format1, format2, ...) (itiso, resolution, firstIteration, uniqverts)} {} ..."
+{
+  .* :: A regex which matches everything
+} ""
+
+INT outer_boundary_cutoff "Number of voxels to cut off the outer boundaries" STEERABLE = ALWAYS
+{
+ 0:* :: "Zero or any positive number"
+} 0
+
+STRING format_str "If the precision specified in the default format string is less than what you need, then insert your own format for floating point numbers here"
+{
+  .* :: A regex which matches everything
+} "%3.3f"
+
+BOOLEAN allow_empty_sends "Allow the Isosurfacer to send zero-length vertex lists to the client application ? Ordinarily these zero-length sends are supressed." STEERABLE = ALWAYS
+{
+} "no"
+
+STRING outdir "Output directory for isosurface data files"
+{
+  .* :: A regex which matches everything
+} "."
+
+#############################################################################
+### import IOUtil parameters
+#############################################################################
+shares: IO
+
+####################
+# verbose flag
+####################
+USES BOOLEAN verbose ""
+{
+}
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..4ef6283
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,10 @@
+# Schedule definitions for thorn IsoSurfacer
+# $Header$
+# schedule IsoSurfacer_Startup at STARTUP before (Webserver_Startup)
+
+if (*isosurfacer != '\0') {
+     schedule IsoSurfacer_Startup at STARTUP
+  {
+    LANG:C
+  } "Startup routine"
+}
diff --git a/src/IsoSurfacer.c b/src/IsoSurfacer.c
new file mode 100644
index 0000000..0ba1c99
--- /dev/null
+++ b/src/IsoSurfacer.c
@@ -0,0 +1,905 @@
+/*@@ 
+   @file      IsoSurfacer.c
+   @date      Wed Oct  8 11:30:45 MET 1997
+   @author    Manuel Panea  <mpd@rzg.mpg.de>
+   @desc 
+   An interface to isosurfacer routines which allows parallel 
+   isosurface finding. The heart of this thorn is @seeroutine
+   FindIsoSurface which takes a GF and returns a vertex* and
+   polygon* for the local (including ownership) chunk of data.
+   The vertices should be based on the x-y-z arrays.<p>
+
+   Note by John Shalf:
+     Currently the lists are collected using a MPI_Gatherv().  I just
+     implemented this for provisional testing.  If the isosurfaces get
+     sufficiently complex, they will flood the master processor with
+     too much data.  In this case we'd need to switch to a more
+     sophisticated method to collect them.  (so the potential problem
+     with the current implementation has been noted and will be dealt
+     with at a later time.)
+   @enddesc 
+ @@*/
+ 
+
+#ifndef TCP
+#define TCP
+#endif
+
+#define Append_Or_EnlargeAndAppend(outs, tmps) \
+q = p; \
+if( (p += strlen(tmps)) >= lfs ) \
+{ \
+    (outs) = ((char *)realloc((void *)(outs), 2 * lfs * sizeof(char))); \
+    lfs *= 2; \
+} \
+strcpy(&(outs)[q], (tmps))
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_DefineThorn.h"
+#include "cctk_Parameters.h"
+
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+
+#include "IsoSurfacerGH.h"
+#include "IsoSurfacerInit.h"
+#include "common.h"
+#include "NuSurfacer.h"
+
+/* Some function forward declarations */
+static void CollectData(cGH *GH, polypatch *perprocessor, polypatch *totals);
+static void WriteBin(cGH *GH, polypatch *totals, par_st *parms,
+                     const char *fullname, CCTK_REAL4 isolevel);
+static void WriteVRML(cGH *GH, polypatch *totals, par_st *parms,
+                      const char *fullname, CCTK_REAL4 isolevel);
+static void WriteASCII(cGH *GH, polypatch *totals, par_st *parms,
+                      const char *fullname, CCTK_REAL4 isolevel);
+static void WriteUCD(cGH *GH, polypatch *totals, par_st *parms,
+                      const char *fullname, CCTK_REAL4 isolevel);
+static void WriteSock(cGH *GH, polypatch *totals, par_st *parms,
+                       const char *fullname, int isoindex, Isolevel&IL);
+static void WriteHDF5(cGH *GH, polypatch *totals, par_st *parms,
+                       const char *fullname, int isoindex, Isolevel&IL);
+
+extern isotimes_st isotimes;
+extern isotimes_st cdtimes;
+
+/***************************************************************************/
+
+/// Ugly global variables
+
+  static polypatch   perprocessor, totals;
+  static GridAll_st  GridAll;
+#ifdef CCTK_MPI
+  static int         print_time=0;
+  static double taccum;
+  double t1;
+#endif
+
+#ifndef true
+#define true 1
+#endif
+#ifndef false
+#define false 0
+#endif
+
+static int GridInit(cGH *GH)
+{
+  static int done = false;
+  int i;
+  if (done) 
+    return true; 
+  
+  done = true; 
+  GridAll.GH = GH;
+  /* get the data pointer to the current timelevel of GF X */
+  i = CCTK_CoordIndex (-1,"x","cart3d");
+  if (i < 0) {
+    CCTK_WARN (1, "Could not get index of 'x' coordinate");
+    return false;
+  }
+  GridAll.timelevel = CCTK_NumTimeLevelsFromVarI (i) - 1;
+  if (GridAll.timelevel > 0)
+      GridAll.timelevel--;
+    if (CCTK_VarTypeI (i) != CCTK_VARIABLE_REAL) {
+      CCTK_WARN (1, "Expected 'x' coordinate to be a REAL variable !");
+      return false;
+    }
+    GridAll.x = (CCTK_REAL *) CCTK_VarDataPtrI (GH, GridAll.timelevel, i);
+
+    /* get the data pointer to the current timelevel of GF Y */
+    i = CCTK_CoordIndex (-1,"y","cart3d");
+    if (i < 0) {
+      CCTK_WARN (1, "Could not get index of 'y' coordinate");
+      return false;
+    }
+    if (CCTK_VarTypeI (i) != CCTK_VARIABLE_REAL) {
+      CCTK_WARN (1, "Expected 'y' coordinate to be a REAL variable !");
+      return false;
+    }
+    GridAll.y = (CCTK_REAL *) CCTK_VarDataPtrI (GH, GridAll.timelevel, i);
+
+    /* get the data pointer to the current timelevel of GF Z */
+    i = CCTK_CoordIndex (-1,"z","cart3d");
+    if (i < 0) {
+      CCTK_WARN (1, "Could not get index of 'z' coordinate");
+      return false;
+    }
+    if (CCTK_VarTypeI (i) != CCTK_VARIABLE_REAL) {
+      CCTK_WARN (1, "Expected 'z' coordinate to be a REAL variable !");
+      return false;
+    }
+    GridAll.z = (CCTK_REAL *) CCTK_VarDataPtrI (GH, GridAll.timelevel, i);
+
+    totals.verts = NULL;
+    totals.polys = NULL;
+    perprocessor.verts = NULL;
+    perprocessor.polys = NULL;
+
+    return true;
+}
+
+static int IsoSurfacerHandleCommands(cGH *GH)
+{
+  isosurfacerGH *myGH = (isosurfacerGH *) GH->extensions [CCTK_GHExtensionHandle ("IsoSurfacer")];
+
+  /* Grab all events that are waiting in queue */
+  /*
+    0)For each command bcast catenated string to all
+    processors
+           
+    Command to change existing isovalue is of format
+    "isosurf",<gridname>[<instance>],<newvalue>
+    1) Select GF
+    2) Scan for tag (if any)
+    If tag does not exist, create one
+    3) change value of existing tag.
+    Command to remove existing surface is of format
+    "isosurf",<gridname>[<instance>],"delete"
+    1) Select GF with scan
+    2) Scan for tag (if any)
+    3) remove existing tag
+    Command to create new isosurface instance
+    "isosurf",<gridname>[<newinstancetag]>],<value>
+    1) Select GF
+    2) Scan for tag (if any)
+    3) If no tag, create new entry
+    Command to create Slice Plane
+    "slice",<gridname>[<newinstancestring>],<coords>
+    Coords are just 2ints + 6 floats
+    2ints are the width and height.  6floats are
+    Just 2 corners of diagonal (assuming ortho)
+    w,h:xyzMin:xyzMax
+    Must be certain to fit in limited space.
+    Byte-normalize the data for performance.
+    1) Currently hardcode for single slice through
+    center plane.  Display as truncated semi
+    transparent surface with texturmap of
+    values (single polygon).  Use mpigatherv
+    to collect the slice.
+    2) Send slice data with command
+    "s:<varname>[<instance>]=<coords>  Range=<min>:<max>"
+    Bcast entry of length 1 if done
+    */
+
+  /*  char catenated[64+64+64] = ""; 
+      if (!http_ReceiveNextCommand(GH, "isosurf", catenated ))
+         break;
+      char *obj = catenated+64,
+           *val = catenated+128;  */
+#if 0
+  HTTPCommand*Cmd = http_ReceiveNextCommand(GH, "isosurf" );
+  
+  if (!Cmd)
+    break;
+  
+  char *obj = Cmd->target,
+    *val = Cmd->value; 
+  
+  char gridname[64]; 
+  int  isoindex = 0; 
+  strcpy(gridname,obj); 
+  
+  /* FUNCTIONALITY: gridfunc[?] means: create a new isosurface */
+  
+  /* gridfunc[4] means: 4th isosurface level  */
+  
+  char   *s = strchr(gridname,'[');
+  if (s){
+    *s='\0'; 
+    s++; 
+    char   *ss=strchr(s,']'); 
+    if (ss){
+      *ss='\0';
+    } 
+    if (*s == '?')
+      isoindex = -1; 
+    else{
+      isoindex=atoi(s); 
+      if (isoindex < 0 )
+	isoindex = 0;
+    }               
+  } 
+  
+  /* Check if the requested gridfunc already has iso parameters */
+  
+  /* param_map::iterator git = myGH->isoparms.find( (const char*)gridname ); */
+  
+  if (git == myGH->isoparms.end() ) { /* not found? --> create new params */
+    /* printf( "No Isoparams for %s\n", gridname ); */ 
+    
+    isoparms_st*params = new isoparms_st(); 
+    params->funcName = strdup( gridname ); 
+    params->formats |= SOCK; 
+    
+    myGH->isoparms[ params->funcName ] = params;
+  } 
+  isoparms_st*params = myGH->isoparms[ gridname ]; 
+  assert( params && params->funcName); 
+  if (isoindex < 0 || (unsigned int) isoindex >= params->isolevels.size() ){
+    params->isolevels.push_back( Isolevel() ); 
+    isoindex = params->isolevels.size() -1 ;
+  } 
+  Isolevel&current = params->isolevels[isoindex]; 
+  current.value = atof(val);  /* Note: Setting of colorinfo not implemented yet.  */
+  delete Cmd;
+#endif 
+  
+  return 0;
+}
+
+static bool doIso(int i, cGH *GH, isosurfacerGH *myGH)
+{
+  /* need to rediscover how to mess with this */
+#if 0
+  char*    fullname = CCTK_FullName (i); 
+  /* see if we are to iterate on this GH.*/
+  param_map::iterator it = myGH->isoparms.find( fullname ); 
+  free (fullname);
+  
+  if (it == myGH->isoparms.end() || !(*it).second)
+    return false; 
+  
+  isoparms_st&myIso = *(*it).second;
+  
+  int firstiter = myIso.firstIteration; 
+  if(GH->cctk_iteration < firstiter)
+    return false; 
+  int itiso = myIso.outfreq; 
+  if (!(itiso && (GH->cctk_iteration - firstiter) % itiso == 0))
+    return false; 
+  return true;
+#endif
+  return false;
+}
+
+
+static void computeIso(int i, cGH *GH, isosurfacerGH *myGH)
+{
+  DECLARE_CCTK_PARAMETERS
+    
+    char*    fullname = CCTK_FullName (i); 
+  
+  isoparms_st*ptrIsofunc = myGH->isoparms[ fullname ]; 
+  
+  if (!ptrIsofunc)
+    return;
+  
+  isoparms_st&isofunc = *ptrIsofunc;
+  
+  /* if the do_iso flag is set */
+  /* then draw an iso for each level of surfaces requested */
+  
+  for (unsigned int j = 0; j < isofunc.isolevels.size(); j++) { 
+    /* Do the bcast for nisolevels here */
+    isotimes.ncalls++;
+#ifndef	_WIN32
+    isotimes.realtime_begin = times(&isotimes.tms_begin);
+#endif
+    
+    GridAll.gfi = i;
+    
+    /* get the current timelevel */
+    GridAll.timelevel = CCTK_NumTimeLevelsFromVarI (i) - 1;
+    if (GridAll.timelevel > 0)
+      GridAll.timelevel--;
+    
+    
+    // is it really required to have this static here??
+    static par_st      parms;
+    
+    parms.outfname              = NULL;
+    parms.outtype              = isofunc.formats; 
+    parms.isovalue              = (CCTK_REAL4)isofunc.isolevels[j].value;
+    parms.step                      = isofunc.resolution;
+    parms.doRemoveDuplicateVertices = isofunc.uniq_verts;
+    parms.useTree = parms.doRemoveDuplicateVertices;
+    parms.doEliminateSmallTriangles = 0;
+    
+#ifdef CCTK_MPI
+    t1=MPI_Wtime();
+#endif
+    if(nusurfacer)           /* use the new isosurfacer */
+      NuFindSurface(GH,i,parms.isovalue,&perprocessor);
+    else                     /* use the old isosurfacer */
+      FindIsoSurface(GH, (void *)&GridAll, &perprocessor, &parms);
+#ifdef CCTK_MPI
+    taccum += (MPI_Wtime()-t1);
+    if(print_time)
+      fprintf(stderr,"       | %s[p=%u] Accumtime=%f\n",
+	      nusurfacer ? "Newsurf" : "Oldsurf", CCTK_MyProc(GH),taccum);
+#endif
+    
+    CollectData(GH, &perprocessor, &totals);
+    
+    if (CCTK_MyProc (GH) == 0){
+      CCTK_INT4 tmppolys[3] = {0,0,0};
+      CCTK_REAL4 tmpverts[3] = {0.0,0.0,0.0};
+      CCTK_INT4 *polybackup=0;
+      CCTK_REAL4 *vertbackup=0;
+      
+      if(allow_empty_sends && totals.nverts<=0){ /* lets kludge it! */
+	totals.nverts=1; /* one vertex */
+	totals.npolys=1; /* one polygon which refers to the same vertex 3 times */
+	polybackup=totals.polys; totals.polys = tmppolys;
+	vertbackup=totals.verts; totals.verts = tmpverts;
+      }
+      
+      if (verbose)
+	printf("       | IsoSurfacer: GF=%s, value %f,",
+	       fullname, parms.isovalue);
+      
+      if(totals.nverts > 0) {
+	
+	if (verbose)
+	  printf("%d vertices, %d triangles\n", totals.nverts, totals.npolys);
+	
+	if(parms.outtype & BIN)
+	  WriteBin(GH, &totals, &parms, fullname, isofunc.isolevels[j].value);
+	
+	if(parms.outtype & ASCII)
+	  WriteASCII(GH, &totals, &parms, fullname,
+		     isofunc.isolevels[j].value);
+	
+	if(parms.outtype & UCD)
+	  WriteUCD(GH, &totals, &parms, fullname, isofunc.isolevels[j].value);
+	
+	if(parms.outtype & VRML)
+	  WriteVRML(GH, &totals, &parms, fullname,isofunc.isolevels[j].value);
+	
+	if(parms.outtype & SOCK)
+	  WriteSock(GH, &totals,&parms,fullname,j, isofunc.isolevels[j]);
+	
+	if(parms.outtype & ISOHDF5)
+	  WriteHDF5(GH, &totals,&parms,fullname,j, isofunc.isolevels[j]);
+	
+      } else
+	printf(" no isosurface\n");
+      
+      
+      if(polybackup || vertbackup){ /* copy back datastructures */
+	totals.polys=polybackup;
+	totals.verts=vertbackup;
+      }
+    }
+#ifndef _WIN32
+    isotimes.realtime_end = times(&isotimes.tms_end);
+    isotimes.realtime_total +=
+      (isotimes.realtime_end - isotimes.realtime_begin);
+    isotimes.tms_total.tms_utime +=
+      (isotimes.tms_end.tms_utime - isotimes.tms_begin.tms_utime);
+    isotimes.tms_total.tms_stime +=
+      (isotimes.tms_end.tms_stime - isotimes.tms_begin.tms_stime);
+#endif
+  } 
+  
+  free (fullname);
+}
+
+int IsoSurfacer(cGH *GH){
+  isosurfacerGH *myGH = (isosurfacerGH *) GH->extensions [CCTK_GHExtensionHandle ("IsoSurfacer")];
+  if (!GridInit(GH))
+    return false; 
+  
+  IsoSurfacerHandleCommands(GH);
+  
+  /* do a check for new isosurfaces */
+  /* Perhaps do a bcast for "changed" flags. 
+     which are embedded in each iso. */
+  for (int i = 0; i<CCTK_NumVars (); i++) { /* for each Grid  */
+    if (!doIso(i, GH, myGH))
+      continue; 
+    
+    if (CCTK_QueryGroupStorageI(GH,CCTK_GroupIndexFromVarI(i))) {
+      computeIso(i, GH, myGH);
+    }
+    else {
+      CCTK_VWarn (9, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "No storage assigned for '%s'", CCTK_VarName (i));
+    }    
+  }
+  return true;
+}
+
+/***************************************************************************/
+
+void CollectData(cGH *GH, polypatch *perprocessor, polypatch *totals) {
+  DECLARE_CCTK_PARAMETERS
+#ifdef CCTK_MPI
+    
+    /* Collect the isosurface data onto processor 0 */
+    /* First collect the numbers of vertices (MPI_Gather)*/
+    /* Then the vertices themselves (MPI_Gatherv)*/
+    /* Then the number of polygons (MPI_Gather)*/
+    /* Then the polygon list. This needs to be offset by the number of
+       vertices ahead of me, so do that after collection. */
+    
+    static  int   a, b;
+  static  int  *vcount = NULL;
+  static  int  *pcount = NULL;
+  static  int  *vdispl = NULL;
+  static  int  *pdispl = NULL;
+  static  int  lastNverts = -1;
+  static  int  lastNpolys = -1;
+  static  CCTK_INT4  my_vpcount [2];
+  static  CCTK_INT4 *vpcount = NULL;
+  pGH     *pughGH;
+  
+  
+  pughGH = PUGH_pGH (GH);
+  
+  cdtimes.ncalls++;
+#ifndef _WIN32
+  cdtimes.realtime_begin = times(&cdtimes.tms_begin);
+#endif
+  
+  /* Allocate temporary arrays for the MPI_Gatherv() operation */
+  if(pughGH->myproc == 0) {
+    if(vpcount == NULL) {
+      vpcount = NEW (2 * pughGH->nprocs, CCTK_INT4);
+      vcount = NEW(4*pughGH->nprocs, int);
+      pcount = &vcount[pughGH->nprocs];
+      vdispl = &vcount[2*pughGH->nprocs];
+      pdispl = &vcount[3*pughGH->nprocs];
+    }
+  }
+  
+  /* **** Gather counts of vertex and polygon lists in one wash */
+  my_vpcount [0] = 3*perprocessor->nverts;
+  my_vpcount [1] = 3*perprocessor->npolys;
+  CACTUS_MPI_ERROR (MPI_Gather (my_vpcount, 2, PUGH_MPI_INT4, vpcount, 2,
+				PUGH_MPI_INT4, 0, pughGH->PUGH_COMM_WORLD));
+  
+  /* Allocate for the vertices */
+  if(pughGH->myproc == 0){
+    totals->nverts = 0;
+    totals->npolys = 0;
+    
+    for (a = 0; a < pughGH->nprocs; a++) {
+      /* now sort out the vcounts and pcounts into the int-arrays
+         used in MPI_Gatherv() */
+      vcount [a] = vpcount [2*a + 0];
+      pcount [a] = vpcount [2*a + 1];
+      
+      /* compute the displacements and total number of elements */
+      vdispl[a] = totals->nverts;
+      pdispl[a] = totals->npolys;
+      totals->nverts += vcount[a];
+      totals->npolys += pcount[a];
+    }
+    
+    totals->nverts /= 3;
+    totals->npolys /= 3;
+    
+    
+    /*
+      Kludge by M. Panea:
+      Allocate 3 more than necessary so as to divert the output of processors
+      with nverts=0 to a dummy area in the upper part of totals->verts.
+      Otherwise the MPI_Gatherv routine gives lots of problems.
+    */
+    if(lastNverts < totals->nverts){
+      REALLOC(totals->verts, 3+3*totals->nverts, CCTK_REAL4);
+      lastNverts = totals->nverts;
+    }
+    if(lastNpolys < totals->npolys){
+      REALLOC(totals->polys, 3+3*totals->npolys, CCTK_INT4);
+      lastNpolys = totals->npolys;
+    }
+  } /* end processor 0 setup */
+  
+  /* Gather the vertex lists from all processors */
+  CACTUS_MPI_ERROR (MPI_Gatherv (perprocessor->verts, perprocessor->nverts*3,
+				 PUGH_MPI_REAL4, totals->verts, vcount, vdispl, PUGH_MPI_REAL4,
+				 0, pughGH->PUGH_COMM_WORLD));
+  
+  /* Gather the polygon lists from all processors */
+  CACTUS_MPI_ERROR (MPI_Gatherv (perprocessor->polys, perprocessor->npolys*3,
+				 PUGH_MPI_INT4, totals->polys, pcount, pdispl, PUGH_MPI_INT4,
+				 0, pughGH->PUGH_COMM_WORLD));
+  /* JMS Comment:
+     Now we need to re-number the connectivity list so that
+     they can be combined in the master list (remember, the
+     connectivity numbering right now is processor local...
+     we need to combine this connectivity info into a single list */
+  if(pughGH->myproc == 0){ /* If I'm the master processor */
+    for(a = 1; a<pughGH->nprocs; a++) { /* for segments gathered from each processor */
+      for(b = 0; b<pcount[a]; b++){ /* for all polys in the segment */
+	/* Add displacement to each connectivity integer */
+	totals->polys[pdispl[a]+b] += vdispl[a]/3;
+      }
+    }
+  }
+  
+#ifndef _WIN32
+  cdtimes.realtime_end = times(&cdtimes.tms_end);
+  cdtimes.realtime_total += (cdtimes.realtime_end - cdtimes.realtime_begin);
+  
+  cdtimes.tms_total.tms_utime +=
+    (cdtimes.tms_end.tms_utime - cdtimes.tms_begin.tms_utime);
+  cdtimes.tms_total.tms_stime +=
+    (cdtimes.tms_end.tms_stime - cdtimes.tms_begin.tms_stime);
+#endif
+  
+#else
+  *totals = *perprocessor;
+#endif /* MPI */
+}
+
+/***************************************************************************/
+
+void WriteSock(cGH *GH, polypatch *totals, par_st *parms,
+	       const char *fullname, int isoindex,
+	       Isolevel&IsoL)
+{
+  // Was: parms->timestep, which was callnumber from above. 
+  sprintf(tmpstring," :%s[%u]=%f,%u range=%f:%f",
+	  full_id,isoindex,parms->isovalue,
+	  GH->cctk_iteration,
+	  parms->minval,parms->maxval);
+  
+  /********Now write vertices**********/
+  *tmpstring='v';
+  verts = totals->nverts > 0 ? totals->verts : tmpvert;
+  DS->sendData(tmpstring,3*totals->nverts,DataType::int2type(REMOTE_IO_FLOAT32),verts);
+  
+  /********Now write polygons**********/
+  *tmpstring='c';
+  polys = totals->npolys > 0 ? totals->polys : tmppoly;
+  DS->sendData(tmpstring,3*totals->npolys,DataType::int2type(REMOTE_IO_INT32),polys);
+}
+
+/***************************************************************************/
+
+static void WriteHDF5(cGH *GH, polypatch *totals, par_st *parms,
+                      const char *fullname, int isoindex, Isolevel&IsoL)
+{
+  DECLARE_CCTK_PARAMETERS
+#ifdef CACTUSPUGHIO_IOHDF5
+  char *filename;
+
+
+  filename = (char *) malloc (strlen (outdir) + strlen (fullname) +
+                              strlen (PUGH_pGH (GH)->identity_string) + 20);
+  sprintf (filename, "%s/%s%s.iso.h5", outdir, fullname,
+           PUGH_pGH (GH)->identity_string);
+  IOHDF5_WriteIsosurface (GH, filename, fullname,
+                          (CCTK_INT) GH->cctk_iteration,
+                          (CCTK_INT) GridAll.timelevel,
+                          (CCTK_REAL) parms->isovalue,
+                          (CCTK_REAL) parms->minval, (CCTK_REAL) parms->maxval,
+                          (int) totals->npolys, totals->polys,
+                          (int) totals->nverts, totals->verts);
+  free (filename);
+
+#else
+
+  CCTK_WARN (1, "No HDF5 isosurface output because CactusPUGHIO/IOHDF5 was not compiled in !");
+
+#endif
+}
+
+/***************************************************************************/
+
+void
+WriteBin(cGH *GH, polypatch *totals, par_st *parms, const char *fullname,
+         CCTK_REAL4 isolevel)
+{
+  /* Write Isosurfaces in Manuel's own format which can
+     then be read by 'isoshow' and 'isoserver' */
+
+  DECLARE_CCTK_PARAMETERS
+  static FILE  *bf;
+  static char   fname[256];
+  static char   tmpstring[100];
+  static CCTK_INT4   tmp;
+  char formatstring[256];
+
+/* Debug printing stuff ---------------------
+  sprintf(fname," :%s[%u]=%f,%u range=%f:%f",
+           fullname,j,parms->isovalue,parms->timestep,
+           parms->minval,parms->maxval);
+  printf("SockString = [%s]\n",fname);
+--------------------------------------------*/
+  if(*format_str=='%')
+    sprintf (formatstring, "%s/%%s%%s_%s_%%d.iso.bin", outdir, format_str);
+  else
+    sprintf (formatstring, "%s/%%s%%s_%%3.3f_%%d.iso.bin", outdir);
+  sprintf (fname, formatstring, fullname, PUGH_pGH (GH)->identity_string,
+           isolevel, GH->cctk_iteration);
+
+  if(!(bf =fopen(fname,"w")))
+    {
+      printf("Attempt to open filename [%s] failed\n",fname);
+      perror("Could not create isosurface bin output file");
+      return;
+    }
+  strcpy(tmpstring, "ISOSURFACE");
+
+  fwrite(tmpstring, 1, strlen(tmpstring)+1, bf);
+  fwrite(fullname, 1, strlen(fullname)+1, bf);
+
+  tmp = (CCTK_INT4)GH->cctk_iteration;
+  fwrite(&tmp, sizeof(CCTK_INT4), 1, bf);
+
+  fwrite(&parms->isovalue, sizeof(CCTK_REAL4), 1, bf);
+
+  fwrite(&parms->step, sizeof(CCTK_INT4), 1, bf);
+  fwrite(&parms->doRemoveDuplicateVertices, sizeof(CCTK_INT4), 1, bf);
+
+  fwrite(&totals->nverts, sizeof(CCTK_INT4), 1, bf);
+  fwrite(&totals->npolys, sizeof(CCTK_INT4), 1, bf);
+
+  fwrite(totals->verts, sizeof(CCTK_REAL4), 3*totals->nverts, bf);
+  if(totals->npolys > 0)
+    fwrite(totals->polys, sizeof(CCTK_INT4), 3*totals->npolys, bf);
+
+  fclose(bf);
+}
+
+/***************************************************************************/
+
+void
+WriteVRML(cGH *GH, polypatch *totals, par_st *parms, const char *fullname,
+          CCTK_REAL4 isolevel)
+{
+  /* Write isosurfaces in VRML 1.0 ASCII format 
+     (can be read by ivview) */
+
+  DECLARE_CCTK_PARAMETERS
+  static int   a, p, q;
+  static FILE *pd;
+  static char  fname[256];
+  static char *outs = NULL;
+  static char  tmps[256];
+  static CCTK_INT4  efs = 0, lfs = 0;
+
+  sprintf(fname, "%s/%s%s_%3.3f_%d.iso.wrl",
+          outdir, fullname, PUGH_pGH (GH)->identity_string, isolevel,
+          GH->cctk_iteration);
+  if (!(pd = fopen(fname, "w"))) {
+      perror("Could not create isosurface vrml output file");
+      return;
+  }
+
+  /* In order to speed up the writing, we put everything in a big string
+   * and write the whole string out in one go.
+   */
+  efs = 30*totals->nverts + 30*totals->npolys; /* Estimated string size */
+  p = 0;
+  if(lfs < efs)
+    {
+      outs = (char *)realloc((void *)outs, efs * sizeof(char));
+      lfs = efs;
+    }
+  outs[0] = '\0';
+
+  sprintf (tmps, "#VRML V1.0 ascii\nSeparator {\n\tPointLight { location 1 1 1 }\n\t  Coordinate3 {\n\t\tpoint [ ");
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+  sprintf (tmps, "%f %f %f",
+            totals->verts[0], totals->verts[1], totals->verts[2]);
+  Append_Or_EnlargeAndAppend(outs, tmps);
+  for (a = 3; a<3*totals->nverts; a+= 3)
+    {
+      sprintf (tmps, ", \n\t\t%f %f %f",
+                totals->verts[a+0], totals->verts[a+1], totals->verts[a+2]);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+  sprintf(tmps, "]\n\t}\n\tIndexedFaceSet {\n\t\tcoordIndex [ ");
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+  sprintf (tmps, "%d, %d, %d, -1", 
+            totals->polys[0], totals->polys[1], totals->polys[2]);
+  Append_Or_EnlargeAndAppend(outs, tmps);
+  for (a = 3; a<3*totals->npolys; a += 3)
+    {
+      sprintf (tmps, ", \n\t\t%d, %d, %d, -1", 
+                totals->polys[a], totals->polys[a+1], totals->polys[a+2]);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+  sprintf(tmps, "\t\t]\n\t}\n}");
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+  fwrite(outs, 1, strlen(outs), pd);
+  fclose(pd);
+}
+
+/***************************************************************************/
+
+void
+WriteASCII(cGH *GH, polypatch *totals, par_st *parms, const char *fullname,
+           CCTK_REAL4 isolevel)
+{
+  /* Write Iso's in Paul's ASCII format */
+
+  DECLARE_CCTK_PARAMETERS
+  static int   a, p, q;
+  static FILE *pd;
+  static char  fname[256];
+  static char *outs = NULL;
+  static char  tmps[256];
+  static CCTK_INT4  efs = 0, lfs = 0;
+
+  sprintf(fname, "%s/%s%s_%3.3f_%d.iso.ascii",
+          outdir, fullname, PUGH_pGH (GH)->identity_string, isolevel,
+          GH->cctk_iteration);
+
+  if (!(pd = fopen(fname, "w"))){
+      perror("Could not create isosurface ascii output file");
+      return;
+    }
+
+  /* In order to speed up the writing, we put everything in a big string
+   * and write the whole string out in one go.
+   */
+  efs = 30*totals->nverts + 20*totals->npolys; /* Estimated string size */
+  p = 0;
+  if(lfs < efs)
+    {
+      outs = (char *)realloc((void *)outs, efs * sizeof(char));
+      lfs = efs;
+    }
+  outs[0] = '\0';
+
+  sprintf (tmps, "%d %d\n", totals->nverts, totals->npolys);
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+  for (a = 0; a<3*totals->nverts; a += 3)
+    {
+      sprintf (tmps, "%f %f %f\n",
+                totals->verts[a+0], totals->verts[a+1], totals->verts[a+2]);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+
+  for (a = 0; a<3*totals->npolys; a += 3)
+    {
+      sprintf (tmps, "%d %d %d\n", 
+                totals->polys[a], totals->polys[a+1], totals->polys[a+2]);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+
+  fwrite(outs, 1, strlen(outs), pd);
+  fclose(pd);
+}
+
+/***************************************************************************/
+
+void
+WriteUCD(cGH *GH, polypatch *totals, par_st *parms, const char *fullname,
+         CCTK_REAL4 isolevel)
+{
+  /* Write Iso's in AVS ascii UCD format */
+
+  DECLARE_CCTK_PARAMETERS
+  static int   a, p, q;
+  static FILE *uf;
+  static char  fname[256];
+  static char *outs = NULL;
+  static char  tmps[256];
+  static CCTK_INT4  efs = 0, lfs = 0;
+
+  sprintf(fname, "%s/%s%s_%3.3f_%d.iso.inp", 
+           outdir, fullname, PUGH_pGH (GH)->identity_string, isolevel,
+           GH->cctk_iteration);
+
+  if (!(uf = fopen(fname, "w"))) {
+      perror("Could not create isosurface ucd output file");
+      return;
+  }
+
+  /* In order to speed up the writing, we put everything in a big string
+   * and write the whole string out in one go.
+   */
+  efs = 60*totals->nverts + 40*totals->npolys; /* Estimated string size */
+  p = 0;
+  if(lfs < efs)
+    {
+      outs = (char *)realloc((void *)outs, efs * sizeof(char));
+      lfs = efs;
+    }
+  outs[0] = '\0';
+
+  sprintf(tmps, "%d %d 1 0 0\n", totals->nverts, totals->npolys);
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+
+  for (a = 0; a<3*totals->nverts; a += 3)
+    {
+      sprintf (tmps, "%d %f %f %f\n", a/3+1,
+                        totals->verts[a+0], totals->verts[a+1], totals->verts[a+2]);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+
+  for (a = 0; a<3*totals->npolys; a += 3)
+    {
+      sprintf (tmps, "%d %d tri %d %d %d\n", a/3+1, a/3+1,
+                totals->polys[a]+1, totals->polys[a+1]+1, totals->polys[a+2]+1);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+
+  sprintf(tmps, "1 1\nisosurface,in\n1 1.20000\n");
+  Append_Or_EnlargeAndAppend(outs, tmps);
+
+  for(a=1; a<totals->nverts; a++)
+    {
+      sprintf(tmps, "%d 1.00000\n", a+1);
+      Append_Or_EnlargeAndAppend(outs, tmps);
+    }
+
+  fwrite(outs, 1, strlen(outs), uf);
+  fclose(uf);
+}
+
+/***************************************************************************/
+
+
+int IsoSurfacer_TimeForOutput(cGH *GH, int i)
+{
+  /* Get the GH extensions for IOUtil and IOASCII */
+  isosurfacerGH *myGH = (isosurfacerGH *)GH->extensions[CCTK_GHExtensionHandle("IsoSurfacer")];
+
+  if (!myGH)
+    return false;
+  
+  /* do a check for new isosurfaces */
+  /* Perhaps do a bcast for "changed" flags. 
+     which are embedded in each iso. */
+  
+  return doIso(i, GH, myGH); 
+
+/*
+
+  if (! myGH->do_iso[i])
+    return 0;
+  
+int  firstiter = myGH->isoparms[i].firstIteration;
+  if(GH->cctk_iteration < firstiter)
+    return 0;
+  
+int  itiso = myGH->isoparms[i].outfreq;
+  if (!(itiso && (GH->cctk_iteration - firstiter) % itiso == 0))
+    return 0;
+  
+  return 1;
+*/
+
+}
+
+int IsoSurfacer_TriggerOutput(cGH *GH,int variable)
+{
+  isosurfacerGH *myGH;
+
+  myGH = (isosurfacerGH *) GH->extensions [CCTK_GHExtensionHandle ("IsoSurfacer")];
+  
+  computeIso(variable,GH,myGH);
+ 
+  /* FIXME MAYBE
+  myGH->out1D_last[variable] = GH->cctk_iteration;
+  */
+  return 0;
+}
diff --git a/src/IsoSurfacerGH.h b/src/IsoSurfacerGH.h
new file mode 100644
index 0000000..b5ef1cc
--- /dev/null
+++ b/src/IsoSurfacerGH.h
@@ -0,0 +1,63 @@
+//   The extensions to the GH structure from IsoSurfacer.
+#ifndef __IsoSurface_GH_h_
+#define __IsoSurface_GH_h_
+struct Isolevel
+{
+  float value;  
+  int	Ncolorinfo;
+  union
+  {
+    struct
+    {
+      CCTK_REAL4 transparency,
+	ForeGround_red, ForeGround_green, ForeGround_blue,
+	BackGround_red, BackGround_green, BackGround_blue;
+    }
+    colors;
+    
+    CCTK_REAL4 colorinfo[7];
+  };
+  
+  Isolevel()
+    : value(0.0)
+    , Ncolorinfo(0)
+  {
+    for(int i=0; i<7; i++)
+      colorinfo[i] = 0.7;
+  }
+};
+
+#define BIN      1
+#define ASCII    2
+#define UCD      4
+#define VRML     8
+#define SOCK     16
+#define ISOHDF5  32
+#define NONE     64
+
+typedef struct isoparms_st
+{
+  char       *funcName;
+  vector<Isolevel> isolevels;
+  short int   formats;
+  short int   outfreq;
+  short int   firstIteration;
+  short int   resolution;
+  short int   uniq_verts;
+  
+  isoparms_st()
+    : funcName(0)
+    , formats(0)
+    , outfreq(1)
+    , firstIteration(1)
+    , resolution(1)
+    , uniq_verts(1)
+  {}
+} isoparms_st;
+
+struct isosurfacerGH
+{
+  isoparms_st isoparms;
+};
+
+#endif // __IsoSurface_GH_h_
diff --git a/src/IsoSurfacerInit.c b/src/IsoSurfacerInit.c
new file mode 100644
index 0000000..5948426
--- /dev/null
+++ b/src/IsoSurfacerInit.c
@@ -0,0 +1,494 @@
+#ifndef TCP /* JMS addition */
+#define TCP
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#ifndef	_WIN32
+#include <time.h>
+#include <unistd.h>
+#endif
+
+#include <ctype.h>
+
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+//#include <cctk_DefineThorn.h>
+
+//#ifdef	CACTUSNET_HTTP
+//#include "CactusNet/http/src/http_GHExtensions.hpp"
+//#endif
+
+#include "http_GHExtensions.hpp"
+
+
+#include "IsoSurfacerGH.h"
+#include "IsoSurfacerInit.h"
+
+
+/* Some function forward declarations */
+CCTK_INT4   ParseIsoString(char *s, param_map&parms);
+CCTK_INT4   NumberOfItems(char *s, char *d);
+char  *GetSubString(char *s, char *d, CCTK_INT4 n);
+
+isotimes_st isotimes;
+isotimes_st cdtimes;
+
+CCTK_INT4   RunIsoSurfacer;
+
+
+#ifdef	CACTUSNET_HTTP
+CactusCommands::Client * IsosurfaceHttpClientCreator(CactusCommands::Reply*msgbuf, CactusCommands::CommandList&cmds );
+#endif
+
+/***************************************************************/
+void *IsoSurfacer_SetupGH (tFleshConfig *config, 
+			   int convergence_level, 
+			   cGH *GH){
+  isosurfacerGH *myGH=(isosurfacerGH*)malloc(sizeof(isosurfacerGH));
+  return myGH;
+}
+
+int IsoSurfacer_InitGH (cGH *GH){
+  DECLARE_CCTK_PARAMETERS
+    int n;
+  isosurfacerGH *myGH;
+  myGH = (isosurfacerGH *) GH->extensions [CCTK_GHExtensionHandle ("IsoSurfacer")];
+  
+  printf("IsoInit\n");
+  /* initialize values */
+  myGH->funcName=0;
+  myGH->formats=0;
+  myGH->outfreq=1;
+  myGH->firstIteration=1;
+  myGH->resolution=1;
+  myGH->uniq_verts=1;
+  myGH->RunIsoSurfacer = 0;
+  
+  /* OK, this is ridiculous.
+     Separate out the isostring components into 
+     separate things and end this craziness */
+  if( (n = ParseIsoString(isosurfacer, myGH->isoparms)) == 0)
+    return 0;
+  
+  RunIsoSurfacer = 1; 
+  /*
+  doIso = 0;
+
+  for (i=0; i<CCTK_NumVars (); i++) {
+    char *fullname = CCTK_FullName (i);
+
+    for(j=0; j<n; j++)
+      if(CCTK_Equals (fullname, allparms[j].funcName) ) {
+        doIso = 1;
+        myGH->do_iso [i] = 1;
+        myGH->isoparms [i] = allparms[j];
+        if((allparms[j]).formats&SOCK){ 
+          printf("usesockets**********\n");
+        }
+        break;
+      }
+    free (fullname);
+  }
+*/ 
+
+  if (CCTK_MyProc (GH) == 0){
+    char *cmd = (char *) malloc (strlen (outdir) + 80); 
+    
+    sprintf (cmd, "mkdir -p %s", outdir); 
+    if (system (cmd) < 0)
+      CCTK_WARN (1, "Problem creating IsoSurfacer output directory"); 
+    free (cmd);
+  }
+  memset(&isotimes, 0, sizeof(isotimes)); 
+  memset(&cdtimes, 0, sizeof(cdtimes));
+  return true;
+}
+
+/************************************************************/
+
+void IsoSurfaceEnd(cGH *GH)
+{
+
+
+  if( RunIsoSurfacer == 0 )
+    return;
+
+#if 0
+  isosurfacerGH *myGH;
+
+  myGH = (isosurfacerGH *) GH->extensions [CCTK_GHExtensionHandle ("IsoSurfacer")];
+
+  /* Close Data Connections (If open) */
+  if(myGH->datasocket) {
+         CCTK_INT4 i=0;
+         SendData(myGH->datasocket,"done",1,INT32,&i);
+         delete static_cast<DataSender*>(myGH->datasocket);
+  }
+  if(myGH->commandsocket) 
+	delete static_cast<CommandReceiver*>(myGH->commandsocket);
+
+#endif	// EXTERNAL_REMOTEIO
+
+  printf("IsoSurfacer timing \n");
+  printf("\n Total number of calls to the IsoSurfacer: %d\n", isotimes.ncalls);
+  printf(" Time spent in IsoSurfacer:\n");
+  printf("  Total   real   time          : %.4f s\n",
+          isotimes.realtime_total * 1.0 / CLK_TCK);
+#ifndef	_WIN32
+  printf("  Total   user   time          : %.4f s\n",
+          isotimes.tms_total.tms_utime * 1.0 / CLK_TCK);
+  printf("  Total   system time          : %.4f s\n",
+          isotimes.tms_total.tms_stime * 1.0 / CLK_TCK);
+  printf("  Average real   time per call : %.4f s\n",
+          isotimes.realtime_total * 1.0 / CLK_TCK / isotimes.ncalls);
+  printf("  Average user   time per call : %.4f s\n",
+          isotimes.tms_total.tms_utime * 1.0 / CLK_TCK / isotimes.ncalls);
+  printf("  Average system time per call : %.4f s\n",
+          isotimes.tms_total.tms_stime * 1.0 / CLK_TCK / isotimes.ncalls);
+#endif
+  printf("\n Total number of calls to CollectData: %d\n", cdtimes.ncalls);
+  printf(" Time spent in CollectData:\n");
+  printf("  Total   real   time          : %.4f s\n",
+          cdtimes.realtime_total * 1.0 / CLK_TCK);
+#ifndef _WIN32
+  printf("  Total   user   time          : %.4f s\n",
+          cdtimes.tms_total.tms_utime * 1.0 / CLK_TCK);
+  printf("  Total   system time          : %.4f s\n",
+          cdtimes.tms_total.tms_stime * 1.0 / CLK_TCK);
+#endif
+  printf("  Average real   time per call : %.4f s\n",
+          cdtimes.realtime_total * 1.0 / CLK_TCK / cdtimes.ncalls);
+#ifndef _WIN32
+  printf("  Average user   time per call : %.4f s\n",
+          cdtimes.tms_total.tms_utime * 1.0 / CLK_TCK / cdtimes.ncalls);
+  printf("  Average system time per call : %.4f s\n",
+          cdtimes.tms_total.tms_stime * 1.0 / CLK_TCK / cdtimes.ncalls);
+  printf("--------------------------------------------------\n");
+#endif
+}
+
+int IsoSurfacer_rfrTraverseGH (cGH *GH, int rfrpoint)
+{
+  return 0;
+}
+
+/*****************************************************************************/
+
+CCTK_INT4 ParseIsoString(char *s, param_map&allparms)
+{
+  CCTK_INT4         i, j;
+  CCTK_INT4         n, m;
+  CCTK_INT4         nlevels, nformats;
+  char        *sub0, *sub1, *sub2;
+
+
+
+  /* First see how many functions we have      */
+  /* by simply counting the number of {} pairs */
+  n = NumberOfItems(s, "{}");
+
+  if(n == 0)
+  {
+    printf("Isosurfacer: empty parameter string '%s'\n", s);
+    return 0;
+  }
+
+  /* Allocate that many isoparms structs */ 
+// *allparms = NEW(n, isoparms_st);
+
+  /* Now fill in each isoparms struct with the info for each function */
+  for(i=0; i<n; i++)
+  {
+    sub0 = GetSubString(s, "{}", i);
+
+    /* Check that we have the right number of () pairs */
+    m = NumberOfItems(sub0, "()");
+
+    if(m != 4)
+    {
+      printf("Isosurfacer: improperly formed parameter string '%s'\n", s);
+      printf("             there should be exactly 4 () pairs\n");
+      free(sub0);
+      return 0;
+    }
+
+    /*- 1 -----------------------------------------------------------------*/
+    /* The first () pair contains the function name */ 
+    char*gridfunc = GetSubString(sub0, "()", 0); 
+
+    isoparms_st*gridparams = new isoparms_st();
+
+	        allparms[gridfunc] = gridparams ;
+
+    	gridparams->funcName = gridfunc;
+
+
+    /*- 2 -----------------------------------------------------------------*/
+    /* The second () pair contains the list of isosurface levels         */
+    sub1 = GetSubString(sub0, "()", 1);
+
+      /* First see how many there are by counting the number of commas + 1 */
+      nlevels = NumberOfItems(sub1, ",");
+
+      /* and get the values */
+      for(j=0; j<nlevels; j++)
+      {      
+         sub2 = GetSubString(sub1, ",", j);
+         gridparams->isolevels.push_back( Isolevel() );
+
+         gridparams->isolevels.back().value = atof(sub2); 
+
+	 gridparams->isolevels.back().Ncolorinfo = 0;
+
+      char*N = strchr(sub2, ':');
+	 if (N)
+	 {
+		N++; 
+	 int	Nvalues = NumberOfItems(N, "/");
+		if (Nvalues > 7 )
+			Nvalues = 7; 
+
+		for(int k=0; k<Nvalues; k++)
+		{
+		char*sub3 = GetSubString(N, "/", k); 
+		char*number = sub3;
+			while(*number && !isdigit(*number)) number++;
+			gridparams->isolevels.back().colorinfo[k] = atof(number);
+			free(sub3);
+		} 
+
+		gridparams->isolevels.back().Ncolorinfo = Nvalues;
+	 }
+
+         free(sub2);
+      }
+
+      free(sub1);
+
+
+    /*- 3 -----------------------------------------------------------------*/
+    /* The third () pair contains the list of output formats */
+    sub1 = GetSubString(sub0, "()", 2);
+
+      /* See how many there are by counting the number of commas + 1 */
+      nformats = NumberOfItems(sub1, ",");
+
+      /* and get the values */
+      gridparams->formats = 0;
+      for(j=0; j<nformats; j++)
+      {
+         sub2 = GetSubString(sub1, ",", j);
+         if( !strcmp(sub2, "BIN") )
+           gridparams->formats |= BIN;
+         else if (!strcmp(sub2, "VRML"))
+           gridparams->formats |= VRML;
+         else if (!strcmp(sub2, "ASCII"))
+           gridparams->formats |= ASCII;
+         else if (!strcmp(sub2, "UCD")){
+#ifdef VERBOSE
+           printf("UCD found in formats string\n");
+#endif
+           gridparams->formats |= UCD;
+        }
+         else if (!strcmp(sub2,"SOCK")){         /* JMS addition */
+#ifdef VERBOSE
+           printf("SOCK found in formats string for %d\n",i);
+#endif
+           gridparams->formats |= SOCK;
+         }
+         else if (!strcmp(sub2,"ISOHDF5")){         /* TR addition */
+#ifdef VERBOSE
+           printf("HDF5 found in formats string for %d\n",i);
+#endif
+           gridparams->formats |= ISOHDF5;
+         }
+         else if(!strcmp(sub2,"NONE")){
+#ifdef VERBOSE
+           printf("NONE string found in formats string\n");
+#endif
+           gridparams->formats |= NONE;
+         }
+         else
+         {
+           printf("Isosurfacer: invalid output format '%s'\n", sub2);
+           printf("             supported formats are: BIN, VRML, ASCII, UCD, SOCK\n");
+           free(sub2);
+           free(sub1);
+           free(sub0);
+           return 0;
+         }
+         free(sub2);
+      }
+
+      free(sub1);
+
+
+    /*- 4 -----------------------------------------------------------------*/
+    /* The fourth () pair contains the other various parameters */
+    sub1 = GetSubString(sub0, "()", 3);
+
+      /* Check that the amount of them is correct (=4) */
+      if( NumberOfItems(sub1, ",") != 4)
+      {
+         printf("Isosurfacer: improperly formed parameter string '%s'\n", s);
+         printf("             there should be exactly 4 parameters in the last () pair\n");
+         free(sub1);
+         free(sub0);
+         return 0;
+      }
+
+      /* The first one is the frequency of output */
+      sub2 = GetSubString(sub1, ",", 0);
+      gridparams->outfreq = (CCTK_INT4)atoi(sub2);
+      free(sub2);
+
+      /* Then the first iteration */
+      sub2 = GetSubString(sub1, ",", 1);
+      gridparams->firstIteration = (CCTK_INT4)atoi(sub2);
+      free(sub2);
+
+      /* Then the resolution */
+      sub2 = GetSubString(sub1, ",", 2);
+      gridparams->resolution = (CCTK_INT4)atoi(sub2);
+      free(sub2);
+
+      /* And finally the "uniq_vertices" flag */
+      sub2 = GetSubString(sub1, ",", 3);
+      gridparams->uniq_verts = (CCTK_INT4)atoi(sub2);
+      free(sub2);
+
+      free(sub1);
+
+    free(sub0);
+  }
+
+  return n;
+}
+
+/*****************************************************************************/
+
+CCTK_INT4
+NumberOfItems(char *s, char *d)
+{
+  int   ld, i;
+  int   result;
+  char *p;
+
+  result = 0;
+  ld = strlen(d);
+
+  switch(ld)
+  {
+    case 1:
+      p = s;
+      while( *p != '\0' )
+      {
+         if( *p++ == d[0] )
+           result++;
+      }
+      if( p != s )
+         result++;
+      return result;
+
+    case 2:
+      p = s;
+      i = 0;
+      while( *p != '\0' )
+      {
+         if ( *p++ == d[i] )
+         {
+           result++;
+           i = (i==0)?1:0;
+         }
+      }
+      if( (result % 2) != 0 )
+         return 0;
+      else
+         return result/2;
+
+    default:
+      return 0;
+  }
+}
+
+/*****************************************************************************/
+
+char *
+GetSubString(char *s, char *d, CCTK_INT4 n)
+{
+  int   ld, count, i;
+  char *result;
+  char *p, *q = NULL;
+
+  result = NULL;
+  ld = strlen(d);
+
+  p = s;
+  count = 0;
+  switch(ld)
+  {
+    case 1:
+      q = p;
+      while( *p != '\0' )
+      {
+         if( *p++ == d[0] )
+         {
+           if(count == n)
+             break;
+           else
+           {
+             count++;
+             q = p;
+           }
+         }
+      }
+      if(*p == '\0')
+         p++;
+      break;
+
+    case 2:
+      i = 0;
+      while( *p != '\0' )
+      {
+         if ( *p++ == d[i] )
+         {
+           if(i==0)
+           {
+             q = p;
+           }
+           else
+           {
+             if(count == n)
+               break;
+             else
+               count++;
+           }
+           i = (i==0)?1:0;
+         }
+      }
+      break;
+
+    default:
+      return NULL;
+  }
+
+  if(count == n)
+  {
+    while( *q == ' ' )
+      q++;
+    p-=2;
+    while( *p == ' ' )
+      p--;
+    p++;
+    result = NEW((p-q+1), char);
+    strncpy(result, q, p-q);
+    result[p-q] = '\0';
+  }
+  return result;
+}
+
+/**********************************************************/
diff --git a/src/IsoSurfacerInit.h b/src/IsoSurfacerInit.h
new file mode 100644
index 0000000..a20b5ac
--- /dev/null
+++ b/src/IsoSurfacerInit.h
@@ -0,0 +1,38 @@
+#ifndef _ISOSURFACERINIT_H
+#define _ISOSURFACERINIT_H
+
+#include <sys/types.h>
+
+#ifndef _WIN32
+#include <sys/times.h>
+#else
+#include <time.h>
+#endif
+
+#include <limits.h>
+#include "common.h"
+
+typedef struct isotimes_str
+{
+#ifndef _WIN32
+  struct tms tms_begin;
+  struct tms tms_end;
+  struct tms tms_total;
+#endif
+
+  clock_t    realtime_begin;
+  clock_t    realtime_end;
+  clock_t    realtime_total;
+
+  CCTK_INT4       ncalls;
+} isotimes_st;
+
+/* prototypes of functions to be registered */
+void *IsoSurfacer_SetupGH (tFleshConfig *config, int convergence_level,cGH *GH);
+int IsoSurfacer_InitGH (cGH *GH);
+int IsoSurfacer (cGH *GH);
+int IsoSurfacer_TriggerOutput (cGH *GH, int);
+int IsoSurfacer_TimeForOutput (cGH *GH, int);
+
+#endif
+
diff --git a/src/NuSurfacer.c b/src/NuSurfacer.c
new file mode 100644
index 0000000..d12aba8
--- /dev/null
+++ b/src/NuSurfacer.c
@@ -0,0 +1,472 @@
+#ifndef TCP
+#define TCP
+#endif
+#include <stdlib.h>
+#include <stdio.h>
+#include "common.h"
+#include "cctk.h"
+
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "NuSurfacer.h"
+
+typedef int EDGE_LIST;
+typedef struct {
+  EDGE_LIST edges[16];
+} TRIANGLE_CASES;
+ 
+static TRIANGLE_CASES triCases[] = { 
+  {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 0 0 */
+  {{ 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 1 1 */
+  {{ 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 2 1 */
+  {{ 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 3 2 */
+  {{ 1, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 4 1 */
+  {{ 0, 3, 8, 1, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 5 3 */
+  {{ 9, 11, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 6 2 */
+  {{ 2, 3, 8, 2, 8, 11, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 7 5 */
+  {{ 3, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 8 1 */
+  {{ 0, 2, 10, 8, 0, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 9 2 */
+  {{ 1, 0, 9, 2, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 10 3 */
+  {{ 1, 2, 10, 1, 10, 9, 9, 10, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 11 5 */
+  {{ 3, 1, 11, 10, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 12 2 */
+  {{ 0, 1, 11, 0, 11, 8, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 13 5 */
+  {{ 3, 0, 9, 3, 9, 10, 10, 9, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 14 5 */
+  {{ 9, 11, 8, 11, 10, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 15 8 */
+  {{ 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 16 1 */
+  {{ 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 17 2 */
+  {{ 0, 9, 1, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 18 3 */
+  {{ 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 19 5 */
+  {{ 1, 11, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 20 4 */
+  {{ 3, 7, 4, 3, 4, 0, 1, 11, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 21 7 */
+  {{ 9, 11, 2, 9, 2, 0, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 22 7 */
+  {{ 2, 9, 11, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}}, /* 23 14 */
+  {{ 8, 7, 4, 3, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 24 3 */
+  {{10, 7, 4, 10, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 25 5 */
+  {{ 9, 1, 0, 8, 7, 4, 2, 10, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 26 6 */
+  {{ 4, 10, 7, 9, 10, 4, 9, 2, 10, 9, 1, 2, -1, -1, -1, -1}}, /* 27 9 */
+  {{ 3, 1, 11, 3, 11, 10, 7, 4, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 28 7 */
+  {{ 1, 11, 10, 1, 10, 4, 1, 4, 0, 7, 4, 10, -1, -1, -1, -1}}, /* 29 11 */
+  {{ 4, 8, 7, 9, 10, 0, 9, 11, 10, 10, 3, 0, -1, -1, -1, -1}}, /* 30 12 */
+  {{ 4, 10, 7, 4, 9, 10, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 31 5 */
+  {{ 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 32 1 */
+  {{ 9, 4, 5, 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 33 3 */
+  {{ 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 34 2 */
+  {{ 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 35 5 */
+  {{ 1, 11, 2, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 36 3 */
+  {{ 3, 8, 0, 1, 11, 2, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 37 6 */
+  {{ 5, 11, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 38 5 */
+  {{ 2, 5, 11, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}}, /* 39 9 */
+  {{ 9, 4, 5, 2, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 40 4 */
+  {{ 0, 2, 10, 0, 10, 8, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 41 7 */
+  {{ 0, 4, 5, 0, 5, 1, 2, 10, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 42 7 */
+  {{ 2, 5, 1, 2, 8, 5, 2, 10, 8, 4, 5, 8, -1, -1, -1, -1}}, /* 43 11 */
+  {{11, 10, 3, 11, 3, 1, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 44 7 */
+  {{ 4, 5, 9, 0, 1, 8, 8, 1, 11, 8, 11, 10, -1, -1, -1, -1}}, /* 45 12 */
+  {{ 5, 0, 4, 5, 10, 0, 5, 11, 10, 10, 3, 0, -1, -1, -1, -1}}, /* 46 14 */
+  {{ 5, 8, 4, 5, 11, 8, 11, 10, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 47 5 */
+  {{ 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 48 2 */
+  {{ 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 49 5 */
+  {{ 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 50 5 */
+  {{ 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 51 8 */
+  {{ 9, 8, 7, 9, 7, 5, 11, 2, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 52 7 */
+  {{11, 2, 1, 9, 0, 5, 5, 0, 3, 5, 3, 7, -1, -1, -1, -1}}, /* 53 12 */
+  {{ 8, 2, 0, 8, 5, 2, 8, 7, 5, 11, 2, 5, -1, -1, -1, -1}}, /* 54 11 */
+  {{ 2, 5, 11, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 55 5 */
+  {{ 7, 5, 9, 7, 9, 8, 3, 2, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 56 7 */
+  {{ 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 10, 7, -1, -1, -1, -1}}, /* 57 14 */
+  {{ 2, 10, 3, 0, 8, 1, 1, 8, 7, 1, 7, 5, -1, -1, -1, -1}}, /* 58 12 */
+  {{10, 1, 2, 10, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 59 5 */
+  {{ 9, 8, 5, 8, 7, 5, 11, 3, 1, 11, 10, 3, -1, -1, -1, -1}}, /* 60 10 */
+  {{ 5, 0, 7, 5, 9, 0, 7, 0, 10, 1, 11, 0, 10, 0, 11, -1}}, /* 61 7 */
+  {{10, 0, 11, 10, 3, 0, 11, 0, 5, 8, 7, 0, 5, 0, 7, -1}}, /* 62 7 */
+  {{10, 5, 11, 7, 5, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 63 2 */
+  {{11, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 64 1 */
+  {{ 0, 3, 8, 5, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 65 4 */
+  {{ 9, 1, 0, 5, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 66 3 */
+  {{ 1, 3, 8, 1, 8, 9, 5, 6, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 67 7 */
+  {{ 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 68 2 */
+  {{ 1, 5, 6, 1, 6, 2, 3, 8, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 69 7 */
+  {{ 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 70 5 */
+  {{ 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}}, /* 71 11 */
+  {{ 2, 10, 3, 11, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 72 3 */
+  {{10, 8, 0, 10, 0, 2, 11, 5, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 73 7 */
+  {{ 0, 9, 1, 2, 10, 3, 5, 6, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 74 6 */
+  {{ 5, 6, 11, 1, 2, 9, 9, 2, 10, 9, 10, 8, -1, -1, -1, -1}}, /* 75 12 */
+  {{ 6, 10, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 76 5 */
+  {{ 0, 10, 8, 0, 5, 10, 0, 1, 5, 5, 6, 10, -1, -1, -1, -1}}, /* 77 14 */
+  {{ 3, 6, 10, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}}, /* 78 9 */
+  {{ 6, 9, 5, 6, 10, 9, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 79 5 */
+  {{ 5, 6, 11, 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 80 3 */
+  {{ 4, 0, 3, 4, 3, 7, 6, 11, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 81 7 */
+  {{ 1, 0, 9, 5, 6, 11, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 82 6 */
+  {{11, 5, 6, 1, 7, 9, 1, 3, 7, 7, 4, 9, -1, -1, -1, -1}}, /* 83 12 */
+  {{ 6, 2, 1, 6, 1, 5, 4, 8, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 84 7 */
+  {{ 1, 5, 2, 5, 6, 2, 3, 4, 0, 3, 7, 4, -1, -1, -1, -1}}, /* 85 10 */
+  {{ 8, 7, 4, 9, 5, 0, 0, 5, 6, 0, 6, 2, -1, -1, -1, -1}}, /* 86 12 */
+  {{ 7, 9, 3, 7, 4, 9, 3, 9, 2, 5, 6, 9, 2, 9, 6, -1}}, /* 87 7 */
+  {{ 3, 2, 10, 7, 4, 8, 11, 5, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 88 6 */
+  {{ 5, 6, 11, 4, 2, 7, 4, 0, 2, 2, 10, 7, -1, -1, -1, -1}}, /* 89 12 */
+  {{ 0, 9, 1, 4, 8, 7, 2, 10, 3, 5, 6, 11, -1, -1, -1, -1}}, /* 90 13 */
+  {{ 9, 1, 2, 9, 2, 10, 9, 10, 4, 7, 4, 10, 5, 6, 11, -1}}, /* 91 6 */
+  {{ 8, 7, 4, 3, 5, 10, 3, 1, 5, 5, 6, 10, -1, -1, -1, -1}}, /* 92 12 */
+  {{ 5, 10, 1, 5, 6, 10, 1, 10, 0, 7, 4, 10, 0, 10, 4, -1}}, /* 93 7 */
+  {{ 0, 9, 5, 0, 5, 6, 0, 6, 3, 10, 3, 6, 8, 7, 4, -1}}, /* 94 6 */
+  {{ 6, 9, 5, 6, 10, 9, 4, 9, 7, 7, 9, 10, -1, -1, -1, -1}}, /* 95 3 */
+  {{11, 9, 4, 6, 11, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 96 2 */
+  {{ 4, 6, 11, 4, 11, 9, 0, 3, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 97 7 */
+  {{11, 1, 0, 11, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 98 5 */
+  {{ 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 11, 1, -1, -1, -1, -1}}, /* 99 14 */
+  {{ 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 100 5 */
+  {{ 3, 8, 0, 1, 9, 2, 2, 9, 4, 2, 4, 6, -1, -1, -1, -1}}, /* 101 12 */
+  {{ 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 102 8 */
+  {{ 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 103 5 */
+  {{11, 9, 4, 11, 4, 6, 10, 3, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 104 7 */
+  {{ 0, 2, 8, 2, 10, 8, 4, 11, 9, 4, 6, 11, -1, -1, -1, -1}}, /* 105 10 */
+  {{ 3, 2, 10, 0, 6, 1, 0, 4, 6, 6, 11, 1, -1, -1, -1, -1}}, /* 106 12 */
+  {{ 6, 1, 4, 6, 11, 1, 4, 1, 8, 2, 10, 1, 8, 1, 10, -1}}, /* 107 7 */
+  {{ 9, 4, 6, 9, 6, 3, 9, 3, 1, 10, 3, 6, -1, -1, -1, -1}}, /* 108 11 */
+  {{ 8, 1, 10, 8, 0, 1, 10, 1, 6, 9, 4, 1, 6, 1, 4, -1}}, /* 109 7 */
+  {{ 3, 6, 10, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 110 5 */
+  {{ 6, 8, 4, 10, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 111 2 */
+  {{ 7, 6, 11, 7, 11, 8, 8, 11, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 112 5 */
+  {{ 0, 3, 7, 0, 7, 11, 0, 11, 9, 6, 11, 7, -1, -1, -1, -1}}, /* 113 11 */
+  {{11, 7, 6, 1, 7, 11, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}}, /* 114 9 */
+  {{11, 7, 6, 11, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 115 5 */
+  {{ 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}}, /* 116 14 */
+  {{ 2, 9, 6, 2, 1, 9, 6, 9, 7, 0, 3, 9, 7, 9, 3, -1}}, /* 117 7 */
+  {{ 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 118 5 */
+  {{ 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 119 2 */
+  {{ 2, 10, 3, 11, 8, 6, 11, 9, 8, 8, 7, 6, -1, -1, -1, -1}}, /* 120 12 */
+  {{ 2, 7, 0, 2, 10, 7, 0, 7, 9, 6, 11, 7, 9, 7, 11, -1}}, /* 121 7 */
+  {{ 1, 0, 8, 1, 8, 7, 1, 7, 11, 6, 11, 7, 2, 10, 3, -1}}, /* 122 6 */
+  {{10, 1, 2, 10, 7, 1, 11, 1, 6, 6, 1, 7, -1, -1, -1, -1}}, /* 123 3 */
+  {{ 8, 6, 9, 8, 7, 6, 9, 6, 1, 10, 3, 6, 1, 6, 3, -1}}, /* 124 7 */
+  {{ 0, 1, 9, 10, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 125 4 */
+  {{ 7, 0, 8, 7, 6, 0, 3, 0, 10, 10, 0, 6, -1, -1, -1, -1}}, /* 126 3 */
+  {{ 7, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 127 1 */
+  {{ 7, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 128 1 */
+  {{ 3, 8, 0, 10, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 129 3 */
+  {{ 0, 9, 1, 10, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 130 4 */
+  {{ 8, 9, 1, 8, 1, 3, 10, 6, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 131 7 */
+  {{11, 2, 1, 6, 7, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 132 3 */
+  {{ 1, 11, 2, 3, 8, 0, 6, 7, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 133 6 */
+  {{ 2, 0, 9, 2, 9, 11, 6, 7, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 134 7 */
+  {{ 6, 7, 10, 2, 3, 11, 11, 3, 8, 11, 8, 9, -1, -1, -1, -1}}, /* 135 12 */
+  {{ 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 136 2 */
+  {{ 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 137 5 */
+  {{ 2, 6, 7, 2, 7, 3, 0, 9, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 138 7 */
+  {{ 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}}, /* 139 14 */
+  {{11, 6, 7, 11, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 140 5 */
+  {{11, 6, 7, 1, 11, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}}, /* 141 9 */
+  {{ 0, 7, 3, 0, 11, 7, 0, 9, 11, 6, 7, 11, -1, -1, -1, -1}}, /* 142 11 */
+  {{ 7, 11, 6, 7, 8, 11, 8, 9, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 143 5 */
+  {{ 6, 4, 8, 10, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 144 2 */
+  {{ 3, 10, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 145 5 */
+  {{ 8, 10, 6, 8, 6, 4, 9, 1, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 146 7 */
+  {{ 9, 6, 4, 9, 3, 6, 9, 1, 3, 10, 6, 3, -1, -1, -1, -1}}, /* 147 11 */
+  {{ 6, 4, 8, 6, 8, 10, 2, 1, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 148 7 */
+  {{ 1, 11, 2, 3, 10, 0, 0, 10, 6, 0, 6, 4, -1, -1, -1, -1}}, /* 149 12 */
+  {{ 4, 8, 10, 4, 10, 6, 0, 9, 2, 2, 9, 11, -1, -1, -1, -1}}, /* 150 10 */
+  {{11, 3, 9, 11, 2, 3, 9, 3, 4, 10, 6, 3, 4, 3, 6, -1}}, /* 151 7 */
+  {{ 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 152 5 */
+  {{ 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 153 8 */
+  {{ 1, 0, 9, 2, 4, 3, 2, 6, 4, 4, 8, 3, -1, -1, -1, -1}}, /* 154 12 */
+  {{ 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 155 5 */
+  {{ 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 11, -1, -1, -1, -1}}, /* 156 14 */
+  {{11, 0, 1, 11, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}}, /* 157 5 */
+  {{ 4, 3, 6, 4, 8, 3, 6, 3, 11, 0, 9, 3, 11, 3, 9, -1}}, /* 158 7 */
+  {{11, 4, 9, 6, 4, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 159 2 */
+  {{ 4, 5, 9, 7, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 160 3 */
+  {{ 0, 3, 8, 4, 5, 9, 10, 6, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 161 6 */
+  {{ 5, 1, 0, 5, 0, 4, 7, 10, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 162 7 */
+  {{10, 6, 7, 8, 4, 3, 3, 4, 5, 3, 5, 1, -1, -1, -1, -1}}, /* 163 12 */
+  {{ 9, 4, 5, 11, 2, 1, 7, 10, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 164 6 */
+  {{ 6, 7, 10, 1, 11, 2, 0, 3, 8, 4, 5, 9, -1, -1, -1, -1}}, /* 165 13 */
+  {{ 7, 10, 6, 5, 11, 4, 4, 11, 2, 4, 2, 0, -1, -1, -1, -1}}, /* 166 12 */
+  {{ 3, 8, 4, 3, 4, 5, 3, 5, 2, 11, 2, 5, 10, 6, 7, -1}}, /* 167 6 */
+  {{ 7, 3, 2, 7, 2, 6, 5, 9, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 168 7 */
+  {{ 9, 4, 5, 0, 6, 8, 0, 2, 6, 6, 7, 8, -1, -1, -1, -1}}, /* 169 12 */
+  {{ 3, 2, 6, 3, 6, 7, 1, 0, 5, 5, 0, 4, -1, -1, -1, -1}}, /* 170 10 */
+  {{ 6, 8, 2, 6, 7, 8, 2, 8, 1, 4, 5, 8, 1, 8, 5, -1}}, /* 171 7 */
+  {{ 9, 4, 5, 11, 6, 1, 1, 6, 7, 1, 7, 3, -1, -1, -1, -1}}, /* 172 12 */
+  {{ 1, 11, 6, 1, 6, 7, 1, 7, 0, 8, 0, 7, 9, 4, 5, -1}}, /* 173 6 */
+  {{ 4, 11, 0, 4, 5, 11, 0, 11, 3, 6, 7, 11, 3, 11, 7, -1}}, /* 174 7 */
+  {{ 7, 11, 6, 7, 8, 11, 5, 11, 4, 4, 11, 8, -1, -1, -1, -1}}, /* 175 3 */
+  {{ 6, 5, 9, 6, 9, 10, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 176 5 */
+  {{ 3, 10, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}}, /* 177 9 */
+  {{ 0, 8, 10, 0, 10, 5, 0, 5, 1, 5, 10, 6, -1, -1, -1, -1}}, /* 178 14 */
+  {{ 6, 3, 10, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 179 5 */
+  {{ 1, 11, 2, 9, 10, 5, 9, 8, 10, 10, 6, 5, -1, -1, -1, -1}}, /* 180 12 */
+  {{ 0, 3, 10, 0, 10, 6, 0, 6, 9, 5, 9, 6, 1, 11, 2, -1}}, /* 181 6 */
+  {{10, 5, 8, 10, 6, 5, 8, 5, 0, 11, 2, 5, 0, 5, 2, -1}}, /* 182 7 */
+  {{ 6, 3, 10, 6, 5, 3, 2, 3, 11, 11, 3, 5, -1, -1, -1, -1}}, /* 183 3 */
+  {{ 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}}, /* 184 11 */
+  {{ 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}}, /* 185 5 */
+  {{ 1, 8, 5, 1, 0, 8, 5, 8, 6, 3, 2, 8, 6, 8, 2, -1}}, /* 186 7 */
+  {{ 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 187 2 */
+  {{ 1, 6, 3, 1, 11, 6, 3, 6, 8, 5, 9, 6, 8, 6, 9, -1}}, /* 188 7 */
+  {{11, 0, 1, 11, 6, 0, 9, 0, 5, 5, 0, 6, -1, -1, -1, -1}}, /* 189 3 */
+  {{ 0, 8, 3, 5, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 190 4 */
+  {{11, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 191 1 */
+  {{10, 11, 5, 7, 10, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 192 2 */
+  {{10, 11, 5, 10, 5, 7, 8, 0, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 193 7 */
+  {{ 5, 7, 10, 5, 10, 11, 1, 0, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 194 7 */
+  {{11, 5, 7, 11, 7, 10, 9, 1, 8, 8, 1, 3, -1, -1, -1, -1}}, /* 195 10 */
+  {{10, 2, 1, 10, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 196 5 */
+  {{ 0, 3, 8, 1, 7, 2, 1, 5, 7, 7, 10, 2, -1, -1, -1, -1}}, /* 197 12 */
+  {{ 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 10, -1, -1, -1, -1}}, /* 198 14 */
+  {{ 7, 2, 5, 7, 10, 2, 5, 2, 9, 3, 8, 2, 9, 2, 8, -1}}, /* 199 7 */
+  {{ 2, 11, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 200 5 */
+  {{ 8, 0, 2, 8, 2, 5, 8, 5, 7, 11, 5, 2, -1, -1, -1, -1}}, /* 201 11 */
+  {{ 9, 1, 0, 5, 3, 11, 5, 7, 3, 3, 2, 11, -1, -1, -1, -1}}, /* 202 12 */
+  {{ 9, 2, 8, 9, 1, 2, 8, 2, 7, 11, 5, 2, 7, 2, 5, -1}}, /* 203 7 */
+  {{ 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 204 8 */
+  {{ 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}}, /* 205 5 */
+  {{ 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}}, /* 206 5 */
+  {{ 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 207 2 */
+  {{ 5, 4, 8, 5, 8, 11, 11, 8, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 208 5 */
+  {{ 5, 4, 0, 5, 0, 10, 5, 10, 11, 10, 0, 3, -1, -1, -1, -1}}, /* 209 14 */
+  {{ 0, 9, 1, 8, 11, 4, 8, 10, 11, 11, 5, 4, -1, -1, -1, -1}}, /* 210 12 */
+  {{11, 4, 10, 11, 5, 4, 10, 4, 3, 9, 1, 4, 3, 4, 1, -1}}, /* 211 7 */
+  {{ 2, 1, 5, 2, 5, 8, 2, 8, 10, 4, 8, 5, -1, -1, -1, -1}}, /* 212 11 */
+  {{ 0, 10, 4, 0, 3, 10, 4, 10, 5, 2, 1, 10, 5, 10, 1, -1}}, /* 213 7 */
+  {{ 0, 5, 2, 0, 9, 5, 2, 5, 10, 4, 8, 5, 10, 5, 8, -1}}, /* 214 7 */
+  {{ 9, 5, 4, 2, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 215 4 */
+  {{ 2, 11, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}}, /* 216 9 */
+  {{ 5, 2, 11, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}}, /* 217 5 */
+  {{ 3, 2, 11, 3, 11, 5, 3, 5, 8, 4, 8, 5, 0, 9, 1, -1}}, /* 218 6 */
+  {{ 5, 2, 11, 5, 4, 2, 1, 2, 9, 9, 2, 4, -1, -1, -1, -1}}, /* 219 3 */
+  {{ 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}}, /* 220 5 */
+  {{ 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 221 2 */
+  {{ 8, 5, 4, 8, 3, 5, 9, 5, 0, 0, 5, 3, -1, -1, -1, -1}}, /* 222 3 */
+  {{ 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 223 1 */
+  {{ 4, 7, 10, 4, 10, 9, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 224 5 */
+  {{ 0, 3, 8, 4, 7, 9, 9, 7, 10, 9, 10, 11, -1, -1, -1, -1}}, /* 225 12 */
+  {{ 1, 10, 11, 1, 4, 10, 1, 0, 4, 7, 10, 4, -1, -1, -1, -1}}, /* 226 11 */
+  {{ 3, 4, 1, 3, 8, 4, 1, 4, 11, 7, 10, 4, 11, 4, 10, -1}}, /* 227 7 */
+  {{ 4, 7, 10, 9, 4, 10, 9, 10, 2, 9, 2, 1, -1, -1, -1, -1}}, /* 228 9 */
+  {{ 9, 4, 7, 9, 7, 10, 9, 10, 1, 2, 1, 10, 0, 3, 8, -1}}, /* 229 6 */
+  {{10, 4, 7, 10, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}}, /* 230 5 */
+  {{10, 4, 7, 10, 2, 4, 8, 4, 3, 3, 4, 2, -1, -1, -1, -1}}, /* 231 3 */
+  {{ 2, 11, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}}, /* 232 14 */
+  {{ 9, 7, 11, 9, 4, 7, 11, 7, 2, 8, 0, 7, 2, 7, 0, -1}}, /* 233 7 */
+  {{ 3, 11, 7, 3, 2, 11, 7, 11, 4, 1, 0, 11, 4, 11, 0, -1}}, /* 234 7 */
+  {{ 1, 2, 11, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 235 4 */
+  {{ 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}}, /* 236 5 */
+  {{ 4, 1, 9, 4, 7, 1, 0, 1, 8, 8, 1, 7, -1, -1, -1, -1}}, /* 237 3 */
+  {{ 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 238 2 */
+  {{ 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 239 1 */
+  {{ 9, 8, 11, 11, 8, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 240 8 */
+  {{ 3, 9, 0, 3, 10, 9, 10, 11, 9, -1, -1, -1, -1, -1, -1, -1}}, /* 241 5 */
+  {{ 0, 11, 1, 0, 8, 11, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}}, /* 242 5 */
+  {{ 3, 11, 1, 10, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 243 2 */
+  {{ 1, 10, 2, 1, 9, 10, 9, 8, 10, -1, -1, -1, -1, -1, -1, -1}}, /* 244 5 */
+  {{ 3, 9, 0, 3, 10, 9, 1, 9, 2, 2, 9, 10, -1, -1, -1, -1}}, /* 245 3 */
+  {{ 0, 10, 2, 8, 10, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 246 2 */
+  {{ 3, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 247 1 */
+  {{ 2, 8, 3, 2, 11, 8, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}}, /* 248 5 */
+  {{ 9, 2, 11, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 249 2 */
+  {{ 2, 8, 3, 2, 11, 8, 0, 8, 1, 1, 8, 11, -1, -1, -1, -1}}, /* 250 3 */
+  {{ 1, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 251 1 */
+  {{ 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 252 2 */
+  {{ 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 253 1 */
+  {{ 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, /* 254 1 */
+  {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}}; /* 255 0 */
+
+#define New(typ,sz) ((typ*)malloc(sizeof(typ)*sz))
+#define Resize(typ,arry,newsz) ((typ*)realloc(arry,sizeof(typ)*newsz))
+#define GrowArray(typ,lst) { lst->maxelem*=2; lst->array=Resize(typ,lst->array,lst->maxelem);}
+#define InitArray(typ,lst) lst->array=New(typ,1024);lst->maxelem=1024;lst->nelem=0
+#define ResetArray(lst) lst->nelem=0
+#define ArraySize(lst) (lst->nelem)
+#define ArrayAddElem(typ,lst,val) { if(lst->maxelem<=lst->nelem) \
+				     GrowArray(typ,lst); \
+				     lst->array[lst->nelem]=val; \
+				     lst->nelem++;}
+
+#define ArrayAddNumElem(typ,lst,vals,nvals) { \
+		     if(lst->maxelem<=(lst->nelem+nvals)) \
+		        GrowArray(typ,lst); \
+		     for(nvals+=lst->nelem;lst->nelem<nvals;lst->nelem++) \
+	                lst->array[lst->nelem] = vals[lst->nelem]; }
+
+#define ComputePoints(pughGH, pindex,p1,p2) { int pointindex=pindex[0]; \
+p1[0]=(isofl)(i+xMapOffset[pointindex])*pughGH->dx0+pughGH->lx0;\
+p1[1]=(isofl)(j+yMapOffset[pointindex])*pughGH->dy0+pughGH->ly0;\
+p1[2]=(isofl)(k+zMapOffset[pointindex])*pughGH->dz0+pughGH->lz0;\
+pointindex=pindex[1];\
+p2[0]=(isofl)(i+xMapOffset[pointindex])*pughGH->dx0+pughGH->lx0;\
+p2[1]=(isofl)(j+yMapOffset[pointindex])*pughGH->dy0+pughGH->ly0;\
+p2[2]=(isofl)(k+zMapOffset[pointindex])*pughGH->dz0+pughGH->lz0;}
+
+#define SetEdgeMaskP(idx_,data_,value_,mask_) { \
+	register int _cmask=1; \
+	mask_=0; \
+        if(data_[idx_+cellOffset[0]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[1]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[2]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[3]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[4]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[5]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[6]]>=value_) mask_|=_cmask; _cmask<<=1; \
+	if(data_[idx_+cellOffset[7]]>=value_) mask_|=_cmask; \
+}
+
+#define AddInterpolatedVertex(isoval,idx,pindex,p1,p2,vertlist) { \
+  isofl _v,_d,_p1val,_p2val; int _i; \
+  _p1val = data[idx+cellOffset[pindex[0]]]; \
+  _p2val = data[idx+cellOffset[pindex[1]]]; \
+  _d=(isoval-_p1val)/(_p2val-_p1val); \
+  for( _i=0;_i<3;_i++){ \
+     _v=(p2[_i]-p1[_i])*_d + p1[_i]; \
+     ArrayAddElem(isofl,vertlist,_v); \
+  } \
+}
+
+#define DeclareArray(typ,typname) typedef struct typname { \
+  int nelem; \
+  int maxelem; \
+  typ *array; } typname;
+
+typedef struct IsoCell {
+  int v[3]; /* index into vertex list for 3 primary verts (if they exist) */
+} IsoCell;
+
+IsoCell *cellslice;
+DeclareArray(isofl,FloatArray); /* declare a list of floats for vertices */
+DeclareArray(CCTK_INT4,IntArray); /* declare a list of intes for triangle connectivity */
+FloatArray *vertexlist;
+IntArray *trianglelist;
+static int edges[12][2] = { {0,1}, {1,2}, {3,2}, {0,3},
+                            {4,5}, {5,6}, {7,6}, {4,7},
+                            {0,4}, {1,5}, {3,7}, {2,6} };
+/*  must also map pairs of vertices to the stored vertex offsets */
+static int xMapOffset[8]={0,1,1,0,0,1,1,0};
+static int yMapOffset[8]={0,0,1,1,0,0,1,1};
+static int zMapOffset[8]={0,0,0,0,1,1,1,1};
+static int cellOffset[8]={0,0,0,0,0,0,0,0};
+/* cell Offsets should be initialized to
+		0 = 0
+		1 = 1
+		2 = 1 + GH->cctk_lsh [0]
+		3 = GH->cctk_lsh [0]
+		4 = GH->cctk_lsh [0]*GH->cctk_lsh [1]
+		5 = 1 + GH->cctk_lsh [0]*GH->cctk_lsh [1]
+		6 = 1+ GH->cctk_lsh [0] + GH->cctk_lsh [0]*GH->cctk_lsh [1]
+		7 = GH->cctk_lsh [0] + GH->cctk_lsh [0]*GH->cctk_lsh [1]
+	*/
+static int isoCellOffset[12]={0,0,0,0, /* -zy,-z,-z,-zx, */
+			      0,0,0,0, /* -y,0,0,-x, */
+			      0,0,0,0};/* -yx,-y,-x,0 */
+static int isoCellVindex[12]={1,0,1,0,
+			      1,0,1,0,
+			      2,2,2,2};
+							
+void InitCellArrays(cGH *GH){
+  int i,j;
+  int nx=GH->cctk_lsh [0],ny=GH->cctk_lsh [1];
+  int ix=1,iy=nx,iz=nx*ny;
+  cellslice = New(IsoCell,2*nx*ny); /* 2 slices */
+  for(i=0;i<nx*ny*2;i++) 
+    for(j=0;j<3;j++) cellslice[i].v[j]=-1; /* init all verts */
+  vertexlist = New(FloatArray,1); InitArray(isofl,vertexlist);
+  trianglelist = New(IntArray,1); InitArray(CCTK_INT4,trianglelist);
+  /* celloffset */
+  cellOffset[0] = 0;
+  cellOffset[1] = ix;
+  cellOffset[2] = ix + iy;
+  cellOffset[3] = iy;
+  cellOffset[4] = iz;
+  cellOffset[5] = ix + iz;
+  cellOffset[6] = ix + iy + iz;
+  cellOffset[7] = iy + iz;
+  /* isocelloffset */
+  isoCellOffset[0]=-(iz+iy); isoCellOffset[1]=-(iz); 
+  isoCellOffset[2]=-(iz);    isoCellOffset[3]=-(iz+ix); /* -zy,-z,-z,-zx, */
+  isoCellOffset[4]=-(iy);    isoCellOffset[5]=0;
+  isoCellOffset[6]=0;        isoCellOffset[7]=-(ix);/* -y,0,0,-x, */
+  isoCellOffset[8]=-(iy+ix); isoCellOffset[9]=-(iy);
+  isoCellOffset[10]=-(ix);   isoCellOffset[11]=0;/* -yx,-y,-x,0*/
+}
+
+void NuFindSurface(cGH *GH, int index,CCTK_REAL isovalue,polypatch *results){
+  static int initialized_nusurf=0;
+  TRIANGLE_CASES *triCase;
+  EDGE_LIST  *edge;
+  int i,j,k,imax,jmax,kmax;
+  int timelevel;
+  int idx,nverts=0; /* vertex and triangle counters for offsets */
+  CCTK_REAL *data;
+  int nx=GH->cctk_lsh [0],ny=GH->cctk_lsh [1],nz=GH->cctk_lsh [2],slicesize,slicemodulo;
+  pGH *pughGH;
+
+  pughGH = PUGH_pGH (GH);
+
+  if(!initialized_nusurf){
+    initialized_nusurf=1;
+    InitCellArrays(GH);
+  }
+  ResetArray(vertexlist);   /* zero the vertex count (but keep storage) */
+  ResetArray(trianglelist); /* zero the triangle count (but keep storage) */
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (index);
+  if (timelevel > 0)
+    timelevel--;
+  data = (CCTK_REAL *) GH->data [index][timelevel];
+
+  /*****Sliceinit */
+  slicesize=nx*ny; slicemodulo=slicesize*2;
+  if(!cellslice) cellslice = New(IsoCell,slicemodulo);
+  /* Initialize the cellslice mapping array */
+  for(idx=0;idx<slicemodulo;idx++)
+    cellslice[idx].v[0]=cellslice[idx].v[1]=cellslice[idx].v[2]=-1;
+  /* idx indexes into data and the cellslice. */
+  for(k=0,kmax=nz-1;k<kmax;k++){
+    int ii;
+    int joff,koff;
+    int isocellidx;
+    int mask;
+    int edgenumber;
+    isofl p1[3],p2[3]; 
+    int isocellnum;
+    int isocellvindex; /* 12 maps to 3 */
+    int elem,*elemp;
+    koff=k*nx*ny;
+    for(j=0,jmax=ny-1;j<jmax;j++){
+      joff=j*nx + koff;
+      for(i=0,imax=nx-1;i<imax;i++){
+	idx=i+joff;
+	isocellidx=idx%slicemodulo;
+	/* Here we clear the cellslice (vertex mapping) array for current cell */
+	cellslice[isocellidx].v[0]=cellslice[isocellidx].v[1]=cellslice[isocellidx].v[2]=-1;
+	SetEdgeMaskP(idx,data,isovalue,mask);
+	if ( mask == 0 || mask == 255 ) continue; /* no surface */
+	triCase = triCases + mask;
+	/* Only 3 edges (must scan for them) and set the values 
+	   So this is adding the verts as well as storing. */
+	for(edge=triCase->edges,ii=0;edge[ii]>=0;ii++){
+	  edgenumber = edge[ii];
+	  /* map to a new isocelloffset based on edge
+	     full 3d offsets into cell array based on edge pair */
+	  isocellnum = (isocellidx +isoCellOffset[edgenumber] + slicemodulo)%slicemodulo;
+	  /* for 12 edges maps to neg offsets 
+	     just a 3 element offset based on edge pair */
+	  isocellvindex = isoCellVindex[edgenumber]; /* 12 maps to 3 */
+	  elemp = &(cellslice[isocellnum].v[isocellvindex]);
+	  elem=*elemp;
+	  if(elem<0){ /* add that vertex */
+	    elem=*elemp=nverts;
+	    ComputePoints(pughGH, edges[edgenumber],p1,p2);
+	    AddInterpolatedVertex(isovalue,idx,edges[edgenumber],p1,p2,vertexlist);
+	    nverts++;
+	  }
+	  ArrayAddElem(CCTK_INT4,trianglelist,elem); /* this is just one vertex! */
+	}
+      }
+    }
+  }
+  results->nverts = nverts; /* actual listlen = 3*nverts */
+  results->verts = vertexlist->array;
+  results->npolys = (trianglelist->nelem)?(trianglelist->nelem/3):(trianglelist->nelem);
+  results->polys = trianglelist->array;
+}	 
+  
diff --git a/src/NuSurfacer.h b/src/NuSurfacer.h
new file mode 100644
index 0000000..8089cdd
--- /dev/null
+++ b/src/NuSurfacer.h
@@ -0,0 +1,8 @@
+#ifndef _NUSURFACER_H_
+#define _NUSURFACER_H_
+
+#include "FindIsoSurface.h"
+void InitCellArrays(cGH *gh);
+void NuFindSurface(cGH *gh, int index,CCTK_REAL isovalue,polypatch *results);
+
+#endif
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..9a93ad9
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,35 @@
+ /*@@
+   @file      Startup.c
+   @date      Fri May 21 1999
+   @author    Thomas Radke
+   @desc 
+   Startup routines for IsoSurfacer.
+   @enddesc 
+   @history
+   @endhistory
+ @@*/
+
+#include <stdio.h>
+
+#include "cctk_Flesh.h"
+#include "cctk_GHExtensions.h"
+#include "cctk_IOMethods.h"
+#include "IsoSurfacerInit.h"
+
+void IsoSurfacer_Startup()
+{
+  int IOMethod;
+  int IsoSurfacer_GHExtension;
+
+  IsoSurfacer_GHExtension = CCTK_RegisterGHExtension ("IsoSurfacer");
+  CCTK_RegisterGHExtensionSetupGH (IsoSurfacer_GHExtension,IsoSurfacer_SetupGH);
+  CCTK_RegisterGHExtensionInitGH (IsoSurfacer_GHExtension, IsoSurfacer_InitGH);
+
+  /* Register the IsoSurfacer routine as output method */
+  IOMethod = CCTK_RegisterIOMethod ("IsoSurfacer");
+  CCTK_RegisterIOMethodOutputGH (IOMethod, IsoSurfacer);
+  CCTK_RegisterIOMethodTimeToOutput (IOMethod, IsoSurfacer_TimeForOutput);
+  CCTK_RegisterIOMethodTriggerOutput (IOMethod, IsoSurfacer_TriggerOutput);
+
+}
+
diff --git a/src/common.h b/src/common.h
new file mode 100644
index 0000000..ef17b67
--- /dev/null
+++ b/src/common.h
@@ -0,0 +1,58 @@
+#ifndef _COMMON_H
+#define _COMMON_H
+
+#include "cctk.h"
+
+#define X 0
+#define Y 1
+#define Z 2
+#define V 3
+
+#define NEW(n,what) ((what*)malloc((n)*sizeof(what)))
+/*
+#define REALLOC(object,n,what) \
+  (object) = (what*)realloc((void *)(object), (n)*sizeof(what)); \
+  if((object) == NULL) printf("IsoSurfacer: realloc failed!\n")
+  */
+#define REALLOC(object,n,what) \
+  if((object) != NULL) free(object); \
+  (object) = NEW(n,what); \
+  if((object) == NULL) printf("IsoSurfacer: REALLOC failed!\n")
+
+#define ABS(x)   ((x) < 0 ? (-(x)) : (x))
+
+typedef CCTK_REAL4  isofl;
+typedef double extfl;
+
+typedef int machint;
+
+#define BIN   1
+#define ASCII 2
+#define UCD   4
+
+typedef struct
+{
+  isofl    *verts;
+  CCTK_INT4      nverts;
+  CCTK_INT4     *polys;
+  CCTK_INT4      npolys;
+} polypatch;
+
+
+typedef struct
+{
+  isofl    isovalue;
+  char    *outfname;
+  CCTK_INT4     outtype;
+  CCTK_INT4     step;
+  CCTK_INT4     doRemoveDuplicateVertices;
+  CCTK_INT4     useTree;
+  CCTK_INT4     doEliminateSmallTriangles;
+       
+//  CCTK_INT4    timestep;  timestep information is really not required here!
+  isofl   minval;
+  isofl   maxval;
+} par_st;
+
+
+#endif
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..4145cb4
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,13 @@
+# Main make.code.defn file for thorn IsoSurfacer
+# $Header$
+
+-include $(CCTK_HOME)/arrangements/external/TCPXX/src/make.code.defn
+
+# Source files in this directory
+SRCS =	Startup.c IsoSurfacerInit.c IsoSurfacer.c \
+       	NuSurfacer.c 
+
+INC_DIRS += $(CCTK_HOME)/arrangements/external/TCPXX/src/   \
+	     $(CCTK_HOME)/arrangements/external/RemoteIO/src 
+
+
diff --git a/src/make.configuration.defn b/src/make.configuration.defn
new file mode 100644
index 0000000..6af43cc
--- /dev/null
+++ b/src/make.configuration.defn
@@ -0,0 +1,8 @@
+
+ifeq ($(findstring CactusPUGH/PUGH,$(THORNS)),)
+.pseudo: MissingPUGHinIsoSurfacer
+MissingPUGHinIsoSurfacer:
+	@echo "IsoSurfacer: requires PUGH"
+	@echo "IsoSurfacer: Please add CactusPUGH/PUGH or remove IsoSurfacer from Thornlist !"
+	exit 2
+endif