diff options
Diffstat (limited to 'Dragon/src/g2s_construire.f90')
| -rw-r--r-- | Dragon/src/g2s_construire.f90 | 1144 |
1 files changed, 1144 insertions, 0 deletions
diff --git a/Dragon/src/g2s_construire.f90 b/Dragon/src/g2s_construire.f90 new file mode 100644 index 0000000..d98d41b --- /dev/null +++ b/Dragon/src/g2s_construire.f90 @@ -0,0 +1,1144 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! Collect the construction functions for different data types used in +! G2S: module. +! +!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: +! - construit_xxxxx : cree un xxxxx, avec xxxxx parmi car2d, carcel, hex, +! hexcel, tri2d ou tube +! +!----------------------------------------------------------------------- +! +module construire + + use celluleBase + use constType + use constUtiles + use segArc + + implicit none + + double precision, parameter :: PI = 3.141592653589793d0 + double precision, parameter :: tg = 0.414213562373095d0 + +contains + + subroutine construit_car2d(cx,cy,sx,sy,turn,inMeshx,inMeshy,splitx,splity, & + inMix,szSA) + double precision,intent(in) :: cx,cy,sx,sy + double precision,dimension(:),intent(in) :: inMeshx,inMeshy !entre 0 et 1 + integer,intent(in) :: turn + integer,dimension(:), intent(in) :: splitx,splity,inMix + integer,intent(inout) :: szSA + + double precision :: xmin,xmax,ymin,ymax,ssx,ssy,dd + double precision,dimension(:),allocatable :: mmeshx,mmeshy,meshx,meshy + integer,dimension(:),allocatable :: mmix,mix,mixDiag + integer :: i,j,dimx,dimy,smix,tmpSzx + integer :: ind,lgx,lgy,tmpSzy + logical :: ldiag + type(t_segArc) :: sg + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + + !prise en compte du split => creation de meshx, meshy et mix a partir des + !donnees d'entree + + lgx = size(inMeshx)-1 ; lgy = size(inMeshy)-1 + ldiag=(lgx == lgy).and.(lgx*(lgx+1)/2 == size(inMix)) + tmpSzx = sum(splitx(1:lgx)) ; tmpSzy = sum(splity(1:lgy)) ; smix = tmpSzx*tmpSzy + allocate(meshx(tmpSzx+1),meshy(tmpSzy+1),mix(smix),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: construit_car2d(1) => allocation pb") + + !sur x + meshx(1) = 0.d0 + ind = 1 + do i = 1,lgx + dd = (inMeshx(i+1)-inMeshx(i)) / splitx(i) + do j = 1,splitx(i)-1 + ind = ind + 1 + meshx(ind) = meshx(ind-1) + dd + end do + ind = ind + 1 + meshx(ind) = inMeshx(i+1) + end do + !sur y + meshy(1) = 0.d0 + ind = 1 + do i = 1,lgy + dd = (inMeshy(i+1)-inMeshy(i)) / splity(i) + do j = 1,splity(i)-1 + ind = ind + 1 + meshy(ind) = meshy(ind-1) + dd + end do + ind = ind + 1 + meshy(ind) = inMeshy(i+1) + end do + !mix + ! correction of plain CAR2D bug by Alain Hebert (May 2016) + if(ldiag) then + allocate(mixDiag(lgx*lgx),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: construit_car2d(2) => allocation pb") + ind = 0 + do j = 1,lgx + do i = j,lgx + ind=ind+1 + if(ind > size(inMix)) call XABORT('overflow1') + if(i+(j-1)*lgx > lgx*lgx) call XABORT('overflow2') + if(j+(i-1)*lgx > lgx*lgx) call XABORT('overflow3') + mixDiag(i+(j-1)*lgx)=inMix(ind) + mixDiag(j+(i-1)*lgx)=inMix(ind) + end do + end do + endif + ind = 1 + do j = 1,lgy + do i = 1,lgx + if(ldiag) then + mix(ind:ind+splitx(i)-1) = mixDiag(i+(j-1)*lgx) + else + mix(ind:ind+splitx(i)-1) = inMix(i+(j-1)*lgx) + endif + ind = ind + splitx(i) + end do + do i = 2,splity(j) + mix(ind:ind+tmpSzx-1) = mix(ind-tmpSzx:ind-1) + ind = ind + tmpSzx + end do + end do + if(ldiag) deallocate(mixDiag) + + !creation des tableaux de travail pour les abscisses, les ordonnees et les + !milieux, en fonction du turn (la construction des segments se faisant + !ensuite toujours de la meme maniere que pour un turn = 1 + if(mod(turn,2) == 0) then + ssx = sy ; ssy = sx ; dimx = tmpSzy ; dimy = tmpSzx + else + ssx = sx ; ssy = sy ; dimx = tmpSzx ; dimy = tmpSzy + end if + allocate(mmeshx(dimx+1),mmeshy(dimy+1),mmix(smix),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: construit_car2d(3) => allocation pb") + select case(turn) + case(1) + mmeshx(:dimx+1) = meshx(:dimx+1) ; mmeshy(:dimy+1) = meshy(:dimy+1) + mmix(:smix) = mix(:smix) + case(2) + mmeshx(:dimx+1) = meshy(:dimx+1) + do i = 1,dimy+1 + mmeshy(i) = 1.-meshx(dimy+2-i) + end do + mmix(:smix) = (/(mix(i:i+(dimx-1)*dimy:dimy),i=dimy,1,-1)/) + case(3) + do i = 1,dimx+1 + mmeshx(i) = 1.-meshx(dimx+2-i) + end do + do i = 1,dimy+1 + mmeshy(i) = 1.-meshy(dimy+2-i) + end do + mmix(:smix) = mix(dimx*dimy:1:-1) + case(4) + do i = 1,dimx+1 + mmeshx(i) = 1.-meshy(dimx+2-i) + end do + mmeshy(:dimy+1) = meshx(:dimy+1) + mmix(:smix) = (/(mix(i:i-(dimx-1)*dimy:-dimy),i=1+(dimx-1)*dimy,dimx*dimy,1)/) + case(5) + do i = 1,dimx+1 + mmeshx(i) = 1.-meshx(dimx+2-i) + end do + mmeshy(:dimy+1) = meshy(:dimy+1) + mmix(:smix) = (/(mix(i:i-(dimx-1):-1),i=dimx,dimx*dimy,dimx)/) + case(6) + mmeshx(:dimx+1) = meshy(:dimx+1) + mmeshy(:dimy+1) = meshx(:dimy+1) + mmix(:smix) = (/(mix(i:i+(dimx-1)*dimy:dimy),i=1,dimy,1)/) + case(7) + mmeshx(:dimx+1) = meshx(:dimx+1) + do i = 1,dimy+1 + mmeshy(i) = 1.-meshy(dimy+2-i) + end do + mmix(:smix) = (/(mix(i:i+(dimx-1):1),i=1+(dimy-1)*dimx,1,-dimx)/) + case(8) + do i = 1,dimx+1 + mmeshx(i) = 1.-meshy(dimx+2-i) + end do + do i = 1,dimy+1 + mmeshy(i) = 1.-meshx(dimy+2-i) + end do + mmix(:smix) = (/(mix(i:i-(dimx-1)*dimy:-dimy),i=dimx*dimy,1+(dimx-1)*dimy,-1)/) + end select + deallocate(meshx,meshy,mix) + + !construction des segments : + !elle se fait dans l'odre suivant: ... n-1 n + ! ----> ---->----> + ! .................... + ! 1 2 ... + ! ---->----> ----> + !puis: |. : |p + ! : |p-1 + ! |2 : + ! |1 : |. + !Remarque: on distingue le cas particulier des bords, pour l'attribution + !du milieux fooMix + xmin = cx-0.5d0*ssx + xmax = xmin+ssx + ymin = cy-0.5d0*ssy + ymax = ymin+ssy + do i=1,dimx + !cas particulier j=1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (1)") + sg=createSeg(xmin+mmeshx(i)*ssx, ymin, xmin+mmeshx(i+1)*ssx, ymin , & + & mmix(i) , fooMix ) + tabSegArc(szSA)=sg + do j=2,dimy + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (2)") + sg=createSeg(xmin+mmeshx(i)*ssx ,ymin+mmeshy(j)*ssy, & + & xmin+mmeshx(i+1)*ssx,ymin+mmeshy(j)*ssy, & + & mmix((j-1)*dimx+i) ,mmix((j-2)*dimx+i) ) + tabSegArc(szSA)=sg + end do + !cas particulier j=dimy+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (3)") + + sg=createSeg(xmin+mmeshx(i)*ssx ,ymax , & + & xmin+mmeshx(i+1)*ssx,ymax , & + & fooMix ,mmix((dimy-1)*dimx+i)) + tabSegArc(szSA)=sg + end do + do j=1,dimy + !cas particulier i=1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (4)") + sg=createSeg(xmin ,ymin+mmeshy(j)*ssy , & + & xmin ,ymin+mmeshy(j+1)*ssy , & + & fooMix ,mmix((j-1)*dimx+1) ) + tabSegArc(szSA)=sg + do i=2,dimx + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (5)") + sg=createSeg(xmin+mmeshx(i)*ssx ,ymin+mmeshy(j)*ssy , & + & xmin+mmeshx(i)*ssx ,ymin+mmeshy(j+1)*ssy, & + & mmix((j-1)*dimx+i-1),mmix((j-1)*dimx+i) ) + tabSegArc(szSA)=sg + end do + !cas particulier i=dimx+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_car2d (6)") + sg=createSeg(xmax ,ymin+mmeshy(j)*ssy , & + & xmax ,ymin+mmeshy(j+1)*ssy, & + & mmix(j*dimx),fooMix ) + tabSegArc(szSA)=sg + end do + deallocate(mmeshx,mmeshy,mmix) + !call PrintTabSegArc(szSA) + end subroutine construit_car2d + + subroutine construit_carcel(cx,cy,sx,sy,turn,radius,offcx,offcy, & + & splitx,splity,mix,sectori,sectorj,cluster,szSA) + double precision,intent(in) :: cx,cy,sx,sy,offcx,offcy + double precision,dimension(:),intent(in) :: radius + integer,intent(in) :: turn,sectori,sectorj + integer,dimension(:),intent(in) :: splitx,splity,mix + type(t_cluster),dimension(:),pointer :: cluster + integer,intent(inout) :: szSA + + integer :: i,j,keepSz,nbseg,tmp,nb + double precision :: ccx,ccy,d,dorig,dextr,tmpx,tmpy + double precision :: tgx, tgy + double precision :: pt1x,pt1y,pt2x,pt2y,ss,ssx,ssy + double precision,dimension(8) :: coef,coef2 + integer,dimension(:),allocatable :: tmpMix + type(t_segArc) :: sg,ar,tmpSg + double precision :: dist1, dist2 + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + keepSz=szSA + select case (turn) + case(1) ; ccx=cx+offcx ; ccy=cy+offcy + case(2) ; ccx=cx+offcy ; ccy=cy-offcx + case(3) ; ccx=cx-offcx ; ccy=cy-offcy + case(4) ; ccx=cx-offcy ; ccy=cy+offcx + case(5) ; ccx=cx-offcx ; ccy=cy+offcy + case(6) ; ccx=cx+offcy ; ccy=cy+offcx + case(7) ; ccx=cx+offcx ; ccy=cy-offcy + case(8) ; ccx=cx-offcy ; ccy=cy-offcx + end select + + if(mod(turn,2)==0) then ; ssx = sy ; ssy = sx + else ; ssx = sx ; ssy = sy ; end if + + !! creation du cadre + allocate(tmpMix(1)) + tmpMix(1) = mix(size(mix)) + call construit_car2d(cx,cy,sx,sy,turn,(/0.d0,1.d0/),(/0.d0,1.d0/), & + splitx,splity,tmpMix,szSA) + deallocate(tmpMix) + nbseg = 2*splitx(1)*splity(1) + splitx(1) + splity(1) + if(nbseg/=4 .and. sectori/=S_not) call XABORT("G2S: mixing SPLIT and SECT& + & not allowed") + + !! gestion des secteurs + ss = min(sx,sy) + if(sectori==s_not) then + !! rien + else if((sectori==S_X_tot)) then + tmp = mix(size(mix)) + coef = 0.5d0 * (/1.d0,-1.d0,-1.d0, 1.d0,& + 1.d0, 1.d0,-1.d0,-1.d0/) * ss + !! 2 1 + !! X + !! 3 4 + if (sectorj==0) then + do i = 1,4 + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (1)") + tmpSg = createSeg(cx,cy,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end do + else + ar = createArc(ccx,ccy,radius(sectorj+1),0.d0,0.d0,tmp,tmp) + do i = 1,4 + sg = createSeg(cx,cy,cx+coef(i),cy+coef(i+4),tmp,tmp) + nb = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y)) + if(nb==1) then + if(.not.(isEqualConst(pt1x,cx+coef(i)).and. & + isEqualConst(pt1y,cy+coef(i+4)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (1)") + tmpSg = createSeg(pt1x,pt1y,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end if + else if(nb==2) then + if(.not.(isEqualConst(pt2x,cx+coef(i)).and. & + isEqualConst(pt2y,cy+coef(i+4)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (3)") + tmpSg = createSeg(pt2x,pt2y,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end if + end if + end do + end if + else if(sectori==S_T_tot) then + tmp=mix(size(mix)) + coef=0.5d0*(/ ssx,0.d0,-ssx,0.d0,& + 0.d0, ssy,0.d0,-ssy/) + !! 2 + !! 3+1 + !! 4 + if (sectorj==0) then + do i = 1,4 + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (4)") + tmpSg = createSeg(cx,cy,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end do + else + ar = createArc(ccx,ccy,radius(sectorj+1),0.d0,0.d0,tmp,tmp) + do i = 1,4 + sg = createSeg(cx,cy,cx+coef(i),cy+coef(i+4),tmp,tmp) + nb = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y)) + if(nb==1) then + if(.not.(isEqualConst(pt1x,cx+coef(i)).and. & + isEqualConst(pt1y,cy+coef(i+4)))) then + nbseg=nbseg+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (5)") + szSA=szSA+1 + tmpSg = createSeg(pt1x,pt1y,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end if + else if(nb==2) then + if(.not.(isEqualConst(pt2x,cx+coef(i)).and. & + isEqualConst(pt2y,cy+coef(i+4)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (6)") + tmpSg = createSeg(pt2x,pt2y,cx+coef(i),cy+coef(i+4),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+2,4)+1 + tabSegArc(szSA) = tmpSg + end if + end if + end do + end if + else if(sectori==S_TX_tot) then + tmp=mix(size(mix)) + coef = 0.5d0 * (/ ssx,ss,0.d0,-ss,-ssx,-ss,0.d0, ss/) + coef2 = 0.5d0 * (/0.d0,ss, ssy, ss,0.d0,-ss,-ssy,-ss/) + !! 432 + !! 5*1 + !! 678 + if (sectorj==0) then + do i = 1,8 + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (7)") + tmpSg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end do + else + ar = createArc(ccx,ccy,radius(sectorj+1),0.d0,0.d0,tmp,tmp) + do i = 1,8 + sg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + nb = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y)) + if(nb==1) then + if(.not.(isEqualConst(pt1x,cx+coef(i)).and. & + isEqualConst(pt1y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (7)") + tmpSg = createSeg(pt1x,pt1y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + else if(nb==2) then + if(.not.(isEqualConst(pt2x,cx+coef(i)).and. & + isEqualConst(pt2y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (8)") + tmpSg = createSeg(pt2x,pt2y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + end if + end do + end if + elseif (sectori == S_TXS_tot) then + tgx = ssx * tg + tgy = ssy * tg + tmp = mix(size(mix)) + coef = 0.5d0 *(/ ssx, tgy, -tgy, -ssx, -ssx, -tgy, tgy, ssx /) + coef2 = 0.5d0 *(/ tgx, ssy, ssy, tgx, -tgx, -ssy, -ssy, -tgx /) + if (sectorj==0) then + do i = 1,8 + nbseg=nbseg+1 ! + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (9)") + tmpSg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end do + else + ar = createArc(ccx,ccy,radius(sectorj+1),0.d0,0.d0,tmp,tmp) + do i = 1,8 + sg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + nb = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y)) + if(nb==1) then + if(.not.(isEqualConst(pt1x,cx+coef(i)).and. & + isEqualConst(pt1y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (10)") + tmpSg = createSeg(pt1x,pt1y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + else if(nb==2) then + if(.not.(isEqualConst(pt2x,cx+coef(i)).and. & + isEqualConst(pt2y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine& + & construit_carcel (11)") + tmpSg = createSeg(pt2x,pt2y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + end if + end do + end if + elseif (sectori == S_WM_tot) then + tgx = ssx * tg + tgy = ssy * tg + tmp = mix(size(mix)) + coef = 0.5d0 *(/ ssx, tgy, -tgy, -ssx, -ssx, -tgy, tgy, ssx /) + coef2 = 0.5d0 *(/ tgx, ssy, ssy, tgx, -tgx, -ssy, -ssy, -tgx /) + if (sectorj==0) then + do i = 1,8 + nbseg=nbseg+1 ! + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (12)") + tmpSg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end do + do i = 1, 4 + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (13)") + tmpSg = createSeg(cx+coef(2*i-1),cy+coef2(2*i-1),cx+coef(2*i),cy+coef2(2*i),tmp,tmp) + tmpSg%sectg=i + tmpSg%sectd= mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + enddo + else + ar = createArc(ccx,ccy,radius(sectorj+1),0.d0,0.d0,tmp,tmp) + do i = 1,8 + sg = createSeg(cx,cy,cx+coef(i),cy+coef2(i),tmp,tmp) + nb = abs(interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y)) + if(nb==1) then + if(.not.(isEqualConst(pt1x,cx+coef(i)).and. & + isEqualConst(pt1y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (14)") + tmpSg = createSeg(pt1x,pt1y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + else if(nb==2) then + if(.not.(isEqualConst(pt2x,cx+coef(i)).and. & + isEqualConst(pt2y,cy+coef2(i)))) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (15)") + tmpSg = createSeg(pt2x,pt2y,cx+coef(i),cy+coef2(i),tmp,tmp) + tmpSg%sectg = i + tmpSg%sectd = mod(i+6,8)+1 + tabSegArc(szSA) = tmpSg + end if + end if + end do + ! creation des 4 segments fermant les ailes + do i = 1, 4 + pt1x = coef(2*i-1)+cx + pt1y = coef2(2*i-1)+cy + pt2x = coef(2*i)+cx + pt2y = coef2(2*i)+cy + dist1 = (pt1x-ccx)*(pt1x-ccx)+(pt1y-ccy)*(pt1y-ccy) + dist2 = (pt2x-ccx)*(pt2x-ccx)+(pt2y-ccy)*(pt2y-ccy) + + if (sqrt(dist1)>radius(sectorj+1).and.sqrt(dist2)>radius(sectorj)) then + nbseg=nbseg+1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_carcel (16)") + tmpSg = createSeg(cx+coef(2*i-1),cy+coef2(2*i-1),cx+coef(2*i),cy+coef2(2*i),tmp,tmp) + !tmpSg%sectg=i + !tmpSg%sectd= mod(i+6,8)+1 + tmpSg%sectg=mod(2*i+5,8)+1 + tmpSg%sectd=2*i + tabSegArc(szSA) = tmpSg + else + call XABORT("G2S : Intersection troubles with WindMill") + endif + enddo + end if + else + call XABORT("G2S : type of sectorisation not recognised") + endif + + !prise en compte des milieux pour un segment completement a l'interieur + !d'une zone annulaire (ie non coupe) + ! CS-IV : La boucle s'applique sur les nbseg derniers segments crees + ! CS-IV : ie sur les segments que l'on vient de faire (nbseg et szSA sont + ! CS-IV : incrementes de la meme facon mais le premier part de 0 alors que + ! CS-IV : part du nombre de SA deja crees. + ! CS-IV : donc le commentaire d'origine est faux, et devrait etre : + ! CS-IV : Prise en compte des milieux pour les segments crees. + do i = 0,nbseg-1 + sg = tabSegArc(szSA-i) + d = distance(ccx,ccy,sg%x,sg%y,sg%dx,sg%dy,tmpx,tmpy) + tmpx=ccx-sg%x ; tmpy=ccy-sg%y ; dorig = sqrt((tmpx*tmpx)+(tmpy*tmpy)) + tmpx=ccx-sg%dx ; tmpy=ccy-sg%dy ; dextr = sqrt((tmpx*tmpx)+(tmpy*tmpy)) + do j = 1,size(radius)-1 + if((min(d,dorig,dextr)+epsilon>radius(j)).and. & + & (max(d,dorig,dextr)-epsilon<radius(j+1))) then + if(sg%mixg/=fooMix) tabSegArc(szSA-i)%mixg=mix(j) + if(sg%mixd/=fooMix) tabSegArc(szSA-i)%mixd=mix(j) + end if + end do + end do + call construit_tube(ccx,ccy,radius,mix,cluster,szSA) + call coupeCercle(keepSz,szSA,nbseg,tRec) + + if(sectori/=S_not) then + call splitSegsForSect(keepSz,szSA) + call majSectori(keepSz,szSA,sectori,tRec,cx,cy) + end if + end subroutine construit_carcel + + subroutine construit_hexhom(cx,cy,sd,mix,szSA,sectori) + double precision,intent(in) :: cx,cy,sd + integer,intent(in) :: mix + integer,intent(inout) :: szSA + integer,intent(in) :: sectori + + double precision,dimension(7) :: xx,yy + double precision :: sqrt3_2S,S_2 + integer :: i + integer,parameter,dimension(7) :: sect=(/1,2,3,4,5,6,1/) + type(t_segArc) :: tmpSg + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + + S_2=5.d-1*sd ; sqrt3_2S=S_2*sqrt(3.d0) + xx(1) = cx + S_2 ; yy(1) = cy + sqrt3_2S + xx(2) = cx - S_2 ; yy(2) = yy(1) + xx(3) = cx - sd ; yy(3) = cy + xx(4) = xx(2) ; yy(4) = cy - sqrt3_2S + xx(5) = xx(1) ; yy(5) = yy(4) + xx(6) = cx + sd ; yy(6) = cy + xx(7) = xx(1) ; yy(7) = yy(1) + do i = 1,6 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_hexhom (1)") + tmpSg = createSeg(xx(i),yy(i),xx(i+1),yy(i+1),mix,fooMix) + tabSegArc(szSA) = tmpSg + end do + if(sectori==S_X_tot) then + do i = 1,6 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_hexhom (2)") + tmpSg = createSeg(cx,cy,xx(i),yy(i),mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end do + else if(sectori/=S_not) then + call XABORT("G2S : type of sectorisation not recognised for & + &homogeneous hexagonal geometry") + endif + end subroutine construit_hexhom + + subroutine construit_hex(cx,cy,sd,mix,szSA,sectori,sectorj,radius,ocx,ocy) + double precision,intent(in) :: cx,cy,sd + integer,intent(in) :: mix + integer,intent(inout) :: szSA + integer,intent(in) :: sectori,sectorj + double precision,dimension(:),intent(in) :: radius + double precision,intent(in) :: ocx,ocy + + double precision,dimension(7) :: xx,yy + double precision :: sqrt3_2S,S_2 + double precision :: pt1x,pt1y,pt2x,pt2y + integer :: i,nb + integer,parameter,dimension(7) :: sect=(/1,2,3,4,5,6,1/) + type(t_segArc) :: sg,ar,tmpSg + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + + S_2=5.d-1*sd ; sqrt3_2S=S_2*sqrt(3.d0) + xx(1) = cx + S_2 ; yy(1) = cy + sqrt3_2S + xx(2) = cx - S_2 ; yy(2) = yy(1) + xx(3) = cx - sd ; yy(3) = cy + xx(4) = xx(2) ; yy(4) = cy - sqrt3_2S + xx(5) = xx(1) ; yy(5) = yy(4) + xx(6) = cx + sd ; yy(6) = cy + xx(7) = xx(1) ; yy(7) = yy(1) + do i = 1,6 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_hex (1)") + tmpSg = createSeg(xx(i),yy(i),xx(i+1),yy(i+1),mix,fooMix) + tabSegArc(szSA) = tmpSg + end do + if (sectori==S_X_tot) then + if (sectorj==0) then + do i = 1,6 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_hex (2)") + tmpSg = createSeg(cx,cy,xx(i),yy(i),mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end do + else + ! creation du + grand cercle impermeable a la sectorisation + ar = createArc(ocx,ocy,radius(sectorj+1),0.d0,0.d0,mix,mix) + do i = 1,6 + sg=createSeg(cx,cy,xx(i),yy(i),mix,mix) + nb=interSgAr(sg,ar,pt1x,pt1y,pt2x,pt2y) + select case(nb) + case(-1) + if(.not.(isEqualConst(pt1x,xx(i)).and. & + isEqualConst(pt1y,yy(i)))) then + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_hex (3)") + tmpSg = createSeg(pt1x,pt1y,xx(i),yy(i),mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end if + case(0) + if(.not.(isEqualConst(cx,xx(i)).and.isEqualConst(cy,yy(i)))) then + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_hex (4)") + tmpSg = createSeg(cx,cy,xx(i),yy(i),mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end if + case(1) + if(.not.(isEqualConst(cx,pt1x).and.isEqualConst(cy,pt1y))) then + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_hex (5)") + tmpSg = createSeg(cx,cy,pt1x,pt1y,mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end if + case(2) + if(.not.(isEqualConst(cx,pt1x).and.isEqualConst(cy,pt1y))) then + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_hex (6)") + tmpSg = createSeg(cx,cy,pt1x,pt1y,mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end if + if(.not.(isEqualConst(pt2x,xx(i)).and. & + isEqualConst(pt2y,yy(i)))) then + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine & + &construit_hex (7)") + tmpSg = createSeg(pt2x,pt2y,xx(i),yy(i),mix,mix) + tmpSg%sectg = sect(i+1) + tmpSg%sectd = sect(i) + tabSegArc(szSA) = tmpSg + end if + end select + end do + endif + else if(sectori/=S_not) then + call XABORT("G2S : type of sectorisation not recognised for & + &hexagonal pincell geometry") + endif + end subroutine construit_hex + + subroutine construit_hexcel(cx,cy,sd,turn,radius,offcx,offcy, & + & mix,secori,secorj,cluster,szSA) + double precision,intent(in) :: cx,cy,sd,offcx,offcy + double precision,dimension(:),intent(in) :: radius + integer,intent(in) :: turn,secori,secorj + integer,dimension(:) :: mix + type(t_cluster),dimension(:),pointer :: cluster + integer,intent(inout) :: szSA + + integer :: i,j,keepSz,nbseg + double precision :: tocx,tocy,tpi_3,costp,sintp + double precision :: d,tmpx,tmpy,dorig,dextr + type(t_segArc) :: sg + + keepSz=szSA + if(turn<=6) then + tpi_3=(turn-1)*pi_3_c ; costp=cos(tpi_3) ; sintp=sin(tpi_3) + tocx = cx + costp*offcx - sintp*offcy + tocy = cy + sintp*offcx + costp*offcy + else + tpi_3=(12-turn)*pi_3_c ; costp=cos(tpi_3) ; sintp=sin(tpi_3) + tocx = cx + costp*offcx - sintp*offcy + tocy = cy - sintp*offcx - costp*offcy + end if + call construit_hex(cx,cy,sd,mix(size(mix)),szSA,secori,secorj,radius,tocx,tocy) + !prise en compte des milieux pour un segement completement a l'interieur + !d'une zone annulaire (ie non coupee) + do i = 0,szSA-keepSz-1 + sg = tabSegArc(szSA-i) + d = distance(tocx,tocy,sg%x,sg%y,sg%dx,sg%dy,tmpx,tmpy) + tmpx=tocx-sg%x ; tmpy=tocy-sg%y ; dorig=sqrt((tmpx*tmpx)+(tmpy*tmpy)) + tmpx=tocx-sg%dx ; tmpy=tocy-sg%dy ; dextr=sqrt((tmpx*tmpx)+(tmpy*tmpy)) + do j = 2,size(radius)-1 + if((min(d,dorig,dextr)>radius(j)).and. & + & (max(d,dorig,dextr)<radius(j+1))) then + if(sg%mixg/=fooMix) tabSegArc(szSA-i)%mixg=mix(j) + if(sg%mixd/=fooMix) tabSegArc(szSA-i)%mixd=mix(j) + end if + end do + end do + nbseg=szSA-keepSz + !creation des anneaux + call construit_tube(tocx,tocy,radius,mix,cluster,szSA) + !creation des anneaux + call coupeCercle(keepSz,szSA,nbseg,tHex) + !creation des anneaux + if(secori/=S_not) call majSectori(keepSz,szSA,secori,tHex,cx,cy) + !creation des anneaux + end subroutine construit_hexcel + + subroutine construit_tri2d(cx,cy,sd,turn,split,mix,szSA) + double precision,intent(in) :: cx,cy,sd + integer,intent(in) :: turn,split + integer,dimension(:),allocatable :: mix + integer,intent(inout) :: szSA + + integer :: i,j,k,dimx,dimy,indMix,deltaIndMix + type(t_segArc) :: sg + double precision :: deltax,deltay + + ! CS-IV : F_C_2 + !double precision,dimension(:),allocatable :: xx,yy + double precision, dimension(2*split+1) :: xx + double precision, dimension(split+1) :: yy + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + + !on cree les triangles sans turn en (0,0), puis on les tourne + !et on les translate + + !dimensionnement et remplissage des coordonnees + ! x + dimx = 2*split+1 + deltax = 5.d-1*sd/split + ! CS-IV : F_C_2 allocate(xx(dimx)) + xx(1) = -5.d-1*sd + do i = 2,dimx + xx(i) = xx(i-1) + deltax + end do + ! y + dimy = split+1 + ! CS-IV : F_C_2 allocate(yy(dimy)) + deltay = deltax*sqrt(3.d0) + yy(1) = -25.d-2*sd*sqrt(3.d0) + do j = 2,dimy + yy(j) = yy(j-1) + deltay + end do + !creation des triangles + !en trois etapes successives : 1=_ ; 2=\ ; 3=/ + + !Premiere ligne + ! Triangle de debut de ligne + ! 1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (1)") + sg = createSeg(xx(1),yy(1),xx(3),yy(1),mix(1),fooMix) + tabSegArc(szSA) = translateAndTurnTriSg(cx,cy,turn,sg) + if(split/=1) then + ! 2 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (2)") + sg = createSeg(xx(3),yy(1),xx(2),yy(2),mix(1),mix(2)) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + end if + ! 3 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (3)") + sg = createSeg(xx(2),yy(2),xx(1),yy(1),mix(1),fooMix) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + ! Triangles de milieu de ligne + do i = 2,split-1 + ! 1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (4)") + sg = createSeg(xx(2*i-1),yy(1),xx(2*i+1),yy(1),mix(2*i-1),fooMix) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + ! 2 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (5)") + sg = createSeg(xx(2*i+1),yy(1),xx(2*i),yy(2),mix(2*i-1),mix(2*i)) + tabSegArc(szSA) = translateAndTurnTriSg(cx,cy,turn,sg) + ! 3 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (6)") + sg = createSeg(xx(2*i),yy(2),xx(2*i-1),yy(1),mix(2*i-1),mix(2*i-2)) + tabSegArc(szSA) = translateAndTurnTriSg(cx,cy,turn,sg) + end do + ! Triangle de fin de ligne + if(split/=1) then + ! 1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (7)") + sg = createSeg(xx(dimx-2),yy(1),xx(dimx),yy(1),mix(2*split-1),fooMix) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + end if + ! 2 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (8)") + sg = createSeg(xx(dimx),yy(1),xx(dimx-1),yy(2),mix(2*split-1),fooMix) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + if(split/=1) then + ! 3 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (9)") + sg = createSeg(xx(dimx-1),yy(2),xx(dimx-2),yy(1),mix(2*split-1), & + & mix(2*split-2)) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + end if + + !Autres lignes + indMix = 2*split+1 + do j = 2,split + k = j + indMix = indMix - 1 + deltaIndMix = 2*(split-j+1) + do i = 1,dimy-j + ! 1 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (10)") + sg = createSeg(xx(k),yy(j),xx(k+2),yy(j), & + & mix(indMix),mix(indMix-deltaIndMix)) + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + ! 2 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (11)") + if(i==dimy-j) then + sg = createSeg(xx(k+2),yy(j),xx(k+1),yy(j+1), & + & mix(indMix),fooMix) + else + sg = createSeg(xx(k+2),yy(j),xx(k+1),yy(j+1), & + & mix(indMix),mix(indMix+1)) + end if + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + ! 3 + szSA=szSA+1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tri2d (12)") + if(i==1) then + sg = createSeg(xx(k+1),yy(j+1),xx(k),yy(j), & + & mix(indMix),fooMix) + else + sg = createSeg(xx(k+1),yy(j+1),xx(k),yy(j), & + & mix(indMix),mix(indMix-1)) + end if + tabSegArc(szSA)=translateAndTurnTriSg(cx,cy,turn,sg) + k = k + 2 + indMix = indMix + 2 + end do + end do + ! CS-IV : F_C_2 deallocate(xx,yy) + end subroutine construit_tri2d + + subroutine construit_tube(cx,cy,radius,mix,cluster,szSA) + double precision,intent(in) :: cx,cy + double precision,dimension(:),intent(in) :: radius + integer,dimension(:), intent(in) :: mix + type(t_cluster),dimension(:),pointer :: cluster + integer,intent(inout) :: szSA + + integer :: i,j,k,sr,sc,nbp,keepSize,nbRing,nbClus,merge + double precision :: ccx,ccy,alpha + type(t_segArc) :: ce + integer,dimension(:),allocatable :: lastMix + + ! programmation defensive + integer :: dimtabSegArc + dimtabSegArc = size(tabSegArc) + + keepSize = szSa + + !anneaux + sr = size(radius) + do i = 2,sr + szSa = szSa + 1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tube (1)") + ce = createArc(cx,cy,radius(i),0.d0,0.d0,mix(i-1),mix(i)) + tabSegArc(szSa) = ce + end do + + nbRing = szSa - keepSize + keepSize = szSa + + merge = size(mix) !servira a identifier les milieux dans les clusters + !clusters + if(.not. associated(cluster)) return !pas de cluster + sc = size(cluster) + if(sc==0) return !pas de cluster + ! calcul du milieux exterieur de chaque cluster + allocate(lastMix(sc)) + do i = 1,sc + lastMix(i) = fooMix + do j = 2,sr + if(cluster(i)%radiusOfPin < radius(j)) then + lastMix(i) = mix(j-1) + exit + end if + end do + end do + ! creation des clusters par rayon croissant (pour que les cercles + ! de rayon max rencontrent les anneaux en premier car la recopie + ! inverse l'ordre) + do i = 1,sc + sr = size(cluster(i)%radius) + nbp = cluster(i)%nbrPin + do j = 1,nbp + ! write(*,*) 'offset cluster ',j,' : ',cluster(i)%angleOfPin + ! alpha = deg2rad*cluster(i)%angleOfPin + (j-1)*dpi_c/nbp + alpha = cluster(i)%angleOfPin + (j-1)*dpi_c/nbp + ccx = cx + cluster(i)%radiusOfPin * cos(alpha) + ccy = cy + cluster(i)%radiusOfPin * sin(alpha) + !tous les cercles sauf le dernier + do k = 1,sr-2 + szSa = szSa + 1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tube (2)") + ce = createArc(ccx,ccy,cluster(i)%radius(k+1),0.d0,0.d0, & + merge+k,merge+k+1) + tabSegArc(szSa) = ce + end do + !cercle le plus grand + szSa = szSa + 1 + ! programmation defensive + if (szSA > dimTabSegArc)& + call XABORT("G2S: memory problem in routine construit_tube (3)") + ce = createArc(ccx,ccy,cluster(i)%radius(sr),0.d0,0.d0, & + merge+sr-1,lastMix(i)) + tabSegArc(szSa) = ce + end do + merge = merge + sr - 1 + end do + deallocate(lastMix) + + nbClus = szSa - keepSize + !coupage des anneaux par les clusters + if(nbRing>0) call cutClusters(nbRing,nbClus,szSa) + end subroutine construit_tube + +end module construire |