Topological partition almost works

eplf.config

@@ -51,6 +51,7 @@ grid_data
  density             real (grid_num_x,grid_num_y,grid_num_z)
  density_grad        real (grid_num_x,grid_num_y,grid_num_z,4)
  density_lapl        real (grid_num_x,grid_num_y,grid_num_z)
+ elf_partition       real (grid_num_x,grid_num_y,grid_num_z)
 
 compute
  nproc               integer

scripts/eplf_edit.py

@@ -214,7 +214,7 @@ def write_main_file(file,inp):
     x = ' '
     exec "if inp.%s: x = 'X'"%(p)
     print >>file, "(%s) %s"%(x,p[8:])
-  print >>file, ""
+
 
 ## Execute temporary input file
 ###############################
@@ -262,21 +262,17 @@ def edit_temp_file(input_file,write_file,read_file,saved_file=None):
 
 
 def build_script(inp):
-  grid = "${EPLF_PATH}/bin/to_cube %s"%(inp.name)
   run_script = "run_"+inp.name
   file = open(run_script,'w')
   buffer = """
   #!/bin/bash
 source ${HOME}/.eplfrc
 $RUNCOMMAND
-$RUNGRID
 """.lstrip()
   command = "${EPLF_PATH}/bin/eplf %s"%(inp.name)
   if inp.has_mpi:
     command = "${MPIRUN} -np %d "%(inp.nproc)+command
-    grid = "${MPIRUN} -np 1 "+command
   buffer = buffer.replace("$RUNCOMMAND",command)
-  buffer = buffer.replace("$RUNGRID",grid)
   print >>file, buffer
   file.close()
   os.chmod(run_script,0744)

scripts/input_wrapper.py

@@ -13,6 +13,7 @@ point_num step_size origin opposite nproc
 compute_elf compute_eplf compute_density
 compute_elf_grad compute_eplf_grad compute_density_grad
 compute_elf_lapl compute_eplf_lapl compute_density_lapl
