! !----------------------------------------------------------------------- ! !Purpose: ! Collect processing functions related to geometric elements. ! !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: ! Une structure "t_segArc" est definie, avec tous les champs necessaires aux ! calculs. On stocke par ailleurs ces elements dans un tableau global de ! pointeurs "tabSegArc". ! \\\\ ! variable globale: ! - tabSegArc : le tableau de travail principal ! \\\\ ! fonctions du module: ! - initializeTabSegArc : mise a zero des tableaux ! - destroyTabSegArc : liberation de la memoire ! - createSeg : constructeur de segment ! - createArc : constructeur de cercle ou d'arc de cercle ! - giveOrigine : retourne les coordonnees de l'origine d'un segArc ! - giveExtremite : retourne les coordonnees de l'extremite d'un segArc ! - extremitesArc : determination des points extremites d'un arc ! - estColieaire : test de la colinearite de 2 segments ! - isAngleInArc : teste l'appartenance d'un angle au domaine d'un arc ! - isIn : teste l'appartenance d'un point, sur le support d'un segment ou ! d'un arc, au segment ou a l'arc ! - isSameWay : teste si deux segments alignes ont meme orientation ! - turnBackSide : retourne le segment donne en entree ! - adjustMix : ajuste les milieux a droite et a gauche d'un segment, par ! rapport a ceux d'un segment de reference le contenant ! - interSgAr : calcule les points d'intersection entre un segment et un arc ! - interSgSg : calcule les points d'intersection entre deux segments ! - interSgSg_V2 : recherche les points d'itersection entre deux segments ! (extremites inclues) ! - interArcArc : calcule les points d'intersection entre deux arcs ! - translateAndTurnTriSg : effectue une translation-rotation d'un cote ! de triangle (fonction tres specifique) ! - giveExtremalsAngles : donne les angles d'incidence des extrmites d'un ! element, par rapport a un point ! - drawSegArc : affiche le tableau des elements dans un fichier postscript ! - coupeCercle : intersecte les couronnes et les clusters ! - splitSegsForSect : coupe l'encadrement d'une cellule rectangulaire en cas ! de sectorisation ! - majSectori : met a jour la donnee sur le secteur des elements ! - addSegsAndClean : ajout de nouveaux segments lors de la mise cote a cote ! des differentes cellules, et elimination des doublons ! - cutClusters : prise en compte des intersections couronnes/cluster, et ! elimination de elements internes aux clusters ! - segWithSameCoord : true if the transmitted segments have the same coordinates ! !----------------------------------------------------------------------- ! module segArc use constUtiles use constType use derivedPSPLOT implicit none interface printSegArc module procedure printSegArc1, printSegArc2 end interface printSegArc type t_segArc integer :: typ !1=segment, 2=arc de cercle double precision :: x,y !origine si 1 , centre si 2 double precision :: dx,dy !extremite si 1 double precision :: r !rayon si 2 double precision :: a,b !angle debut et fin si 2 (entre -pi et pi) integer :: mixg =0!numero de couronne a gauche ou interieur integer :: mixd =0!numero de couronne a droite ou exterieur integer :: nodeg =0!node a gauche ou interieur integer :: noded =0!node a droite ou exterieur integer :: indCellPg=0 !cellulePlaced d'origine a gauche integer :: indCellPd=0 !cellulePlaced d'origine a droite integer :: sectg=0 !numero du secteur a gauche dans la cellule integer :: sectd=0 !numero du secteur a droite dans la cellule integer :: neutronicMixg=0 !milieux neutronique a gauche integer :: neutronicMixd=0 !milieux neutronique a droite logical :: clusg=.false. ! circle belong to a cluster logical :: clusd=.false. ! circle belong to a cluster end type t_segArc integer,parameter :: tseg=1 , tcer=2 , tarc=3 , fooMix=-99 , fooNode=-999 , & & tRec=1 , tHex=2 , tTri=3 !type of geometry integer :: typgeo !variable globale de type tableau d'elements geometriques type(t_segArc), dimension(:), pointer :: tabSegArc ! storage of cutting straights and segments type(t_segArc), dimension(:), allocatable :: tabStrCut,tabSegCut type segArcArrayBis type(t_segArc) :: sa logical :: keep end type segArcArrayBis type segArcArrayTer type(t_segArc) :: sa logical :: keep logical :: cl !condition limite ou cluster selon les cas end type segArcArrayTer ! Programmation defensive integer :: alloc_ok contains function createSeg(ox,oy,ex,ey,mg,md) double precision,intent(in) :: ox,oy,ex,ey integer,intent(in) :: mg,md type(t_segArc) :: createSeg createSeg%typ = tseg createSeg%x = ox ; createSeg%y = oy createSeg%dx = ex ; createSeg%dy = ey createSeg%r = 0.d0 createSeg%a = 0.d0 ; createSeg%b = 0.d0 createSeg%mixg = mg ; createSeg%mixd = md createSeg%nodeg = fooNode ; createSeg%noded = fooNode createSeg%indCellPg = 0 ; createSeg%indCellPd = 0 createSeg%sectg = 0 ; createSeg%sectd = 0 createSeg%neutronicMixg = 0 ; createSeg%neutronicMixd = 0 createSeg%clusg = .false. ; createSeg%clusd = .false. end function createSeg function copySegExceptOrigin(ox,oy,Seg) double precision,intent(in) :: ox,oy type(t_segArc),intent(in) :: Seg type(t_segArc) :: copySegExceptOrigin if (Seg%typ /= tseg) call XABORT("copySegExceptOrigin : inconsistency") copySegExceptOrigin = Seg copySegExceptOrigin%x = ox ; copySegExceptOrigin%y = oy end function copySegExceptOrigin function copySegExceptEnd(ex,ey,Seg) double precision,intent(in) :: ex,ey type(t_segArc),intent(in) :: Seg type(t_segArc) :: copySegExceptEnd if (Seg%typ /= tseg) call XABORT("copySegExceptEnd : inconsistency") copySegExceptEnd = Seg copySegExceptEnd%dx = ex ; copySegExceptEnd%dy = ey end function copySegExceptEnd function copySegWithNewExtrems(ox,oy,ex,ey,seg) double precision, intent(in) :: ox, oy, ex, ey type(t_segArc),intent(in) :: Seg type(t_segArc) :: copySegWithNewExtrems if (Seg%typ /= tseg) call XABORT("copySegWithNewExtrems: inconsistency") copySegWithNewExtrems = Seg copySegWithNewExtrems%x = ox ; copySegWithNewExtrems%y = oy copySegWithNewExtrems%dx = ex ; copySegWithNewExtrems%dy = ey end function copySegWithNewExtrems function copyUTurnSeg(seg) type(t_segArc), intent(in) :: seg type(t_segArc) :: copyUTurnSeg if (Seg%typ /= tseg) call XABORT("copyUTurnSeg : inconsistency") copyUTurnSeg%typ = tseg copyUTurnSeg%x = seg%dx ; copyUTurnSeg%y = seg%dy copyUTurnSeg%dx = seg%x ; copyUTurnSeg%dy = seg%y copyUTurnSeg%r = seg%r ; copyUTurnSeg%a = seg%a ; copyUTurnSeg%b = seg%b copyUTurnSeg%mixg = seg%mixd ; copyUTurnSeg%mixd = seg%mixg copyUTurnSeg%nodeg = seg%noded ; copyUTurnSeg%noded = seg%nodeg copyUTurnSeg%indCellPg = seg%indCellPd ; copyUTurnSeg%indCellPd = seg%indCellPg copyUTurnSeg%sectg = seg%sectd ; copyUTurnSeg%sectd = seg%sectg copyUTurnSeg%neutronicMixg = seg%neutronicMixd copyUTurnSeg%neutronicMixd = seg%neutronicMixg end function copyUTurnSeg function createArc(cx,cy,r,a,b,mi,me) double precision,intent(in) :: cx,cy,r,a,b integer,intent(in) :: mi,me type(t_segArc) :: createArc createArc%a = a ; createArc%b = b if (isEqual(createArc%a,createArc%b)) then ; createArc%typ = tcer else ; createArc%typ = tarc end if createArc%x = cx ; createArc%y = cy createArc%dx = 0.d0 ; createArc%dy = 0.d0 createArc%r = r createArc%mixg = mi ; createArc%mixd = me createArc%nodeg = fooNode ; createArc%noded = fooNode createArc%indCellPg = 0 ; createArc%indCellPd = 0 createArc%sectg = 0 ; createArc%sectd = 0 createArc%neutronicMixg = 0 ; createArc%neutronicMixd = 0 end function createArc function copyArcExceptOrigin(o,Arc) double precision,intent(in) :: o type(t_segArc),intent(in) :: Arc type(t_segArc) :: copyArcExceptOrigin if (Arc%typ /= tarc) call XABORT("copyArcExceptOrigin : inconsistency") copyArcExceptOrigin = Arc copyArcExceptOrigin%a = o end function copyArcExceptOrigin function copyArcExceptEnd(e,Arc) double precision,intent(in) :: e type(t_segArc),intent(in) :: Arc type(t_segArc) :: copyArcExceptEnd if (Arc%typ /= tarc) call XABORT("copyArcExceptEnd : inconsistency") copyArcExceptEnd = Arc copyArcExceptEnd%b = e end function copyArcExceptEnd function copyArcWithNewAngles(o,e,Arc) double precision,intent(in) :: o, e type(t_segArc),intent(in) :: Arc type(t_segArc) :: copyArcWithNewAngles if (Arc%typ /= tarc) call XABORT("copyArcWithNewAngles : inconsistency") copyArcWithNewAngles = Arc copyArcWithNewAngles%a = o copyArcWithNewAngles%b = e end function copyArcWithNewAngles function createArcFromCircle(o,e,Circ) double precision,intent(in) :: o, e type(t_segArc),intent(in) :: Circ type(t_segArc) :: createArcFromCircle if (Circ%typ /= tcer) call XABORT("createArcFromCircle : inconsistency") createArcFromCircle = Circ createArcFromCircle%typ = tarc createArcFromCircle%a = o createArcFromCircle%b = e end function createArcFromCircle subroutine giveOrigine(sa,xx,yy) type(t_segArc),intent(in) :: sa double precision,intent(out) :: xx,yy select case(sa%typ) case(tseg) ; xx = sa%x ; yy = sa%y case(tarc) ; xx = sa%x+sa%r*cos(sa%a) ; yy=sa%y+sa%r*sin(sa%a) case(tcer) ; xx = sa%x+sa%r ; yy = sa%y end select end subroutine giveOrigine subroutine giveExtremite(sa,xx,yy) type(t_segArc),intent(in) :: sa double precision,intent(out) :: xx,yy select case(sa%typ) case(tseg) ; xx = sa%dx ; yy = sa%dy case(tarc) ; xx = sa%x+sa%r*cos(sa%b) ; yy=sa%y+sa%r*sin(sa%b) case(tcer) ; xx = sa%x+sa%r ; yy = sa%y end select end subroutine giveExtremite subroutine giveFourPointsOnCircle(sa,P1,P2,P3,P4) type(t_segArc), intent(in) :: sa type(t_point), intent(out) :: P1,P2,P3,P4 if (sa%typ/=tcer) call XABORT("giveFourPointsOnCircle : bad use") P1 = t_point(sa%x+sa%r,sa%y) P2 = t_point(sa%x-sa%r,sa%y) P3 = t_point(sa%x,sa%y+sa%r) P4 = t_point(sa%x,sa%y-sa%r) end subroutine giveFourPointsOnCircle subroutine extremitesArc(ar,pt1x,pt1y,pt2x,pt2y) type(t_segArc),intent(in) :: ar double precision,intent(out) :: pt1x,pt1y,pt2x,pt2y pt1x=ar%x+ar%r*cos(ar%a) ; pt1y=ar%y+ar%r*sin(ar%a) pt2x=ar%x+ar%r*cos(ar%b) ; pt2y=ar%y+ar%r*sin(ar%b) end subroutine extremitesArc subroutine MediumPointOnArc(Ar,P) type(t_segArc), intent(in) :: Ar type(t_point), intent(out) :: P double precision :: angle if (Ar%a < Ar %b) then angle = angleNormal((Ar%a + Ar%b) * 0.5) else angle = angleNormal((Ar%a + Ar%b ) *0.5 + pi_c) endif p%x = Ar%x + Ar%r * cos(angle) p%y = Ar%y + Ar%r * sin(angle) end subroutine MediumPointOnArc function estColineaire(sa1,sa2) type(t_segArc),intent(in) :: sa1,sa2 logical :: estColineaire if (sa1%typ/=tseg .or. sa2%typ/=tseg) then estColineaire = .false. else estcolineaire = estcoli(sa1%dx-sa1%x,sa1%dy-sa1%y,sa2%dx-sa2%x,sa2%dy-sa2%y) end if end function estColineaire function isAngleInArc(a,arc) double precision,intent(in) :: a type(t_segArc),intent(in) :: arc integer :: isAngleInArc double precision :: angl !dit si un angle est sur l'intervalle d'un arc : 0 -> pas dessus, ! 1 -> c'est l'origine, 2 -> entre les deux, 3 -> c'est l'extremite if (arc%typ==tcer) then isAngleInArc = 2 return end if angl = angleNormal(a) if (isEqualConst(arc%a,angl)) then ; isAngleInArc = 1 else if (isEqualConst(arc%b,angl)) then ; isAngleInArc = 3 else if (arc%a pas dessus, ! 1 -> c'est l'origine, 2 -> entre les deux, 3 -> c'est l'extremite !ATTENTION :le point doit etre deja sur la droite ou le cercle d'appui if (sa%typ==tseg) then if (isEqualConst(sa%x,ptx) .and. isEqualConst(sa%y,pty)) then isIn = 1 else if (isEqualConst(sa%dx,ptx) .and. isEqualConst(sa%dy,pty)) then isIn = 3 else if (((sa%x-ptx)*(sa%dx-ptx)+(sa%y-pty)*(sa%dy-pty))<0) then isIn = 2 else isIn = 0 end if else angl = calculeAngle(sa%x,sa%y,ptx,pty) isIn = isAngleInArc(angl,sa) end if end function isIn function isIn_V2(ptx,pty,sa) double precision,intent(in) :: ptx,pty type(t_segArc),intent(in) :: sa integer :: isIn_V2 double precision :: angl, tmp !dit si un point est sur un segment ou un arc : 0 -> pas dessus, ! 1 -> c'est l'origine, 2 -> entre les deux, 3 -> c'est l'extremite !ATTENTION :le point doit etre deja sur la droite ou le cercle d'appui if (sa%typ==tseg) then if (isEqualConst(sa%x,ptx) .and. isEqualConst(sa%y,pty)) then isIn_V2 = 1 else if (isEqualConst(sa%dx,ptx) .and. isEqualConst(sa%dy,pty)) then isIn_V2 = 3 else tmp = (sa%x-ptx)*(sa%dx-ptx)+(sa%y-pty)*(sa%dy-pty) if (tmp < dp_0) then isIn_V2 = 2 else if (isEqualConst(tmp,dp_0)) then isIn_V2 = 4 else isIn_V2 = 0 end if end if elseif (sa%typ == tarc) then angl = calculeAngle(sa%x,sa%y,ptx,pty) isIn_V2 = isAngleInArc(angl,sa) else tmp = sqrt((sa%x-ptx)*(sa%x-ptx)+(sa%y-pty)*(sa%y-pty)) if (isEqualConst(tmp,sa%r)) then isIn_V2 = 2 else isIn_V2 = 0 endif end if end function isIn_V2 function isSameWay(sg1,sg2) type(t_segArc),intent(in) :: sg1,sg2 logical :: isSameWay !dit si 2 segments alignes ont meme orientation isSameWay=(((sg1%dx-sg1%x)*(sg2%dx-sg2%x)+(sg1%dy-sg1%y)*(sg2%dy-sg2%y))>0) end function isSameWay subroutine AddSA(SA,szSA) type(t_segArc), intent(in) :: SA integer, intent(inout) :: szSA if (szSA == size(tabSegArc)) call XABORT("AddSA : memory problem") szSA = szSA + 1 tabSegArc(szSA) = SA end subroutine AddSA subroutine Add2SA(SA1,SA2,szSA) type(t_segArc), intent(in) :: SA1,SA2 integer, intent(inout) :: szSA if (szSA == size(tabSegArc)) call XABORT("Add2SA : memory problem") szSA = szSA + 2 tabSegArc(szSA-1) = SA1 tabSegArc(szSA) = SA2 end subroutine Add2SA function ReplaceSA(newSA,iSA) type(t_segArc), intent(in) :: newSA integer, intent(in) :: iSA type(t_segArc) :: ReplaceSA ReplaceSA = newSA tabSegArc(iSA) = newSA end function ReplaceSA subroutine SlideSA(iSA,szSA) integer, intent(in) :: iSA integer, intent(inout) :: szSA tabSegArc(iSA:szSA-1) = tabSegArc(iSA+1:szSA) szSA = szSA - 1 end subroutine SlideSA subroutine AddSC(SC,szSC) type(t_segArc), intent(in) :: SC integer, intent(inout) :: szSC if (szSC == size(tabSegCut)) then write(6,*) szSC,size(tabSegCut) call XABORT("AddSC : memory problem") endif szSC = szSC + 1 tabSegCut(szSC) = SC end subroutine AddSC subroutine Add2SC(SC1,sC2,szSC) type(t_segArc), intent(in) :: SC1,SC2 integer, intent(inout) :: szSC if (szSC == size(tabSegCut)) call XABORT("Add2SC : memory problem") szSC = szSC + 2 tabSegCut(szSC-1) = SC1 tabSegCut(szSC) = SC2 end subroutine Add2SC function ReplaceSC(newSC,iSC) type(t_segArc), intent(in) :: newSC integer, intent(in) :: iSC type(t_segArc) :: ReplaceSC ReplaceSC = newSC tabSegCut(iSC) = newSC end function ReplaceSC subroutine OverlaidSegmentsManagement(n1, n2, iSC, szSC, SC, iSA, szSA, SA) integer, intent(in) :: n1, n2 integer, intent(inout) :: iSC, szSC, iSA, szSA type(t_SegArc), intent(inout) :: SC, SA type(t_segArc) :: SC1, SC2, SC3, SA1, SA2 select case(n1) case(11) ! SC and SA have same origin select case(n2) case(23) SA1 = CopySegExceptOrigin(SC%dx,SC%dy,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(32) ! end of SA internal to SC SC1 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SC2 = CopySegExceptEnd(SA%dx,SA%dy,SC) call AddSC(SC1,szSC) SC = ReplaceSC(SC2,iSC) call SlideSA(iSA,szSA) case(33) ! SC and SA have same end call SlideSA(iSA,szSA) case default call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(12) ! the origin of SC is internal to SA select case (n2) case(21) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SA1 = CopySegExceptOrigin(SC%x,SC%y,SA) SC = ReplaceSC(SC1,iSC) SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSA = iSA + 1 case(23) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SC3 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SA1 = CopySegExceptOrigin(SC%dx,SC%dy,SA) SC = ReplaceSC(SC3,iSC) SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSa = iSA + 1 case(31) SC1 = CopySegExceptOrigin(SA%x,SA%y,SC) SC2 = CopySegWithNewExtrems(SA%dx,SA%dy,SA%x,SA%y,SC) SC = ReplaceSC(SC2,iSC) call AddSC(SC1,szSC) call SlideSA(iSA,szSA) case(32) if (isSameWay(SC,SA)) then SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SC3 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SC = ReplaceSC(SC3,iSC) else SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SC3 = CopySegWithNewExtrems(SA%dx,SA%dy,SA%x,SA%y,SC) SC = ReplaceSC(SC3,iSC) endif call Add2SC(SC1,SC2,szSC) call SlideSA(iSA,szSA) case(33) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SC = ReplaceSC(SC2,iSC) call AddSC(SC1,szSC) call SlideSA(iSA,szSA) case default call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(13) ! the origin of SC is the end of SA select case(n2) case(21) SA1 = CopySegExceptOrigin(SC%x,SC%y,SA) call AddSA(SA1,szSA) iSA = iSA + 1 case(31) SC1 = CopySegWithNewExtrems(SA%dx,SA%dy,SA%x,SA%y,SC) SC = ReplaceSC(SC1,iSC) call SlideSA(iSA,szSA) case(32) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dx,SC) SC = ReplaceSC(SC2,iSC) call AddSC(SC1,szSC) call SlideSA(iSa,szSA) case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(21) select case(n2) case(12) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SA1 = CopySegExceptOrigin(SC%x,SC%y,SA) SC = ReplaceSC(SC1,iSC) call AddSC(SC2,szSC) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(13) SA1 = CopySegExceptOrigin(SC%x,SC%y,SC) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(23) if (isSameWay(SC,SA)) then SA1 = CopySegExceptEnd(SC%x,SC%y,SA) SA2 = CopySegExceptOrigin(SC%dx,SC%dy,SA) else SA1 = CopySegExceptEnd(SC%Dx,SC%Dy,SA) SA2 = CopySegExceptOrigin(SC%x,SC%y,SA) endif SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSA = iSA + 1 case(32) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SA1 = CopySegExceptEnd(SC%x,SC%y,SA) SC = ReplaceSC(SC1,iSC) call AddSC(SC2,szSC) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(33) SA1 = CopySegExceptEnd(SC%x,SC%y,SC) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(23) select case (n2) case(11) SA1 = CopySegExceptOrigin(SC%dx,SC%dy,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(12) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SC = ReplaceSC(SC1,iSC) call AddSC(SC2,szSC) SA1 = CopySegExceptOrigin(SC%dx,SC%dy,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(21) if (isSameWay(SC,SA)) then SA1 = CopySegExceptEnd(SC%x,SC%y,SA) SA2 = CopySegExceptOrigin(SC%dx,SC%dy,SA) else SA1 = CopySegExceptEnd(SC%dx,SC%dy,SA) SA2 = CopySegExceptOrigin(SC%x,SC%y,SA) endif SA = ReplaceSA(SA1,iSA) call AddSA(SA2,szSA) iSA = iSA + 1 case(31) SA1 = CopySegExceptEnd(SC%dx,SC%dy,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case(32) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SA1 = CopySegExceptEnd(SC%dx,SC%dy,SA) SC = ReplaceSC(SC1,iSC) SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSA = iSA + 1 case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(31) select case (n2) case(12) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SC = ReplaceSC(SC1,iSC) call AddSc(SC2,szSC) call SlideSA(iSA,szSA) case(13) SC1 = CopySegWithNewExtrems(SA%dx,SA%dy,SA%x,SA%y,SC) SC = ReplaceSC(SC1,iSC) call SlideSA(iSA,szSA) case(23) SA1 = CopySegExceptEnd(SC%dx,SC%dy,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(32) select case (n2) case(11) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SC = ReplaceSC(SC1,iSC) call AddSC(SC2,szSC) iSA = iSA + 1 case(12) if (isSameWay(SA,SC)) then SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SC3 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SC = ReplaceSC(SC3,iSC) else SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%x,SA%y,SC) SC3 = CopySegWithNewExtrems(SA%dx,SA%dy,SA%x,SA%y,SC) SC = ReplaceSC(SC3,iSC) end if call Add2SC(SC1,SC2,szSC) call SlideSA(iSA,szSA) case(13) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SC = ReplaceSC(SC1,iSC) call AddSC(SC2,szSC) call SlideSA(iSA,szSA) case(21) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SA1 = CopySegExceptOrigin(SC%x,SC%y,SA) SC = ReplaceSC(SC1,iSC) SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSA = iSA + 1 case(23) SC1 = CopySegExceptEnd(SA%dx,SA%dy,SC) SC2 = CopySegExceptOrigin(SA%dx,SA%dy,SC) SA1 = CopySegExceptOrigin(SC%dx,SC%dy,SA) SC = ReplaceSC(SC1,iSC) SA = ReplaceSA(SA1,iSA) call AddSC(SC2,szSC) iSA = iSA + 1 case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case(33) select case (n2) case(11) SC1 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SC = ReplaceSC(SC1,iSC) call SlideSA(iSA,szSA) case(12) SC1 = CopySegExceptEnd(SA%x,SA%y,SC) SC2 = CopySegWithNewExtrems(SA%x,SA%y,SA%dx,SA%dy,SC) SC = ReplaceSC(SC2,iSC) call AddSC(SC1,szSC) call SlideSA(iSA,szSA) case(21) SA1 = CopySegExceptEnd(SC%x,SC%y,SA) SA = ReplaceSA(SA1,iSA) iSA = iSA + 1 case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManagement : impossible circumstance") end select case default write(6,*) "Pb n1,n2",n1,N2 call XABORT("OverlaidSegmentsManaement : ??") end select end subroutine OverlaidSegmentsManagement function turnBackSide(sg) type(t_segArc),intent(inout) :: sg type(t_segArc) :: turnBackSide !tourne le sens du segment en entree turnBackSide=createSeg(sg%dx,sg%dy,sg%x,sg%y,sg%mixd,sg%mixg) turnBackSide%nodeg=sg%noded turnBackSide%noded=sg%nodeg turnBackSide%indCellPg=sg%indCellPd turnBackSide%indCellPd=sg%indCellPg turnBackSide%sectg=sg%sectd turnBackSide%sectd=sg%sectg turnBackSide%neutronicMixg=sg%neutronicMixd turnBackSide%neutronicMixd=sg%neutronicMixg sg=turnBackSide end function turnBackSide subroutine adjustMix(sgToAdjust,sgOfRef) type(t_segArc),intent(inout) :: sgToAdjust type(t_segArc),intent(in) :: sgOfRef !ajuste les milieux a droite et a gauche de sgToAdjust !par rapport a ceux de sgOfRef (sgToAdjust doit etre inclus dans sgOfRef) if (sgOfRef%mixg/=fooMix) then sgToAdjust%mixg = sgOfRef%mixg sgToAdjust%indCellPg = sgOfRef%indCellPg sgToAdjust%sectg = sgOfRef%sectg sgToAdjust%neutronicMixg = sgOfRef%neutronicMixg end if if (sgOfRef%mixd/=fooMix) then sgToAdjust%mixd = sgOfRef%mixd sgToAdjust%indCellPd = sgOfRef%indCellPd sgToAdjust%sectd = sgOfRef%sectd sgToAdjust%neutronicMixd = sgOfRef%neutronicMixd end if end subroutine adjustMix logical function interSgAr_V2(sg,ar, n1, n2, P1, P2) ! Return true only if there is at least an intersection between sg and ar. ! If there is an point of intersection or an tangency point, P1 is this point. ! If there are two intersection points, P1 and P2 are these points. ! N1 and n2 let us pinpoint the point of intersection on both segment and arc. ! If there is no point of intersection, the values for n1 and n2 are 0. type(t_segArc),intent(in) :: sg,ar integer, intent(out):: n1, n2 type (t_point), intent(out) :: P1, P2 double precision :: sox, soy, sfx, sfy, cx, cy, ray, aox, aoy double precision :: prjx, prjy, d, u, v, norm, tmp, Ix, Iy, Jx, Jy integer :: IinS,IinA,JinS,JinA ! Initialization n1 = 0 ; n2 = 0 P1 = t_point(dp_0,dp_0) ; P2 = t_point(dp_0,dp_0) InterSgAr_V2 = .false. ! Coordinates sox = sg%x ; sfx = sg%dx ; soy = sg%y ; sfy = sg%dy cx = ar%x ; cy = ar%y ; ray = ar%r ! Special case : ar is a point if ((ar%typ == tarc.and.isEqualConst(angleNormal(ar%a),angleNormal(ar%b)))) then call giveOrigine(ar,aox,aoy) if (isIn(aox,aoy,sg) /= 0) then interSgAr_V2 =.true. ; n1 = isIn(aox,aoy,sg)+10 P1 = t_point(aox,aoy) endif goto 10 endif ! Special case : sg is a point if (isEqualConst(sox,sfx).and.isEqualConst(soy,sfy)) then d = sqrt((sox-cx)*(sox-cx) + (soy-cy)*(soy-cy)) if (isEqualConst(d,ray)) then interSgAr_V2 =.true. ; n1 = 1+isIn(sox,soy,ar)*10 P1 = t_point(sox,soy) endif goto 10 endif ! Distance Segment-Arc d = distance(cx,cy,sox,soy,sfx,sfy,prjx,prjy) ! Distance analysis if (.not. isEqualConst(d,ar%r)) then if (d < ray) then ! no one, one or two distinct intersection points ! P1 before and P2 after on the straight for segment u = sfx - sox ; v = sfy - soy ; norm = u*u + v*v if (.not.isEqualConst(sqrt(norm),dp_0)) then tmp = sqrt((ray*ray - d*d)/norm) u = tmp*u ; v = tmp*v Ix=prjx-u ; Iy=prjy-v ; IinS=isIn_V2(Ix,Iy,sg) ; IinA=isIn(Ix,Iy,ar) Jx=prjx+u ; Jy=prjy+v ; JinS=isIn_V2(Jx,Jy,sg) ; JinA=isIn(Jx,Jy,ar) if (IinS/=0.and.IinA/=0) then interSgAr_V2 = .true. ; n1 = IinS + 10*IinA ; P1 = t_point(Ix,Iy) if (JinS/=0.and.JinA/=0) then n2 = JinS + 10*JinA ; P2 = t_point(Jx,Jy) goto 10 endif else if (JinS/=0.and.JinA/=0) then interSgAr_V2 = .true. ; n1 = JinS + 10*JinA ; P1 = t_point(Jx,Jy) goto 10 end if else call XABORT("interSgAr_V2 : seg was not detected as a point but its norm is zero") end if end if else ! Tangency point Ix=prjx ; Iy=prjy ; IinS=isIn(Ix,Iy,sg) ; IinA=isIn(Ix,Iy,ar) if ((IinS/=0.and.IinA/=0)) then interSgAr_V2=.true. ; n1 = 100 + 10 * IinA + IinS ; P1 = t_point(Ix,Iy) end if end if ! 10 continue return end function interSgAr_V2 function interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y) type(t_segArc),intent(in) :: sg,ar double precision,intent(out) :: pt1x,pt1y,pt2x,pt2y integer :: interSgAr double precision :: prjx,prjy,d double precision :: u,v,tmp,Ix,Iy,Jx,Jy integer :: IinS,IinA,JinS,JinA !renvoie le nombre de points d'intersection entre sg et ar, !et leurs coordonnees eventuelles !Si un seul point d'intersection, renvoi 1 si pt==I et -1 si pt==J !Si l'arc et le segment sont tangents, renvoie 2, et pt1=pt2 interSgAr = 0 pt1x=0.d0 ; pt1y=0.d0 ; pt2x=1.d0 ; pt2y=1.d0 d = distance(ar%x,ar%y,sg%x,sg%y,sg%dx,sg%dy,prjx,prjy) if (.not. isEqualConst(d,ar%r)) then if (d>ar%r) then interSgAr = 0 else !2 points d'intersection distincts I et J avec !I avant J sur la droite d'appui du segment u=sg%dx-sg%x ; v=sg%dy-sg%y tmp=sqrt((ar%r*ar%r-d*d)/(u*u+v*v)) u=tmp*u ; v=tmp*v Ix=prjx-u ; Iy=prjy-v ; IinS=isIn(Ix,Iy,sg) ; IinA=isIn(Ix,Iy,ar) Jx=prjx+u ; Jy=prjy+v ; JinS=isIn(Jx,Jy,sg) ; JinA=isIn(Jx,Jy,ar) if ((IinS/=0.and.IinA/=0).and.(IinS==2.or.IinA==2)) then interSgAr=1 ; pt1x=Ix ; pt1y=Iy end if if ((JinS/=0.and.JinA/=0).and.(JinS==2.or.JinA==2)) then if (interSgAr==1) then interSgAr=2 ; pt2x=Jx ; pt2y=Jy else interSgAr=-1 ; pt1x=Jx ; pt1y=Jy end if end if end if else Ix=prjx ; Iy=prjy ; IinS=isIn(Ix,Iy,sg) ; IinA=isIn(Ix,Iy,ar) if ((IinS/=0.and.IinA/=0).and.(IinS==2.or.IinA==2)) then interSgAr=2 ; pt1x=Ix ; pt1y=Iy ; pt2x=Ix ; pt2y=Iy end if end if end function interSgAr logical function InterSgSg_V2(sg1,sg2,n1, n2, P1,P2) type(t_segArc),intent(in) :: sg1,sg2 integer, intent(inout) :: n1, n2 type(t_point), intent(out) :: P1, P2 double precision :: xo1, xf1, xo2, xf2, yo1, yf1, yo2, yf2 double precision :: min1, min2, max1, max2, m, prodvect double precision :: a1, a2, b1, b2, y ! Initialization n1 = 0 ; n2 = 0 P1 = t_point(dp_0,dp_0) ; P2 = t_point(dp_0,dp_0) InterSgSg_V2 = .false. ! Coordinates xo1=sg1%x ; xf1 =sg1%dx ; yo1 = sg1%y ; yf1 = sg1%dy xo2=sg2%x ; xf2 =sg2%dx ; yo2 = sg2%y ; yf2 = sg2%dy ! Very Special case : point instead of segment min1 = sqrt((xf1-xo1)*(xf1-xo1)+(yf1-yo1)*(yf1-yo1)) min2 = sqrt((xf2-xo2)*(xf2-xo2)+(yf2-yo2)*(yf2-yo2)) if (isEqualConst(min1,dp_0)) then prodvect = (xf2-xo2)*(yo1-yo2)-(yf2-yo2)*(xo1-xo2) if (isEqualConst(prodvect,dp_0) .and. (isIn(xo1,yo1,Sg2)/=0)) then InterSgSg_V2 = .true. ; n1 = 1 ; P1 = t_point(xo1,yo1) endif goto 10 elseif(isEqualConst(min2,dp_0)) then prodvect = (xf1-xo1)*(yo2-yo1)-(yf1-yo1)*(xo2-xo1) if (isEqualConst(prodvect,dp_0) .and. (isIn(xo2,yo2,Sg1)/=0)) then InterSgSg_V2 = .true. ; n1 = 1 ; P1 = t_point(xo2,yo2) endif goto 10 else ! Special case : two vertical segments if (IsEqualConst(xo1,xf1).and.IsEqualConst(xo2,xf2)) then if (IsEqualConst(xo1,xo2)) then min1 = min(yo1,yf1) ; max1 = max(yo1,yf1) min2 = min(yo2,yf2) ; max2 = max(yo2,yf2) if (isEqualConst(min1,max2).or.isEqualConst(min2,max1)) then InterSgSg_V2 = .true. ; n1 = 1 if (isEqualConst(min1,max2)) P1 = t_point(xo1,min1) if (isEqualConst(max1,min2)) P1 = t_point(xo2,min2) elseif (max2 < min1 .or. max1 < min2) then InterSgSg_V2 = .false. else InterSgSg_V2 = .true. n1 = 1 ; P1 = t_point(xo1,max(min1,min2)) n2 = 1 ; P2 = t_point(xo1,min(max1,max2)) endif endif else ! Special case : only one vertical segment if (isEqualConst(xo1,xf1).and..not.isEqualConst(xo2,xf2)) then a2 = (yf2 - yo2) / (xf2 - xo2) ; b2 = yo2 - a2 * xo2 if ((isEqualConst(xo1,min(xo2,xf2)).or.xo1>min(xo2,xf2)) .and.& (isEqualConst(xo2,max(xo2,xf2)).or.xo1min(yo1,yf1)).and.& (isEqualConst(y,max(yo1,yf1)).or.ymin(yo2,yf2)).and.& (isEqualConst(y,max(yo2,yf2)).or.ymin(xo1,xf1) ).and.& (isEqualConst(xo2,max(xo1,xf1)).or.xo2min(yo1,yf1)).and.& (isEqualConst(y,max(yo1,yf1)).or.ymin(yo2,yf2)).and.& (isEqualConst(y,max(yo2,yf2)).or.ymin1) .and. & (isEqualConst(m,max1).or.mmin2) .and. & (isEqualConst(m,max2).or.m(d2+rc2)).or.(rc2>(d2+rr2)) if (d>(maxD+epsilon).or.d<(minD-epsilon)) return !pas d'intersection if (isEqualConst(d,maxD).or.isEqualConst(d,minD)) then !tangence entre les deux cercles d'appuis aClus = calculeAngle(clus%x,clus%y,ring%x,ring%y) if (revert) aClus = angleNormal(aClus - pi_c) inRing = isIn(clus%x+clus%r*cos(aClus),clus%y+clus%r*sin(aClus),ring) aRing = calculeAngle(ring%x,ring%y,clus%x,clus%y) inClus = isIn(ring%x+ring%r*cos(aRing),ring%y+ring%r*sin(aRing),clus) if (inRing==2 .and. inClus==2) interArcArc = 1 !dans les arcs return end if !il y a intersection reelle des deux cercles d'appuis ! l'angle a1 adjacent a r1 est determinable a l'aide de la formule ! cos(a1) = (d\B2+r1\B2-r2\B2)/(2*d*r1) (idem mutatis-mutandis pour a2) aC = acos((d2+rc2-rr2)/(2*d*clus%r)) aClus = angleNormal(calculeAngle(clus%x,clus%y,ring%x,ring%y) - aC) bClus = angleNormal(calculeAngle(clus%x,clus%y,ring%x,ring%y) + aC) aR = acos((d2+rr2-rc2)/(2*d*ring%r)) aRing = angleNormal(calculeAngle(ring%x,ring%y,clus%x,clus%y) - aR) bRing = angleNormal(calculeAngle(ring%x,ring%y,clus%x,clus%y) + aR) aIn = (isAngleInArc(aClus,clus)==2) .and. (isAngleInArc(aRing,ring)==2) bIn = (isAngleInArc(bClus,clus)==2) .and. (isAngleInArc(bRing,ring)==2) if (aIn .and. bIn) then interArcArc = 2 else if (aIn) then interArcArc = 1 else if (bIn) then interArcArc = 1 aClus = bClus aRing = bRing end if end function interArcArc function giveExtremalsAngles(sa,cx,cy,ao,ae) type(t_segArc),intent(in) :: sa double precision,intent(in) :: cx,cy double precision,intent(out) :: ao,ae logical :: giveExtremalsAngles double precision :: ox,oy,ex,ey select case(sa%typ) case(tcer) ! on sort avec une valeur false ao = 0.d0 ; ae=0.d0 ; giveExtremalsAngles = .false. return case(tarc) call extremitesArc(sa,ox,oy,ex,ey) case(tseg) ox = sa%x ; oy = sa%y ; ex = sa%dx ; ey = sa%dy end select ao = angleNormal(calculeAngle(cx,cy,ox,oy)) ae = angleNormal(calculeAngle(cx,cy,ex,ey)) giveExtremalsAngles = .not.( (isEqualConst(cx,ox).and.isEqualConst(cy,oy))& .or.(isEqualConst(cx,ex).and.isEqualConst(cy,ey)) ) end function giveExtremalsAngles subroutine drawSegArc(fileNbr,szSA,drawMix,zoomx,zoomy) integer,intent(in) :: fileNbr,szSA,drawMix real,intent(in) :: zoomx(2),zoomy(2) type(t_segArc) :: sa integer :: i real :: cx,cy,angl,lx,ly,tailleNbr,delx,dely g_psp_bBoxXmin = 1.e10 g_psp_bBoxYmin = 1.e10 g_psp_bBoxXmax = -1.e10 g_psp_bBoxYmax = -1.e10 !recuperation des donnees permettants le centrage do i = 1,szSA sa = tabSegArc(i) if (sa%typ==tseg) then g_psp_bBoxXmin = min(real(sa%x),real(sa%dx),g_psp_bBoxXmin) g_psp_bBoxYmin = min(real(sa%y),real(sa%dy),g_psp_bBoxYmin) g_psp_bBoxXmax = max(real(sa%x),real(sa%dx),g_psp_bBoxXmax) g_psp_bBoxYmax = max(real(sa%y),real(sa%dy),g_psp_bBoxYmax) else g_psp_bBoxXmin = min(real(sa%x-sa%r),g_psp_bBoxXmin) g_psp_bBoxYmin = min(real(sa%y-sa%r),g_psp_bBoxYmin) g_psp_bBoxXmax = max(real(sa%x+sa%r),g_psp_bBoxXmax) g_psp_bBoxYmax = max(real(sa%y+sa%r),g_psp_bBoxYmax) end if end do delx = g_psp_bBoxXmax - g_psp_bBoxXmin dely = g_psp_bBoxYmax - g_psp_bBoxYmin g_psp_bBoxXmax = g_psp_bBoxXmin + zoomx(2)*delx g_psp_bBoxXmin = g_psp_bBoxXmin + zoomx(1)*delx g_psp_bBoxYmax = g_psp_bBoxYmin + zoomy(2)*dely g_psp_bBoxYmin = g_psp_bBoxYmin + zoomy(1)*dely lx = 0.1 * ( g_psp_bBoxXmax - g_psp_bBoxXmin ) ly = 0.1 * ( g_psp_bBoxYmax - g_psp_bBoxYmin ) write(*,10) g_psp_bBoxXmin,g_psp_bBoxXmax,g_psp_bBoxYmin,g_psp_bBoxYmax 10 format(' g2s_segArc: plot domain=(',1p,e14.7,':',e14.7,',',e14.7,':',e14.7,')') g_psp_bBoxXmin = g_psp_bBoxXmin - lx g_psp_bBoxYmin = g_psp_bBoxYmin - ly g_psp_bBoxXmax = g_psp_bBoxXmax + lx g_psp_bBoxYmax = g_psp_bBoxYmax + ly tailleNbr = min(595./(g_psp_bBoxXmax-g_psp_bBoxXmin), & 842./(g_psp_bBoxYmax-g_psp_bBoxYmin)) tailleNbr = 3.6 / tailleNbr !impression call psinit(fileNbr,.true.) do i = 1,szSA sa = tabSegArc(i) if (sa%typ==tseg) then call line(sa%x,sa%y,sa%dx,sa%dy) cx=real((sa%dx+sa%x)*0.5d0) ; cy=real((sa%dy+sa%y)*0.5d0) angl = real(calculeAngle(sa%x,sa%y,sa%dx,sa%dy)*rad2deg-90.d0) if (drawMix==1) then call keknum(cx,cy,tailleNbr,real(sa%nodeg),angl,-1,2) call keknum(cx,cy,tailleNbr,real(sa%noded),angl,-1,0) else if (drawMix==2) then call keknum(cx,cy,tailleNbr,real(sa%neutronicMixg),angl,-1,2) call keknum(cx,cy,tailleNbr,real(sa%neutronicMixd),angl,-1,0) else if (drawMix==3) then call keknum(cx,cy,tailleNbr,real(i),angl,-1,2) call keknum(cx,cy,tailleNbr,real(i),angl,-1,0) end if else if (sa%typ==tcer) then call arc(sa%x,sa%y,sa%r,0.d0,180.d0) call arc(sa%x,sa%y,sa%r,180.d0,360.d0) cx=real(sa%x+sa%r) ; cy=real(sa%y) if (drawMix==1) then call keknum(cx,cy,tailleNbr,real(sa%nodeg),0.,-1,2) call keknum(cx,cy,tailleNbr,real(sa%noded),0.,-1,0) else if (drawMix==2) then call keknum(cx,cy,tailleNbr,real(sa%neutronicMixg),0.,-1,2) call keknum(cx,cy,tailleNbr,real(sa%neutronicMixd),0.,-1,0) else if (drawMix==3) then call keknum(cx,cy,tailleNbr,real(i),0.,-1,2) call keknum(cx,cy,tailleNbr,real(i),0.,-1,0) end if else call arc(sa%x,sa%y,sa%r,sa%a*rad2deg,sa%b*rad2deg) if (sa%b>sa%a) then angl=real((sa%b+sa%a)*0.5d0) else angl=real((sa%b+sa%a)*0.5d0+pi_c) end if cx=real(sa%x+cos(angl)*sa%r) ; cy=real(sa%y+sin(angl)*sa%r) angl=real(angl*rad2deg) if (drawMix==1) then call keknum(cx,cy,tailleNbr,real(sa%nodeg),angl,-1,2) call keknum(cx,cy,tailleNbr,real(sa%noded),angl,-1,0) else if (drawMix==2) then call keknum(cx,cy,tailleNbr,real(sa%neutronicMixg),angl,-1,2) call keknum(cx,cy,tailleNbr,real(sa%neutronicMixd),angl,-1,0) else if (drawMix==3) then call keknum(cx,cy,tailleNbr,real(i),angl,-1,2) call keknum(cx,cy,tailleNbr,real(i),angl,-1,0) end if end if end do call plotnd() end subroutine drawSegArc subroutine coupeCercle(indDeb,szSA,nbSeg,typ) integer,intent(in) :: indDeb,nbSeg,typ integer,intent(inout) :: szSA type(t_segArc) :: sg,ar double precision :: pt1x,pt1y,pt2x,pt2y,angl1,angl2,midx,midy,tmp integer :: i,j,k,tailleAv,tailleAp,nbPtInter type(segArcArrayBis),dimension(:),allocatable :: tmpTabSegArc ! pour programmation defensive integer :: taille_table_temp !RQ: indDeb est l'indice du dernier element du tableau qui ! ne doit pas etre pris en compte. nbSeg est le nombre ! de segments paralleles aux axes, et typ correspond au type (3,4,6,ou...) ! CS-IV : Correction des commentaires d'origine ! CS-IV : indDeb : indice du dernier element a ne pas prendre en compte ! CS-IV : (point de depart) ! CS-IV : nbSeg : nombre de SA a prendre en charge ! CS-IV : szSA : nombre total de SA dans le tableau !preparation du tableau de travail tailleAp = szSA-indDeb ! CS-IV : pourquoi ce 25 ?(on considere que la decoupe ne va pas ! CS-IV : accroitre le nombre de segments de plus de 25 fois ! CS-IV : ce coef devrait etre reflechi et dependre du nbr de cercles et ! CS-IV : et de segments ! pour programmation defensive taille_table_temp = 25*tailleAp allocate(tmpTabSegArc(1:taille_table_temp),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: coupeCercle (1) => allocation pb") !copie des segments et arcs construits dans le tableau !temporaire ! CS-IV : et suppression dans le tableau initial tmpTabSegArc(1:tailleAp)%sa = tabSegArc(indDeb+1:indDeb+tailleAp) tmpTabSegArc(1:tailleAp)%keep = .true. szSa = indDeb ! on coupe tout les elements possible presents, en ne gardant pas ceux qui ! sont a couper, mais seulement les morceaux obtenus ! pour i->segment, pour j->arc ou cercle do tailleAv = tailleAp i = 0 do i = i+1 if (.not. tmpTabSegArc(i)%keep) then if (i/=tailleAv ) cycle exit end if sg = tmpTabSegArc(i)%sa if (sg%typ/=tseg) then if (i/=tailleAv ) cycle exit end if ! CS-IV : le SA tmpTabSegArc(i) est un segment do j = 1,tailleAv if ((j==i).or.(.not. tmpTabSegArc(j)%keep)) cycle ar = tmpTabSegArc(j)%sa if (ar%typ==tseg) cycle ! CS-IV : le SA tmpTabSegArc(j) est un arc nbPtInter = interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y) if (nbPtInter==2) then !!deux points d'intersection : ! CS-IV : on obtient 3 segments et 2 arcs tmpTabSegArc(i)%keep=.false. ; tmpTabSegArc(j)%keep=.false. ! => 3 segments (qui gardent les caracteristiques sect*) do k = 1,3 tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (1)") tmpTabSegArc(tailleAp)%sa=sg end do tmpTabSegArc(tailleAp-2)%sa%dx=pt1x tmpTabSegArc(tailleAp-2)%sa%dy=pt1y tmpTabSegArc(tailleAp-1)%sa%x=pt1x tmpTabSegArc(tailleAp-1)%sa%y=pt1y tmpTabSegArc(tailleAp-1)%sa%dx=pt2x tmpTabSegArc(tailleAp-1)%sa%dy=pt2y tmpTabSegArc(tailleAp)%sa%x=pt2x tmpTabSegArc(tailleAp)%sa%y=pt2y tmpTabSegArc(tailleAp-2)%sa%mixg=ar%mixd tmpTabSegArc(tailleAp-2)%sa%mixd=ar%mixd tmpTabSegArc(tailleAp-1)%sa%mixg=ar%mixg tmpTabSegArc(tailleAp-1)%sa%mixd=ar%mixg tmpTabSegArc(tailleAp)%sa%mixg=ar%mixd tmpTabSegArc(tailleAp)%sa%mixd=ar%mixd ! elimination d'un eventuel segment nul do k = 0,2 tmpTabSegArc(tailleAp-k)%keep= .not. & & ( (isEqual( tmpTabSegArc(tailleAp-k)%sa%x , & & tmpTabSegArc(tailleAp-k)%sa%dx) ) .and. & & (isEqual( tmpTabSegArc(tailleAp-k)%sa%y , & & tmpTabSegArc(tailleAp-k)%sa%dy) ) ) end do angl1=calculeAngle(ar%x,ar%y,pt1x,pt1y) angl2=calculeAngle(ar%x,ar%y,pt2x,pt2y) !travail sur les arcs if (ar%typ==tcer) then !c'etait un cercle complet if (isEqualAngl(angl1,angl2)) then !tangence tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (2)") tmpTabSegArc(tailleAp)%sa=ar tmpTabSegArc(tailleAp)%keep=.true. tmpTabSegArc(tailleAp)%sa%typ=tarc tmpTabSegArc(tailleAp)%sa%a=angl1 tmpTabSegArc(tailleAp)%sa%b=angl1 !positionnement des sect* if (estAGauche(sg%x,sg%y,sg%dx,sg%dy,ar%x,ar%y)) then !cercle completement a gauche de sg tmpTabSegArc(tailleAp)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp)%sa%sectd=sg%sectg else !cercle completement a droite de sg tmpTabSegArc(tailleAp)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp)%sa%sectd=sg%sectd end if else !pas tangence !un arc de chaque cote de sg (1\B0 a drt, 2\B0 a gch) do k = 1,2 tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (3)") tmpTabSegArc(tailleAp)%sa=ar tmpTabSegArc(tailleAp)%keep=.true. tmpTabSegArc(tailleAp)%sa%typ=tarc end do tmpTabSegArc(tailleAp-1)%sa%a=angl1 tmpTabSegArc(tailleAp-1)%sa%b=angl2 tmpTabSegArc(tailleAp-1)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp-1)%sa%sectd=sg%sectd tmpTabSegArc(tailleAp)%sa%a=angl2 tmpTabSegArc(tailleAp)%sa%b=angl1 tmpTabSegArc(tailleAp)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp)%sa%sectd=sg%sectg end if else !c'etait deja un arc de cercle do k = 1,3 tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (4)") tmpTabSegArc(tailleAp)%sa=ar end do tmpTabSegArc(tailleAp-2)%sa%b=angl1 !on test si on est bien dans le bon sens if (isIn(pt2x,pt2y,tmpTabSegArc(tailleAp-2)%sa)==2) then !on est a l'envers (ouverture a droite) tmpTabSegArc(tailleAp-2)%sa%b=angl2 tmpTabSegArc(tailleAp-1)%sa%a=angl2 tmpTabSegArc(tailleAp-1)%sa%b=angl1 tmpTabSegArc(tailleAp)%sa%a=angl1 !ajustement des sect* tmpTabSegArc(tailleAp-2)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp-2)%sa%sectd=sg%sectd tmpTabSegArc(tailleAp-1)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp-1)%sa%sectd=sg%sectg tmpTabSegArc(tailleAp)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp)%sa%sectd=sg%sectd else !on est a l'endroit (ouverture a gauche) tmpTabSegArc(tailleAp-1)%sa%a=angl1 tmpTabSegArc(tailleAp-1)%sa%b=angl2 tmpTabSegArc(tailleAp)%sa%a=angl2 !ajustement des sect* tmpTabSegArc(tailleAp-2)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp-2)%sa%sectd=sg%sectg tmpTabSegArc(tailleAp-1)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp-1)%sa%sectd=sg%sectd tmpTabSegArc(tailleAp)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp)%sa%sectd=sg%sectg end if do k = 0,2 tmpTabSegArc(tailleAp-k)%keep=.not. & & ( isEqualAngl( tmpTabSegArc(tailleAp-k)%sa%a , & & tmpTabSegArc(tailleAp-k)%sa%b ) ) end do end if i=0 exit else if (abs(nbPtInter)==1) then !! un seul point d'intersection ! CS-IV : on obtient 2 segments & 2 arcs tmpTabSegArc(i)%keep=.false. tmpTabSegArc(j)%keep=.false. do k = 1,2 tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (5)") tmpTabSegArc(tailleAp)%sa=sg end do tmpTabSegArc(tailleAp-1)%sa%dx=pt1x tmpTabSegArc(tailleAp-1)%sa%dy=pt1y tmpTabSegArc(tailleAp)%sa%x=pt1x tmpTabSegArc(tailleAp)%sa%y=pt1y if (nbPtInter==1) then !intersection en I (origne de sg a l'exterieur de ar) tmpTabSegArc(tailleAp-1)%sa%mixg=ar%mixd tmpTabSegArc(tailleAp-1)%sa%mixd=ar%mixd tmpTabSegArc(tailleAp)%sa%mixg=ar%mixg tmpTabSegArc(tailleAp)%sa%mixd=ar%mixg else !intersection en J (origne de sg a l'interieur de ar) tmpTabSegArc(tailleAp-1)%sa%mixg=ar%mixg tmpTabSegArc(tailleAp-1)%sa%mixd=ar%mixg tmpTabSegArc(tailleAp)%sa%mixg=ar%mixd tmpTabSegArc(tailleAp)%sa%mixd=ar%mixd end if do k = 0,1 tmpTabSegArc(tailleAp-k)%keep=.not. & & ( ( isEqual( tmpTabSegArc(tailleAp-k)%sa%x , & & tmpTabSegArc(tailleAp-k)%sa%dx ) ) .and.& & ( isEqual( tmpTabSegArc(tailleAp-k)%sa%y , & & tmpTabSegArc(tailleAp-k)%sa%dy ) ) ) end do angl1=calculeAngle(ar%x,ar%y,pt1x,pt1y) if (ar%typ==tcer) then !ar etait un cercle tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (6)") tmpTabSegArc(tailleAp)%sa=ar tmpTabSegArc(tailleAp)%sa%typ=tarc tmpTabSegArc(tailleAp)%sa%a=angl1 tmpTabSegArc(tailleAp)%sa%b=angl1 tmpTabSegArc(tailleAp)%keep=.true. !pas d'ajustement des sect*, car l'arc n'est ni d'un cote !ni de l'autre else do k = 1,2 tailleAp=tailleAp+1 ! pour programmation defensive if (tailleAP > taille_table_temp) & call XABORT("G2S : memory problem in routine CoupeCercle (7)") tmpTabSegArc(tailleAp)%sa=ar end do tmpTabSegArc(tailleAp-1)%sa%b=angl1 tmpTabSegArc(tailleAp)%sa%a=angl1 !ajustement des sect* if (nbPtInter==1) then !intersection en I (origne de sg a l'exterieur de ar) tmpTabSegArc(tailleAp-1)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp-1)%sa%sectd=sg%sectg tmpTabSegArc(tailleAp)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp)%sa%sectd=sg%sectd else !intersection en J (origne de sg a l'interieur de ar) tmpTabSegArc(tailleAp-1)%sa%sectg=sg%sectd tmpTabSegArc(tailleAp-1)%sa%sectd=sg%sectd tmpTabSegArc(tailleAp)%sa%sectg=sg%sectg tmpTabSegArc(tailleAp)%sa%sectd=sg%sectg end if do k = 0,1 tmpTabSegArc(tailleAp-k)%keep=.not. & & ( isEqualAngl( tmpTabSegArc(tailleAp-k)%sa%a , & & tmpTabSegArc(tailleAp-k)%sa%b ) ) end do end if i=0 exit end if end do if (i==tailleAv) exit end do if ((tailleAp==tailleAv) .and. (i==tailleAv)) exit end do !on elimine les elements en dehors des limites de la cellule, et on remet le !milieu exterieur a la cellule a fooMix !Remarque : les elements tmpTabSegArc(1:nbSeg)%sa sont les segments !initiaux paralleles aux bords du domaine. Ils peuvent etre a jeter !(keep==.false.), mais ils servent tout de meme de reference. Pour savoir !si ils sont au bord de la cellule, on teste si le milieu est fooMix. !Si oui, on elimine tous les elements qui sont de ce cote. if (typ==tRec .or. typ==tHex) then do i = 1,nbSeg sg = tmpTabSegArc(i)%sa do j = 1,tailleAp ar = tmpTabSegArc(j)%sa if (ar%typ==tseg) then if (.not. & (estColi(sg%dx-sg%x,sg%dy-sg%y,ar%dx-ar%x,ar%dy-ar%y) & .and. pointsAlignes(sg%x,sg%y,sg%dx,sg%dy,ar%x,ar%y)) & ) cycle if (sg%mixd==fooMix) then tmpTabSegArc(j)%sa%mixd=fooMix else if (sg%mixg==fooMix) then tmpTabSegArc(j)%sa%mixg=fooMix end if cycle end if if ((j==i) .or. (.not. tmpTabSegArc(j)%keep)) cycle if (ar%b>ar%a) then angl1=(ar%b+ar%a)*0.5d0 else angl1=(ar%b+ar%a)*0.5d0+pi_c end if !point milieu de l'arc midx=ar%x+ar%r*cos(angl1) midy=ar%y+ar%r*sin(angl1) !sinus de l'angle (OE,OM) tmp = sin( calculeAngle(sg%x,sg%y,midx,midy) - & & calculeAngle(sg%x,sg%y,sg%dx,sg%dy) ) if (sg%mixd==fooMix) then !on garde a gauche tmpTabSegArc(j)%keep=(tmp>0.).or.(isEqualAngl(ar%a,ar%b).and. & & estAGauche(sg%x,sg%y,sg%dx,sg%dy,ar%x,ar%y)) else if (sg%mixg==fooMix) then !on garde a droite tmpTabSegArc(j)%keep=(tmp<0.).or.(isEqualAngl(ar%a,ar%b).and. & & estADroite(sg%x,sg%y,sg%dx,sg%dy,ar%x,ar%y)) end if end do end do else call XABORT("G2S : type of geometrie not supported") end if !elimination des segments nuls eventuellement restants do i = 1,tailleAp sg = tmpTabSegArc(i)%sa if (sg%typ==tseg) tmpTabSegArc(i)%keep=tmpTabSegArc(i)%keep .and. & .not.(isEqual(tmpTabSegArc(i)%sa%x,tmpTabSegArc(i)%sa%dx) .and. & isEqual(tmpTabSegArc(i)%sa%y,tmpTabSegArc(i)%sa%dy)) end do !on recupere dans le tableau general les elements restants do i = 1,tailleAp if (.not. tmpTabSegArc(i)%keep) cycle szSa = szSa+1 tabSegArc(szSa) = tmpTabSegArc(i)%sa end do !on nettoie le tableau temporaire deallocate(tmpTabSegArc) end subroutine coupeCercle subroutine splitSegsForSect(indDeb,szSA) integer,intent(in) :: indDeb integer,intent(inout) :: szSA type(t_segArc) :: sgi,sgj integer :: i,j ! pour programmation defensive integer :: taille_table_tabSegArc !! permet de spliter les segments en bordure de domaine, !! losqu'ils sont interceptes par les segments delimitants !! la sectorisation. Si on a intersection, c'est toujours a !! l'interieur d'un segment de bordure, et au bout d'un segment !! de secteur. -->^-----> <-----^<-- !! Dans ce cas, on split ainsi: (1) i | szSA ou: (2) szSA | i ! pour programmation defensive taille_table_tabSegArc = size(tabSegArc) i = indDeb do i = i + 1 if (i>szSA) exit sgi=tabSegArc(i) if (sgi%typ/=tseg) cycle do j = indDeb+1,szSA if (j==i) cycle sgj=tabSegArc(j) if (sgj%typ/=tseg) cycle if (PointsAlignes(sgi%x,sgi%y,sgi%dx,sgi%dy,sgj%dx,sgj%dy).and.& (isIn(sgj%dx,sgj%dy,sgi)==2).and. & .not.PointsAlignes(sgi%x,sgi%y,sgi%dx,sgi%dy,sgj%x,sgj%y)) then szSA=szSA+1 ! pour programmation defensive if (szSA > taille_table_tabSegArc) & call XABORT("G2S : memory problem in routine splitSegsForSect") tabSegArc(szSA) = sgi tabSegArc(i)%dx = sgj%dx ; tabSegArc(i)%dy = sgj%dy tabSegArc(szSA)%x = sgj%dx ; tabSegArc(szSA)%y = sgj%dy if (estAGaucheStrict(sgj%x,sgj%y,sgj%dx,sgj%dy,sgi%x,sgi%y)) then !! cas (1) tabSegArc(i)%sectg = sgj%sectg tabSegArc(i)%sectd = sgj%sectg tabSegArc(szSA)%sectg = sgj%sectd tabSegArc(szSA)%sectd = sgj%sectd else !! cas (2) tabSegArc(i)%sectg = sgj%sectd tabSegArc(i)%sectd = sgj%sectd tabSegArc(szSA)%sectg = sgj%sectg tabSegArc(szSA)%sectd = sgj%sectg end if sgi = tabSegArc(i) end if end do end do end subroutine splitSegsForSect subroutine majSectori(indDeb,szSA,sectori,typGeo,cx,cy) double precision, parameter :: PI = 3.141592653589793d0 integer,intent(in) :: indDeb,szSA,sectori,typGeo double precision,intent(in) :: cx,cy !met a jour les secteurs pour les elements double precision :: ao,ae !angles origine et extremite integer :: i,j,nbSect,oIn,eIn logical :: goodAngles double precision,dimension(9) :: limit if ((typGeo/=tRec).and.(typGeo/=tHex)) & call XABORT("G2S: internal error with typGeo in subroutine majSectori") !creation des zones angulaires nbSect = 0 select case(sectori) case(S_not) return !rien a faire case(S_X_tot) if (typGeo == tRec) then nbSect = 4 limit(1) = pi_2_c*5.d-1 ; limit(2) = pi_2_c*1.5d0 limit(3) = -limit(2) ; limit(4) = -limit(1) limit(5) = limit(1) else nbSect = 6 limit(1) = 0.d0 ; limit(2) = pi_3_c ; limit(3) = pi_3_c*2.d0 limit(4) = pi_c ; limit(5) = -limit(3) ; limit(6) = -limit(2) limit(7) = limit(1) end if case(S_T_tot) nbSect = 4 limit(1) = 0.d0 ; limit(2) = pi_2_c limit(3) = pi_c ; limit(4) = -limit(2) limit(5) = limit(1) case(S_TX_tot) nbSect = 8 limit(1) = 0.d0 ; limit(2) = pi_2_c*5.d-1 limit(3) = pi_2_c ; limit(4) = pi_2_c*1.5d0 limit(5) = pi_c ; limit(6) = -limit(4) limit(7) = -limit(3) ; limit(8) = -limit(2) limit(9) = limit(1) case(S_TXS_tot) nbsect = 8 limit(1) = PI/8. do i=1, 7 limit(i+1) = limit(i) + PI/4. enddo limit(9) = limit(1) case(S_WM_tot) nbsect = 8 limit(1) = PI/8. do i=1, 7 limit(i+1) = limit(i) + PI/4. enddo limit(9) = limit(1) case default call XABORT("G2S: internal error with sectori in subroutine majSectori") end select !traitement des elements do i = indDeb+1,szSA if (tabSegArc(i)%typ == tcer) cycle !on ne s'interesse pas aux cercles !on recupere les angles des bouts goodAngles = giveExtremalsAngles(tabSegArc(i),cx,cy,ao,ae) if (.not.goodAngles) cycle !une des extremites est le centre ! => on ne s'en occupe pas do j = 1,nbSect !boucle sur les secteurs oIn = isAngleInInterval(ao,limit(j),limit(j+1)) eIn = isAngleInInterval(ae,limit(j),limit(j+1)) if ((oIn==0) .or. (eIn==0)) then !pas dans le secteur cycle !on passe au secteur suivant else if ((oIn==2) .or. (eIn==2)) then !on est completement dans le secteur tabSegArc(i)%sectg = j tabSegArc(i)%sectd = j exit !on sort de la boucle en j else if (oIn/=eIn) then !element occupant le secteur angulaire total => idem cas precedent tabSegArc(i)%sectg = j tabSegArc(i)%sectd = j exit !on sort de la boucle en j else !les 2 extremites sont sur le meme separateur ! => on ne s'en occupe pas exit !on sort de la boucle en j end if end do end do end subroutine majSectori subroutine addSegsAndClean(sizeSA) integer,intent(inout) :: sizeSA integer :: i,j,sizeTmp integer :: oj,ej !test l'appartenance a sgi des bouts de sgj type(t_segArc) :: sgi,sgj type(segArcArrayBis),dimension(:),allocatable :: tmpTab ! pour programmation defensive integer :: taille_table_temp, taille_table_tabSegArc ! copie dans un tableau temporaire des segments et nettoyage du tableau ! global ! pour programmation defensive taille_table_temp = sizeSA*10 taille_table_tabSegArc = size(tabSegArc) allocate(tmpTab(taille_table_temp),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: addSegsAndClean(1) => allocation pb") tmpTab(1:sizeSA)%sa=tabSegArc(1:sizeSA) tmpTab(1:sizeSA)%keep=.true. sizeTmp=sizeSA i=0 do i = i+1 if (i>sizeTmp) exit if (.not. tmpTab(i)%keep) cycle if (tmpTab(i)%sa%typ/=tseg) cycle sgi=tmpTab(i)%sa !c'est bien un segment a garder pour le moment j = 0 do j = j+1 if (j>sizeTmp) exit if (j==i) cycle if (.not. tmpTab(j)%keep) cycle if (tmpTab(j)%sa%typ/=tseg) cycle sgj=tmpTab(j)%sa if (.not. estColineaire(sgi,sgj)) cycle if (.not. pointsAlignes(sgi%x,sgi%y,sgi%dx,sgi%dy,sgj%x,sgj%y)) cycle !on renverse le sens de sgj si besoin if (.not. isSameWay(sgi,sgj)) sgj=turnBackSide(tmpTab(j)%sa) oj=isIn(sgj%x,sgj%y,sgi) ; ej=isIn(sgj%dx,sgj%dy,sgi) select case(ej) !seuls cas 0,2,3 utiles pour ej case(0) !seuls cas 0,1,2 utiles pour oj select case(oj) case(0) !si isIn(sgi%x,sgi%y,sgj)/=2 rien sinon, sgi dans sgj if (isIn(sgi%x,sgi%y,sgj)/=2) cycle call adjustMix(tmpTab(i)%sa,sgj) sizeTmp=sizeTmp+1 ! pour programmation defensive if (sizeTmp > taille_table_temp) & call XABORT("G2S : memory problem in addSegsAndClean (1)") tmpTab(sizeTmp)%keep=.true. ; tmpTab(sizeTmp)%sa=sgj tmpTab(j)%sa%dx=sgi%x ; tmpTab(j)%sa%dy=sgi%y tmpTab(sizeTmp)%sa%x=sgi%dx ; tmpTab(sizeTmp)%sa%y=sgi%dy case(1) !origines confondues et sgi dans sgj call adjustMix(tmpTab(i)%sa,sgj) tmpTab(j)%sa%x=sgi%dx ; tmpTab(j)%sa%y=sgi%dy case(2) !sgi et sgj s'intersectent reellement (sgj au dessus) sizeTmp=sizeTmp+1 ! pour programmation defensive if (sizeTmp > taille_table_temp) & call XABORT("G2S : memory problem in addSegsAndClean (2)") tmpTab(sizeTmp)%keep=.true. ; tmpTab(sizeTmp)%sa=sgi tmpTab(i)%sa%dx=sgj%x ; tmpTab(i)%sa%dy=sgj%y tmpTab(sizeTmp)%sa%x=sgj%x ; tmpTab(sizeTmp)%sa%y=sgj%y tmpTab(sizeTmp)%sa%dx=sgi%dx ; tmpTab(sizeTmp)%sa%dy=sgi%dy call adjustMix(tmpTab(sizeTmp)%sa,sgj) tmpTab(j)%sa%x=sgi%dx ; tmpTab(j)%sa%y=sgi%dy end select case(2) !seuls cas 0,1,2 utiles pour oj select case(oj) case(0) !sgi et sgj s'intersectent reellement (sgj en dessous) sizeTmp=sizeTmp+1 ! pour programmation defensive if (sizeTmp > taille_table_temp) & call XABORT("G2S : memory problem in addSegsAndClean (3)") tmpTab(sizeTmp)%keep=.true. ; tmpTab(sizeTmp)%sa=sgi tmpTab(i)%sa%x=sgj%dx ; tmpTab(i)%sa%y=sgj%dy tmpTab(sizeTmp)%sa%x=sgi%x ; tmpTab(sizeTmp)%sa%y=sgi%y tmpTab(sizeTmp)%sa%dx=sgj%dx ; tmpTab(sizeTmp)%sa%dy=sgj%dy call adjustMix(tmpTab(sizeTmp)%sa,sgj) tmpTab(j)%sa%dx=sgi%x ; tmpTab(j)%sa%dy=sgi%y case(1) !origines confondues et sgj dans sgi call adjustMix(tmpTab(j)%sa,sgi) tmpTab(i)%sa%x=sgj%dx ; tmpTab(i)%sa%y=sgj%dy case(2) !sgj dans sgi call adjustMix(tmpTab(j)%sa,sgi) sizeTmp=sizeTmp+1 ! pour programmation defensive if (sizeTmp > taille_table_temp) & call XABORT("G2S : memory problem in addSegsAndClean (4)") tmpTab(sizeTmp)%keep=.true. ; tmpTab(sizeTmp)%sa=sgi tmpTab(i)%sa%dx=sgj%x ; tmpTab(i)%sa%dy=sgj%y tmpTab(sizeTmp)%sa%x=sgj%dx ; tmpTab(sizeTmp)%sa%y=sgj%dy end select case(3) !seuls cas 0,1,2 utiles pour oj select case(oj) case(0) !extremites confondues et sgi dans sgj call adjustMix(tmpTab(i)%sa,sgj) tmpTab(j)%sa%dx=sgi%x ; tmpTab(j)%sa%dy=sgi%y case(1) !sgi et sgj confondus call adjustMix(tmpTab(i)%sa,sgj) tmpTab(j)%keep=.false. case(2) !extremites confondues et sgj dans sgi call adjustMix(tmpTab(j)%sa,sgi) tmpTab(i)%sa%dx=sgj%x ; tmpTab(i)%sa%dy=sgj%y end select end select end do end do do i = 1,sizeTmp if (tmpTab(i)%keep) then if (tmpTab(i)%sa%typ==tseg) then if (isEqual(tmpTab(i)%sa%x,tmpTab(i)%sa%dx).and. & isEqual(tmpTab(i)%sa%y,tmpTab(i)%sa%dy)) tmpTab(i)%keep=.false. end if end if end do !recopie dans le tableau global des elements a conserver sizeSA=0 do i = 1,sizeTmp if (tmpTab(i)%keep) then sizeSA=sizeSA+1 ! pour programmation defensive if (sizeSA > taille_table_tabSegArc) & call XABORT("G2S : memory problem in addSegsAndClean (5)") tabSegArc(sizeSA) = tmpTab(i)%sa end if end do deallocate(tmpTab) end subroutine addSegsAndClean subroutine cutClusters(nbRing,nbClus,szSa) integer,intent(in) :: nbRing,nbClus integer,intent(inout) :: szSa integer :: i,j,szTmp,nbInter double precision :: aClus,bClus,aRing,bRing logical :: revert,eqA,eqB type(t_segArc) :: ring,clus type(segArcArrayTer),dimension(:),allocatable :: tmpTabAr ! pour programmation defensive integer :: taille_table_temp taille_table_temp = nbRing*nbClus*4 !ca me semble suffisant (note de l'auteur ) allocate(tmpTabAr(taille_table_temp),stat=alloc_ok) if (alloc_ok /= 0) call XABORT("G2S: cutCluster(1) => allocation pb") szTmp = 0 !initialisation du tableau temporaire !recopie des cluster et des anneaux dans l'ordre inverse de leur !creation => du plus grand vers le plus petit do i = 1,nbClus szTmp = szTmp + 1 tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .true. tmpTabAr(szTmp)%sa = tabSegArc(szSa) szSa = szSa - 1 end do do i = 1,nbRing szTmp = szTmp + 1 tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .false. tmpTabAr(szTmp)%sa = tabSegArc(szSa) szSa = szSa - 1 end do i=0 do i = i+1 if (i>szTmp) exit if (.not. (tmpTabAr(i)%keep .and. tmpTabAr(i)%cl)) cycle clus = tmpTabAr(i)%sa j = 0 do j = j+1 if (j>szTmp) exit if ((.not. tmpTabAr(j)%keep) .or. tmpTabAr(j)%cl) cycle ring = tmpTabAr(j)%sa nbInter = interArcArc(clus,ring,aClus,bClus,aRing,bRing,revert) select case(nbInter) case(2) if (clus%typ==tcer) then !demi cercle interieur tmpTabAr(i)%sa%typ = tarc tmpTabAr(i)%sa%mixd = ring%mixg tmpTabAr(i)%sa%a = aClus tmpTabAr(i)%sa%b = bClus clus = tmpTabAr(i)%sa !demi cercle exterieur szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (1)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .true. tmpTabAr(szTmp)%sa = clus tmpTabAr(szTmp)%sa%mixd = ring%mixd tmpTabAr(szTmp)%sa%a = bClus tmpTabAr(szTmp)%sa%b = aClus else eqA = isEqualAngl(clus%a,aClus) eqB = isEqualAngl(clus%b,bClus) !debut if (.not.eqA) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (2)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .true. tmpTabAr(szTmp)%sa = clus tmpTabAr(szTmp)%sa%mixd = ring%mixd tmpTabAr(szTmp)%sa%b = aClus end if !fin if (.not.eqB) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (3)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .true. tmpTabAr(szTmp)%sa = clus tmpTabAr(szTmp)%sa%mixd = ring%mixd tmpTabAr(szTmp)%sa%a = bClus end if !milieu if ((.not.eqA).or.(.not.eqB)) then tmpTabAr(i)%sa%a = aClus tmpTabAr(i)%sa%b = bClus tmpTabAr(i)%sa%mixd = ring%mixg clus = tmpTabAr(i)%sa end if end if if (ring%typ==tcer) then tmpTabAr(j)%sa%typ = tarc tmpTabAr(j)%sa%a = bRing tmpTabAr(j)%sa%b = aRing ring = tmpTabAr(j)%sa else tmpTabAr(j)%keep = .false. eqA = isEqualAngl(ring%a,aRing) eqB = isEqualAngl(ring%b,bRing) !premier bout d'arc if (.not.eqA) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (4)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .false. tmpTabAr(szTmp)%sa = ring tmpTabAr(szTmp)%sa%b = aRing end if !second bout d'arc if (.not.eqB) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (5)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .false. tmpTabAr(szTmp)%sa = ring tmpTabAr(szTmp)%sa%a = bRing end if end if case(1) !a priori, tel que l'algo est fait, le cas ne se produit que si !le cluster est tangent a un anneau !travail sur l'anneau if (ring%typ==tcer) then tmpTabAr(j)%sa%typ = tarc tmpTabAr(j)%sa%a = aRing tmpTabAr(j)%sa%b = aRing ring = tmpTabAr(j)%sa else tmpTabAr(j)%keep = .false. eqA = isEqualAngl(ring%a,aRing) eqB = isEqualAngl(ring%b,aRing) !premier bout d'arc if (.not.eqA) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (6)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .false. tmpTabAr(szTmp)%sa = ring tmpTabAr(szTmp)%sa%b = aRing end if !second bout d'arc if (.not.eqB) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (7)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .false. tmpTabAr(szTmp)%sa = ring tmpTabAr(szTmp)%sa%a = aRing end if end if !travail sur le cluster if (clus%typ==tcer) then tmpTabAr(i)%sa%typ = tarc tmpTabAr(i)%sa%a = aClus tmpTabAr(i)%sa%b = aClus clus = tmpTabAr(i)%sa if (revert) then !cluster a l'exterieur de l'anneau tmpTabAr(i)%sa%mixd = ring%mixd else !cluster a l'interieur de l'anneau tmpTabAr(i)%sa%mixd = ring%mixg end if else ! le cluster a deja ete coupe => tangence exterieur ! et derniere utilisation eqA = isEqualAngl(clus%a,aClus) eqB = isEqualAngl(clus%b,aClus) !debut if (.not.eqA) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (8)") tmpTabAr(szTmp)%keep = .true. tmpTabAr(szTmp)%cl = .true. tmpTabAr(szTmp)%sa = clus tmpTabAr(szTmp)%sa%mixd = ring%mixd tmpTabAr(szTmp)%sa%b = aClus end if !fin if (.not.eqB) then szTmp = szTmp + 1 ! programmation defensive if (sztmp > taille_table_temp) & call XABORT("G2S : memory problem in routine cutClusters (9)") tmpTabAr(szTmp)%sa%a = aClus tmpTabAr(szTmp)%sa = clus tmpTabAr(szTmp)%sa%mixd = ring%mixd tmpTabAr(szTmp)%sa%a = aClus end if !milieu if ((.not.eqA).or.(.not.eqB)) then tmpTabAr(i)%keep = .false. i = 0 end if end if end select end do end do !recopie des resultats do i = 1,szTmp if (.not. tmpTabAr(i)%keep) cycle szSa = szSa + 1 tabSegArc(szSa) = tmpTabAr(i)%sa end do deallocate(tmpTabAr) end subroutine cutClusters logical function segWithSameCoord(S1,S2,n) type(t_segArc), intent(in) :: S1, S2 integer, intent(out) :: n segWithSameCoord = .false. ; n=0 if (isEqualConst(S1%x,S2%x)) then if (isEqualConst(S1%dx,S2%dx)) then if (isEqualConst(S1%y,S2%y)) then if (isEqualConst(S1%dy,S2%dy)) then segWithSameCoord = .true. ; n=1 endif endif endif elseif (isEqualConst(S1%x,S2%dx)) then if (isEqualConst(S1%dx,S2%x)) then if (isEqualConst(S1%y,S2%dy)) then if (isEqualConst(S1%dy,S2%y)) then segWithSameCoord = .true. ; n=-1 endif endif endif endif end function segWithSameCoord subroutine PrintTabSegArc(szSA) integer, intent(in) :: szSA type(t_segArc) :: sa integer :: i write(*,10) szSA do i = 1, szSA sa = tabSegArc(i) if (sa%typ == tseg) then write(*,20) i, sa%x, sa%y, sa%dx, sa%dy elseif (sa%typ == tarc) then write(*,30) i, sa%x, sa%y, sa%r, sa%a, sa%b elseif (sa%typ == tcer) then write(*,35) i, sa%x, sa%y, sa%r else write(*,38) i, sa%typ, sa%x, sa%y endif write(*,40) sa%sectg, sa%sectd write(*,50) sa%mixg, sa%mixd write(*,60) sa%nodeg, sa%noded write(*,70) sa%indcellpg, sa%indcellpd write(*,80) sa%neutronicMixg,sa%neutronicMixd enddo 10 format("Number of SegArc :", i4) 20 format("N. SegArc ", i6," segment type : ",/, & "|-> origin/extrem : (", f7.4, ";", f7.4,")/(", f7.4, ";", f7.4,")") 30 format("N. SegArc ", i6," arc type:",/, & "|-> center/radius/angles : (", f7.4,";",f7.4,")/", f7.4,"/",f7.4,";", f7.4) 35 format("N. SegArc ", i6," circle type :",/, & "|-> center/radius : (", f7.4,";",f7.4,")/", f7.4) 38 format("N. SegArc ", i6," unknown type :",/, & "|-> type/origin/ : (", i6,"/",f7.4,";",f7.4,")") 40 format( "|-> sectg sectd : ", i6,2x,i6) 50 format( "|-> mixg mixd : ", i6,2x,i6) 60 format( "|-> nodeg noded : ", i6,2x,i6) 70 format( "|-> IndCellPg IndCellPd : ", i6,2x,i6) 80 format( "|-> neutronicMixg neutronicMixd : ", i6,2x,i6) end subroutine PrintTabSegArc subroutine PrintSegArc1(i) integer, intent(in) ::i type(t_segArc) :: sa sa = tabSegArc(i) if (sa%typ == tseg) then write(*,20) i, sa%x, sa%y, sa%dx, sa%dy elseif (sa%typ == tarc) then write(*,30) i, sa%x, sa%y, sa%r, sa%a, sa%b elseif (sa%typ == tcer) then write(*,35) i, sa%x, sa%y, sa%r else write(*,38) i, sa%typ, sa%x, sa%y endif write(*,40) sa%sectg, sa%sectd write(*,50) sa%mixg, sa%mixd write(*,60) sa%nodeg, sa%noded write(*,70) sa%indcellpg, sa%indcellpd write(*,80) sa%neutronicMixg,sa%neutronicMixd 20 format("N. SegArc ", i6," segment type : ",/, & "|-> origin/extrem : (", f7.4, ";", f7.4,")/(", f7.4, ";", f7.4,")") 30 format("N. SegArc ", i6," arc type:",/, & "|-> center/radius/angles : (", f7.4,";",f7.4,")/", f7.4,"/",f7.4,";", f7.4) 35 format("N. SegArc ", i6," circle type :",/, & "|-> center/radius : (", f7.4,";",f7.4,")/", f7.4) 38 format("N. SegArc ", i6," unknown type :",/, & "|-> type/origin/ : (", i6,"/",f7.4,";",f7.4,")") 40 format( "|-> sectg sectd : ", i6,2x,i6) 50 format( "|-> mixg mixd : ", i6,2x,i6) 60 format( "|-> nodeg noded : ", i6,2x,i6) 70 format( "|-> IndCellPg IndCellPd : ", i6,2x,i6) 80 format( "|-> neutronicMixg neutronicMixd : ", i6,2x,i6) end subroutine PrintSegArc1 subroutine PrintSegArc2(sa) type(t_segArc), intent(in) :: sa if (sa%typ == tseg) then write(*,20) sa%x, sa%y, sa%dx, sa%dy elseif (sa%typ == tarc) then write(*,30) sa%x, sa%y, sa%r, sa%a, sa%b elseif (sa%typ == tcer) then write(*,35) sa%x, sa%y, sa%r else write(*,38) sa%typ, sa%x, sa%y endif write(*,40) sa%sectg, sa%sectd write(*,50) sa%mixg, sa%mixd write(*,60) sa%nodeg, sa%noded write(*,70) sa%indcellpg, sa%indcellpd write(*,80) sa%neutronicMixg,sa%neutronicMixd 20 format(" Segment type : ",/, & "|-> origin/extrem : (", f7.4, ";", f7.4,")/(", f7.4, ";", f7.4,")") 30 format(" Arc type:",/, & "|-> center/radius/angles : (", f7.4,";",f7.4,")/", f7.4,"/",f7.4,";", f7.4) 35 format("Circle type :",/, & "|-> center/radius : (", f7.4,";",f7.4,")/", f7.4) 38 format("unknown type :",/, & "|-> type/origin/ : (", i6,"/",f7.4,";",f7.4,")") 40 format( "|-> sectg sectd : ", i6,2x,i6) 50 format( "|-> mixg mixd : ", i6,2x,i6) 60 format( "|-> nodeg noded : ", i6,2x,i6) 70 format( "|-> IndCellPg IndCellPd : ", i6,2x,i6) 80 format( "|-> neutronicMixg neutronicMixd : ", i6,2x,i6) end subroutine PrintSegArc2 end module segArc