The beginning of adding Florin's isosurfacer to calculate triangulations

of the EH and calculate it's area. Don't use it yet. This is still in
a development phase.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@183 2a26948c-0e4f-0410-aee8-f1d3e353619c

6 files changed, 716 insertions, 144 deletions

diff --git a/param.ccl b/param.ccl
index efa14b0..0bb8e7e 100644
--- a/param.ccl
+++ b/param.ccl
@@ -264,6 +264,12 @@ CCTK_INT area_interpolation_order "What order should be used for the area interp
   1:*           :: "A valid positive interpolation order"
 } 3
 
+KEYWORD area_calculation_method "How should the areas be calculated"
+{
+  "standard"    :: "Using a angular coordinate system on the surface"
+  "isosurface"  :: "Using an isosurface triangulation"
+} "standard"
+
 STRING generator_interpolator "What interpolator should be used for the generators"
 {
   ".+"          :: "A valid interpolator name"
diff --git a/schedule.ccl b/schedule.ccl
index 9707f63..1550011 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -118,154 +118,172 @@ schedule EHFinder_InitWeights at CCTK_POST_RECOVER_VARIABLES
   LANG: Fortran
 } "Setup weights for Simpson integration"
 
-schedule GROUP EHFinder_Level_Sets at CCTK_ANALYSIS
+if (CCTK_EQUALS(area_calculation_method,"standard"))
 {
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Loop over the level set functions"
-
-schedule EHFinder_LevelSetLoopInit in EHFinder_Level_Sets
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Initialize the loop counter over the level set functions"
-
-schedule GROUP EHFinder_Surfaces in EHFinder_Level_Sets while ehfinder::more_levelsets
-{
-  STORAGE: surface_integers
-  STORAGE: surface_reals
-  STORAGE: surface_index
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Count the number of surfaces and integrate over them"
-
-
-schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Initialize while loop control"
-
-schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces
-{
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Counting and mark surfaces"
-
-schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-  SYNC: surface_index
-} "Check if there are more surfaces"
-
-schedule GROUP EHFinder_MarkPoints in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points
-{
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Marking surfaces"
-
-schedule EHFinder_MarkSurfaces in EHFinder_MarkPoints
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-  SYNC: surface_index
-} "Mark points inside the current surface"
-
-schedule EHFinder_ApplySymSC in EHFinder_MarkPoints after EHFinder_MarkSurfaces
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Select the surface counter grid function for boundary conditions"
-
-schedule GROUP ApplyBCs as EHFinderSC_ApplyBSc in EHFinder_MarkPoints after EHFinder_ApplySymSC
-{
-} "Apply boundary conditions (symmetries)"
-
-schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Output info about found surfaces"
-
-schedule group EHFinder_Integration in EHFinder_Surfaces after EHFinder_InfoSurfaces while ehfinder::integrate_counter
-{
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Find and integrate over surfaces"
-
-schedule EHFinder_FindSurface in EHFinder_Integration
-{
-  LANG: Fortran
-  STORAGE: surface_tmp_arrays, surface_int_array
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-  SYNC: surface_arrays
-} "Find Surface"
-
-schedule EHFinder_FindSurfaceElement in EHFinder_Integration after EHFinder_FindSurface
-{
-  LANG: Fortran
-  STORAGE: surface_tmp_arrays, interp_metric_arrays 
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Find Surface Area Element"
-
-schedule EHFinder_IntegrateArea in EHFinder_Integration after EHFinder_FindSurfaceElement
-{
-  LANG: Fortran
-  STORAGE: integrate_tmp_array 
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Calculate area integrals"
-
-schedule EHFinder_IntegrateCentroid in EHFinder_Integration after EHFinder_IntegrateArea
-{
-  LANG: Fortran
-  STORAGE: surface_tmp_arrays, integrate_tmp_array
-  TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
-} "Calculate centroid integrals"
-
-schedule EHFinder_IntegrateCircumference in EHFinder_Integration after EHFinder_IntegrateArea
-{
-  LANG: Fortran
-  STORAGE: surface_tmp_arrays, integrate_tmp_array
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Calculate circumferences"
-
-schedule EHFinder_UpdateCounter in EHFinder_Surfaces
-{
-  LANG: Fortran
-  TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Update the loop variables"
-
-schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Level_Sets
-{
-  LANG: Fortran
-  STORAGE: eh_area
-  TRIGGERS: eh_area
-} "Copy areas to output variable"
+  schedule GROUP EHFinder_Level_Sets at CCTK_ANALYSIS
+  {
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Loop over the level set functions"
+  
+  schedule EHFinder_LevelSetLoopInit in EHFinder_Level_Sets
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Initialize the loop counter over the level set functions"
+  
+  schedule GROUP EHFinder_Surfaces in EHFinder_Level_Sets while ehfinder::more_levelsets
+  {
+    STORAGE: surface_integers
+    STORAGE: surface_reals
+    STORAGE: surface_index
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Count the number of surfaces and integrate over them"
+  
+  
+  schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Initialize while loop control"
+  
+  schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces
+  {
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Counting and mark surfaces"
+  
+  schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+    SYNC: surface_index
+  } "Check if there are more surfaces"
+  
+  schedule GROUP EHFinder_MarkPoints in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points
+  {
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Marking surfaces"
+  
+  schedule EHFinder_MarkSurfaces in EHFinder_MarkPoints
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+    SYNC: surface_index
+  } "Mark points inside the current surface"
+  
+  schedule EHFinder_ApplySymSC in EHFinder_MarkPoints after EHFinder_MarkSurfaces
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Select the surface counter grid function for boundary conditions"
+  
+  schedule GROUP ApplyBCs as EHFinderSC_ApplyBSc in EHFinder_MarkPoints after EHFinder_ApplySymSC
+  {
+  } "Apply boundary conditions (symmetries)"
+  
+  schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Output info about found surfaces"
+  
+  schedule group EHFinder_Integration in EHFinder_Surfaces after EHFinder_InfoSurfaces while ehfinder::integrate_counter
+  {
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Find and integrate over surfaces"
+  
+  schedule EHFinder_FindSurface in EHFinder_Integration
+  {
+    LANG: Fortran
+    STORAGE: surface_tmp_arrays, surface_int_array
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+    SYNC: surface_arrays
+  } "Find Surface"
+  
+  schedule EHFinder_FindSurfaceElement in EHFinder_Integration after EHFinder_FindSurface
+  {
+    LANG: Fortran
+    STORAGE: surface_tmp_arrays, interp_metric_arrays 
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Find Surface Area Element"
+  
+  schedule EHFinder_IntegrateArea in EHFinder_Integration after EHFinder_FindSurfaceElement
+  {
+    LANG: Fortran
+    STORAGE: integrate_tmp_array 
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Calculate area integrals"
+  
+  schedule EHFinder_IntegrateCentroid in EHFinder_Integration after EHFinder_IntegrateArea
+  {
+    LANG: Fortran
+    STORAGE: surface_tmp_arrays, integrate_tmp_array
+    TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
+  } "Calculate centroid integrals"
+  
+  schedule EHFinder_IntegrateCircumference in EHFinder_Integration after EHFinder_IntegrateArea
+  {
+    LANG: Fortran
+    STORAGE: surface_tmp_arrays, integrate_tmp_array
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Calculate circumferences"
+  
+  schedule EHFinder_UpdateCounter in EHFinder_Surfaces
+  {
+    LANG: Fortran
+    TRIGGERS: eh_area, eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Update the loop variables"
+  
+  schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Level_Sets
+  {
+    LANG: Fortran
+    STORAGE: eh_area
+    TRIGGERS: eh_area
+  } "Copy areas to output variable"
+  
+  schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Level_Sets
+  {
+    LANG: Fortran
+    STORAGE: eh_centroid_x, eh_centroid_y, eh_centroid_z
+    TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
+  } "Copy centroids to output variable"
+  
+  schedule EHFinder_CopyCircumference at CCTK_ANALYSIS after EHFinder_Level_Sets
+  {
+    LANG: Fortran
+    STORAGE: eh_circ_eq, eh_circ_pol
+    TRIGGERS: eh_circ_eq, eh_circ_pol
+  } "Copy circumferences to output variable"
+}
 
-schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Level_Sets
+if (CCTK_EQUALS(area_calculation_method,"isosurface"))
 {
-  LANG: Fortran
-  STORAGE: eh_centroid_x, eh_centroid_y, eh_centroid_z
-  TRIGGERS: eh_centroid_x, eh_centroid_y, eh_centroid_z
-} "Copy centroids to output variable"
+  schedule GROUP EHFinder_IsoSurfaceArea at CCTK_ANALYSIS
+  {
+    TRIGGERS: eh_area
+  } "Find isosurfaces and calculate the area"
 
-schedule EHFinder_CopyCircumference at CCTK_ANALYSIS after EHFinder_Level_Sets
-{
-  LANG: Fortran
-  STORAGE: eh_circ_eq, eh_circ_pol
-  TRIGGERS: eh_circ_eq, eh_circ_pol
-} "Copy circumferences to output variable"
+  schedule EHFinder_IsoSurface in EHFinder_IsoSurfaceArea
+  {
+    LANG: Fortran
+    STORAGE: eh_area
+    TRIGGERS: eh_area
+  } "Find isosurfaces"
+}
 
 # Read in the data used in reconstructing the 4-metric if necessary
 
diff --git a/src/EHFinder_IsoSurface.F90 b/src/EHFinder_IsoSurface.F90
new file mode 100644
index 0000000..a5ab9d8
--- /dev/null
+++ b/src/EHFinder_IsoSurface.F90
@@ -0,0 +1,531 @@
+! Finding isosurface(s).
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+  use IsoSurface_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: n_coords, n_triangles
+  CCTK_REAL, dimension(:), allocatable :: x_coord, y_coord, z_coord
+  CCTK_INT, dimension(:), allocatable :: t_index
+  CCTK_REAL, parameter :: iso_value = zero
+  CCTK_INT :: l
+
+! Find index ranges for points excluding ghost zones and symmetry points for
+! the 3D grid functions.
+#include "include/physical_part.h"
+
+  call CCTK_INFO ('Entered IsoSurface')
+  
+  print*,nx,ny,nz
+  print*,size(f)
+  print*,Xf0,Xf1,Xf2,f_length
+  print*,size(x)
+!  allocate(x_coord(5),y_coord(5),z_coord(5),t_index(20))
+  do l = 1, eh_number_level_sets
+!    call triangularization ( n_coords, n_triangles, nx, ny, nz, x_coord, &
+!                             y_coord, z_coord, t_index, f(:,:,:,l), x, y, z )
+    call triangularization ( n_coords, n_triangles, nx, ny, nz, f(:,:,:,l), x, y, z )
+
+  print*,'The surface consists of ',n_triangles/3,' triangles and ',n_coords,' coords'
+  end do
+
+end subroutine EHFinder_IsoSurface
+
+
+!subroutine triangularization(l,contor,Nx,Ny,Nz,coord_x,coord_z,coord_y,triangle_index,grid_value,grid_x,grid_y,grid_z)
+subroutine triangularization(l,contor,Nx,Ny,Nz,grid_value,grid_x,grid_y,grid_z)
+
+  use IsoSurface_mod
+
+  implicit none
+
+!declaring variables that will be used in the program
+  CCTK_INT,intent(in)::Nx,Ny,Nz
+  CCTK_REAL:: u,ABx,ABy,ABz,BCx,BCy,BCz,dx
+  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_value
+  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_x
+  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_y
+  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_z
+  CCTK_INT::i,j,k,jj,kk,index,sum,nsize
+  CCTK_INT,intent(out)::l,contor
+!  CCTK_REAL,dimension(:),allocatable,intent(out):: coord_x,coord_y,coord_z
+!  CCTK_INT,dimension (:),allocatable,intent(out)::triangle_index
+!  CCTK_INT,dimension(:),allocatable::resize_triangle_index
+!  CCTK_REAL,dimension(:),allocatable::resize_coord_x,resize_coord_y,resize_coord_z
+!we store table with triangles in a 2 dimensional constant array using reshape function; Fortran stores column-like and we need row-like
+  CCTK_INT, dimension (256,16):: a
+  CCTK_INT, dimension (16,256):: b(16,256)=reshape&
+((/-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1,&
+3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1,&
+3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1,&
+3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1,&
+9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1,&
+9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1,&
+2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1,&
+8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1,&
+9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1,&
+4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1,&
+3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1,&
+1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1,&
+4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1,&
+4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1,&
+9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1,&
+5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1,&
+2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1,&
+9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1,&
+0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1,&
+2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1,&
+10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1,&
+4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1,&
+5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1,&
+5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1,&
+9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1,&
+0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1,&
+1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1,&
+10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1,&
+8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1,&
+2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1,&
+7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1,&
+9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1,&
+2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1,&
+11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1,&
+9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1,&
+5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1,&
+11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1,&
+11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1,&
+1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1,&
+9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1,&
+5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1,&
+2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1,&
+0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1,&
+5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1,&
+6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1,&
+3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1,&
+6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1,&
+5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1,&
+1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1,&
+10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1,&
+6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1,&
+8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1,&
+7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1,&
+3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1,&
+5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1,&
+0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1,&
+9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1,&
+8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1,&
+5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1,&
+0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1,&
+6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1,&
+10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1,&
+10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1,&
+8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1,&
+1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1,&
+3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1,&
+0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1,&
+10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1,&
+3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1,&
+6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1,&
+9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1,&
+8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1,&
+3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1,&
+6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1,&
+0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1,&
+10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1,&
+10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1,&
+2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1,&
+7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1,&
+7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1,&
+2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1,&
+1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1,&
+11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1,&
+8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1,&
+0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1,&
+7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1,&
+10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1,&
+2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1,&
+6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1,&
+7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1,&
+2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1,&
+1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1,&
+10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1,&
+10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1,&
+0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1,&
+7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1,&
+6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1,&
+8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1,&
+9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1,&
+6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1,&
+4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1,&
+10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1,&
+8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1,&
+0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1,&
+1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1,&
+8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1,&
+10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1,&
+4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1,&
+10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1,&
+5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1,&
+11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1,&
+9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1,&
+6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1,&
+7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1,&
+3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1,&
+7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1,&
+9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1,&
+3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1,&
+6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1,&
+9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1,&
+1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1,&
+4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1,&
+7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1,&
+6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1,&
+3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1,&
+0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1,&
+6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1,&
+0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1,&
+11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1,&
+6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1,&
+5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1,&
+9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1,&
+1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1,&
+1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1,&
+10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1,&
+0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1,&
+5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1,&
+10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1,&
+11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1,&
+9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1,&
+7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1,&
+2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1,&
+8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1,&
+9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1,&
+9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1,&
+1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1,&
+9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1,&
+9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1,&
+5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1,&
+0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1,&
+10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1,&
+2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1,&
+0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1,&
+0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1,&
+9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1,&
+5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1,&
+3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1,&
+5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1,&
+8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1,&
+0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1,&
+9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1,&
+0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1,&
+1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1,&
+3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1,&
+4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1,&
+9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1,&
+11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1,&
+11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1,&
+2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1,&
+9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1,&
+3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1,&
+1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1,&
+4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1,&
+4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1,&
+0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1,&
+3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1,&
+3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1,&
+0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1,&
+9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1,&
+1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+1, 3, 8, 9, 1, 8, -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,&
+0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
+-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/),(/16,256/))
+
+  print*,'1'
+  do j=1,256
+    do i=1,16
+      a(j,i)=b(i,j)
+    end do 
+  end do
+
+  !allocating initial space for coordinate and triangle index arrays.These arrays will be extended if needed during the execution of the program using allocate and deallocate functions.
+  nsize=5
+  allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize))
+
+  !marching through the grid and computing "sum" for each individual cube. "l" is the contor for the vectors that keeps the coordinates of intersection points. we will also produce a vector called triangle_index that tells which points build which triangles.
+  l=0
+  contor=0
+  do k=1,Nz-1
+    do j=1,Ny-1
+      do i=1,Nx-1                                 
+        sum=1
+          if (grid_value(i,j,k)<0) then 
+            sum=sum+2**0
+          end if
+          if (grid_value(i+1,j,k)<0) then 
+            sum=sum+2**3
+          end if
+          if (grid_value(i+1,j+1,k)<0) then 
+            sum=sum+2**2
+          end if
+          if (grid_value(i,j+1,k)<0) then 
+            sum=sum+2**1
+          end if
+          if (grid_value(i+1,j+1,k+1)<0) then 
+            sum=sum+2**6
+          end if
+          if (grid_value(i+1,j,k+1)<0) then 
+            sum=sum+2**7
+          end if
+          if (grid_value(i,j,k+1)<0) then 
+            sum=sum+2**4
+          end if
+          if (grid_value(i,j+1,k+1)<0) then 
+            sum=sum+2**5
+          end if
+
+          if (a(sum,1)>-1) then
+            do jj=1,15,3
+              if (a(sum,jj)>-1) then
+                do kk=jj,jj+2
+
+                  if (a(sum,kk)==0) then 
+                    u=abs(grid_value(i,j,k))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i,j,k))+abs(grid_value(i,j+1,k)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)
+                    coord_y(l)=grid_y(i,j,k)+u
+                    coord_z(l)=grid_z(i,j,k)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==1) then 
+                    u=abs(grid_value(i,j+1,k))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j+1,k))+abs(grid_value(i+1,j+1,k)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)+u
+                    coord_y(l)=grid_y(i,j+1,k)
+                    coord_z(l)=grid_z(i,j,k)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==2) then 
+                    u=abs(grid_value(i+1,j,k))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i+1,j+1,k))+abs(grid_value(i+1,j,k)))
+                    l=l+1
+                    coord_x(l)=grid_x(i+1,j,k)
+                    coord_y(l)=grid_y(i,j,k)+u
+                    coord_z(l)=grid_z(i,j,k)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==3) then
+                    u=abs(grid_value(i,j,k))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i+1,j,k))+abs(grid_value(i,j,k)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)+u
+                    coord_y(l)=grid_y(i,j,k)
+                    coord_z(l)=grid_z(i,j,k)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==4) then 
+                    u=abs(grid_value(i,j,k+1))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i,j,k+1))+abs(grid_value(i,j+1,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)
+                    coord_y(l)=grid_y(i,j,k)+u
+                    coord_z(l)=grid_z(i,j,k+1)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==5) then 
+                    u=abs(grid_value(i,j+1,k+1))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j+1,k+1))+abs(grid_value(i+1,j+1,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)+u
+                    coord_y(l)=grid_y(i,j+1,k)
+                    coord_z(l)=grid_z(i,j,k+1)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==6) then 
+                    u=abs(grid_value(i+1,j,k+1))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i+1,j+1,k+1))+abs(grid_value(i+1,j,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i+1,j,k)
+                    coord_y(l)=grid_y(i,j,k)+u
+                    coord_z(l)=grid_z(i,j,k+1)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+ 
+                  if (a(sum,kk)==7) then
+                    u=abs(grid_value(i,j,k+1))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j,k+1))+abs(grid_value(i+1,j,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)+u
+                    coord_y(l)=grid_y(i,j,k)
+                    coord_z(l)=grid_z(i,j,k+1)
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==8) then 
+                    u=abs(grid_value(i,j,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i,j,k))+abs(grid_value(i,j,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)
+                    coord_y(l)=grid_y(i,j,k)
+                    coord_z(l)=grid_z(i,j,k)+u
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==9) then 
+                    u=abs(grid_value(i,j+1,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i,j+1,k))+abs(grid_value(i,j+1,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i,j,k)
+                    coord_y(l)=grid_y(i,j+1,k)
+                    coord_z(l)=grid_z(i,j,k)+u
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==10) then 
+                    u=abs(grid_value(i+1,j+1,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i+1,j+1,k))+abs(grid_value(i+1,j+1,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i+1,j,k)
+                    coord_y(l)=grid_y(i,j+1,k)
+                    coord_z(l)=grid_z(i,j,k)+u
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+
+                  if (a(sum,kk)==11) then 
+                    u=abs(grid_value(i+1,j,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i+1,j,k))+abs(grid_value(i+1,j,k+1)))
+                    l=l+1
+                    coord_x(l)=grid_x(i+1,j,k)
+                    coord_y(l)=grid_y(i,j,k)
+                    coord_z(l)=grid_z(i,j,k)+u
+                    contor=contor+1
+                    triangle_index(contor)=l-1
+                  end if
+                end do
+
+                  if (l>3) then
+                    ABx=coord_x(l)-coord_x(l-1)
+                    BCx=coord_x(l-1)-coord_x(l-2)
+                    ABy=coord_y(l)-coord_y(l-1)
+                    BCy=coord_y(l-1)-coord_y(l-2)
+                    ABz=coord_z(l)-coord_z(l-1)
+                    BCz=coord_z(l-1)-coord_z(l-2)
+                  if (((ABx*BCy)==(ABy*BCx)).and.((ABx*BCz)==(ABz*BCx)).and.((ABy*BCz)==(ABz*BCy))) then
+                    l=l-3 
+                    contor=contor-3
+                  end if
+                end if
+              end if
+
+              if ((nsize-l)<10) then
+                allocate(resize_coord_x(nsize),resize_coord_y(nsize),resize_coord_z(nsize),resize_triangle_index(nsize))
+                resize_coord_x(1:l)=coord_x(1:l)
+                resize_coord_y(1:l)=coord_y(1:l)
+                resize_coord_z(1:l)=coord_z(1:l)
+                resize_triangle_index(1:l)=triangle_index(1:l)
+                deallocate(coord_x,coord_y,coord_z,triangle_index)
+                nsize=nsize*2
+                allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize))
+                coord_x(1:l)=resize_coord_x(1:l)
+                coord_y(1:l)=resize_coord_y(1:l)
+                coord_z(1:l)=resize_coord_z(1:l)
+                triangle_index(1:l)=resize_triangle_index(1:l)
+                deallocate(resize_coord_x,resize_coord_y,resize_coord_z,resize_triangle_index)
+            end if
+          end do
+        end if
+      end do
+    end do
+  end do 
+end subroutine triangularization
diff --git a/src/IsoSurface_mod.F90 b/src/IsoSurface_mod.F90
new file mode 100644
index 0000000..b2c29c1
--- /dev/null
+++ b/src/IsoSurface_mod.F90
@@ -0,0 +1,13 @@
+! Module with various common variable.
+! $Header$
+#include "cctk.h"
+
+module IsoSurface_mod
+
+  use EHFinder_Constants
+
+  CCTK_REAL,dimension(:),allocatable:: coord_x,coord_y,coord_z
+  CCTK_INT,dimension (:),allocatable::triangle_index
+  CCTK_INT,dimension(:),allocatable::resize_triangle_index
+  CCTK_REAL,dimension(:),allocatable::resize_coord_x,resize_coord_y,resize_coord_z
+end module IsoSurface_mod
diff --git a/src/make.code.defn b/src/make.code.defn
index d68fe1a..6a92f55 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -3,6 +3,7 @@
 
 # Source files in this directory
 SRCS = EHFinder_mod.F90 \
