! !----------------------------------------------------------------------- ! !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