path: root/Dragon/src/g2s_boundCond.f90

!
!-----------------------------------------------------------------------
!
!Purpose:
! Process data relative to boundary conditions.
!
!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:
! Deux structures differentes ont ete definies en variables globales:
!  - bCData : pour stocker les donnees en entree du code
!  - SALbCData : pour stocker les donnees correspondantes a fournir dans
!                le jeux de donnees SAL genere
! Le fonction definies sont
!  - appliBoundariConditions : calcule les valeurs des champs des structures
!                              en appelant les fonctions specifiques a chaque
!                              type de geometrie
! \\\\
!  - initializebCData : mise a zero des stuctures
!  - destroybCData : liberation de la memoire
!  - appliBoundariConditionsForXXX : fonction specifique au type XXX avec
!                                    XXX = Rec, Hex, Tri ou Tub
!  - setBoundSide : applique une condition limite de contact aux elements
!  - setBoundCut : applique une condition limite avec reduction du domaine
!                  aux elements
!  - prepareSALBCData : prepare les donnees de conditions limites pour SAL
!
!-----------------------------------------------------------------------
!
module boundCond
  use cellulePlaced
  use constType
  use constUtiles
  use GANLIB
  use segArc

  implicit none  

  type t_bCData
     double precision,dimension(2) :: sidexy !longueur en xy de la gigogne
     !ext (pour rec et tri ST60 ou COMPLETE), ou cote d'un hexagone et cote
     !de l'assemblage (pour hexa et tri S30 et SA60)
     double precision,dimension(4) :: minmaxXY
     integer,dimension(6)          :: bc
     double precision,dimension(6) :: albedo
     integer,dimension(6)          :: albInd
     double precision,dimension(2) :: toOrig_xy !vecteur de translation du
     !centre vers le coins inferieur gauche
     integer                       :: iHex !type de geo si hexagone
     integer                       :: iTri !type de geo si triangle
  end type t_bCData

  type(t_bCData),save  :: bCData

  type t_SALbCData
     integer                      :: sunsetType
     integer                      :: SALtype
     integer                      :: nber
     integer,dimension(:),allocatable :: elemNb
     real                         :: albedo
     real                         :: tx,ty
     real                         :: cx,cy,angle
  end type t_SALbCData

  type(t_SALbCData),dimension(6),save  :: SALbCDataTab

