From f1cfca91b8d2e196c2c5d96ce9832446aed27aa4 Mon Sep 17 00:00:00 2001 From: diener Date: Fri, 20 Aug 2004 19:02:20 +0000 Subject: 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 --- param.ccl | 6 + schedule.ccl | 306 +++++++++++++------------ src/EHFinder_IsoSurface.F90 | 531 ++++++++++++++++++++++++++++++++++++++++++++ src/IsoSurface_mod.F90 | 13 ++ src/make.code.defn | 2 + src/make.code.deps | 2 + 6 files changed, 716 insertions(+), 144 deletions(-) create mode 100644 src/EHFinder_IsoSurface.F90 create mode 100644 src/IsoSurface_mod.F90 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 -- cgit v1.2.3