summaryrefslogtreecommitdiff
path: root/Dragon/src/g2s_boundCond.f90
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/g2s_boundCond.f90')
-rw-r--r--Dragon/src/g2s_boundCond.f901738
1 files changed, 1738 insertions, 0 deletions
diff --git a/Dragon/src/g2s_boundCond.f90 b/Dragon/src/g2s_boundCond.f90
new file mode 100644
index 0000000..a7223c7
--- /dev/null
+++ b/Dragon/src/g2s_boundCond.f90
@@ -0,0 +1,1738 @@
+!
+!-----------------------------------------------------------------------
+!
+!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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:02 +0000