+compute_elf_partition
 """.split()
 rw_data_full = rw_data + geom_data
 
@@ -26,6 +27,9 @@ DEFAULT_COMPUTE_DENSITY_GRAD  = False
 DEFAULT_COMPUTE_EPLF_LAPL     = False
 DEFAULT_COMPUTE_ELF_LAPL      = False
 DEFAULT_COMPUTE_DENSITY_LAPL  = False
+DEFAULT_COMPUTE_EPLF_PART     = False
+DEFAULT_COMPUTE_ELF_PART      = False
+DEFAULT_COMPUTE_DENSITY_PART  = False
 DEFAULT_POINT_NUM             = [80,80,80]
 
 ######################################################################
@@ -322,6 +326,17 @@ class InputFile(object):
       raise InputFileExn("Wrong type")
     ezfio.set_compute_density_lapl(value)
 
+  def get_compute_elf_partition(self):
+    if ezfio.has_compute_elf_partition():
+      return ezfio.get_compute_elf_partition()
+    else:
+      return DEFAULT_COMPUTE_ELF_PART
+
+  def set_compute_elf_partition(self,value):
+    if not isinstance(value,bool):
+      raise InputFileExn("Wrong type")
+    ezfio.set_compute_elf_partition(value)
+
   # Build the corresponding properties
   for i in ro_data:
     exec "%s = property(fget=get_%s,fset=None)"%(i,i)

src/ezfio_interface.irp.f

@@ -40,6 +40,7 @@ data = [ \
 ("grid_data_density"               , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ("grid_data_density_grad"          , "real"          , "(grid_x_num,grid_y_num,grid_z_num,4)" ),
 ("grid_data_density_lapl"          , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
+("grid_data_elf_partition"         , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ]
 
 data_no_set = [\

src/grid.irp.f

@@ -249,13 +249,149 @@ SUBST [ X ]
 
 END_TEMPLATE
 
-IRP_IF caca
-BEGIN_PROVIDER [ integer, grid_partition_$X, (grid_x_num,grid_y_num,grid_z_num) ]
+BEGIN_TEMPLATE
+
+BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ]
  implicit none
  BEGIN_DOC  
 ! Create the topological partition of $X
  END_DOC
  
-END_PROVIDER
-IRP_ENDIF
+ integer :: ix, iy, iz
+ real, parameter :: UNDEFINED=123456789.123456789
+
+ if (mpi_master) then
+
+   grid_$X_partition(1,1,1) = UNDEFINED
+   call get_grid_data_$X_partition(grid_$X_partition)
+   if (grid_$X_partition(1,1,1) /= UNDEFINED) then
+     return
+   endif
+
+   do iz=1,grid_z_num 
+    do iy=1,grid_y_num 
+     do ix=1,grid_x_num 
+      grid_$X_partition(ix,iy,iz) = 0.
+     enddo
+    enddo
+   enddo
+
+   do iz=1,grid_z_num 
+    do iy=1,grid_y_num 
+     do ix=1,grid_x_num 
+      call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition)
+     enddo
+    enddo
+   enddo
+
+   call set_grid_data_$X_partition(grid_$X_partition)
+ endif
+
+END_PROVIDER
+
+SUBST [ X ]
+ elf;;
+
+END_TEMPLATE
+
+BEGIN_PROVIDER [ real, current_partition_index ]
+ implicit none
+ BEGIN_DOC  
+! Current number for the space partitioning
+ END_DOC
+ current_partition_index = 1.
+END_PROVIDER
+
+recursive subroutine find_attribution(ix,iy,iz,g,a)
+ implicit none
+
+ integer, intent(in) :: ix, iy, iz
+ real, intent(in) :: g(grid_x_num,grid_y_num,grid_z_num)
+ real, intent(inout) :: a(grid_x_num,grid_y_num,grid_z_num)
+
+ real :: buffer(-1:1,-1:1,-1:1)
+
+ integer :: i,j,k
+
+ if (a(ix,iy,iz) /= 0) then
+  return
+ else if (ix == 1) then
+  call find_attribution(2,iy,iz,g,a)
+  a(ix,iy,iz) = a(2,iy,iz)
+ else if (ix == grid_x_num) then
+  call find_attribution(ix-1,iy,iz,g,a)
+  a(ix,iy,iz) = a(ix-1,iy,iz)
+ else if (iy == 1) then
+  call find_attribution(ix,2,iz,g,a)
+  a(ix,iy,iz) = a(ix,2,iz)
+ else if (iy == grid_y_num) then
+  call find_attribution(ix,iy-1,iz,g,a)
+  a(ix,iy,iz) = a(ix,iy-1,iz)
+ else if (iz == 1) then
+  call find_attribution(ix,iy,2,g,a)
+  a(ix,iy,iz) = a(ix,iy,2)
+ else if (iz == grid_z_num) then
+  call find_attribution(ix,iy,iz-1,g,a)
+  a(ix,iy,iz) = a(ix,iy,iz-1)
+
+ else
+
+   ! Copy the environment of the current point
+   do i=-1,1
+    do j=-1,1
+     buffer(-1,j,i) = g(ix-1, iy+j, iz+i)
+     buffer( 0,j,i) = g(ix  , iy+j, iz+i)
+     buffer( 1,j,i) = g(ix+1, iy+j, iz+i)
+    enddo
+   enddo 
+
+   ! Avoid strict equality between points
+   buffer(0,0,0) *= 1.000001
+
+   ! Find the index of the maximum value
+   real :: m,mm
+   integer :: im(3)
+   mm = buffer(-1,-1,-1)
+   m  = buffer(-1,-1,-1)
+   do i=-1,1
+    do j=-1,1
+     do k=-1,1
+      if (buffer(k,j,i) > m) then
+       m = buffer(k,j,i)
+       im(1) = k
+       im(2) = j
+       im(3) = i
+      endif
+      mm = min(buffer(k,j,i),mm)
+     enddo
+    enddo
+   enddo 
+
+
+   if ( (im(1) == 0).and.(im(2) == 0).and.(im(3) == 0) ) then
+     ! If the maximum is at the center, a new maximum is found
+     if ( (m-mm)/(m+mm+1.e-12) > 0.01) then
+       do i=-1,1
+        do j=-1,1
+         do k=-1,1
+          if (g(ix+k,iy+j,iz+i) == g(ix,iy,iz)) then
+            a(ix+k,iy+j,iz+i) = current_partition_index
+          endif
+         enddo
+        enddo
+       enddo
+       current_partition_index += 1.
+       TOUCH current_partition_index 
+     endif
+   else
+     ! Otherwise, get the partition index of the max value
+     call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a)
+     if (a(ix,iy,iz) == 0) then
+       a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3))
+     endif
+   endif
+
+ endif
+
+end
 

src/main.irp.f

@@ -17,6 +17,10 @@ program eplf
    PROVIDE grid_elf
    FREE grid_elf
  endif
+ if (comp_elf_partition) then
+   PROVIDE grid_elf_partition
+   FREE grid_elf_partition
+ endif
  if (comp_elf_grad.or.comp_elf_lapl) then
    PROVIDE grid_elf_grad
    PROVIDE grid_elf_lapl

src/to_cube.irp.f

@@ -52,6 +52,7 @@ program to_cube
    density_grad;4;,4;;
    elf_grad;4;,4;;
    eplf_grad;4;,4;;
+   elf_partition;3; ;;
   END_TEMPLATE
   10   format (2X,I3,3(2X,F10.6))
   11   format (2X,I3,4(2X,F10.6))