+       IsoSurface_mod.F90 \
        EHFinder_ParamCheck.F90 \
        EHFinder_Init.F90 \
        MoL_Init.F90 \
@@ -13,6 +14,7 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_Check.F90 \
        EHFinder_Integrate.F90 \
        EHFinder_FindSurface.F90 \
+       EHFinder_IsoSurface.F90 \
        EHFinder_Generator_Sources.F90 \
        EHFinder_Generator_Sources2.F90 \
        EHFinder_ReInitialize.F90 \
diff --git a/src/make.code.deps b/src/make.code.deps
index 383006f..0b34034 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -2,6 +2,7 @@
 # $Header$
 
 EHFinder_mod.F90.o:			EHFinder_Constants.F90.o
+EHFinder_IsoSurface_mod.F90.o:		EHFinder_Constants.F90.o
 EHFinder_Init.F90.o:			EHFinder_mod.F90.o
 EHFinder_Sources.F90.o:			EHFinder_mod.F90.o
 EHFinder_Generator_Sources.F90.o:	EHFinder_mod.F90.o
@@ -14,3 +15,4 @@ EHFinder_ReadData.F90.o:		EHFinder_mod.F90.o
 EHFinder_Check.F90.o:			EHFinder_mod.F90.o
 EHFinder_Integrate.F90.o:		EHFinder_mod.F90.o
 EHFinder_FindSurface.F90.o:		EHFinder_mod.F90.o
+EHFinder_IsoSurface.F90.o:		EHFinder_mod.F90.o IsoSurface_mod.F90.o