Initial commit from Polytechnique Montreal

1 files changed, 679 insertions, 0 deletions

diff --git a/Dragon/src/g2s_celluleBase.f90 b/Dragon/src/g2s_celluleBase.f90
new file mode 100644
index 0000000..e5b89d0
--- /dev/null
+++ b/Dragon/src/g2s_celluleBase.f90
@@ -0,0 +1,679 @@
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Creation of an array of type(t_celluleBase) structures.
+!
+!Copyright:
+! Copyright (C) 2001 Ecole Polytechnique de Montreal.
+! This library is free software; you can redistribute it and/or
+! modify it under the terms of the GNU Lesser General Public
+! License as published by the Free Software Foundation; either
+! version 2.1 of the License, or (at your option) any later version.
+!
+!Author(s):
+! G. Civario (CS-SI)
+!
+!Comments:
+! Il s'agit de toutes les cellules terminales de la geometrie (=celles qui
+! ne comportent pas de sous-cellules).
+! La creation se fait a l'aide de la routine recursive buildCellsBase.
+! Elle explore toute l'arboressence de l'objet PyLCM donne en entree, et appelle
+!  a chaque niveau la routine writeCellBase. 
+! Celle-ci teste si la cellule est terminale, et cree une entree dans le tableau
+! des cellules de base si c'est le cas.
+! La cellule creee est ensuite completee si besoin, et sa coherence est
+! verifiee
+! \\\\
+! variable globale
+!  - tabCelluleBase : tableau des cellules de base
+! \\\\
+! fonctions du module
+!  - initializeTabCelluleBase : mise a zero du tableau
+!  - destoyTabCelluleBase : liberation de la memoire
+!  - decale : decale les valeurs d'un tableau pour les faire demarer a 0
+!  - createCB : constructeur d'une cellule de base
+!  - destroyCB : destructeur d'une cellule de base
+!  - createCluster : constructeur d'un cluster
+!  - sortClusterTab : trie un tableau de clusters
+!  - verrifieCB : verification de la coherence d'une cellule de base
+!  - exploiteSplit : exploitation de la donnee split dans la cellule
+!  - buildCellsBase : fonction recursive de creation du tableau de cellules
+!  - writeCellBase : creation d'une cellule
+!
+!-----------------------------------------------------------------------
+!
+module celluleBase
+
+  use constType
+  use GANLIB
+  use segArc, only : alloc_ok
+
+  implicit none
+
+  !cluster
+  type t_cluster
+     character*12                          :: name        !nom
+     integer                               :: nbrPin      !nombre de crayons
+     double precision                      :: radiusOfPin !rayon de la couronne
+     double precision                      :: angleOfPin  !angle du 1er crayon
+     double precision,dimension(:),pointer :: radius      !rayons des anneaux
+     integer,dimension(:),pointer          :: mix         !milieux des crayons
+  end type t_cluster
+
+  !cellule generique de plus bas niveau hierarchique dans les gigognes
+  type t_celluleBase
+     character*12                          :: name      !nom
+     integer,dimension(40)                 :: sv        !state vector
+     integer,dimension(:),allocatable          :: mix       !milieux
+     integer,dimension(:),allocatable          :: merge     !regroupements
+     double precision,dimension(:),allocatable :: radius    !rayons
+     double precision,dimension(3)         :: offcenter !x , y et z
+     double precision,dimension(:),allocatable :: meshx     !en commancant a 0.0
+     double precision,dimension(:),allocatable :: meshy     !en commancant a 0.0
+     double precision                      :: side      !pour tri et hex
+     integer,dimension(:),allocatable          :: splitr    !>0  rayon; <0  surface
+     integer,dimension(:),allocatable          :: splitx    !>0
+     integer,dimension(:),allocatable          :: splity    !>0
+     type(t_cluster),dimension(:),pointer  :: cluster   !les clusters
+     !tableau de bool qui donne la presence ou non de chacun des 10 champs
+     !avec en plus l'indice 0 disant si la cellule est viable (complete)
+     logical,dimension(0:12)               :: ok
+  end type t_celluleBase
+
+  !parametres pour le vecteur ok
+  integer,parameter :: n_sv=1 , n_mix=2 , n_radius=3 , n_offcenter=4 , &
+                     & n_meshx=5 , n_meshy=6 , n_side=7 , n_splitr=8 , &
+                     & n_splitx=9 , n_splity=10 , n_cluster=11 , n_merge=12 , &
+                     & n_tot=0
+
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!  declaration d'une variable tableau globale
+  !!
+  type(t_celluleBase),dimension(:),allocatable :: tabCelluleBase
+  !!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+contains
+
+  subroutine initializeTabCelluleBase(dimTabCelluleBase)
+    integer,intent(in) :: dimTabCelluleBase
+
+    allocate(tabCelluleBase(dimTabCelluleBase),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: initializeTabCelluleBase(1) => allocation pb")
+  end subroutine initializeTabCelluleBase
+  
+  subroutine destroyTabCelluleBase(szB)
+    integer,intent(in) :: szB
+    integer :: i
+
+    do i = 1,szB
+       call destroyCB(tabCelluleBase(i))
+    end do
+    deallocate(tabCelluleBase)
+   end subroutine destroyTabCelluleBase
+
+  !fait demarer un tableau de reels a 0.0
+  subroutine decale(tab)
+!    double precision,dimension(:),pointer :: tab
+    double precision,dimension(:) :: tab
+    integer :: i, lg
+    
+    lg = size(tab)
+    do i = lg,1,-1
+       tab(i) = tab(i) - tab(1)
+    end do
+  end subroutine decale
+  
+  !remplisage d'une cellule de base
+  subroutine createCB(cell,name,ip)
+    type(t_celluleBase),intent(out)  :: cell
+    character*12,intent(in)          :: name
+    type(c_ptr),intent(in)           :: ip
+
+    integer :: lg,typ,i,lgm
+    real,dimension(:),allocatable         :: tmpTabReal
+    character*12,dimension(:),allocatable :: clusterName
+
+    cell%name = name
+    
+    call LCMGET(ip,'STATE-VECTOR',cell%sv)
+    cell%ok(n_sv) = .true.
+
+    call LCMLEN(ip,'MIX         ',lgm,typ)
+    cell%ok(n_mix) = (lgm/=0)
+    if(cell%ok(n_mix)) then
+       allocate(cell%mix(lgm),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(1) => allocation pb")
+       call LCMGET(ip,'MIX         ',cell%mix)
+    else
+       call XABORT("G2S: no mix in the cellule " // name)
+    end if
+
+    call LCMLEN(ip,'MERGE       ',lg,typ)
+    cell%ok(n_merge) = (lg/=0)
+    if(cell%ok(n_merge)) then
+       if(lg/=lgm) call XABORT("G2S: bad dimension for merge in the &
+            &cellule " // name)
+       allocate(cell%merge(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(2) => allocation pb")
+       call LCMGET(ip,'MERGE       ',cell%merge)
+    else
+       allocate(cell%merge(lgm),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(3) => allocation pb")
+       cell%merge(:lgm) = (/(i,i=1,lgm)/)
+    end if
+
+    call LCMLEN(ip,'RADIUS      ',lg,typ)
+    cell%ok(n_radius) = (lg/=0)
+    if(cell%ok(n_radius)) then
+       allocate(cell%radius(lg),tmpTabReal(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(4) => allocation pb")
+       call LCMGET(ip,'RADIUS      ',tmpTabReal)
+       cell%radius(:lg)=tmpTabReal(:lg)
+       deallocate(tmpTabReal)
+    else if(cell%sv(1)==G_Carcel.or.cell%sv(1)==G_Hexcel) then
+       !cas carcel 0 ou hexcel 0
+       allocate(cell%radius(1))
+       cell%radius(1)=0.d0
+    else 
+       !nullify(cell%radius)
+    end if
+    
+    call LCMLEN(ip,'OFFCENTER   ',lg,typ)
+    cell%ok(n_offcenter) = (lg/=0)
+    if(cell%ok(n_offcenter)) then
+       allocate(tmpTabReal(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(5) => allocation pb")
+       call LCMGET(ip,'OFFCENTER   ',tmpTabReal)
+       cell%offcenter=tmpTabReal
+       deallocate(tmpTabReal)
+    else
+       cell%offcenter=0.d0
+    end if
+    
+    call LCMLEN(ip,'MESHX       ',lg,typ)
+    cell%ok(n_meshx) = (lg/=0)
+    if(cell%ok(n_meshx)) then
+       allocate(cell%meshx(lg),tmpTabReal(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(6) => allocation pb")
+       call LCMGET(ip,'MESHX       ',tmpTabReal)
+       cell%meshx(:lg)=tmpTabReal(:lg)
+       deallocate(tmpTabReal)
+       call decale(cell%meshx) !decalage pour demarer a 0.0
+    else
+       !nullify(cell%meshx)
+    end if
+    
+    call LCMLEN(ip,'MESHY       ',lg,typ)
+    cell%ok(n_meshy) = (lg/=0)
+    if(cell%ok(n_meshy)) then
+       allocate(cell%meshy(lg),tmpTabReal(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(7) => allocation pb")
+       call LCMGET(ip,'MESHY       ',tmpTabReal)
+       cell%meshy(:lg)=tmpTabReal(:lg)
+       deallocate(tmpTabReal)
+       call decale(cell%meshy) !decalage pour demarer a 0.0
+    else
+       !nullify(cell%meshy)
+    end if
+
+    call LCMLEN(ip,'SIDE        ',lg,typ)
+    cell%ok(n_side) = (lg/=0)
+    if(cell%ok(n_side)) then
+       allocate(tmpTabReal(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(8) => allocation pb")
+       call LCMGET(ip,'SIDE        ',tmpTabReal)
+       cell%side=tmpTabReal(1)
+       deallocate(tmpTabReal)
+    else
+       cell%side=0.d0
+    end if
+
+    call LCMLEN(ip,'SPLITR      ',lg,typ)
+    cell%ok(n_splitr) = (lg/=0)
+    if(cell%ok(n_splitr)) then
+       allocate(cell%splitr(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(9) => allocation pb")
+       call LCMGET(ip,'SPLITR      ',cell%splitr)
+    else
+       !nullify(cell%splitr)
+    end if
+
+    call LCMLEN(ip,'SPLITX      ',lg,typ)
+    cell%ok(n_splitx) = (lg/=0)
+    if(cell%ok(n_splitx)) then
+       allocate(cell%splitx(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(10) => allocation pb")
+       call LCMGET(ip,'SPLITX      ',cell%splitx)
+    else if(cell%sv(1)==G_Car2d.or.cell%sv(1)==G_Carcel) then
+       allocate(cell%splitx(cell%sv(3)),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(11) => allocation pb")
+       cell%splitx = 1
+    else
+       allocate(cell%splitx(1))
+       cell%splitx(1) = 1
+    end if
+
+    call LCMLEN(ip,'SPLITY      ',lg,typ)
+    cell%ok(n_splity) = (lg/=0)
+    if(cell%ok(n_splity)) then
+       allocate(cell%splity(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(12) => allocation pb")
+       call LCMGET(ip,'SPLITY      ',cell%splity)
+    else if(cell%sv(1)==G_Car2d.or.cell%sv(1)==G_Carcel) then
+       allocate(cell%splity(cell%sv(4)),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(13) => allocation pb")
+       cell%splity = 1
+    else
+       allocate(cell%splity(1))
+       cell%splity(1) = 1
+    end if
+
+    call LCMLEN(ip,'CLUSTER     ',lg,typ)
+    cell%ok(n_cluster) = (lg/=0)
+    if(cell%ok(n_cluster)) then
+       lg = lg/3
+       allocate(cell%cluster(lg),clusterName(lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: createCB(14) => allocation pb")
+       call LCMGTC(ip,'CLUSTER     ',12,lg,clusterName)
+       do i = 1,lg
+          cell%cluster(i) = createCluster(ip,clusterName(i))
+       end do
+       !on trie le tableau des cluster par rayon croissant
+       call sortClusterTab(cell%cluster)
+       deallocate(clusterName)
+    else
+       nullify(cell%cluster)
+    end if
+
+    call verrifieCB(cell)
+    call exploiteSplit(cell)
+  end subroutine createCB
+
+  subroutine destroyCB(cell)
+    type(t_celluleBase),intent(inout) :: cell
+    integer :: i
+
+    if(allocated(cell%mix)) deallocate(cell%mix)
+    if(allocated(cell%merge)) deallocate(cell%merge)
+    if(allocated(cell%radius)) deallocate(cell%radius)
+    if(allocated(cell%meshx)) deallocate(cell%meshx)
+    if(allocated(cell%meshy)) deallocate(cell%meshy)
+    if(allocated(cell%splitr)) deallocate(cell%splitr)
+    if(allocated(cell%splitx)) deallocate(cell%splitx)
+    if(allocated(cell%splity)) deallocate(cell%splity)
+    if(associated(cell%cluster)) then
+       do i = 1,size(cell%cluster)
+          deallocate(cell%cluster(i)%mix,cell%cluster(i)%radius)
+          !nullify(cluster%mix,cell%cluster(i)%radius)
+       end do
+       deallocate(cell%cluster)
+    end if
+  end subroutine destroyCB
+
+  function createCluster(cellBIp,clusterName)
+    type(c_ptr),intent(in)      :: cellBIp
+    character*12,intent(in) :: clusterName
+    type(t_cluster)         :: createCluster
+
+    type(c_ptr) :: ip
+    integer :: lg,lgx,lgy,ty
+    real    :: rpin,apin,cx,cy
+    real,dimension(:),allocatable :: radius
+    real, parameter :: pi = 3.141592653589793
+
+    ip = cellBIp
+    call LCMSIX(ip,clusterName,1)
+
+    createCluster%name = clusterName
+    
+    call LCMLEN(ip,'CPINX       ',lgx,ty)
+    call LCMLEN(ip,'CPINY       ',lgy,ty)
+    if ((lgx==1).and.(lgy==1)) then
+
+      createCluster%nbrPin = 1
+      call LCMGET(ip,'CPINX       ',cx)
+      call LCMGET(ip,'CPINY       ',cy)
+      rpin = sqrt(cx**2+cy**2)
+      if((cy == 0.0).and.(cy == 0.0)) then
+        rpin = 0.0
+        apin = 0.0
+      else if(cy >= 0.0) then
+        apin = acos(cx/rpin)
+      else
+        apin = 2.0*pi-acos(cx/rpin)
+      endif
+
+    else
+
+      call LCMGET(ip,'NPIN        ',createCluster%nbrPin)
+      call LCMGET(ip,'RPIN        ',rpin)
+      call LCMGET(ip,'APIN        ',apin)
+
+    endif
+    createCluster%radiusOfPin = rpin
+    createCluster%angleOfPin = apin
+
+    call LCMLEN(ip,'MIX         ',lg,ty)
+    allocate(radius(lg+1),createCluster%radius(lg+1),createCluster%mix(lg),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: createCluster => allocation pb")
+    call LCMGET(ip,'RADIUS      ',radius)
+    createCluster%radius(:lg+1) = radius(:lg+1)
+    deallocate(radius)
+
+    call LCMGET(ip,'MIX         ',createCluster%mix)
+  end function createCluster
+
+  subroutine sortClusterTab(clusterPrt)
+    type(t_cluster),dimension(:),pointer :: clusterPrt
+    !trie le tableau des clusters par rayon croissant des couronnes
+    type(t_cluster) :: tmpCluster
+    integer         :: sz,indMax,i
+    logical         :: trie
+ 
+    if(.not. associated(clusterPrt)) return !pas de cluster
+    sz = size(clusterPrt)
+    do indMax = sz,2,-1
+       trie = .true.
+       do i = 1,indMax-1
+          if(clusterPrt(i+1)%radiusOfPin<clusterPrt(i)%radiusOfPin) then
+             tmpCluster=clusterPrt(i+1)
+             clusterPrt(i+1)=clusterPrt(i)
+             clusterPrt(i)=tmpCluster
+             trie = .false.
+          end if
+       end do
+       if(trie) exit
+    end do
+  end subroutine sortClusterTab
+
+  !verification de la validite d'une cellule en fonction de son vecteur ok
+  !et mise a jour du ok(0)
+  subroutine verrifieCB(cell)
+    type(t_celluleBase),intent(inout) :: cell
+
+    integer :: typCell,sectori,sectorj,nx,ny
+
+    typCell = cell%sv(1)
+    sectori = cell%sv(14)
+    sectorj = cell%sv(15)
+    select case(typCell)
+    case(G_Virtual)
+       call XABORT("G2S: Type of geometry VIRTUAL not supported")
+    case(G_Car1d)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_meshx) &
+            & .and. .not. &
+            & ( cell%ok(n_radius) .or. cell%ok(n_offcenter) &
+            & .or. cell%ok(n_meshy) .or. cell%ok(n_side) &
+            & .or. cell%ok(n_splitr) .or. cell%ok(n_splity) )
+    case(G_Tube)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_radius) &
+            & .and. .not. &
+            & ( cell%ok(n_offcenter) .or. cell%ok(n_meshx) &
+            & .or. cell%ok(n_meshy) .or. cell%ok(n_side) &
+            & .or. cell%ok(n_splitx) .or. cell%ok(n_splity) )
+       if(sectori/=S_not) call XABORT("G2S: SECT not allowed for TUBE")
+    case(G_Car2d)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_meshx) .and. cell%ok(n_meshy) &
+            & .and. .not. &
+            & ( cell%ok(n_radius) .or. cell%ok(n_offcenter) &
+            & .or. cell%ok(n_side) .or. cell%ok(n_splitr) )
+       if(sectori/=S_not) call XABORT("G2S: SECT not allowed for CAR2D")
+    case(G_Hex)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_side) &
+            & .and. .not. &
+            & ( cell%ok(n_meshx) .or. cell%ok(n_meshy) &
+            & .or. cell%ok(n_splitr) .or. cell%ok(n_splitx) &
+            & .or. cell%ok(n_splity) .or. cell%ok(n_radius) )
+       if((sectori/=S_not).and.(size(cell%mix)/=6)) &
+            call XABORT("G2S: wrong len for MIX with SECT")
+    case(G_Tri)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_side) &
+            & .and. .not. cell%ok(n_radius)
+       if(sectori/=S_not) call XABORT("G2S: SECT not allowed for TRI2D")
+    case(G_Carcel)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_radius) .and. cell%ok(n_meshx) &
+            & .and. cell%ok(n_meshy) &
+            & .and. .not. &
+            & ( cell%ok(n_side) )
+       select case(sectori)
+       case(S_not)
+          nx = size(cell%meshx)-1
+          ny = size(cell%meshy)-1
+          if(size(cell%mix)/=nx*ny*size(cell%radius)) &
+               call XABORT("G2S: wrong len for MIX")
+       case(S_X_tot,S_T_tot)
+          if(size(cell%mix)/=4*size(cell%radius)-3*sectorj) &
+               call XABORT("G2S: wrong len for MIX with SECT")
+       case(S_TX_tot,S_TXS_tot)
+          if(size(cell%mix)/=8*size(cell%radius)-7*sectorj) &
+               call XABORT("G2S: wrong len for MIX with SECT")
+       case(S_WM_tot)
+          if(size(cell%mix)/=8*size(cell%radius)-7*sectorj+4) &
+               call XABORT("G2S: wrong len for MIX with SECT")
+       end select
+    case(G_Hexcel)
+       cell%ok(n_tot) = cell%ok(n_sv) .and. cell%ok(n_mix) &
+            & .and. cell%ok(n_radius) .and. cell%ok(n_side) &
+            & .and. .not. &
+            & ( cell%ok(n_meshx) .or. cell%ok(n_meshy) &
+            & .or. cell%ok(n_splitx) .or. cell%ok(n_splity) )
+       select case(sectori)
+       case(S_not)
+           if(size(cell%mix)/=size(cell%radius)) &
+               call XABORT("G2S: wrong len for MIX")
+       case(S_X_tot)
+          if(size(cell%mix)/=6*size(cell%radius)-5*sectorj) &
+               call XABORT("G2S: wrong len for MIX with SECT")
+       end select
+    case default
+       call XABORT("G2S: the type of geometry is not allowed")
+    end select
+  end subroutine verrifieCB
+
+  subroutine exploiteSplit(cell)
+    type(t_celluleBase),intent(inout) :: cell
+    !modification des vecteurs radius, mix et merge si splitr active
+
+    double precision,dimension(:),allocatable :: tmpRad
+    integer,dimension(:),allocatable          :: tmpMix,tmpMrg
+    double precision                          :: interval
+    integer                                   :: longueur,i,j,k,sect,valS,sectori,sectorj,typecel,longueur2, &
+                                                 k2,nsectext,nsectint,ip,nx,ny
+
+    if(.not. cell%ok(n_splitr)) return
+    if((size(cell%splitr)+1)/=size(cell%radius)) call XABORT("G2S : wrong length for SPLITR")
+    longueur=0
+    sectori = cell%sv(14)
+    sectorj = cell%sv(15)
+    typecel = cell%sv(1)
+    do i = 1,size(cell%splitr)
+       longueur=longueur+abs(cell%splitr(i))
+    end do
+    nx=0
+    ny=0
+    nsectint=0
+    nsectext=0
+    select case(typecel)
+    case(G_Carcel)
+       nx=size(cell%meshx)-1
+       ny=size(cell%meshy)-1
+       if((nx*ny.gt.1).and.sectori.ne.S_not) call XABORT("G2S : CARCEL lr lx>1 ly>1 is incompatinle with SECT")
+       select case(sectori)
+       case(S_not)
+          nsectint=1
+          nsectext=1
+       case(S_X_tot,S_T_tot)
+          nsectint=4
+          if(sectorj /= 0) nsectint=1
+          nsectext=4
+       case(S_TX_tot,S_TXS_tot)
+          nsectint=8
+          if(sectorj /= 0) nsectint=1
+          nsectext=8
+       case(S_WM_tot)
+          nsectint=12
+          if(sectorj /= 0) nsectint=1
+          nsectext=12
+       end select
+    case(G_Hexcel)
+       select case(sectori)
+       case(S_not)
+          nsectint=1
+          nsectext=1
+       case(S_X_tot)
+          nsectint=6
+          if(sectorj /= 0) nsectint=1
+          nsectext=6
+       end select
+    end select
+    if(nx*ny.gt.1) then
+       longueur2=nx*ny*(longueur+1)
+    else
+       longueur2=nsectint*longueur+nsectext
+    endif
+    allocate(tmpRad(longueur+1),tmpMix(longueur2),tmpMrg(longueur2),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: exploitSplit(1) => allocation pb")
+    tmpRad(1)=0.d0
+    k = 1
+    k2= 1
+    do i = 1,size(cell%splitr)
+       valS = cell%splitr(i)
+       if(valS==0) then
+          call XABORT("G2S : SPLITR may not be null")
+       else if(valS>0) then
+          !on coupe en sous-rayons egaux -> moyenne arthmetique des rayons
+          interval = (cell%radius(i+1)-cell%radius(i))/valS
+          do j = 1,valS
+             if(nx*ny.gt.1) then
+                do ip=1,nx*ny
+                   tmpMix((ip-1)*(longueur+1)+k) = cell%mix((ip-1)*(size(cell%radius))+i)
+                   tmpMrg((ip-1)*(longueur+1)+k) = (ip-1)*(longueur+1)+k
+                end do
+             else
+                do sect=1,nsectint
+                   tmpMix(k2) = cell%mix(nsectint*(i-1)+sect)
+                   !! BEFORE tmpMrg(k) = cell%merge(i)
+                   tmpMrg(k2) = k2
+                   k2= k2+ 1
+                end do
+             endif
+             k = k + 1
+             tmpRad(k) = tmpRad(k-1) + interval
+          end do
+       else
+          !on coupe en sous-surfaces egales -> moyenne geometrique des rayons
+          valS = abs(valS)
+          interval = (cell%radius(i+1)**2-cell%radius(i)**2)/valS
+          do j = 1,valS
+             if(nx*ny.gt.1) then
+                do ip=1,nx*ny
+                   tmpMix((ip-1)*(longueur+1)+k) = cell%mix((ip-1)*(size(cell%radius))+i)
+                   tmpMrg(k) = (ip-1)*(longueur+1)+k
+                end do
+             else
+                do sect=1,nsectint
+                   tmpMix(k2) = cell%mix(nsectint*(i-1)+sect)
+                   !! BEFORE tmpMrg(k) = cell%merge(i)
+                   tmpMrg(k2) = k2
+                   k2= k2+ 1
+                end do
+             endif
+             k = k + 1
+             tmpRad(k) = sqrt(tmpRad(k-1)**2 + interval)
+          end do
+       end if
+    end do
+    if(nx*ny.gt.1) then
+       do ip=1,nx*ny
+          tmpMix(ip*(longueur+1)) = cell%mix(ip*(size(cell%radius)))
+          tmpMrg(ip*(longueur+1)) = ip*(longueur+1)
+       end do
+    else
+       do sect=1,nsectext
+          tmpMix(k2) = cell%mix(nsectint*size(cell%splitr)+sect)
+          !!BEFORE tmpMrg(longueur+1) = cell%merge(size(cell%merge))
+          tmpMrg(k2) = k2
+          k2= k2+ 1
+       end do
+    endif
+    deallocate(cell%radius) ; allocate(cell%radius(longueur+1),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: exploiteSplit(2) => allocation pb")
+    deallocate(cell%mix) ; allocate(cell%mix(longueur2),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: exploiteSplit(3) => allocation pb")
+    deallocate(cell%merge) ; allocate(cell%merge(longueur2),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: exploiteSplit(4) => allocation pb")
+    cell%radius(:longueur+1) = tmpRad(:longueur+1)
+    cell%mix(:longueur2) = tmpMix(:longueur2) ; cell%merge(:longueur2) = tmpMrg(:longueur2)
+    deallocate(tmpRad,tmpMix,tmpMrg)
+  end subroutine exploiteSplit
+  
+  recursive subroutine buildCellsBase(ip,sz,dirname)
+    type(c_ptr),intent(inout):: ip
+    integer,intent(inout)    :: sz
+    character*12 ,intent(in) :: dirname
+
+    character*12             :: namp, savename, subdirname
+    integer                  :: type, long
+
+    integer,parameter        :: dir=0
+
+    call writeCellBase(ip,sz,dirname)
+    namp = ' '
+    call LCMNXT(ip,namp)
+    savename = namp
+    do 
+       if(namp == ' ') exit
+       call LCMLEN(ip,namp,long,type)
+       if(type == dir) then
+          if(namp /= 'BIHET') then
+             subdirname = namp
+             call LCMSIX(ip,namp,1)
+             call buildCellsBase(ip,sz,subdirname)
+             call LCMSIX(ip,namp,2)
+          end if
+       end if
+       call LCMNXT(ip,namp)
+       if(namp == savename) exit
+    end do
+  end subroutine buildCellsBase
+
+  subroutine writeCellBase(ip,sz,dirname)
+    type(c_ptr),intent(in)   :: ip
+    integer,intent(inout)    :: sz
+    character*12 ,intent(in) :: dirname
+
+    integer,dimension(40)    :: st !state vector
+    integer                  :: type, long, i
+    logical                  :: toCreate
+
+    call LCMLEN(ip,'STATE-VECTOR',long,type)
+    if(long==0) then
+       !on est dans la partie nouvelle (donnees ajoutees par pretaitement
+       !python => on sort de la subroutine
+       return
+    end if
+    call LCMGET(ip,'STATE-VECTOR',st)
+    call LCMLEN(ip,'NPIN        ',long,type) ! pour tester si c'est un cluster
+    if(st(8)==0 .and. long==0) then !pas de sous-cellules
+       toCreate = .true.
+       do i = 1,sz
+          if(tabCelluleBase(i)%name==dirname) toCreate = .false.
+       end do
+       if(toCreate) then
+          sz = sz + 1
+!!$          allocate(tabCelluleBase(sz)%p,stat=alloc_ok)
+!!$          if (alloc_ok /= 0) call XABORT("G2S: writeCellBase => allocation pb")
+          call createCB(tabCelluleBase(sz),dirname,ip)
+       end if
+    end if
+  end subroutine writeCellBase
+
+end module celluleBase