contains

  subroutine initializebCData()
    SALbCDataTab(:6)%nber = 0
  end subroutine initializebCData

  subroutine destroybCData()
    integer :: i
    do i = 1,6
       if (allocated(SALbCDataTab(i)%elemNb)) then
          deallocate(SALbCDataTab(i)%elemNb)
       end if
    end do
  end subroutine destroybCData

  subroutine appliBoundariConditions(ip,szSA,nbCLP)
    type(c_ptr),intent(in):: ip
    integer,intent(inout) :: szSA
    integer,intent(out)   :: nbCLP

    select case(geomTyp)
    case(RecTyp)
       !       write(*,*) 'entering appliBoundariConditionsForRec'
       call appliBoundariConditionsForRec(ip,szSA,nbCLP)
    case(HexTyp)
       !       write(*,*) 'entering appliBoundariConditionsForHex'
       call appliBoundariConditionsForHex(ip,szSA,nbCLP)
    case(TriaTyp)
       !       write(*,*) 'entering appliBoundariConditionsForTri'
       call appliBoundariConditionsForTri(ip,szSA,nbCLP)
    case(TubeTyp)
       !       write(*,*) 'entering appliBoundariConditionsForTub'
       call appliBoundariConditionsForTub(ip,nbCLP)
    end select
  end subroutine appliBoundariConditions

  subroutine appliBoundariConditionsForRec(geoIp,szSA,nbCLP)
    type(c_ptr),intent(in)    :: geoIp
    integer,intent(inout) :: szSA
    integer,intent(out)   :: nbCLP

    double precision              :: rminx,rminy,rmaxx,rmaxy
    double precision,dimension(4) :: x,y,xx,yy,cx,cy,cxx,cyy
    type(c_ptr)                   :: ip
    integer                       :: i,nbCut
    real,dimension(2)             :: tmpTab2
    real,dimension(4)             :: tmpTab4
    real,dimension(6)             :: tmpTab6
    type(t_segArc)                :: sg

    ! programmation defensive
    integer :: dimTabSegArc

    dimTabSegArc = size(tabSegArc)

    ip = geoIp
    !recuperations des donnees sur les conditions aux limites
    call LCMGET(ip,'NCODE       ',bCData%bc)
    call LCMGET(ip,'ZCODE       ',tmpTab6) ; bCData%albedo=tmpTab6
    call LCMGET(ip,'ICODE       ',bCData%albInd)
    !    write(*,*) 'NCODE :',bCData%bc
    !    write(*,*) 'ZCODE :',bCData%albedo
    !    write(*,*) 'ICODE :',bCData%albInd
    call LCMSIX(ip,'NEW-DATA    ',1)
    call LCMSIX(ip,'BOUND-DATA  ',1)
    call LCMGET(ip,'SIDEXY      ',tmpTab2) ; bCData%sidexy=tmpTab2
    call LCMGET(ip,'MINMAXXY    ',tmpTab4) ; bCData%minmaxXY=tmpTab4
    call LCMSIX(ip,'BOUND-DATA  ',2)
    call LCMSIX(ip,'NEW-DATA    ',2)

    !exploitation des donnees
    rmaxx = 0.5d0*bCData%sidexy(1) ; rminx = -rmaxx
    rmaxy = 0.5d0*bCData%sidexy(2) ; rminy = -rmaxy
    x(1) = rminx ; y(1) = rmaxy ; xx(1) = rminx ; yy(1) = rminy
    x(2) = rmaxx ; y(2) = rminy ; xx(2) = rmaxx ; yy(2) = rmaxy
    x(3) = rminx ; y(3) = rminy ; xx(3) = rmaxx ; yy(3) = rminy
    x(4) = rmaxx ; y(4) = rmaxy ; xx(4) = rminx ; yy(4) = rmaxy 
    cx(1)  = bCData%minmaxXY(1) ; cy(1)  = rmaxy
    cxx(1) = bCData%minmaxXY(1) ; cyy(1) = rminy
    cx(2)  = bCData%minmaxXY(3) ; cy(2)  = rminy
    cxx(2) = bCData%minmaxXY(3) ; cyy(2) = rmaxy
    cx(3)  = rminx              ; cy(3)  = bCData%minmaxXY(2)
    cxx(3) = rmaxx              ; cyy(3) = bCData%minmaxXY(2)
    cx(4)  = rmaxx              ; cy(4)  = bCData%minmaxXY(4)
    cxx(4) = rminx              ; cyy(4) = bCData%minmaxXY(4)

    !creation des conditions aux limites (pour une geometrie carre)
    !lorsqu'on coupe, on garde le cote gauche (sens trigo)
    nbCLP=0
    nbCut=0
    do i = 1,4
       select case(bCData%bc(i))
       case(B_Void,B_Refl,B_Ssym,B_Albe,B_Zero,B_Tran,B_Pi_2,B_Pi)
          !write(*,*) 'reflective boundary condition ',i
          if (bCData%bc(i)/=B_Void) nbCLP=nbCLP+1
          call setBoundSide(x(i),y(i),xx(i),yy(i),bCData%bc(i)+100*i,szSA)
       case(B_Diag)
          !write(*,*) 'diagonal symmetry ',i,"seg ",szsa+1
          if (i>=3) cycle
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1 
          ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForRec (1)")
          if (i==1) then
             sg=createSeg(x(4),y(4),x(3),y(3),fooMix,-(bCData%bc(i)+100*i))
          else
             sg=createSeg(x(3),y(3),x(4),y(4),fooMix,-(bCData%bc(i)+100*i))
          end if
          tabSegArc(szSA)=sg
       case(B_Syme)
          !write(*,*) 'syme ', i,"seg ",szsa+1
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1  
          ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForRec (1)")
          sg=createSeg(cx(i),cy(i),cxx(i),cyy(i),fooMix,-(bCData%bc(i)+100*i))
          tabSegArc(szSA)=sg
       end select
    end do
    if (nbCut/=0) call setBoundCut_V2(nbCut,szSA)
  end subroutine appliBoundariConditionsForRec

  subroutine setBoundSide(x,y,xx,yy,nbCL,szSA)
    double precision,intent(in) :: x,y,xx,yy
    integer,intent(in)          :: nbCL,szSA

    type(t_segArc) :: sa
    integer        :: i

    do i = 1,szSA
       sa = tabSegArc(i)
       if (sa%typ/=tseg) cycle
       if (.not. ( estColi(xx-x,yy-y,sa%dx-sa%x,sa%dy-sa%y) .and. &
            &   pointsAlignes(x,y,xx,yy,sa%x,sa%y) ) ) cycle
       if (sa%mixd==fooMix) then
          tabSegArc(i)%mixd=-nbCL
       else if (sa%mixg==fooMix) then
          tabSegArc(i)%mixg=-nbCL
       end if
    end do
  end subroutine setBoundSide

  subroutine setBoundCut_V2(nbCut,szSA)
    integer, intent(in)    :: nbCut
    integer, intent(inout) :: szSA

    ! segments management
    integer :: iSA, iSC, iStr, thisSA, thisOtherSA, szSC
    type(t_SegArc) :: SA, SA1, SA2, SA3, SC, SC1, SC2, SC3
    double precision :: interangle, interAngle2, RefAngle, RefOtherAngle
    double precision :: thisAngle, thisOtherAngle, angle1,angle2
    ! follow up for intersections
    logical :: L_Inter, test1, test2, test3, test4       
    type(t_point) :: P1, P2, P3, P4
    integer :: n1, n2
    double precision :: aox,aoy,afx,afy
    ! sort
    ! find minimal angle
    double precision :: SCAngle, SA_x,SA_y,SA_dx,SA_dy

    ! Part 1 :  Backup of cutting straights & updating
    allocate(tabStrCut(nbCut), tabSegCut(szSA),stat=alloc_ok)
    if (alloc_ok /= 0) call XABORT("G2S: setBoundCut_V2 => allocation pb")
    tabStrCut(1:nbCut) = tabSegArc(szSA-nbCut+1:szSA)
    szSA = szSA - nbCut
    !
    ! Part2 : Adaptation for cutting Straights 
    B_A1: do iStr = 1, nbCut
       SC = tabStrCut(iStr)
       B_B1: do iSA = iStr+1, nbCut
          SA = tabStrCut(iSA)
          L_Inter = InterSgSg_V2(SC,SA,n1,n2,P1,P2)
          if (L_Inter.and.n2==0) then
             SA1 = copySegExceptEnd(P1%x,P1%y,SA)
             if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA1%x,SA1%y)) then
                tabStrCut(iSA) = SA1
             else
                SA1 = copySegExceptOrigin(P1%x,P1%y,SA)
                tabStrCut(iSA) = SA1
             endif
          end if
       enddo B_B1
    enddo B_A1

    ! Part 3 : process for the current straight
    B_A2: do iStr = 1, nbCut  
       szSC = 1 
       tabSegCut(szSC) = tabStrCut(iStr)
       iSC = szSC 

       B_B2: do
          if (iSC > szSC)  exit B_B2
          SC = tabSegCut(iSC)
          iSA = 1
          B_C2: do
             if (iSA > szSA) exit B_C2
             SA = tabSegArc(iSA)
             if (SA%typ == tseg) then  
                ! SA is a segment
                L_Inter = InterSgSg_V2(SC,SA,n1,n2,P1,P2)
                if (.not.L_Inter) then
                   ! no intersection
                   if  (estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%x,SA%y).or.&
                        estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%dx,SA%dy)) then
                      iSA = iSA + 1
                   else
                      call SlideSA(iSA,szSA)
                   endif
                   cycle B_C2 
                elseif (L_Inter.and.(n2==0)) then
                   ! only one intersection point
                   if (n1 == 22) then
                      SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                      SC2 = createSeg(P1%x,P1%y,SC%dx,SC%dy,SC%mixg,SC%mixd)
                      SC = ReplaceSC(SC1,iSC)
                      call AddSC(SC2,szSC)
                      SA1 = copySegExceptEnd(P1%x,P1%y,SA)
                      if (.not.estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA1%x,SA1%y)) &
                         SA1 = copySegExceptOrigin(P1%x,P1%y,SA)
                      SA = ReplaceSA(SA1,iSA) ; iSA = iSA + 1

                   elseif ((n1 == 21).or.(n1 == 23)) then
                      ! the intersection point 
                      SA1 = copySegExceptEnd(P1%x,P1%y,SA)
                      if (.not.estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA1%x,SA1%y))&
                           SA1 = copySegExceptOrigin(P1%x,P1%y,SA)
                      SA = ReplaceSA(SA1,iSA) ; iSA = iSA + 1

                   elseif ((n1==11).or.(n1==13).or.(n1==31).or.(n1==33)) then
                      if ( ((n1==11).and.&
                           estAGAuche(SC%x,SC%y,SC%dx,SC%dy,SA%dx,SA%dy)) .or. &
                           ((n1==13).and.&
                           estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%dx,SA%dy)) .or. &
                           ((n1==31).and.&
                           estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%x,SA%y  )) .or. &
                           ((n1==33).and.&
                           estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%x, SA%y ))) then
                         iSA = iSA + 1
                      else
                         call SlideSA(iSA,szSA)
                      endif

                   elseif ((n1 == 12).or.(n1 == 32)) then
                      SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                      SC2 = createSeg(P1%x,P1%y,SC%dx,SC%dy,SC%mixg,SC%mixd)
                      SC = ReplaceSC(SC1,iSC)
                      call AddSC(SC2,szSC)
                      if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,SA%x,SA%y)) then
                         iSA = iSA + 1
                      else
                         call SlideSA(iSA,szSA)
                      endif
                   else
                      call XABORT("SetBoundCut_V2 : inconsistency")
                   endif
                   cycle B_C2

                else  ! (L_Inter.and.n2 /= 0 )
                  call OverlaidSegmentsManagement(n1,n2,iSC,szSC,SC,iSA,szSA,SA)
                  cycle B_C2    
                endif
             else      
                ! Sa is an arc or a circle
                L_Inter = interSgAr_V2(SC,SA,n1,n2,P1,P2)
                if (SA%typ == tarc) then
                   call giveOrigine(SA,Aox,Aoy) ; call giveExtremite(SA,Afx,Afy)
                endif
                if (.not.L_Inter) then
                   if (SA%typ /= tcer) then
                      if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,aox,aoy).or. &
                           estAGauche(SC%x,SC%y,SC%dx,SC%dy,afx,afy)) then
                         iSA = iSA + 1
                      else
                         call SlideSA(iSA,szSA)
                      endif
                   else 
                      ! SA is a circle
                      call giveFourPointsOnCircle(SA,P1,P2,P3,P4)
                      test1 = estAGauche(SC%x,SC%y,SC%dx,SC%dy,P1%x,P1%y)
                      test2 = estAGauche(SC%x,SC%y,SC%dx,SC%dy,P2%x,P2%y)
                      test3 = estAGauche(SC%x,SC%y,SC%dx,SC%dy,P3%x,P3%y)
                      test4 = estAGauche(SC%x,SC%y,SC%dx,SC%dy,P4%x,P4%y)
                      if  (test1.or.test2.or.test3.or.test4) then
                         ! circle will be cut by another SC
                         iSA = iSA + 1
                      else
                         call SlideSA(iSA,szSA)
                      endif
                   endif
                   cycle B_C2
                elseif (n2 == 0) then
                   if (n1>100) then
                      ! one intersection point which is a tangency point
                      if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,aox,aoy)) then
                         SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                         SC2 = createSeg(P1%x,P1%y,SC%x,SC%y,SC%mixg,SC%mixd)
                         SC = ReplaceSC(SC1,iSC)
                         call AddSC(SC2,szSC)
                         interAngle = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                         if (SA%typ == tarc) then
                            SA1 = copyArcExceptEnd(interAngle,SA)
                            SA2 = copyArcExceptOrigin(interAngle,SA)
                            SA = ReplaceSA(SA1,iSA)
                            call AddSA(SA2,szSA)
                         else
                            interangle2 = anglenormal(interangle+pi_c)
                            SA1 = createArcFromCircle(interAngle,interAngle2,SA)
                            SA2 = createArcFromCircle(interAngle2,interAngle,SA)
                            SA = ReplaceSA(SA1,iSA)
                            call AddSA(SA2,szSA)
                         end if
                         iSA = iSA + 1
                      else
                         call SlideSA(iSA,szSA)
                      end if
                   elseif (n1<100) then
                      ! one intersection point P1
                      if (mod(n1,10) == 2) then
                         ! P1 is inside segment
                         SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                         SC2 = createSeg(P1%x,P1%y,SC%dx,SC%dy,SC%mixg,SC%mixd)
                         SC = ReplaceSC(SC1,iSC)
                         call AddSC(SC2,szSC)
                         interAngle = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                         if (n1/10 == 2) then
                            if (SA%typ == tarc) then
                               SA1 = copyArcExceptEnd(interAngle,SA)
                               if (.not.estAGauche &
                                    (SC%x,SC%y,SC%dx,SC%dy,aox,aoy)) &
                                   SA1 = copyArcExceptOrigin(interAngle,SA)
                               SA = ReplaceSA(SA1,iSA)
                            else
                               interangle2 = anglenormal(interangle+pi_c)
                               SA1 = createArcFromCircle &
                                    (interAngle,interAngle2,SA) 
                               call MediumPointOnArc(SA1,P1)
                               if (.not.estAGauche &
                                    (SC%x,SC%y,SC%dx,SC%dy,P1%x,P1%y)) &
                                    SA1 = createArcFromCircle &
                                          (interAngle2,interAngle,SA) 
                               SA = ReplaceSA(SA1,iSA)
                            endif
                         endif
                      else
                         ! P1 is the begin or the end of SC
                         if (n1/10 == 2) then
                            ! P1 is inside SA
                            interAngle = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                            if (SA%typ == tarc) then
                               SA1 = copyArcExceptEnd(interAngle,SA)
                               if (.not.estAGauche &
                                        (SC%x,SC%y,SC%dx,SC%dy,aox,aoy)) &
                                  SA1 = copyArcExceptOrigin(interAngle,SA)
                            else
                               interangle2 = anglenormal(interangle+pi_c)
                               SA1 = createArcFromCircle &
                                     (interAngle,interAngle2,SA) 
                               call MediumPointOnArc(SA1,P1)
                               if (.not.estAGauche &
                                        (SC%x,SC%y,SC%dx,SC%dy,P1%x,P1%y)) &
                                  SA1 = createArcFromCircle &
                                        (interAngle2,interAngle,SA) 
                            endif
                            SA = ReplaceSA(SA1,iSA)
                            SC = ReplaceSC(SC1,iSC)
                            iSA = iSA + 1
                         else
                            ! P1 is at the begin or the end of SA
                            if (((n1/10==1).and.                               &
                                 estAGauche(SC%x,SC%y,SC%dx,SC%dy,afx,afy))    &
                                 .or.                                          &
                                 ((n1/10==3).and.                              &
                                 estAGauche(SC%x,SC%y,SC%dx,SC%dy,aox,aoy))) then
                               iSA =iSA + 1
                            else
                               call SlideSA(iSA,szSA)
                            endif
                         endif
                      endif
                   endif
                else 
                   ! two points of intersection
                   if ((n1==21.and.n2==23) .or.(n1==23.and.n2==21)) then
                      ! the intersection points are the extremities of SA
                      interangle  = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                      interangle2 = calculeAngle(SA%x,SA%y,P2%x,P2%y)
                      angle1 = min(interangle,interangle2)
                      angle2 = max(interangle,interangle2)
                      if (SA%a > SA%b) then
                         SA1 = CopyArcExceptEnd(angle1,SA)
                         SA2 = CopyArcWithNewAngles(angle1,angle2,SA)
                         SA3 = CopyArcExceptOrigin(angle2,SA)
                      else
                         SA1 = CopyArcExceptEnd(angle2,SA)
                         SA2 = CopyArcWithNewAngles(angle2,angle1,SA)
                         SA3 = CopyArcExceptOrigin(angle1,SA)
                      endif
                      call MediumPointOnArc(SA2,P1)
                      if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,P1%x,P1%y)) then
                         SA = ReplaceSA(SA2,iSA)
                      else
                         SA = ReplaceSA(SA1,iSA)
                         call AddSA(SA3,szSA)
                      end if
                      iSA = iSA + 1
                   elseif((n1==12.and.n2==32).or.(n1==32.and.n2==12)) then
                      SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                      SC2 = createSeg(P1%x,P1%y,P2%x,P2%y,SC%mixg,SC%mixd)
                      SC3 = createSeg(P2%x,P2%y,SC%dx,SC%dy,SC%mixg,SC%mixd)
                      SC = ReplaceSC(SC1,iSC)
                      call Add2SC(SC2,SC3,szSC)
                      iSA = iSA + 1
                   elseif ((n1==11.and.n2==33).or. &
                           (n1==31.and.n2==13).or. &
                           (n1==13.and.n2==31)) then
                      ! P1 and P2 are at the begin or the end of SA and SC
                      ! nothing to do
                      iSA = iSA + 1
                   else
                      SC1 = createSeg(SC%x,SC%y,P1%x,P1%y,SC%mixg,SC%mixd)
                      SC2 = createSeg(P1%x,P1%y,P2%x,P2%y,SC%mixg,SC%mixd)
                      SC3 = createSeg(P2%x,P2%y,SC%dx,SC%dy,SC%mixg,SC%mixd)
                      if (SA%typ == tarc) then
                         interAngle = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                         SA1 = copyArcExceptEnd(interAngle,SA)
                         interAngle = calculeAngle(SA%x,SA%y,P2%x,P2%y)
                         SA2 = copyArcWithNewAngles(SA1%b,interAngle,SA)
                         SA3 = copyArcExceptOrigin(SA2%a,SA)
                         ! if SA1 is on left side, SA3 is on left side too
                         ! if SA1 isn't on left side, SA2 is on left side
                         call giveOrigine(SA1,aox,aoy)
                         if (estAGauche(SC%x,SC%y,SC%dx,SC%dy,aox,aoy)) then
                            SA = ReplaceSA(SA1,iSA)
                            call AddSA(SA3,szSA)
                         else
                            SA = ReplaceSA(SA2,iSA)
                         endif
                      else
                         interAngle = calculeAngle(SA%x,SA%y,P1%x,P1%y)
                         interAngle2 = calculeAngle(SA%x,SA%y,P2%x,P2%y)
                         SA1 = createArcFromCircle(interAngle2,interAngle,SA)
                         call MediumPointOnArc(SA1,P3)
                         if (.not.estAGauche(SC%x,SC%y,SC%dx,SC%dy,P3%x,P3%y)) &
                              SA1 = createArcFromCircle &
                                    (interAngle,interAngle2,SA)
                         SA = ReplaceSA(SA1,iSA)
                      endif
                      SC = ReplaceSC(SC1,iSC)
                      call Add2SC(SC2,SC3,szSC)
                       iSA = iSA + 1
                   endif
                endif
                cycle B_C2
             endif
          enddo B_C2
          iSC = iSC + 1
       enddo B_B2
       
!!$       ! Part 3 : Sort
!!$       B_B3:do iSC = 1,szSC-1
!!$          P1%x = tabSegCut(iSC)%dx ; P1%y = tabSegCut(iSC)%dy
!!$          B_C3:do jSC = iSC+1,szSC
!!$             if (tabSegCut(jSC)%x == P1%x .and.tabSegCut(jSC)%y == P1%y) then
!!$                if(jSC==iSC+1)  then ; cycle B_B3
!!$                else
!!$                   SCBuffer = tabSegCut(iSC+1)
!!$                   tabSegCut(iSC+1) = tabSegCut(jSC)
!!$                   tabSegCut(jSC) = SCBuffer
!!$                endif
!!$             end if
!!$          end do B_C3
!!$       end do B_B3
       
       ! Part 4 : search for neighbourhood informations
       B_B4: do iSC = 1, szSC
          SC = tabSegCut(iSC)
          thisSA  = 0         ; thisOtherSA = 0
          RefAngle = infinity ; RefOtherAngle = infinity
          SCAngle = calculeAngle(SC%x,SC%y,SC%dx,SC%dy)
          ! Mix awarding
          B_C4: do iSA = 1, szSA
             SA = tabSegArc(iSA)
             if (SA%typ == tseg) then
                if (IsEqualConst(SC%x,SA%x).and.IsEqualConst(SC%y,SA%y)) then
                   thisAngle = calculeAngle(SA%x,SA%y,SA%dx,SA%dy) - SCAngle
                   if (thisAngle < RefAngle) then
                      thisSA = iSA ; RefAngle = thisAngle
                   endif
                elseif (IsEqualConst(SC%x,SA%dx).and.IsEqualConst(SC%y,SA%dy)) then
                   thisAngle = calculeAngle(SA%x,SA%y,SA%dx,SA%dy) + SCAngle
                   if (thisAngle < RefAngle) then
                      thisSA = iSA ; RefAngle = thisAngle
                   endif
                endif
             else if (SA%typ == tarc) then
                call giveOrigine(SA,SA_x,SA_y)
                call giveExtremite(SA,SA_dx,SA_dy)
                if (IsEqualConst(SC%x,SA_x) .and. IsEqualConst(SC%y,SA_y)) then
                   thisAngle = SA%a + pi_2_c - SCAngle
                   if (thisAngle < RefAngle) then
                      thisSA = iSA ; RefAngle = thisAngle
                   endif
                elseif (IsEqualConst(SC%x,SA_dx).and.IsEqualConst(SC%y,SA_dy)) then
                   thisAngle = SA%b - pi_2_c + SCAngle
                   if (thisAngle < RefAngle) then
                      thisSA = iSA ; RefAngle = thisAngle
                   endif
                endif
             endif
          enddo B_C4
          if (thisSA == 0) then
            call XABORT('unable to find element thisSA')
          endif
          if (tabSegArc(thisSA)%typ == tseg) then
             if (isEqualConst(tabSegArc(ThisSA)%x,SC%x).and.           &
                  isEqualConst(tabSegArc(ThisSA)%y,SC%y)) then
                SC%mixg      = tabSegArc(thisSA)%mixd
                SC%nodeg     = tabSegArc(thisSA)%noded
                SC%IndCellPg = tabSegArc(thisSA)%indCellPd 
             elseif(isEqualConst(tabSegArc(ThisSA)%dx,SC%x).and.       &
                  isEqualConst(tabSegArc(ThisSA)%dy,SC%y)) then
                SC%mixg = tabSegArc(thisSA)%mixg
                SC%nodeg     = tabSegArc(thisSA)%nodeg
                SC%IndCellPg = tabSegArc(thisSA)%indCellPg
              else 
                call XABORT("SetBoundCut_V2 : mix error (1)")
             endif
          else
             call giveOrigine(tabSegArc(thisSA),SA_x,SA_y)
             call giveExtremite(tabSegArc(thisSA),SA_dx,SA_dy)
             if (isEqualConst(SA_x,SC%x).and.isEqualConst(SA_y,SC%y)) then
                SC%mixg = tabSegArc(thisSA)%mixd
                SC%nodeg     = tabSegArc(thisSA)%noded
                SC%IndCellPg = tabSegArc(thisSA)%indCellPd 
              elseif (isEqualConst(SA_dx,SC%x).and.isEqualConst(SA_dy,SC%y)) then
                SC%mixg = tabSegArc(thisSA)%mixg
                SC%nodeg     = tabSegArc(thisSA)%nodeg
                SC%IndCellPg = tabSegArc(thisSA)%indCellPg
              else
                call XABORT("SetBoundCut_V2 : mix error (2)")
             endif
          endif

          ! Mix control
          SCAngle = calculeAngle(SC%dx,SC%dy,SC%x,SC%y)
          B_C5: do iSA = 1,szSA
             SA = tabSegArc(iSA)
             if (SA%typ == tseg) then
                if (IsEqualConst(SC%dx,SA%dx).and.IsEqualConst(SC%dy,SA%dy)) then
                   thisOtherAngle = calculeAngle(SA%x,SA%y,SA%dx,SA%dy) + SCAngle
                   if (thisOtherAngle < RefOtherAngle) then
                      thisOtherSA = iSA ; RefOtherAngle = thisOtherAngle
                   endif
                elseif (IsEqualConst(SC%dx,SA%x).and.IsEqualConst(SC%dy,SA%y)) then
                   thisOtherAngle = calculeAngle(SA%x,SA%y,SA%dx,SA%dy) - SCAngle
                   if (thisOtherAngle < RefOtherAngle) then
                      thisOtherSA = iSA ; RefOtherAngle = thisOtherAngle
                   endif
                endif
             elseif (SA%typ == tarc) then
                call giveOrigine(SA,SA_x,SA_y)
                call giveExtremite(SA,SA_dx,SA_dy)
                if (IsEqualConst(SC%dx,SA_dx).and.IsEqualConst(SC%dy,SA_dy)) then
                   thisOtherAngle = SA%B + pi_2_c - SCAngle
                   if (thisOtherAngle < RefOtherAngle) then
                      thisOtherSA = iSA ; RefOtherAngle = thisOtherAngle
                   endif
                elseif (IsEqualConst(SC%dx,SA_x).and.IsEqualConst(SC%dy,SA_y)) then
                   thisOtherAngle = SA%a - pi_2_c + SCAngle
                   if (thisOtherAngle < RefOtherAngle) then
                      thisOtherSA = iSA ; RefOtherAngle = thisOtherAngle
                   endif
                endif
             endif
          enddo B_C5
          if (thisOtherSA == 0) then
             ! we do not find a neigbour for SC
             call XABORT("SetBoundCut_V2 : incredible situation too")
          endif
          if (tabSegArc(thisOtherSA)%typ == tseg) then
             if (IsEqualConst(SC%dx,tabSegArc(thisOtherSA)%x).and. &
                  IsEqualConst(SC%dy,tabSegArc(thisOtherSA)%y)) then
                if (SC%mixg /= tabSegArc(thisOtherSA)%mixg) then
                   call XABORT("SetBoundCut_V2 : mix error (3)")
                endif
             elseif (IsEqualConst(SC%dx,tabSegArc(thisOtherSA)%dx).and. &
                  IsEqualConst(SC%dy,tabSegArc(thisOtherSA)%dy)) then
                if (SC%mixg /= tabSegArc(thisOtherSA)%mixd)  then
                   call XABORT("SetBoundCut_V2 : mix error (4)")
                end if
             else
                call XABORT("SetBoundCut_V2 : mix error (5)")
             endif
          else
             call giveOrigine(tabSegArc(thisOtherSA),SA_x,SA_y)
             call giveExtremite(tabSegArc(thisOtherSA),SA_dx,SA_dy)
             if (isEqualConst(SC%dx,SA_x).and.isEqualConst(SC%dy,SA_y)) then
                if (SC%mixg /= tabSegArc(thisOtherSA)%mixg) then
                   call XABORT("SetBoundCut_V2 : mix error (6)")
                end if
             elseif(isEqualCOnst(SC%dx,SA_dx).and. &
                    isEqualConst(SC%dy,SA_dy)) then
                if (SC%mixg /= tabSegArc(thisOtherSA)%mixd) then
                    call XABORT("SetBoundCut_V2 : mix error (7)")
                endif
             else
                call XABORT("SetBoundCut_V2 : mix error (8)")
             endif
          endif
          tabSegCut(iSC) = SC
       enddo B_B4
       if (szSA+szSC > size(TabSegArc)) call XABORT("InterSgSg_V2 : memory problem (1)")
       ! Transfer
       tabSegArc(szSA+1:szSA+szSC) = tabSegCut(1:szSC)
       szSA = szSA + szSC
       szSC = 0
    enddo B_A2

    deallocate(tabSegCut, tabStrCut)
  end subroutine setBoundCut_V2

  subroutine setBoundCut(nbCut,szSA)
    integer,intent(in)    :: nbCut
    integer,intent(inout) :: szSA

    type(t_segArc)    :: sa,sgi,sg,ar
    type(segArcArrayTer),dimension(:),allocatable :: tmpTab
    double precision  :: intx,inty,pt1x,pt1y,pt2x,pt2y,angl
    integer           :: i,j,sizeTmp,nbPtInter,sztmp
    logical           :: oInSa,eInSa,oInSgi,eInSgi

    ! programmation defensive
    integer :: taille_table_tmpTab

    !on coupe le domaine de maniere symetrique
    !  et on garde les elements a gauche de l'axe de coupe
    !copie dans un tableau temporaire des segments
    !et nettoyage du tableau global

    sztmp = szSA*10
    ! programmation defensive 
    taille_table_tmpTab = sztmp
    allocate(tmpTab(sztmp),stat=alloc_ok)
    if (alloc_ok /= 0) call XABORT("G2S: setBoundCut(1) => allocation pb")
    do i = 1,nbCut
       tmpTab(i)%sa=tabSegArc(szSA-nbCut+i)
       tmpTab(i)%keep=.true. ; tmpTab(i)%cl=.true.
    end do
    do i = 1,szSA-nbCut
       tmpTab(i+nbCut)%sa=tabSegArc(i)
       tmpTab(i+nbCut)%keep=.true. ; tmpTab(i+nbCut)%cl=.false.
    end do
    sizeTmp=szSA

    !coupage des segments   
    j=0
    do
       j = j+1
       if (j>sizeTmp) exit
       if (.not. tmpTab(j)%cl .or. .not. tmpTab(j)%keep) cycle
       sa = tmpTab(j)%sa
       if (isEqualConst(sa%x,sa%dx).and.isEqualConst(sa%y,sa%dy)) then
          tmpTab(j)%keep=.false.
          cycle
       end if
       i=0
       do
          i = i+1
          if (i>sizeTmp) exit
          if ((.not. tmpTab(i)%keep) .or. i==j) cycle
          if (tmpTab(i)%sa%typ/=tseg) cycle
          sgi=tmpTab(i)%sa
          !cas ou les segments de coupe rencontrent d'autres segments sans les
          ! couper. exemple: les coins des carres pour une symetrie diagonnale
          oInSa=pointsAlignes(sgi%x,sgi%y,sa%x,sa%y,sa%dx,sa%dy).and. &
               & (isIn(sgi%x,sgi%y,sa)==2)
          eInSa=pointsAlignes(sgi%dx,sgi%dy,sa%x,sa%y,sa%dx,sa%dy).and. &
               & (isIn(sgi%dx,sgi%dy,sa)==2)
          if (oInSa.and.(.not.eInSa)) then
             tmpTab(j)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%dx,sgi%dy)) then
                if (sgi%mixd/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixd
             else
                if (sgi%mixg/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixg
             end if
             tmpTab(sizeTmp)%sa%x = sgi%x ; tmpTab(sizeTmp)%sa%y = sgi%y
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%dx,sgi%dy)) then
                if (sgi%mixg/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixg
             else
                if (sgi%mixd/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixd
             end if
             tmpTab(sizeTmp)%sa%dx = sgi%x ; tmpTab(sizeTmp)%sa%dy = sgi%y
             exit
          else if ((.not.oInSa).and.eInSa) then
             tmpTab(j)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%x,sgi%y)) then
                !write(*,*) "zaza sbc 5", i, j, sizetmp, sgi%mixG
                if (sgi%mixg/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixg
             else
                !write(*,*) "zaza sbc 6" ,i, j, sizetmp, sgi%mixd
                if (sgi%mixd/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixd
             end if
             tmpTab(sizeTmp)%sa%x = sgi%dx ; tmpTab(sizeTmp)%sa%y = sgi%dy
              sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%x,sgi%y)) then
                !write(*,*) "zaza sbc 7", i, j, sizetmp, sgi%mixd
                if (sgi%mixd/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixd
             else
                !write(*,*) "zaza sbc 8", i, j, sizetmp, sgi%mixG
                if (sgi%mixg/=fooMix) tmpTab(sizeTmp)%sa%mixg = sgi%mixg
             end if
             tmpTab(sizeTmp)%sa%dx = sgi%dx ; tmpTab(sizeTmp)%sa%dy = sgi%dy
             exit
          end if
          !cas ou les segments a couper interceptent le segment de coupe a
          ! une de ses extremites. Exemple: bord du domaine pour la symetrie
          ! non diagonale sur un carre
          ! Remarque: pas d'exit , mais un cycle a la fin pour que le segment
          ! de coupe puisse etre encore utilise
          oInSgi=pointsAlignes(sa%x,sa%y,sgi%x,sgi%y,sgi%dx,sgi%dy).and. &
               & (isIn(sa%x,sa%y,sgi)==2)
          eInSgi=pointsAlignes(sa%dx,sa%dy,sgi%x,sgi%y,sgi%dx,sgi%dy).and. &
               & (isIn(sa%dx,sa%dy,sgi)==2)
          if (oInSgi.and.(.not.eInSgi)) then
             tmpTab(i)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sgi ; tmpTab(sizeTmp)%cl=tmpTab(i)%cl
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%x,sgi%y)) then
                tmpTab(sizeTmp)%sa%dx=sa%x ;  tmpTab(sizeTmp)%sa%dy=sa%y
                !write(*,*) "zaza sbc 9", i, j, sizetmp, sgi%mixg
                if (sgi%mixg/=fooMix) tmpTab(j)%sa%mixg=sgi%mixg
             else
                tmpTab(sizeTmp)%sa%x=sa%x ;  tmpTab(sizeTmp)%sa%y=sa%y
                !write(*,*) "zaza sbc 10", i, j, sizetmp, sgi%mixd
                if (sgi%mixd/=fooMix) tmpTab(j)%sa%mixg=sgi%mixd
             end if
             tmpTab(sizeTmp)%keep=.true.
             cycle
          else if (eInSgi.and.(.not.oInSgi)) then
             tmpTab(i)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sgi ; tmpTab(sizeTmp)%cl=tmpTab(i)%cl
             if (estAGauche(sa%x,sa%y,sa%dx,sa%dy,sgi%x,sgi%y)) then
                tmpTab(sizeTmp)%sa%dx=sa%dx ;  tmpTab(sizeTmp)%sa%dy=sa%dy
                !write(*,*) "zaza sbc 11", i, j, sizetmp, sgi%mixd
                if (sgi%mixd/=fooMix) tmpTab(j)%sa%mixg=sgi%mixd
             else
                tmpTab(sizeTmp)%sa%x=sa%dx ;  tmpTab(sizeTmp)%sa%y=sa%dy
                !write(*,*) "zaza sbc 12", i, j, sizetmp, sgi%mixg
                if (sgi%mixg/=fooMix) tmpTab(j)%sa%mixg=sgi%mixg
             end if
             tmpTab(sizeTmp)%keep=.true.
             cycle
          end if
          !cas ou on coupe vraiment
          if ( interSgSg(sa,sgi,intx,inty)) then
             tmpTab(j)%keep=.false. ; tmpTab(i)%keep=.false.
             sizeTmp=sizeTmp+1 ; 
             tmpTab(sizeTmp)%sa=sgi ; tmpTab(sizeTmp)%cl=tmpTab(i)%cl
             tmpTab(sizeTmp)%sa%x=intx ; tmpTab(sizeTmp)%sa%y=inty
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")

             tmpTab(sizeTmp)%sa=sgi ; tmpTab(sizeTmp)%cl=tmpTab(i)%cl
             tmpTab(sizeTmp)%sa%dx=intx ; tmpTab(sizeTmp)%sa%dy=inty
             tmpTab(sizeTmp-1)%keep = estAGauche(sa%x,sa%y,sa%dx,sa%dy, &
                  & sgi%dx,sgi%dy)
             tmpTab(sizeTmp)%keep = .not.tmpTab(sizeTmp-1)%keep
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%sa%x=intx ; tmpTab(sizeTmp)%sa%y=inty
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")

             tmpTab(sizeTmp)%sa=sa ; tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%sa%dx=intx ; tmpTab(sizeTmp)%sa%dy=inty
             if (tmpTab(sizeTmp-2)%keep) then
                tmpTab(sizeTmp)%sa%mixg=sgi%mixd
                tmpTab(sizeTmp-1)%sa%mixg=sgi%mixg
                !write(*,*) "zaza sbc 13", i, j, sizetmp, sgi%mixd,sgi%mixg
             else
                tmpTab(sizeTmp)%sa%mixg=sgi%mixg
                tmpTab(sizeTmp-1)%sa%mixg=sgi%mixd
                !write(*,*) "zaza sbc 14" ,i, j, sizetmp, sgi%mixd,sgi%mixg
             end if
             exit
          end if
       end do
    end do
    do j = 1,nbCut
       sa = tmpTab(j)%sa
       do i = nbCut+1,sizeTmp
          if (.not. tmpTab(i)%keep) cycle
          sgi=tmpTab(i)%sa
          if (sgi%typ/=tseg) then
             !elimination des arcs et cercles du mauvais cote
             ! le centre de l'arc est a droite du segment de coupe
             ! => on l'elimine
             tmpTab(i)%keep=.not.estADroiteStrict(sa%x,sa%y,sa%dx,sa%dy, &
                  & sgi%x,sgi%y)
          else !elimination des segments du mauvais cote
             tmpTab(i)%keep=.not.(estADroiteStrict(sa%x,sa%y,sa%dx,sa%dy, &
                  & sgi%x,sgi%y) .or. &
                  & estADroiteStrict(sa%x,sa%y,sa%dx,sa%dy,sgi%dx,sgi%dy) )
          end if
       end do
    end do
    !coupage des arcs
    j=0
    do
       j = j+1
       if (j>sizeTmp) exit
       if ((.not. tmpTab(j)%cl).or.(.not. tmpTab(j)%keep)) cycle
       sg = tmpTab(j)%sa
       if (isEqual(tmpTab(j)%sa%x,tmpTab(j)%sa%dx) .and. &
            &  isEqual(tmpTab(j)%sa%y,tmpTab(j)%sa%dy)) then
          tmpTab(j)%keep=.false.
          cycle
       end if
       i=0
       do
          i = i+1
          if (i>sizeTmp) exit
          if (.not. tmpTab(i)%keep) cycle
          if (tmpTab(i)%sa%typ==tseg) cycle
          ar=tmpTab(i)%sa
          nbPtInter = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y))
          if (nbPtInter==2) then
             !on coupe le segment en trois morceaux
             tmpTab(j)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             !morceau 1
             tmpTab(sizeTmp)%sa=sg ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%sa%dx=pt1x ; tmpTab(sizeTmp)%sa%dy=pt1y
             tmpTab(sizeTmp)%sa%mixg=ar%mixd
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             !morceau 2
             tmpTab(sizeTmp)%sa=sg ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%sa%x=pt1x ; tmpTab(sizeTmp)%sa%y=pt1y
             tmpTab(sizeTmp)%sa%dx=pt2x ; tmpTab(sizeTmp)%sa%dy=pt2y
             tmpTab(sizeTmp)%sa%mixg=ar%mixg
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             !morceau 3
             tmpTab(sizeTmp)%sa=sg ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%sa%x=pt2x ; tmpTab(sizeTmp)%sa%y=pt2y
             tmpTab(sizeTmp)%sa%mixg=ar%mixd
             !on coupe l'arc en 2 morceaux et on n'en garde qu'un
             tmpTab(i)%keep=.false.
             if (ar%typ==tcer) then
                sizeTmp=sizeTmp+1
                tmpTab(sizeTmp)%sa=ar ; tmpTab(sizeTmp)%cl=.false.
                tmpTab(sizeTmp)%keep=.true.
                tmpTab(sizeTmp)%sa%typ=tarc
                tmpTab(sizeTmp)%sa%a=calculeAngle(ar%x,ar%y,pt2x,pt2y)
                tmpTab(sizeTmp)%sa%b=calculeAngle(ar%x,ar%y,pt1x,pt1y)
             else
                ! il s'agit ici d'un arc centre sur le segment et le coupant 2
                ! fois, dans une geometrie symetrique... => a priori,
                !  cela ne doit jamais arriver
                call XABORT("G2S : Error, not symetrical geometry")
             end if
             exit
          else if (nbPtInter==1) then
             !on coupe le segment en deux morceaux
             tmpTab(j)%keep=.false.
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut") 
             !morceau 1
             tmpTab(sizeTmp)%sa=sg ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%sa%dx=pt1x ; tmpTab(sizeTmp)%sa%dy=pt1y
             sizeTmp=sizeTmp+1 
             ! programmation defensive
             if (sizeTmp > taille_table_tmpTab) &
                  call XABORT("G2S : memory problem in routine SetBoundCut")
             !morceau 2
             tmpTab(sizeTmp)%sa=sg ; tmpTab(sizeTmp)%cl=.true.
             tmpTab(sizeTmp)%keep=.true.
             tmpTab(sizeTmp)%sa%x=pt1x ; tmpTab(sizeTmp)%sa%y=pt1y
             !determination des milieux pour les segments selon l'ordre
             !relatif des points Osg, InterSgAr, Car... 
             !(en fait selon le sens relatif des segments OI et CI)
             if (((pt1x-sg%x)*(pt1x-ar%x)+(pt1y-sg%y)*(pt1y-ar%y))>0.) then
                !dans le meme sens => premier segment a l'interieur
                tmpTab(sizeTmp-1)%sa%mixg=ar%mixg
                tmpTab(sizeTmp)%sa%mixg=ar%mixd
             else !=> premier segment a l'exterieur
                tmpTab(sizeTmp-1)%sa%mixg=ar%mixd
                tmpTab(sizeTmp)%sa%mixg=ar%mixg
             end if
             !on coupe l'arc en 1 ou 2 morceaux et on n'en garde qu'un
             tmpTab(i)%keep=.false.
             if (ar%typ==tcer) then
                sizeTmp=sizeTmp+1 
                ! programmation defensive
                if (sizeTmp > taille_table_tmpTab) &
                     call XABORT("G2S : memory problem in routine SetBoundCut")
                tmpTab(sizeTmp)%sa=ar ; tmpTab(sizeTmp)%cl=.false.
                tmpTab(sizeTmp)%keep=.true.
                tmpTab(sizeTmp)%sa%typ=tarc
                angl=calculeAngle(ar%x,ar%y,pt1x,pt1y)
                if (isPI(angl)) then
                   tmpTab(sizeTmp)%sa%a=-pi_c
                   tmpTab(sizeTmp)%sa%b=pi_c
                else
                   tmpTab(sizeTmp)%sa%a=angl
                   tmpTab(sizeTmp)%sa%b=angl
                end if
             else
                sizeTmp=sizeTmp+1 
                tmpTab(sizeTmp)%sa=ar ; tmpTab(sizeTmp)%cl=.false.
                tmpTab(sizeTmp)%keep=.true.
                !pour savoir quelle partie on garde, on se sert du meme test
                !que pour les milieux
                if (((pt1x-sg%x)*(pt1x-ar%x)+(pt1y-sg%y)*(pt1y-ar%y))>0.) then
                   !dans le meme sens => on garde la deuxieme moitie de l'arc
                   tmpTab(sizeTmp)%sa%a=calculeAngle(ar%x,ar%y,pt1x,pt1y)
                else !=> on garde la premiere moitie de l'arc
                   tmpTab(sizeTmp)%sa%b=calculeAngle(ar%x,ar%y,pt1x,pt1y)
                end if
             end if
             exit
          else !pas d'intersection mais l'arc est peut-etre a eliminer
             if (ar%typ/=tarc) cycle
             call extremitesArc(ar,pt1x,pt1y,pt2x,pt2y)
             if ( ( estADroiteStrict(sg%x,sg%y,sg%dx,sg%dy,pt1x,pt1y) .or.   &
                  & estADroiteStrict(sg%x,sg%y,sg%dx,sg%dy,pt2x,pt2y) )      &
                  &                        .and.     .not.                   &
                  & ( estAGaucheStrict(sg%x,sg%y,sg%dx,sg%dy,pt1x,pt1y) .or. &
                  &   estAGaucheStrict(sg%x,sg%y,sg%dx,sg%dy,pt2x,pt2y) ) )  &
                  & tmpTab(i)%keep=.false.
          end if
       end do
    end do
    !elimination des segments nuls eventuellement restants
    do i = 1,sizeTmp
       sa = tmpTab(i)%sa
       if (sa%typ==tseg) tmpTab(i)%keep=tmpTab(i)%keep .and. &
            & .not.(isEqual(tmpTab(i)%sa%x,tmpTab(i)%sa%dx) .and.&
            &       isEqual(tmpTab(i)%sa%y,tmpTab(i)%sa%dy))
    end do
    !elimination des segments non de coupe se trouvant sur les segments de
    ! coupe (ca arrive pour les sectorisations)
    do j = 1,sizeTmp
       if (.not.(tmpTab(j)%cl .and. tmpTab(j)%keep)) cycle
       sa = tmpTab(j)%sa
       if (sa%typ/=tseg) cycle
       do i = 1,sizeTmp
          if (.not. tmpTab(i)%keep .or. tmpTab(i)%cl ) cycle
          sgi=tmpTab(i)%sa
          if (sgi%typ/=tseg) cycle
          if (.not. estColineaire(sa,sgi)) cycle
          if (.not. pointsAlignes(sa%x,sa%y,sa%dx,sa%dy,sgi%x,sgi%y)) cycle
          !on renverse le sens de sgi si besoin
          if (.not. isSameWay(sa,sgi)) sgi=turnBackSide(tmpTab(i)%sa)
          if (isEqualConst(sa%x,sgi%x).and.isEqualConst(sa%dx,sgi%dx).and. &
               & isEqualConst(sa%y,sgi%y).and.isEqualConst(sa%dy,sgi%dy)) then
             tmpTab(j)%sa%mixg=sgi%mixg
             tmpTab(i)%keep=.false.
             exit
          end if
       end do
    end do
    !recopie dans le tableau global des elements a conserver
    szSA=0
    ! programmation defensive
    if (sizeTmp > size(tabSegArc)) &
         call XABORT("G2S : memory problem in routine SetBoundCut (2)")
    do i = 1,sizeTmp
       if (tmpTab(i)%keep) then
          szSA=szSA+1
          tabSegArc(szSA)=tmpTab(i)%sa
       end if
    end do
    deallocate(tmpTab)
  end subroutine setBoundCut

  subroutine prepareSALBCData(szSA,nbCLP)
    integer,intent(in) :: szSA,nbCLP

    integer        :: i,j,k,sunsetType,tmpTabSize,indice
    logical        :: newBC
    double precision :: basex
    type(t_segArc) :: sa
    integer,dimension(:),allocatable :: tmpTabElemNb

    !enregistrement des numeros des differentes conditions limites
    do i = 1,nbCLP
       SALbCDataTab(i)%sunsetType = -1
       sunsetType = 0
       do j = 1,szSA
          sa = tabSegArc(j)
          if (sa%mixg==fooMix .or. sa%mixd==fooMix) then
             cycle
          else if (sa%mixg<0) then
             sunsetType = -sa%mixg
          else if (sa%mixd<0) then
             sunsetType = -sa%mixd
          else
             cycle
          end if
          newBC = .true.
          do k = 1,i
             newBC = (newBC .and. (sunsetType /= SALbCDataTab(k)%sunsetType) &
                  & .and. (mod(sunsetType,100) /= B_Void) )
          end do
          if (.not. newBC) cycle
          !on a un type de condition limite non encore pris en compte
          SALbCDataTab(i)%sunsetType = sunsetType
          exit
       end do
    end do
    !creation des donnees SAL pour chaque CL
    do i = 1,nbCLP
       SALbCDataTab(i)%nber = 0
       sunsetType = mod(SALbCDataTab(i)%sunsetType,100)
       indice = SALbCDataTab(i)%sunsetType / 100
       select case(geomTyp)
       case(RecTyp)
          select case(sunsetType)
          case(B_Refl,B_Ssym)
             SALbCDataTab(i)%albedo = 1.
             !SALbCDataTab(i)%SALtype = 1 !ne marche pas pour le moment
             ! => on passe par une symetrie, mais exterieure
             SALbCDataTab(i)%SALtype = 4
             select case(indice)
             case(1)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 90.
             case(2)
                SALbCDataTab(i)%cx = real(0.5*bCData%sidexy(1)-bCData%toOrig_xy(1))
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 90.
             case(3)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 0.
             case(4)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = real(0.5*bCData%sidexy(2)-bCData%toOrig_xy(2))
                SALbCDataTab(i)%angle = 0.
             end select
          case(B_Tran)
             SALbCDataTab(i)%SALtype = 2
             select case(indice)
             case(1)
                SALbCDataTab(i)%tx = real(bCData%sidexy(1))
                SALbCDataTab(i)%ty = 0.
             case(2)
                SALbCDataTab(i)%tx = real(-bCData%sidexy(1))
                SALbCDataTab(i)%ty = 0.
             case(3)
                SALbCDataTab(i)%tx = 0.
                SALbCDataTab(i)%ty = real(bCData%sidexy(2))
             case(4)
                SALbCDataTab(i)%tx = 0.
                SALbCDataTab(i)%ty = real(-bCData%sidexy(2))
             end select
          case(B_Diag)
             SALbCDataTab(i)%SALtype = 4
             SALbCDataTab(i)%cx = 0.
             SALbCDataTab(i)%cy = 0.
             SALbCDataTab(i)%angle = 45.          
          case(B_Syme)
             SALbCDataTab(i)%SALtype = 4
             select case(indice)
             case(1)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 90.
             case(2)
                SALbCDataTab(i)%cx = real(bCData%minmaxXY(3)-bCData%toOrig_xy(1))
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 90.
             case(3)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = 0.
                SALbCDataTab(i)%angle = 0.
             case(4)
                SALbCDataTab(i)%cx = 0.
                SALbCDataTab(i)%cy = real(bCData%minmaxXY(4)-bCData%toOrig_xy(2))
                SALbCDataTab(i)%angle = 0.
             end select
          case(B_Albe)
             SALbCDataTab(i)%SALtype = 0
             SALbCDataTab(i)%albedo = real(bCData%albedo(indice))
          case(B_Zero)
             SALbCDataTab(i)%SALtype = 0
             SALbCDataTab(i)%albedo = 0.
          case(B_Pi_2)
             SALbCDataTab(i)%SALtype = 3
             SALbCDataTab(i)%cx = 0.
             SALbCDataTab(i)%cy = 0.
             if (indice==1) then
                SALbCDataTab(i)%angle = -90.
             else if (indice==3) then
                SALbCDataTab(i)%angle = 90.
             end if
          case(B_Pi)
             SALbCDataTab(i)%SALtype = 3
             SALbCDataTab(i)%cx = real(-bCData%toOrig_xy(1))
             SALbCDataTab(i)%cy = real(-bCData%toOrig_xy(2))
             SALbCDataTab(i)%angle = 180.
          end select
       case(HexTyp)
          SALbCDataTab(i)%SALtype = 4
          if (bCData%iHex==H_Complete) SALbCDataTab(i)%SALtype = 2
          basex = 0.5d0*bCData%sidexy(1)
          select case(indice)
          case(1)
             if (bCData%iHex==H_S30) then
                SALbCDataTab(i)%cx    = 0.
                SALbCDataTab(i)%cy    = 0.
                SALbCDataTab(i)%angle = 0.
             else if (bCData%iHex==H_SA60) then
                SALbCDataTab(i)%cx    = real(-bCData%toOrig_xy(1))
                SALbCDataTab(i)%cy    = real(-bCData%toOrig_xy(2))
                SALbCDataTab(i)%angle = -30.
             else if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = 0.0
                SALbCDataTab(i)%ty    = real(2.0d0*sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = 0.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(1)")
             end if
          case(2)
             SALbCDataTab(i)%angle = 30.
             if (bCData%iHex==H_S30) then
                SALbCDataTab(i)%cx    = real(0.5d0*basex)
                SALbCDataTab(i)%cy    = 0.
             else if (bCData%iHex==H_SA60) then
                SALbCDataTab(i)%cx    = real(-bCData%toOrig_xy(1))
                SALbCDataTab(i)%cy    = real(-bCData%toOrig_xy(2))
             else if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = real(1.5d0*basex)
                SALbCDataTab(i)%ty    = real(sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = -60.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(2)")
             end if
          case(3)
             SALbCDataTab(i)%cy    = 0.
             SALbCDataTab(i)%angle = 90.
             if (bCData%iHex==H_S30) then
                SALbCDataTab(i)%cx = real(.5*sqrt(3.)*bCData%sidexy(2))
             else if (bCData%iHex==H_SA60) then
                SALbCDataTab(i)%cx = real(.5*sqrt(3.)*bCData%sidexy(2) &
                     -bCData%toOrig_xy(1))
             else if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = real(1.5d0*basex)
                SALbCDataTab(i)%ty    = real(-sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = -120.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(3)")
             end if
          case(4)
             if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = 0.0
                SALbCDataTab(i)%ty    = real(-2.0d0*sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = 180.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(4)")
             end if
          case(5)
             if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = real(-1.5d0*basex)
                SALbCDataTab(i)%ty    = real(-sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = 120.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(5)")
             end if
          case(6)
             if (bCData%iHex==H_Complete) then
                SALbCDataTab(i)%tx    = real(-1.5d0*basex)
                SALbCDataTab(i)%ty    = real(sqrt(0.75d0)*basex)
                SALbCDataTab(i)%angle = 60.0
             else
                call XABORT("G2S: internal error in prepareSALBCData(6)")
             end if
          end select
       case(TriaTyp)
          select case(sunsetType)
          case(B_Refl)
             SALbCDataTab(i)%SALtype = 1
             SALbCDataTab(i)%albedo = 1.
          case(B_Albe)
             SALbCDataTab(i)%SALtype = 0
             SALbCDataTab(i)%albedo = real(bCData%albedo(indice))          
          case(B_Zero)
             SALbCDataTab(i)%SALtype = 0
             SALbCDataTab(i)%albedo = 0.
          case(B_Syme)
             SALbCDataTab(i)%SALtype = 4
             select case(indice)
             case(1)
                if (bCData%iTri==T_S30) then
                   SALbCDataTab(i)%cx    = 0.
                   SALbCDataTab(i)%cy    = 0.
                   SALbCDataTab(i)%angle = 0.
                else if (bCData%iTri==T_SA60) then
                   SALbCDataTab(i)%cx    = real(-bCData%toOrig_xy(1))
                   SALbCDataTab(i)%cy    = real(-bCData%toOrig_xy(2))
                   SALbCDataTab(i)%angle = -30.
                else
                   call XABORT("G2S: internal error in prepareSALBCData(7)")
                end if
             case(2)
                if (bCData%iTri==T_S30) then
                   SALbCDataTab(i)%cx    = 0.
                   SALbCDataTab(i)%cy    = 0.
                   SALbCDataTab(i)%angle = 30.
                else if (bCData%iTri==T_SA60) then
                   SALbCDataTab(i)%cx    = real(-bCData%toOrig_xy(1))
                   SALbCDataTab(i)%cy    = real(-bCData%toOrig_xy(2))
                   SALbCDataTab(i)%angle = 30.
                else
                   call XABORT("G2S: internal error in prepareSALBCData(8)")
                end if
             case(3)
                SALbCDataTab(i)%cx    = real(bCData%minmaxXY(3)-bCData%toOrig_xy(1))
                SALbCDataTab(i)%cy    = 0.
                SALbCDataTab(i)%angle = 90.
             end select
          end select
       end select
    end do
    !remplissage du tableau des indices des segArc concernes par la condition
    do i = 1,nbCLP
       sunsetType = SALbCDataTab(i)%sunsetType
       allocate(tmpTabElemNb(szSA),stat=alloc_ok)
       if (alloc_ok /= 0) call XABORT("G2S: prepareSALBCData(1) => allocation pb")
       tmpTabSize = 0
       do j = 1,szSA
          sa = tabSegArc(j)
          if (sa%mixg==-sunsetType .or. sa%mixd==-sunsetType) then
             tmpTabSize = tmpTabSize + 1
             tmpTabElemNb(tmpTabSize) = j
          end if
       end do
       SALbCDataTab(i)%nber = tmpTabSize
       allocate(SALbCDataTab(i)%elemNb(tmpTabSize),stat=alloc_ok)
       if (alloc_ok /= 0) call XABORT("G2S: prepareSALBCData(2) => allocation pb")
       do j = 1,tmpTabSize
          SALbCDataTab(i)%elemNb(j) = tmpTabElemNb(j)
       end do
       deallocate(tmpTabElemNb)
    end do
  end subroutine prepareSALBCData

  subroutine appliBoundariConditionsForHex(geoIp,szSA,nbCLP)
    type(c_ptr),intent(in):: geoIp
    integer,intent(inout) :: szSA
    integer,intent(out)   :: nbCLP

    real,dimension(2) :: tmpTab2
    real,dimension(6) :: tmpTab6
    type(c_ptr)       :: ip
    integer           :: i,nbCut
    double precision  :: fx,fy,tx,ty,sqrt3_2,x,y,xx,yy,basex
    double precision  :: rayIntCell,rayIntCore,rayExtCell,rayExtCore
    type(t_segArc)    :: sg

    ! programmation defensive
    integer :: dimTabSegArc
    !
    dimTabSegArc = size(tabSegArc)

    ip = geoIp
    call LCMGET(ip,'NCODE       ',bCData%bc)
    call LCMGET(ip,'ZCODE       ',tmpTab6) ; bCData%albedo=tmpTab6
    call LCMGET(ip,'ICODE       ',bCData%albInd)
    call LCMGET(ip,'IHEX        ',bCData%iHex)
    bCData%sidexy(:)=0.0d0
    sqrt3_2 = 5.d-1*sqrt(3.d0)
    nbCut = 0 ; nbCLP = 0
    rayExtCell = bCData%sidexy(1)
    rayIntCell = sqrt3_2*rayExtCell
    rayExtCore = bCData%sidexy(2)
    rayIntCore = sqrt3_2*rayExtCore
    select case(bCData%iHex)
    case(H_S30)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1           
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = -rayExtCell           ; fy = 0.d0
       tx = rayIntCore+rayExtCell ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = rayIntCore+sqrt3_2*rayIntCell
       fy = 5.d-1*(rayExtCore+rayIntCell)
       tx = 0.d0 ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg
       if (bCData%bc(1)==B_Syme) then
          !3
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1            
          ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForHex")
          fx = rayIntCore ; fy = 0.d0
          tx = fx         ; ty = 5.d-1*rayExtCore
          sg = createSeg(fx,fy,tx,ty,fooMix,-300-bCData%bc(1))
          tabSegArc(szSA) = sg
       end if
    case(H_SA60)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = -sqrt3_2*rayIntCell
       fy = 5.d-1*rayIntCell
       tx = rayIntCore+sqrt3_2*rayIntCell
       ty = -5.d-1*(rayExtCore+rayIntCell)
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = rayIntCore+sqrt3_2*rayIntCell
       fy = 5.d-1*(rayExtCore+rayIntCell)
       tx = 0.d0 ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg
       if (bCData%bc(1)==B_Syme) then
          !3
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1            
          ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForHex")
          fx = rayIntCore ; fy = -5.d-1*rayExtCore
          tx = fx         ; ty = -fy
          sg = createSeg(fx,fy,tx,ty,fooMix,-300-bCData%bc(1))
          tabSegArc(szSA) = sg
       end if
    case(H_SB60)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = -rayExtCell           ; fy = 0.d0
       tx = rayIntCore+rayExtCell ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = 5.d-1*(rayIntCore+rayExtCell)
       fy = sqrt3_2*(rayIntCore+rayExtCell) 
       tx = 0.d0 ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg       
    case(H_S90)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = -rayExtCell           ; fy = 0.d0
       tx = rayIntCore+rayExtCell ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = 0.d0 ; fy = rayExtCore+rayIntCell
       tx = 0.d0 ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg              
    case(H_R120)
       !!!appliquer les rotations adequat (pas simple)
    case(H_R180)
       !!!idem mais plus simple
    case(H_SA180)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = 0.d0 ; fy = rayExtCore+rayIntCell
       tx = 0.d0 ; ty = -fy
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg                    
    case(H_SB180)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForHex")
       fx = -rayExtCore-rayIntCell ; fy = 0.d0
       tx = -fx                    ; ty = 0.d0
       sg = createSeg(fx,fy,tx,ty,fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg       
    case(H_Complete)
       do i=1,6
         if (bCData%bc(1) == B_Tran) then
           nbCLP=nbCLP+1
           call LCMGET(ip,'SIDEXY      ',tmpTab2) ; bCData%sidexy=tmpTab2
           y=0.0d0 ; yy=0.0d0
           basex = 0.5d0*bCData%sidexy(1)
           select case(i)
           case(1)
             x = -0.5d0*basex
             y = -sqrt(0.75d0)*basex
             xx = 0.5d0*basex
             yy = y
           case(2)
             x = -basex
             xx = -0.5d0*basex
             yy = -sqrt(0.75d0)*basex
           case(3)
             x = -0.5d0*basex
             y = sqrt(0.75d0)*basex
             xx = -basex
           case(4)
             x = 0.5d0*basex
             y = sqrt(0.75d0)*basex
             xx = -0.5d0*basex
             yy = y
           case(5)
             x = basex
             xx = 0.5d0*basex
             yy = sqrt(0.75d0)*basex
           case(6)
             x = 0.5d0*basex
             y = -sqrt(0.75d0)*basex
             xx = basex
           end select
           call setBoundSide(x,y,xx,yy,100*i+B_Tran,szSA)
         endif
       enddo
    end select
    if (nbCut/=0) call setBoundCut(nbCut,szSA)
  end subroutine appliBoundariConditionsForHex

  subroutine appliBoundariConditionsForTri(geoIp,szSA,nbCLP)
    type(c_ptr),intent(in):: geoIp
    integer,intent(inout) :: szSA
    integer,intent(out)   :: nbCLP

    real,dimension(2) :: tmpTab2
    real,dimension(4) :: tmpTab4
    real,dimension(6) :: tmpTab6
    type(c_ptr)       :: ip
    integer           :: nbCut,i
    double precision  :: dy
    type(t_segArc)    :: sg
    integer,dimension(40)         :: sv
    double precision,dimension(4) :: x,y
    
    ! programmation defensive
    integer :: dimTabSegArc 
    dimTabSegArc = size(tabSegArc)
    
    ip = geoIp
    call LCMGET(ip,'STATE-VECTOR',sv)
    call LCMGET(ip,'NCODE       ',bCData%bc)
    call LCMGET(ip,'ZCODE       ',tmpTab6) ; bCData%albedo=tmpTab6
    call LCMGET(ip,'ICODE       ',bCData%albInd)
    call LCMGET(ip,'ITRI        ',bCData%iTri)
    call LCMSIX(ip,'NEW-DATA    ',1)
    call LCMSIX(ip,'BOUND-DATA  ',1)
    call LCMGET(ip,'SIDEXY      ',tmpTab2) ; bCData%sidexy=tmpTab2
    call LCMGET(ip,'MINMAXXY    ',tmpTab4) ; bCData%minmaxXY=tmpTab4

    nbCut = 0
    nbCLP = 0
    select case(bCData%iTri)
    case(T_S30)
       !1 (on ne coupe pas car inutile)
       nbCLP=nbCLP+1
       call setBoundSide(0.d0,0.d0,bCData%sidexy(2),0.d0,100+B_Syme,szSA)
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1            
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
            &appliBoundaryConditionsForTri")
       dy = bCData%sidexy(1) * sqrt3_2d
       if (estPaire(sv(4))) then
          x(1) = bCData%sidexy(2)
          y(1) = dy * sv(4)
       else
          x(1) = bCData%sidexy(2) - 2.5d-1*bCData%sidexy(1)
          y(1) = dy * (sv(4) - 5.d-1)
       end if
       sg = createSeg(x(1),y(1),0.d0,0.d0,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg
       if (bCData%bc(1)==B_Syme) then
          !3
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1        ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForTri")

          if (estPaire(sv(4))) then
             x(2) = bCData%sidexy(2)
             y(2) = dy * sv(4)
          else
             x(2) = bCData%sidexy(2) - 5.d-1*bCData%sidexy(1)
             y(2) = dy * (sv(4) - 2.d0/3)
          end if
          bCData%minmaxXY(3)=x(2) !stockage de l'abscisse de l'axe de coupe
          sg = createSeg(x(2),0.d0,x(2),y(2),fooMix,-300-B_Syme)
          tabSegArc(szSA) = sg
       end if
    case(T_SA60)
       !1
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1        
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForTri")
       dy = bCData%sidexy(1) * sqrt3_2d
       if (estPaire(sv(4))) then
          x(1) = bCData%sidexy(2)
          y(1) = -dy * sv(4)
       else
          x(1) = bCData%sidexy(2) - 2.5d-1*bCData%sidexy(1)
          y(1) = -dy * (sv(4) - 5.d-1)
       end if
       sg = createSeg(0.d0,0.d0,x(1),y(1),fooMix,-100-B_Syme)
       tabSegArc(szSA) = sg
       !2
       nbCut=nbCut+1 ; nbCLP=nbCLP+1
       szSA=szSA+1        
       ! programmation defensive
       if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForTri")
       sg = createSeg(x(1),-y(1),0.d0,0.d0,fooMix,-200-B_Syme)
       tabSegArc(szSA) = sg
       if (bCData%bc(1)==B_Syme) then
          !3
          nbCut=nbCut+1 ; nbCLP=nbCLP+1
          szSA=szSA+1        ! programmation defensive
          if (szSA > dimTabSegArc) call XABORT("G2S: memory problem in routine &
               &appliBoundaryConditionsForTri")
          if (estPaire(sv(4))) then
             x(2) = bCData%sidexy(2)
             y(2) = dy * sv(4)
          else
             x(2) = bCData%sidexy(2) - 5.d-1*bCData%sidexy(1)
             y(2) = dy * (sv(4) - 2.d0/3)
          end if
          bCData%minmaxXY(3)=x(2) !stockage de l'abscisse de l'axe de coupe
          sg = createSeg(x(2),-y(2),x(2),y(2),fooMix,-300-B_Syme)
          tabSegArc(szSA) = sg
       end if
    case(T_ST60)
       x(1) = -5.d-1*bCData%sidexy(1) ; x(2) = 0. ; x(3) = -x(1) ; x(4) = x(1)
       y(1) = sqrt3_2d*x(1) ; y(2) = -y(1) ; y(3) = y(1) ; y(4) = y(1)
       do i = 1,3
          if (bCData%bc(i)/= B_Void) nbCLP=nbCLP+1
          call setBoundSide(x(i+1),y(i+1),x(i),y(i),bCData%bc(i)+100*i,szSA)
       end do
    case(T_Complete)
       y(1) = -5.d-1*bCData%sidexy(2) ; y(2) =  y(1)
       y(3) = -y(1)                   ; y(4) =  y(1)
       x(1) = -5.d-1*bCData%sidexy(1) ; x(2) = -x(1)
       x(3) =  x(2) + y(3)            ; x(4) =  x(1) + y(3)
       call setBoundSide(x(4),y(4),x(1),y(1),bCData%bc(1)+100,szSA)
       call setBoundSide(x(2),y(2),x(3),y(3),bCData%bc(2)+200,szSA)
       call setBoundSide(x(1),y(1),x(2),y(2),bCData%bc(3)+300,szSA)
       call setBoundSide(x(3),y(3),x(4),y(4),bCData%bc(4)+400,szSA)
    case default
       call XABORT("G2S: internal error in subroutine &
            &appliBoundariConditionsForTri")
    end select
    if (nbCut/=0) call setBoundCut(nbCut,szSA)
  end subroutine appliBoundariConditionsForTri

  subroutine appliBoundariConditionsForTub(geoIp,nbCLP)
    type(c_ptr),intent(in):: geoIp
    integer,intent(out)   :: nbCLP

    real,dimension(6) :: tmpTab6

    nbCLP = 0
    call LCMGET(geoIp,'NCODE       ',bCData%bc)
    call LCMGET(geoIp,'ZCODE       ',tmpTab6) ; bCData%albedo=tmpTab6
    call LCMGET(geoIp,'ICODE       ',bCData%albInd)
  end subroutine appliBoundariConditionsForTub

end module boundCond