diff options
Diffstat (limited to 'Dragon/src/g2s_pretraitement.f90')
| -rw-r--r-- | Dragon/src/g2s_pretraitement.f90 | 2243 |
1 files changed, 2243 insertions, 0 deletions
diff --git a/Dragon/src/g2s_pretraitement.f90 b/Dragon/src/g2s_pretraitement.f90 new file mode 100644 index 0000000..baf20ad --- /dev/null +++ b/Dragon/src/g2s_pretraitement.f90 @@ -0,0 +1,2243 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! Preprocessing of geometric data recovered from a geometry LCM data +! structure. +! +!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: +! Ce fichier regroupe toutes les actions effectuees en amont du code, pour +! verifier et completer le jeux de donnees geometriques. Il utilise de maniere +! intensive l'api de PyLCM. +! Les fonctions de ce fichier ont ete initialement developpees en python, puis +! traduite en fortran90. C'est pourquoi la forme des algorithmes peut etre +! parfois deroutante... +! \\\\ +! Le module pretraitement definit quelques constantes, dont en particulier la +! precision admise en entree du code (pour la coherence des donnees). +! De plus, les matrices "composeRotRec" et "composeRotTri" permettent +! d'effectuer la composition de deux transformations definies selon le mode +! Dragon. +! \\\\ +! Le seul pretraitement qui soit reellement delicat, est celui qui correspond +! a une geometrie de type gigogne rectangulaire. Il est effectue par la routine +! "creerEtChargerNewDataRec", et est documente dans le source de cette routine. +! +!----------------------------------------------------------------------- +! +module pretraitement + use boundcond + use cast + use celluleBase + use cellulePlaced + use constType + use GANLIB + + implicit none + + real,parameter :: geomPrec = 1.e-4 + real,parameter :: sqrt_3f = 1.732050807568 + + !composeRotXXX(t1,t2)=t2ot1 + integer,dimension(8,8),parameter :: composeRotRec = & + & reshape((/1,2,3,4,5,6,7,8, & + & 2,3,4,1,8,5,6,7, & + & 3,4,1,2,7,8,5,6, & + & 4,1,2,3,6,7,8,5, & + & 5,6,7,8,1,2,3,4, & + & 6,7,8,5,4,1,2,3, & + & 7,8,5,6,3,4,1,2, & + & 8,5,6,7,2,3,4,1/) , (/8,8/)) + + integer,dimension(12,12),parameter :: composeRotTri = & + & reshape((/ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, & + & 2, 3, 4, 5, 6, 1, 8, 9,10,11,12, 7, & + & 3, 4, 5, 6, 1, 2, 9,10,11,12, 7, 8, & + & 4, 5, 6, 1, 2, 3,10,11,12, 7, 8, 9, & + & 5, 6, 1, 2, 3, 4,11,12, 7, 8, 9,10, & + & 6, 1, 2, 3, 4, 5,12, 7, 8, 9,10,11, & + & 7,12,11,10, 9, 8, 1, 6, 5, 4, 3, 2, & + & 8, 7,12,11,10, 9, 2, 1, 6, 5, 4, 3, & + & 9, 8, 7,12,11,10, 3, 2, 1, 6, 5, 4, & + & 10, 9, 8, 7,12,11, 4, 3, 2, 1, 6, 5, & + & 11,10, 9, 8, 7,12, 5, 4, 3, 2, 1, 6, & + & 12,11,10, 9, 8, 7, 6, 5, 4, 3, 2, 1/) , (/12,12/)) + +contains + subroutine prepareData(geoIp,sizeB,sizeP,lgMaxGig) + type(c_ptr),intent(in):: geoIp + integer,intent(inout) :: sizeB,sizeP,lgMaxGig + + integer,dimension(40) :: st + type(c_ptr) :: ip + integer :: nbNewMix,i + + !premier completion de la structure geometrique pour transformer + !les constructions melangeant MIX et CELL, en des constructions + !uniformes de ce point de vue + + nbNewMix = 0 + call separateMixAndCell(geoIp,nbNewMix) + + ip = geoIp + !creation des cellules de bases + call buildCellsBase(ip,sizeB,'/ ') + ip = geoIp + + call LCMGET(ip,'STATE-VECTOR',st) + select case(st(1)) + case(G_Car2d,G_Carcel) + geomTyp = RecTyp + call creerEtChargerNewDataRec(ip,sizeB,sizeP) + case(G_Hex,G_Hexcel) + geomTyp = HexTyp + call creerEtChargerNewDataHex(ip,sizeB,sizeP) + case(G_Tri) + geomTyp = TriaTyp + call creerEtChargerNewDataTri(ip,sizeB,sizeP) + case(G_Tube) + geomTyp = TubeTyp + if (sizeB/=1) call XABORT("G2S: more than one cellule in a cylindrical& + & geometry") + sizeP = 1 + tabCellulePlaced(1)%indice = 1 + tabCellulePlaced(1)%xcenter = 0.d0 + tabCellulePlaced(1)%ycenter = 0.d0 + tabCellulePlaced(1)%turn = 1 + allocate (tabCellulePlaced(1)%gig(1),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: prepareData(2) => allocation pb") + tabCellulePlaced(1)%gig = 1 + allocate (tabCellulePlaced(1)%mrg(1),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: prepareData(3) => allocation pb") + tabCellulePlaced(1)%mrg = 1 + case default + call XABORT("G2S: Type of geometry not supported") + end select + + lgMaxGig = 0 + do i = 1,sizeP + lgMaxGig = max(lgMaxGig,size(tabCellulePlaced(i)%gig)) + end do + end subroutine prepareData + + recursive subroutine separateMixAndCell(ipIn,nbNewMix) + type(c_ptr),intent(in):: ipIn + integer,intent(inout) :: nbNewMix + + type(c_ptr) :: ip,jp + integer :: i,lgC,lg,ty,indOfCellToAdd,mixToPut,nbCellToAdd,lgC2 + character*12 :: newCellName,text12 + integer,dimension(40) :: sv,newSv + integer,dimension(:),allocatable :: mix + character*12,dimension(:),allocatable :: cell + + ip = ipIn + call LCMLEN(ip,'CELL ',lgC,ty) + lgC = lgC/3 + if (lgC/=0) then !il y a des sous-cellules + allocate(cell(lgC),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: separateMixAndCell(1) => allocation pb") + call LCMGTC(ip,'CELL ',12,lgC,cell) + do i = 1,lgC + !traitement recursif sur les sous cellules de la cellule courante + call LCMLEN(ip,cell(i),lg,ty) + if (ty==0) then !la sous cellule est bien definie ici + jp=LCMDID(ip,cell(i)) + call separateMixAndCell(jp,nbNewMix) + end if + end do + deallocate(cell) + !test sur la mixite de la cellule courante + call LCMLEN(ip,'MIX ',lg,ty) + allocate(mix(lg),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: prepareData(2) => allocation pb") + call LCMGET(ip,'MIX ',mix) + nbCellToAdd = count(mix(:lg)>0) + if (nbCellToAdd/=0) then !la cellule courante est mixte MIX/CELL + call LCMGET(ip,'STATE-VECTOR',sv) + newSv(:) = 0 + newSv(1) = sv(1) + newSv(2) = 0 + newSv(3) = 1 + newSv(4) = 1 + newSv(5) = 0 + newSv(6) = 1 + newSv(7) = 1 + allocate(cell(lg),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: prepareData(3) => allocation pb") + cell(:lg)=' ' + call LCMLEN(ip,'CELL ',lgC2,ty) + lgC2 = lgC2/3 + call LCMGTC(ip,'CELL ',12,lgC2,cell) + indOfCellToAdd = lgC + do i = 1,lg + if (mix(i)>0) then + !creation d'une nouvelle cellule + nbNewMix = nbNewMix+1 + indOfCellToAdd = indOfCellToAdd+1 + mixToPut = mix(i) + mix(i) = -indOfCellToAdd + text12 = i2s(nbNewMix) + newCellName = 'Nmix'//text12(:8) + cell(indOfCellToAdd) = newCellName + jp=LCMDID(ip,newCellName) + call LCMPUT(jp,'MIX ',1,1,mixToPut) + newSv(7) = mixToPut + call LCMPUT(jp,'STATE-VECTOR',40,1,newSv) + call LCMPUT(jp,'NCODE ',6,1,(/0,0,0,0,0,0/)) + call LCMPUT(jp,'ICODE ',6,1,(/0,0,0,0,0,0/)) + text12='L_GEOM' + call LCMPTC(jp,'SIGNATURE',12,text12) + call LCMPUT(jp,'ZCODE ',6,2,(/0.,0.,0.,0.,0.,0./)) + end if + end do + !modification du STATE-VECTOR, du MIX et du CELL de la + !cellule courante + sv(9) = sv(9)+nbCellToAdd + call LCMPUT(ip,'STATE-VECTOR',40,1,sv) + call LCMPUT(ip,'MIX ',lg,1,mix) + call LCMPTC(ip,'CELL ',12,lg,cell) + deallocate(cell) + end if + deallocate(mix) + end if + end subroutine separateMixAndCell + + subroutine creerEtChargerNewDataRec(geoIp,szB,szP) + type(c_ptr),intent(in):: geoIp + integer,intent(inout) :: szB,szP + + !effectue le pretraitement initialement ecrit en python. + !celui-ci utilisait intensivement les dictionnaires python. + !Ceux-ci seront donc remplaces par des objets LCM + type(c_ptr) :: ip,sidesIp,lCentreIp,lMinMaxIp,centreCalculesIp,lcIp,dicoIp + integer :: szLc + + ip = geoIp + !modification de la geometrie dans le cas diagonal, par ajout des + !cellules manquantes dans la gigogne exterieure (ceci pour permettre + !la mise en oeuvre du traitement general aux geometries rectangulaires) + call traiteConditionDiagonale(ip) + !creation de la liste des carres + call LCMOP(lcIp,'lc ',0,1,0) + call creeListeCellules(lcIp,ip,szLc) + !creation du dictionnaire de verification des donnees geometriques + call LCMOP(dicoIp,'dico ',0,1,0) + call creeDico(dicoIp,lcIp,szLc) + + !creation de la liste des longueurs des cotes et resolution du systeme + call LCMOP(sidesIp,'sides ',0,1,0) + call resoudDico(dicoIp,sidesIp) + !cree et prepare la liste des centres des cellules de base, + !avec leur gigogne d'origine + !et destruction du dictionnaire + call LCMOP(lCentreIp,'centres ',0,1,0) + call prepareCentres(lCentreIp,lcIp,sidesIp,dicoIp) + call LCMCL(dicoIp,2) + !preparation des longeurs de la gigogne exterieure, et des coordonnees + !min et max de ses sous cellules + call LCMOP(lMinMaxIp,'minAndMaxXY ',0,1,0) + call prepareMinMax(lMinMaxIp,lCentreIp) + !resolution du systeme donnant les coordonnees des centres des cellules + !de base et stockage dans une nouvelle structure + call LCMOP(centreCalculesIp,'resCentres ',0,1,0) + call calculeCentre(centreCalculesIp,lCentreIp,0) + !compilation des resultats et ajout dans la geometrie des nouvelles donnees + call compileResutats(ip,centreCalculesIp,sidesIp,lMinMaxIp) + call LCMCL(lcIp,2) + call LCMCL(centreCalculesIp,2) + call LCMCL(sidesIp,2) + call LCMCL(lMinMaxIp,2) + !utilisation des donnees (on accroche le traitement classique developpe + !pour python) + call chargerNewData(geoIp,szB,szP) + end subroutine creerEtChargerNewDataRec + + subroutine traiteConditionDiagonale(geoIp) + type(c_ptr),intent(in) :: geoIp + + integer :: i,j,k,pos,n,lgAv,lgAp,lg,ty + integer,dimension(6) :: ncode + integer,dimension(40) :: sv + integer,dimension(:),allocatable :: mixAv,mixAp,turnAv,turnAp + integer,dimension(:),allocatable :: mergeAv,mergeAp + + call LCMGET(geoIp,'NCODE ',ncode) + if ((ncode(1)/=B_Diag) .and. (ncode(2)/=B_Diag)) & + return !pas de condition DIAG + call LCMGET(geoIp,'STATE-VECTOR',sv) + if (sv(8)==0) return !pas de sous cellules + call LCMLEN(geoIp,'MIX ',lgAv,ty) + n = (nint(sqrt(1.+8.*lgAv)) - 1) / 2 + lgAp = n*n + + allocate(mixAv(lgAv),turnAv(lgAv),mergeAv(lgAv),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: traiteConditionDiagonale(1) => allocation pb") + allocate(mixAp(lgAp),turnAp(lgAp),mergeAp(lgAp),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: traiteConditionDiagonale(2) => allocation pb") + call LCMGET(geoIp,'MIX ',mixAv) + call LCMLEN(geoIp,'TURN ',lg,ty) + if (lg == 0) then + turnAv(:lgAv) = 1 + else + call LCMGET(geoIp,'TURN ',turnAv) + end if + call LCMLEN(geoIp,'MERGE ',lg,ty) + if (lg == 0) then + mergeAV(:lgAv) = (/(i,i=1,lgAv)/) + else + call LCMGET(geoIp,'MERGE ',mergeAV) + end if + + pos = 0 + if ((ncode(1)==B_Diag) .and. (ncode(4)==B_Diag)) then !triangle inf + do i = 1,n + do j = 1,i-1 + k = i + ((j-1)*(2*n-j))/2 + pos = pos + 1 + mixAp(pos) = mixAv(k) + mergeAP(pos) = mergeAV(k) + turnAp(pos) = composeRotRec(turnAv(k),6) + end do + do j = i,n + k = j + ((i-1)*(2*n-i))/2 + pos = pos + 1 + mixAp(pos) = mixAv(k) + mergeAP(pos) = mergeAV(k) + turnAp(pos) = turnAv(k) + end do + end do + else if ((ncode(2)==B_Diag) .and. (ncode(3)==B_Diag)) then !triangle sup + do i = 1,n + do j = 1,i + k = j + (i*(i-1))/2 + pos = pos + 1 + mixAp(pos) = mixAv(k) + mergeAP(pos) = mergeAV(k) + turnAp(pos) = turnAv(k) + end do + do j = i+1,n + k = i + (j*(j-1))/2 + pos = pos + 1 + mixAp(pos) = mixAv(k) + mergeAP(pos) = mergeAV(k) + turnAp(pos) = composeRotRec(turnAv(k),6) + end do + end do + else + call XABORT("G2S: internal error in routine traiteConditionDiagonale") + end if + + call LCMPUT(geoIp,'MIX ',pos,1,mixAp) + call LCMPUT(geoIp,'TURN ',pos,1,turnAp) + call LCMPUT(geoIp,'MERGE ',pos,1,mergeAp) + deallocate(mixAv,turnAv,mixAp,turnAp,mergeAv,mergeAp) + end subroutine traiteConditionDiagonale + + subroutine creeListeCellules(lcIp,geoIp,nbCarre) + type(c_ptr),intent(in) :: lcIp,geoIp + integer,intent(out) :: nbCarre + + integer :: nbAV + + nbCarre = 0 + call rempliListeCellules(lcIp,geoIp,nbCarre) + do + nbAV = nbCarre + call rempliListeCellules(lcIp,geoIp,nbCarre) + if (nbCarre==nbAV) exit + end do + end subroutine creeListeCellules + + subroutine rempliListeCellules(lcIp,geoIp,nbCarre) + type(c_ptr),intent(in):: lcIp,geoIp + integer,intent(inout) :: nbCarre + + character*12 :: subName + type(c_ptr) :: subIp,jplist + integer :: i,j,long,typ,nbCell,long2 + character*12,dimension(:),allocatable :: dataCh12 + + if (nbCarre==0) then + nbCarre = 1 + jplist=LCMDID(lcIp,'root ') + call LCMEQU(geoIp,jplist) + end if + subName = ' ' + do i = 1,nbCarre + call LCMNXT(lcIp,subName) + subIp=LCMGID(lcIp,subName) + call LCMLEN(subIp,'CELL ',long,typ) + if (long/=0) then + nbCell = long/3 + allocate(dataCh12(nbCell),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: rempliListeCellules => allocation pb") + call LCMGTC(subIp,'CELL ',12,nbCell,dataCh12) + do j = 1,nbCell + call LCMLEN(subIp,dataCh12(j),long,typ) + call LCMLEN(lcIp,dataCh12(j),long2,typ) + if (long/=0.and.long2==0) then + call LCMSIX(subIp,dataCh12(j),1) + nbCarre = nbCarre + 1 + jplist=LCMDID(lcIp,dataCh12(j)) + call LCMEQU(subIp,jplist) + call LCMSIX(subIp,dataCh12(j),2) + end if + end do + deallocate(dataCh12) + end if + end do + end subroutine rempliListeCellules + + subroutine creeDico(dicoIp,lcIp,szLc) + type(c_ptr),intent(in) :: dicoIp,lcIp + integer,intent(in) :: szLc + + character*12 :: carreName,sideName,tmpStr,arrayName,nbStr + type(c_ptr) :: cIp,dIp,dJp + integer :: i,j,k,long,typ,lenMix,ind,lignNbr,lenCell + real :: side + integer,dimension(40) :: sv + integer,dimension(:),allocatable :: mix,turn,turnStr,merg,mergStr + character*12,dimension(:),allocatable :: cell,tabCh + + !initialisation + dIp = dicoIp + carreName = ' ' + do i = 1,szLc + call LCMNXT(lcIp,carreName) + cIp=LCMGID(lcIp,carreName) + sideName = 'x' // carreName(1:11) + do j = 1,2 + call getSquareSide(cIp,j,side) + dJp=LCMDID(dIp,sideName) + call LCMPUT(dJp,'value',1,2,side) + sideName(1:1)='y' + end do + end do + !remplissage + dIp = dicoIp + carreName = ' ' + lignNbr = 0 + do k = 1,szLc + !tavail sur la kieme cellule + call LCMNXT(lcIp,carreName) + cIp=LCMGID(lcIp,carreName) + call LCMGET(cIp,'STATE-VECTOR',sv) + lenMix = sv(3)*sv(4) + if ((sv(1)/=G_Car2d).or.(sv(9)==0)) cycle !pas d'info a retirer + call LCMLEN(cIp,'CELL ',lenCell,typ) + lenCell=lenCell/3 + allocate(mix(lenMix),turn(lenMix),cell(lenCell),merg(lenMix),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: creeDico(1) => allocation pb") + call LCMLEN(cIp,'TURN ',long,typ) + if (long==0) then + turn(:lenMix) = 1 + else + call LCMGET(cIp,'TURN ',turn) + end if + call LCMLEN(cIp,'MERGE ',long,typ) + if (long==0) then + merg(:lenMix) = (/(i,i=1,lenMix)/) + else + call LCMGET(cIp,'MERGE ',merg) + end if + call LCMGET(cIp,'MIX ',mix) + !on teste si tous les milieux sont bien negatifs + do i = 1,lenMix + if (mix(i)<0) cycle + call XABORT("G2S: error, meltig MIX and CELL not supported") + end do + call LCMGTC(cIp,'CELL ',12,lenCell,cell) + !travail sur l'axe des x + allocate(tabCh(2*sv(3)-1),turnStr(sv(3)),mergStr(sv(3)),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: creeDico(2) => allocation pb") + tabCh = '+' + sideName = 'x' // carreName(1:11) + dJp=LCMDID(dIp,sideName) + do j = 1,sv(4) + do i = 1,sv(3) + ind = i+(j-1)*sv(3) + tmpStr = tourne('x',turn(ind)) // cell(-mix(ind))(1:11) + tabCh(2*i-1) = tmpStr + turnStr(i) = turn(ind) + mergStr(i) = merg(ind) + end do + lignNbr = lignNbr + 1 + nbStr = i2s(lignNbr) + arrayName = "array" // nbStr(:7) + call LCMPTC(dJp,arrayName,12,2*sv(3)-1,tabCh) + arrayName = "turns" // nbStr(:7) + call LCMPUT(dJp,arrayName,sv(3),1,turnStr) + arrayName = "merge" // nbStr(:7) + call LCMPUT(dJp,arrayName,sv(3),1,mergStr) + call exploiteStr(tabCh,cIp,sideName,dicoIp,lignNbr) + end do + deallocate(tabCh,turnStr,mergStr) + !travail sur l'axe des y + allocate(tabCh(2*sv(4)-1),turnStr(sv(4)),mergStr(sv(4)),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: creeDico(3) => allocation pb") + tabCh = '+' + sideName(1:1) = 'y' + dJp=LCMDID(dIp,sideName) + do i = 1,sv(3) + do j = 1,sv(4) + ind = i+(j-1)*sv(3) + tmpStr = tourne('y',turn(ind)) // cell(-mix(ind))(1:11) + tabCh(2*j-1) = tmpStr + turnStr(j) = turn(ind) + mergStr(j) = merg(ind) + end do + lignNbr = lignNbr + 1 + nbStr = i2s(lignNbr) + arrayName = "array" // nbStr(:7) + call LCMPTC(dJp,arrayName,12,2*sv(4)-1,tabCh) + arrayName = "turns" // nbStr(:7) + call LCMPUT(dJp,arrayName,sv(4),1,turnStr) + arrayName = "merge" // nbStr(:7) + call LCMPUT(dJp,arrayName,sv(4),1,mergStr) + call exploiteStr(tabCh,cIp,sideName,dicoIp,lignNbr) + end do + deallocate(tabCh,turnStr,mergStr) + !fin du travail sur la cellule + deallocate(mix,turn,cell,merg) + end do + end subroutine creeDico + + function tourne(axe,turn) + character,intent(in) :: axe + integer,intent(in) :: turn + character :: tourne + + tourne = axe + if (mod(turn,2)==1) return + if (axe=='x') then + tourne = 'y' + else + tourne = 'x' + end if + end function tourne + + subroutine exploiteStr(inStr,carreIp,sideName,dicoIp,lignNbr) + character*12,dimension(:),intent(in) :: inStr + type(c_ptr),intent(in) :: carreIp,dicoIp + character*12,intent(in) :: sideName + integer,intent(inout) :: lignNbr + + character*12,dimension((size(inStr)+1)/2) :: tmpStr,workStr + character*12,dimension(:),allocatable :: resStr + real,dimension(:),allocatable :: mesh + character*12 :: name2 + type(c_ptr) :: workIp + integer :: i,j,k,nbOcc,szTS,szIS,nbRest,long,typ + character*12 :: arrayName,name1,meshName,number,text12 + real :: value,res + + workIp = dicoIp + szIS = size(inStr) + szTS = (szIS+1)/2 + !on enleve les "+" dans le tableau + do i = 1,szTS + tmpStr(i)=inStr(2*i-1) + enddo + !creation a partir de r=a+b+c+b+a, de a=r-b-c-b/2, + ! b=r-a-c-a/2 et c=r-a-b-b-a + do i = 1,szTS + nbOcc = 0 + do j = 1,szTS + if (tmpStr(j)==tmpStr(i)) then + nbOcc = nbOcc + 1 + else if (j > nbOcc) then + workStr(j-nbOcc) = tmpStr(j) + end if + end do + nbRest = szTS - nbOcc + if (nbOcc==1) then + allocate(resStr(2*nbRest+1),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: exploiteStr(2) => allocation pb") + else + allocate(resStr(2*nbRest+3),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: exploiteStr(3) => allocation pb") + end if + resStr(1) = sideName + do k = 1,nbRest + resStr(2*k) = '-' + resStr(2*k+1) = workStr(k) + end do + long=2*nbRest+1 + if (nbOcc/=1) then + resStr(2*nbRest+2) = '/' + number = i2s(nbOcc) + resStr(2*nbRest+3) = number + long=2*nbRest+3 + end if + lignNbr = lignNbr + 1 + text12 = i2s(lignNbr) + arrayName = "expl" // text12(:8) + call LCMSIX(workIp,tmpStr(i),1) + call LCMPTC(workIp,arrayName,12,long,resStr) + call LCMSIX(workIp,tmpStr(i),2) + deallocate(resStr) + !exploitation orthogonale de la chaine ('xa + xb' => ya=yb) + name1 = tourne(tmpStr(i)(1:1),2) // tmpStr(i)(2:12) + do k = 1,nbRest + name2 = tourne(workStr(k)(1:1),2) // workStr(k)(2:12) + lignNbr = lignNbr + 1 + arrayName = "orth" // text12(:8) + call LCMSIX(workIp,name1,1) + call LCMPTC(workIp,arrayName,12,name2) + call LCMSIX(workIp,name1,2) + end do + end do + !exploitation eventuelle du maillage en entree + if (sideName(1:1)=='x') then + meshName = 'MESHX' + else if (sideName(1:1)=='y') then + meshName = 'MESHY' + else + call XABORT("G2S: internal error in subroutine exploiteStr") + end if + call LCMLEN(carreIp,meshName,long,typ) + if (long/=0) then + allocate(mesh(long),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: exploiteStr(4) => allocation pb") + call LCMGET(carreIp,meshName,mesh) + do i = 1,szTS + call getValueOf(tmpStr(i),workIp,value) + res = mesh(i+1)-mesh(i) + if ((value>0.) .and. (abs(value-res)>geomPrec)) then + write(6,*) " ",meshName," -> ",tmpStr(i) + call XABORT("G2S: error, incoherent geometry(1)") + endif + call setValueOf(tmpStr(i),workIp,res) + end do + deallocate(mesh) + end if + end subroutine exploiteStr + + subroutine getSquareSide(cIp,axis,side) + type(c_ptr),intent(in) :: cIp + integer,intent(in) :: axis + real,intent(out) :: side + + integer :: typ,long + real,dimension(:),allocatable :: tr + + side = -1. + if (axis==1) then !axe des x + call LCMLEN(cIp,'MESHX ',long,typ) + if (long/=0) then + allocate(tr(long),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: getSquareSide(1) => allocation pb") + call LCMGET(cIp,'MESHX ',tr) + side = tr(long) - tr(1) + deallocate(tr) + end if + else if (axis==2) then !axe des y + call LCMLEN(cIp,'MESHY ',long,typ) + if (long/=0) then + allocate(tr(long),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: getSquareSide(2) => allocation pb") + call LCMGET(cIp,'MESHY ',tr) + side = tr(long) - tr(1) + deallocate(tr) + end if + else + call XABORT("G2S: internal error in subroutine getSquareSide") + end if + end subroutine getSquareSide + + subroutine getValueOf(sideName,dicoIp,value) + character*12,intent(in) :: sideName + type(c_ptr),intent(in) :: dicoIp + real,intent(out) :: value + + type(c_ptr) :: ip,jp + character*131 :: hsmg + integer :: long,typ + + ip = dicoIp + jp = LCMGID(ip,sideName) + call LCMLEN(jp,'value',long,typ) + if(long == 0) then + write(hsmg,'(36hG2S: value key missing for sideName=,a,5h (1).)') trim(sideName) + call LCMLIB(ip) + call LCMLIB(jp) + call XABORT(hsmg) + endif + call LCMGET(jp,'value',value) + end subroutine getValueOf + + subroutine setValueOf(sideName,dicoIp,value) + character*12,intent(in) :: sideName + type(c_ptr),intent(in) :: dicoIp + real,intent(in) :: value + + type(c_ptr) :: ip,jp + + ip = dicoIp + jp = LCMDID(ip,sideName) + call LCMPUT(jp,'value',1,2,value) + end subroutine setValueOf + + subroutine resoudDico(dicoIp,resIp) + type(c_ptr),intent(in) :: dicoIp,resIp + + logical :: newVal,flag1,flag2,empty,lcm + character*12 :: sideName,saveSideName,namlcm,nammy + character*131:: hsmg + type(c_ptr) :: ip,jp + integer :: ilong,long,typ + real :: val + + !la resolution se fait en travaillant sur toutes les variables une a une + !on resout les unes apres les autres toutes les equations regissant + !chaque variable (si possible) + ! on teste a chaque passe si on a trouve une nouvelle valeur + ! si oui on continue, sinon (apres une nouvelle passe infructueuse), on + ! arrete le traitement et on teste si tous les resultats sont bons + ! oui -> on sort ; non -> on abandonne car manque de donnees + + flag1 = .true. + do + flag2 = flag1 !garde la valeur de flag1 a l'essai precedent + flag1 = .false. + sideName = ' ' + call LCMINF(dicoIp,namlcm,nammy,empty,ilong,lcm) + if (.not. empty) then + call LCMNXT(dicoIp,sideName) + saveSideName = sideName + do !boucle sur toutes les donnees du dico + newVal = evaluateEquation(dicoIp,sideName) + !test si newVal a ete vrai au moins une fois + if (newVal) flag1 = .true. + call LCMNXT(dicoIp,sideName) + if (sideName == saveSideName) exit + end do + end if + if (.not.(flag1 .or. flag2)) then + !il ne s'est rien passe au cours des 2 tentatives + !precedentes => on abandonne si une valeur est non resolue + saveSideName = sideName + do !boucle sur toutes les donnees du dico + call getValueOf(sideName,dicoIp,val) + if (val<0) then !une valeur est non resolue -> ABORT + write(6,*) "Curent values of geometry :" + saveSideName = sideName + do !boucle sur toutes les donnees du dico + call getValueOf(sideName,dicoIp,val) + write(6,*) " ",sidename," -> ",val + call LCMNXT(dicoIp,sideName) + if (sideName == saveSideName) exit + end do + call XABORT("G2S: not enought data in the geometry(1)") + end if + call LCMNXT(dicoIp,sideName) + if (sideName == saveSideName) exit + end do + exit !toutes les valeurs sont bonnes + end if + end do + !recopie des resultats + ip = dicoIp + sideName = ' ' + call LCMINF(ip,namlcm,nammy,empty,ilong,lcm) + if (empty) return + call LCMNXT(ip,sideName) + saveSideName = sideName + do + jp=LCMGID(ip,sideName) + call LCMLEN(jp,'value',long,typ) + if(long == 0) then + write(hsmg,'(36hG2S: value key missing for sideName=,a,5h (2).)') trim(sideName) + call LCMLIB(ip) + call LCMLIB(jp) + call XABORT(hsmg) + endif + call LCMGET(jp,'value',val) + call LCMPUT(resIp,sideName,1,2,val) + call LCMNXT(ip,sideName) + if (sideName == saveSideName) exit + end do + end subroutine resoudDico + + function evaluateEquation(dicoIp,sideName) + type(c_ptr),intent(in) :: dicoIp + character*12,intent(in) :: sideName + logical :: evaluateEquation + !evaluation une a une des equations de dico relatives a sideName + !retourne .true. si nouvelle affectation et .false. sinon + + real :: res,val + double precision :: dres + type(c_ptr) :: sideIp + integer :: i,long,typ,ilong + logical :: goodLine,empty,lcm + character*12 :: eqName,saveEqName,namlcm,nammy + character*12,dimension(:),allocatable :: str + + evaluateEquation = .false. + sideIp = dicoIp + call LCMSIX(sideIp,sideName,1) + eqName = ' ' + call LCMINF(sideIp,namlcm,nammy,empty,ilong,lcm) + if (empty) call XABORT("G2S: intenal error in data structure in & + &function evaluateEquation") + call LCMNXT(sideIp,eqName) + saveEqName = eqName + do !recherche de la premiere occurence d'une arrayXXX + call LCMNXT(sideIp,eqName) + if (eqName(1:5)=='array') exit + if (eqName==saveEqName) then + !pas de donnee interessante => on quitte + return + end if + end do + saveEqName = eqName + do !boucle sur toutes les equations de sideName + if (eqName(1:5)/='array') then + call LCMNXT(sideIp,eqName) + if (eqName==saveEqName) exit + cycle + end if + !on travaille sur une arrayXXX + goodLine = .true. + call LCMLEN(sideIp,eqName,long,typ) + long = long/3 + allocate(str(long),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: evaluateEquation => allocation pb") + call LCMGTC(sideIp,eqName,12,long,str) + call getValueOf(str(1),dicoIp,res) + if (res<0) then + goodLine = .false. + end if + dres = res + do i = 2,long,2 + if (str(i)=='+') then + call getValueOf(str(i+1),dicoIp,val) + if (val<0) then + goodLine = .false. + exit + end if + dres = dres + val + else if (str(i)=='-') then + call getValueOf(str(i+1),dicoIp,val) + if (val<0) then + goodLine = .false. + exit + end if + dres = dres - val + else if (str(i)=='/') then + val = s2i(str(i+1)) + dres = dres / val + else + call XABORT("G2S: internal error in function evaluateEquation") + end if + end do + res = real(dres) + !si on sort normalement de la boucle, on teste la coherence + !ou on assigne le resultat obtenu + if (goodLine) then + call getValueOf(sideName,dicoIp,val) + if (val<0) then + call setValueOf(sideName,dicoIp,res) + evaluateEquation = .true. !on a resolu une nouvelle equation + else if (abs(val-res)>geomPrec) then + call XABORT("G2S: error, incoherent geometry(2)") + end if + end if + deallocate(str) + call LCMNXT(sideIp,eqName) + if (eqName==saveEqName) exit + end do + end function evaluateEquation + + subroutine prepareCentres(lCentreIpIn,lcIp,sidesIp,dicoIpIn) + type(c_ptr),intent(in) :: lCentreIpIn,lcIp,sidesIp,dicoIpIn + + type(c_ptr) :: tmpIp,coordCentIp,ip,dicoIp,lCentreIp + integer :: i,j,lg,lgc,ty,dimValToPut,nbArray,gigNum,ilong + logical :: empty,lcm + character*12 :: namlcm,nammy + character*12 :: carreName,saveCarreName,tmpName,eqName,saveEqName + character*12 :: sideName,saveSideName,arrayName,saveArrayName + character*12 :: turnName,mergeName + real :: saveVal,val + real,dimension(:),allocatable :: valToPut,xx,yy + integer,dimension(40) :: sv + integer,dimension(:),allocatable :: turns,merg + character*12,dimension(:),allocatable :: eqStr,arrayNameStr + + lCentreIp = lCentreIpIn + dicoIp = dicoIpIn + !initialisation de la liste et creation et initialisation + !d'une liste temporaire a partir du dictionnaire utilise a + !l'etape precedente, et d'un autre temporaire appelle coordCent + call LCMOP(tmpIp,'lCtemp ',0,1,0) + carreName = ' ' + call LCMINF(lcIp,namlcm,nammy,empty,ilong,lcm) + if (empty) call XABORT("G2S: intenal error in data structure in & + &subroutine prepareCentres") + call LCMNXT(lcIp,carreName) + saveCarreName = carreName + do + !travail sur lCtemp + call LCMSIX(tmpIp,carreName,1) + tmpName = 'x' // carreName(1:11) + call LCMSIX(dicoIp,tmpName,1) + eqName = ' ' + call LCMNXT(dicoIp,eqName) + saveEqName = eqName + do + if (eqName(1:5)=='array') then !l'equation est a conserver + call putEquationIn(dicoIp,eqName,tmpIp) + end if + call LCMNXT(dicoIp,eqName) + if (eqName==saveEqName) exit + end do + call LCMSIX(dicoIp,tmpName,2) + call LCMSIX(tmpIp,carreName,2) + !travail sur lCentre + call LCMSIX(lCentreIp,carreName,1) + call LCMSIX(lCentreIp,carreName,2) + call LCMNXT(lcIp,carreName) + if (carreName==saveCarreName) exit + end do + + !creation d'un autre temporaire appelle coordCent + call LCMOP(coordCentIp,'lCoordCent ',0,1,0) + sideName = ' ' + call LCMNXT(dicoIp,sideName) + saveSideName = sideName + do + call LCMSIX(dicoIp,sideName,1) + eqName = ' ' + call LCMNXT(dicoIp,eqName) + saveEqName = eqName + do + if (eqName(1:5)=='array') then !l'equation est a exploiter + call LCMLEN(dicoIp,eqName,lg,ty) + lg = lg / 3 + dimValToPut = (lg+1) / 2 + allocate(eqStr(lg),valToPut(dimValToPut),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT("G2S: prepareCentres(1) => allocation pb") + call LCMGTC(dicoIp,eqName,12,lg,eqStr) + !calcul des coordonnees des sous cellules de la cellule + !consideree et stockage de la valeur obtenu dans coordCent + call LCMGET(sidesIp,sideName,val) + valToPut(1) = -0.5 * val + call LCMGET(sidesIp,eqStr(1),val) + valToPut(1) = valToPut(1) + 0.5 * val + do i = 2,dimValToPut + saveVal = val + call LCMGET(sidesIp,eqStr(2*i-1),val) + valToPut(i) = valToPut(i-1) + 0.5 * (val+saveVal) + end do + call LCMPUT(coordCentIp,sideName,dimValToPut,2,valToPut) + deallocate(eqStr,valToPut) + exit !on sort car on n'a pas a verifier la coherence + end if + call LCMNXT(dicoIp,eqName) + if (eqName==saveEqName) exit + end do + call LCMSIX(dicoIp,sideName,2) + call LCMNXT(dicoIp,sideName) + if (sideName==saveSideName) exit + end do + + carreName = ' ' + call LCMNXT(lcIp,carreName) + saveCarreName = carreName + do + ip=LCMGID(lcIp,carreName) + call LCMGET(ip,'STATE-VECTOR',sv) + !recuperation des noms des equations utiles + call LCMSIX(tmpIp,carreName,1) + call LCMINF(tmpIp,namlcm,nammy,empty,ilong,lcm) + if (.not. empty) then + allocate(arrayNameStr(max(sv(3),sv(4)))) + nbArray = 0 + arrayName = ' ' + call LCMNXT(tmpIp,arrayName) + saveArrayName = arrayName + do + if (arrayName(1:5)/='array') then + call LCMNXT(tmpIp,arrayName) + if (arrayName == saveArrayName) exit + cycle + end if + nbArray = nbArray + 1 + arrayNameStr(nbArray) = arrayName + call LCMNXT(tmpIp,arrayName) + if (arrayName == saveArrayName) exit + end do + !trie des noms des equations + call sortTabCh12(arrayNameStr,nbArray) + !travail sur les equations + gigNum = 0 + do i = 1,nbArray + arrayName = arrayNameStr(i) + call LCMLEN(tmpIp,arrayName,lgc,ty) + lgc = lgc / 3 + allocate(eqStr(lgc)) + lg = (lgc+1) / 2 + allocate(turns(lg),merg(lg),xx(lg),yy(nbArray)) + call LCMGTC(tmpIp,arrayName,12,lgc,eqStr) + turnName = 'turns' // arrayName(6:12) + call LCMGET(tmpIp,turnName,turns) + mergeName = 'merge' // arrayName(6:12) + call LCMGET(tmpIp,mergeName,merg) + do j = 1,lg + gigNum = gigNum + 1 + sideName = 'x' // carreName(1:11) + call LCMGET(coordCentIp,sideName,xx) + sideName(1:1) = 'y' + call LCMGET(coordCentIp,sideName,yy) + sideName=eqStr(2*j-1)(2:12) + !ajout du resultat dans lCentre + call addInRes(lCentreIp,sideName,carreName, & + & xx(j),yy(i),turns(j),gigNum,merg(j)) + end do + deallocate(eqStr,turns,merg,xx,yy) + end do + deallocate(arrayNameStr) + end if + + call LCMSIX(tmpIp,carreName,2) + call LCMNXT(lcIp,carreName) + if (carreName==saveCarreName) exit + end do + !destruction du repertoire 'root' dans lCentre + call LCMDEL(lCentreIp,'root ') + !fermeture des temporaires + call LCMCL(coordCentIp,2) + call LCMCL(tmpIp,2) + end subroutine prepareCentres + + subroutine addInRes(ipIn,extName,intName,xx,yy,turn,gigNum,merg) + type(c_ptr),intent(in) :: ipIn + integer,intent(in) :: turn,gigNum,merg + character*12,intent(in) :: extName,intName + real,intent(in) :: xx,yy + + real,dimension(:),allocatable :: tabR + integer,dimension(:),allocatable :: tabI + character*12 :: reg,gigName + integer :: lg,ty + type(c_ptr) :: ip + + ip = ipIn + call LCMSIX(ip,extName,1) + call LCMSIX(ip,intName,1) + reg = 'x' + call LCMLEN(ip,reg,lg,ty) + allocate(tabR(lg+1)) + !ajout de xx + if (lg /= 0) call LCMGET(ip,reg,tabR) + tabR(lg+1) = xx + call LCMPUT(ip,reg,lg+1,2,tabR) + !ajout de yy + reg = 'y' + if (lg /= 0) call LCMGET(ip,reg,tabR) + tabR(lg+1) = yy + call LCMPUT(ip,reg,lg+1,2,tabR) + deallocate(tabR) + !ajout de turn + reg = 't' + allocate(tabI(lg+1)) + if (lg /= 0) call LCMGET(ip,reg,tabI) + tabI(lg+1) = turn + call LCMPUT(ip,reg,lg+1,1,tabI) + deallocate(tabI) + !ajout de la gigogne et du merge + gigName = i2s(lg+1) + call LCMSIX(ip,gigName,1) + reg = 'gig' + call LCMPUT(ip,reg,1,1,gigNum) + reg = 'mrg' + call LCMPUT(ip,reg,1,1,merg) + call LCMSIX(ip,gigName,2) + call LCMSIX(ip,intName,2) + call LCMSIX(ip,extName,2) + end subroutine addInRes + + subroutine sortTabCh12(tab,sz) + character*12,dimension(:),intent(inout) :: tab + integer,intent(in) :: sz + + integer :: indMax,i,s1,s2 + logical :: trie + character*12 :: tmp + + do indMax = sz,2,-1 + trie = .true. + do i = 1,indMax-1 + s1 = s2i(tab(i)(6:12)) + s2 = s2i(tab(i+1)(6:12)) + if (s2<s1) then + tmp=tab(i+1) ; tab(i+1)=tab(i) ; tab(i)=tmp + trie = .false. + end if + end do + if (trie) exit + end do + end subroutine sortTabCh12 + + subroutine putEquationIn(dicoIp,eqName,lCentreIp) + type(c_ptr),intent(in) :: dicoIp,lCentreIp + character*12,intent(in) :: eqName + + character*12,dimension(:),allocatable :: str12 + integer,dimension(:),allocatable :: turns,merg + character*12 :: turnsName,mergeName + integer :: lg,ty,lgt + + call LCMLEN(dicoIp,eqName,lg,ty) + if (lg/=0) then + turnsName = 'turns' // eqName(6:12) + mergeName = 'merge' // eqName(6:12) + lg = lg/3 + lgt = (lg+1)/2 + allocate(str12(lg),turns(lgt),merg(lgt)) + call LCMGTC(dicoIp,eqName,12,lg,str12) + call LCMGET(dicoIp,turnsName,turns) + call LCMGET(dicoIp,mergeName,merg) + call LCMPTC(lCentreIp,eqName,12,lg,str12) + call LCMPUT(lCentreIp,turnsName,lgt,1,turns) + call LCMPUT(lCentreIp,mergeName,lgt,1,merg) + deallocate(str12,turns,merg) + else + call XABORT("G2S: internal error in subroutine putEquationIn") + end if + end subroutine putEquationIn + + subroutine prepareMinMax(lMinMaxIp,lCentreIp) + type(c_ptr),intent(in) :: lMinMaxIp,lCentreIp + + character*12 :: carreName,saveCarreName,namlcm,nammy + type(c_ptr) :: ipLC + integer :: lg,ty,ilong + logical :: empty,lcm + real :: minX,maxX,minY,maxY + real,dimension(:),allocatable :: xx,yy + + ipLC = lCentreIp + minX = 0. ; maxX = 0. ; minY = 0. ; maxY = 0. + carreName = ' ' + call LCMINF(ipLC,namlcm,nammy,empty,ilong,lcm) + if (.not. empty) then !plus d'une seule cellule + call LCMNXT(ipLC,carreName) + saveCarreName = carreName + do + call LCMSIX(ipLC,carreName,1) + call LCMLEN(ipLC,'root ',lg,ty) + if (lg/=0) then + call LCMSIX(ipLC,'root ',1) + call LCMLEN(ipLC,'x ',lg,ty) + allocate(xx(lg)) + call LCMGET(ipLC,'x ',xx) + call LCMLEN(ipLC,'y ',lg,ty) + allocate(yy(lg)) + call LCMGET(ipLC,'y ',yy) + minX = min(minX,minval(xx)) + maxX = max(maxX,maxval(xx)) + minY = min(minY,minval(yy)) + maxY = max(maxY,maxval(yy)) + deallocate(xx,yy) + call LCMSIX(ipLC,'root ',2) + end if + call LCMSIX(ipLC,carreName,2) + call LCMNXT(ipLC,carreName) + if (carreName==saveCarreName) exit + end do + end if + call LCMPUT(lMinMaxIp,'minX ',1,2,minX) + call LCMPUT(lMinMaxIp,'maxX ',1,2,maxX) + call LCMPUT(lMinMaxIp,'minY ',1,2,minY) + call LCMPUT(lMinMaxIp,'maxY ',1,2,maxY) + end subroutine prepareMinMax + + recursive subroutine calculeCentre(resIp,lcIp,deep) + type(c_ptr),intent(inout) :: resIp + type(c_ptr),intent(in) :: lcIp + integer,intent(in) :: deep + + character*12 :: dirName,saveDirName,subDirName,saveSubDirName + character*12 :: namlcm,nammy + type(c_ptr) :: valuesIp,tmpLcIp,cpResIp,tmpDic + integer :: i,lg,ty,nbGig,ilong + logical :: empty,lcm + + call LCMEQU(lcIp,resIp) + tmpLcIp = lcIp + dirName = ' ' + call LCMINF(tmpLcIp,namlcm,nammy,empty,ilong,lcm) + if (empty) then + !une seule cellule => traitement particulier + call LCMSIX(resIp,'root ',1) + call LCMSIX(resIp,'root ',1) + call LCMPUT(resIp,'x ',1,2,0.) + call LCMPUT(resIp,'y ',1,2,0.) + call LCMPUT(resIp,'t ',1,1,1) + call LCMSIX(resIp,'1 ',1) + call LCMPUT(resIp,'gig ',1,1,1) + call LCMPUT(resIp,'mrg ',1,1,1) + call LCMSIX(resIp,'1 ',2) + call LCMSIX(resIp,'root ',2) + call LCMSIX(resIp,'root ',2) + return + end if + call LCMOP(valuesIp,'temporariObj',0,1,0) + call LCMNXT(tmpLcIp,dirName) + saveDirName = dirName + nbGig = 0 + do + call LCMSIX(tmpLcIp,dirName,1) + subDirName = ' ' + call LCMINF(tmpLcIp,namlcm,nammy,empty,ilong,lcm) + if (.not. empty) then + call LCMNXT(tmpLcIp,subDirName) + saveSubDirName = subDirName + do + call LCMLEN(valuesIp,subDirName,lg,ty) + if (lg==0) then + nbGig = nbGig + 1 + call LCMSIX(valuesIp,subDirName,1) + call LCMSIX(valuesIp,subDirName,2) + end if + call LCMNXT(tmpLcIp,subDirName) + if (subDirName==saveSubDirName) exit + end do + end if + call LCMSIX(tmpLcIp,dirName,2) + call LCMNXT(tmpLcIp,dirName) + if (dirName==saveDirName) exit + end do + !test de sortie de la recurrence + if (nbGig <= 1) then + call LCMEQU(lcIp,resIp) + call LCMCL(valuesIp,2) + return + end if + !liste des cellules a garder + dirName = ' ' + call LCMNXT(lcIp,dirName) + saveDirName = dirName + do + call LCMLEN(valuesIp,dirName,lg,ty) + if (lg==0) call LCMDEL(resIp,dirName) + call LCMNXT(lcIp,dirName) + if (dirName==saveDirName) exit + end do + + dirName = ' ' + call LCMNXT(lcIp,dirName) + saveDirName = dirName + do + call LCMLEN(valuesIp,dirName,lg,ty) + if (lg==0) then + call LCMSIX(tmpLcIp,dirName,1) + subDirName = ' ' + call LCMNXT(tmpLcIp,subDirName) + saveSubDirName = subDirName + do + call LCMSIX(tmpLcIp,subDirName,1) + if (subDirName/='root') then + call LCMLEN(tmpLcIp,'t ',lg,ty) + do i = 1,lg + call LCMOP(tmpDic,'tmpDic ',0,1,0) + call translateDic(tmpDic,lcIp,dirName,subDirName,i) + call mergeDic(resIp,dirName,tmpDic) + call LCMCL(tmpDic,2) + end do + else + call LCMOP(tmpDic,'tmpDic ',0,1,0) + call LCMSIX(tmpDic,subDirName,1) + call LCMEQU(tmpLcIp,tmpDic) !copie uniquement de la branche + !interresante (lc[dirname]['root']) + call LCMSIX(tmpDic,subDirName,2) + call mergeDic(resIp,dirName,tmpDic) + call LCMCL(tmpDic,2) + end if + call LCMSIX(tmpLcIp,subDirName,2) + call LCMNXT(tmpLcIp,subDirName) + if (subDirName==saveSubDirName) exit + end do + call LCMSIX(tmpLcIp,dirName,2) + end if + call LCMNXT(lcIp,dirName) + if (dirName==saveDirName) exit + end do + !iteration + call LCMOP(cpResIp,'resCentres ',0,1,0) + call LCMEQU(resIp,cpResIp) + call calculeCentre(resIp,cpResIp,deep+1) + call LCMCL(cpResIp,2) + !elimination des cellules non terminales si deep==0 + if (deep==0) then + dirName = ' ' + call LCMNXT(valuesIp,dirName) + saveDirName = dirName + do + call LCMLEN(resIp,dirName,lg,ty) + if (lg/=0) call LCMDEL(resIp,dirName) + call LCMNXT(valuesIp,dirName) + if (dirName==saveDirName) exit + end do + end if + call LCMCL(valuesIp,2) + end subroutine calculeCentre + + subroutine translateDic(resIp,dicIp,dName,sDName,ind) + type(c_ptr),intent(in) :: resIp,dicIp + integer,intent(in) :: ind + character*12,intent(in) :: dName,sDName + !effectue la translation dans le repere de la cellule englobante + !de la liste des cellules de dicIp[sDName], suivant les vecteurs de + !dicIp[dName][sDName] (qui comprennent x(i),y(i),t(i),[`i`][gig] i=1,sz) + + character*12 :: dirName,saveDirName,gigName + type(c_ptr) :: tmpRes,tmpDic + integer :: i,lg,ty,tval,lggigI,lggigval + real :: xval,yval + integer,dimension(:),allocatable :: tt,gigI,gigval,mrgI,mrgval + real,dimension(:),allocatable :: xx,yy + + tmpRes = resIp ; tmpDic = dicIp + !recuperation des donnees de translation + call LCMSIX(tmpDic,dName,1) + call LCMSIX(tmpDic,sDName,1) + call LCMLEN(tmpDic,'t ',lg,ty) + allocate(xx(lg),yy(lg),tt(lg)) + call LCMGET(tmpDic,'t ',tt) ; tval = tt(ind) + call LCMGET(tmpDic,'x ',xx) ; xval = xx(ind) + call LCMGET(tmpDic,'y ',yy) ; yval = yy(ind) + deallocate(xx,yy,tt) + gigName = i2s(ind) + call LCMSIX(tmpDic,gigName,1) + call LCMLEN(tmpDic,'gig ',lggigval,ty) + allocate(gigval(lggigval),mrgval(lggigval)) + call LCMGET(tmpDic,'gig ',gigval) + call LCMGET(tmpDic,'mrg ',mrgval) + tmpDic = dicIp + !translation + call LCMSIX(tmpDic,sDName,1) + dirName = ' ' + call LCMNXT(tmpDic,dirName) + saveDirName = dirName + do + call LCMSIX(tmpDic,dirName,1) + call LCMSIX(tmpRes,dirName,1) + call LCMLEN(tmpDic,'t ',lg,ty) + allocate(xx(lg),yy(lg),tt(lg)) + call LCMGET(tmpDic,'t ',tt) + call LCMGET(tmpDic,'x ',xx) + call LCMGET(tmpDic,'y ',yy) + do i = 1,lg + gigName = i2s(i) + call LCMSIX(tmpDic,gigName,1) + call LCMSIX(tmpRes,gigName,1) + call LCMLEN(tmpDic,'gig ',lggigI,ty) + allocate(gigI(lggigI+lggigval),mrgI(lggigI+lggigval)) + call LCMGET(tmpDic,'gig ',gigI) + call LCMGET(tmpDic,'mrg ',mrgI) + call translateWithTurn(xx(i),yy(i),tt(i),xval,yval,tval) + gigI(lggigI+1:lggigI+lggigval) = gigval + call LCMPUT(tmpRes,'gig ',lggigI+lggigval,1,gigI) + mrgI(lggigI+1:lggigI+lggigval) = mrgval + call LCMPUT(tmpRes,'mrg ',lggigI+lggigval,1,mrgI) + deallocate(gigI,mrgI) + call LCMSIX(tmpRes,gigName,2) + call LCMSIX(tmpDic,gigName,2) + end do + call LCMPUT(tmpRes,'t ',lg,1,tt) + call LCMPUT(tmpRes,'x ',lg,2,xx) + call LCMPUT(tmpRes,'y ',lg,2,yy) + deallocate(xx,yy,tt) + call LCMSIX(tmpRes,dirName,2) + call LCMSIX(tmpDic,dirName,2) + call LCMNXT(tmpDic,dirName) + if (dirName==saveDirName) exit + end do + deallocate(gigval,mrgval) + end subroutine translateDic + + subroutine translateWithTurn(xx,yy,tt,xval,yval,tval) + real,intent(inout) :: xx,yy + integer,intent(inout) :: tt + real,intent(in) :: xval,yval + integer,intent(in) :: tval + + select case(geomTyp) + case(RecTyp) + select case (tt) + case(1) + xx=xx+xval ; yy=yy+yval + case(2) + xx=xx+yval ; yy=yy-xval + case(3) + xx=xx-xval ; yy=yy-yval + case(4) + xx=xx-yval ; yy=yy+xval + case(5) + xx=xx-xval ; yy=yy+yval + case(6) + xx=xx+yval ; yy=yy+xval + case(7) + xx=xx+xval ; yy=yy-yval + case(8) + xx=xx-yval ; yy=yy-xval + end select + tt = composeRotRec(tval,tt) + case(TriaTyp) + select case (tt) + case(1) + xx=xx+xval ; yy=yy+yval + case(2) + xx=xx+0.5*(xval-sqrt_3f*yval) + yy=yy+0.5*(yval+sqrt_3f*xval) + case(3) + xx=xx+0.5*(-xval-sqrt_3f*yval) + yy=yy+0.5*(-yval+sqrt_3f*xval) + case(4) + xx=xx-xval ; yy=yy-yval + case(5) + xx=xx+0.5*(-xval+sqrt_3f*yval) + yy=yy+0.5*(-yval-sqrt_3f*xval) + case(6) + xx=xx+0.5*(xval+sqrt_3f*yval) + yy=yy+0.5*(yval-sqrt_3f*xval) + case(7) + xx=xx+xval ; yy=yy-yval + case(8) + xx=xx+0.5*(xval-sqrt_3f*yval) + yy=yy-0.5*(yval+sqrt_3f*xval) + case(9) + xx=xx+0.5*(-xval-sqrt_3f*yval) + yy=yy-0.5*(-yval+sqrt_3f*xval) + case(10) + xx=xx-xval ; yy=yy+yval + case(11) + xx=xx+0.5*(-xval+sqrt_3f*yval) + yy=yy+0.5*(yval+sqrt_3f*xval) + case(12) + xx=xx+0.5*(xval+sqrt_3f*yval) + yy=yy-0.5*(yval-sqrt_3f*xval) + end select + tt = composeRotTri(tval,tt) + case default + call XABORT("G2S: internal error in subroutine translateWithTurn") + end select + end subroutine translateWithTurn + + subroutine mergeDic(resIp,dirName,otherDicIp) + type(c_ptr),intent(in) :: resIp,otherDicIp + character*12,intent(in) :: dirName + !concatene deux dictionnaires (res[dirName] et otherDic) + !contenants des positions (x,y,t,(1:gig),..,(size(t):gig)) + + character*12 :: sdName,nGig,namlcm,nammy + type(c_ptr) :: tmpResIp,tmpOtherIp + integer :: i,lg,ty,lgN,lgO,ilong + logical :: empty,lcm + real,dimension(:),allocatable :: nra,ora + integer,dimension(:),allocatable :: nia,oia,mrg + + tmpResIp = resIp ; tmpOtherIp = otherDicIp + call LCMSIX(tmpResIp,dirName,1) + call LCMINF(tmpOtherIp,namlcm,nammy,empty,ilong,lcm) + if (empty) return !### + sdName = ' ' + call LCMNXT(tmpOtherIp,sdName) + call LCMSIX(tmpOtherIp,sdName,1) + call LCMSIX(tmpResIp,sdName,1) + call LCMLEN(tmpResIp,'t ',lgN,ty) + if (lgN/=0) then !le repertoire existait deja => ajout de donnees + call LCMLEN(tmpOtherIp,'t ',lgO,ty) + allocate(nra(lgO+lgN),ora(lgO),nia(lgO+lgN),oia(lgO)) + !t + call LCMGET(tmpResIp,'t ',nia) + call LCMGET(tmpOtherIp,'t ',oia) + nia(lgN+1:lgO+lgN) = oia + call LCMPUT(tmpResIp,'t ',lgO+lgN,1,nia) + !x + call LCMGET(tmpResIp,'x ',nra) + call LCMGET(tmpOtherIp,'x ',ora) + nra(lgN+1:lgO+lgN) = ora + call LCMPUT(tmpResIp,'x ',lgO+lgN,2,nra) + !y + call LCMGET(tmpResIp,'y ',nra) + call LCMGET(tmpOtherIp,'y ',ora) + nra(lgN+1:lgO+lgN) = ora + call LCMPUT(tmpResIp,'y ',lgO+lgN,2,nra) + deallocate(nra,ora,nia,oia) + lg = lgO + do i = 1,lg + nGig = i2s(i) + call LCMSIX(tmpOtherIp,nGig,1) + call LCMLEN(tmpOtherIp,'gig ',lgO,ty) + allocate(oia(lgO),mrg(lgO)) + call LCMGET(tmpOtherIp,'gig ',oia) + call LCMGET(tmpOtherIp,'mrg ',mrg) + call LCMSIX(tmpOtherIp,nGig,2) + nGig = i2s(i+lgN) + call LCMSIX(tmpResIp,nGig,1) + call LCMPUT(tmpResIp,'gig ',lgO,1,oia) + call LCMPUT(tmpResIp,'mrg ',lgO,1,mrg) + call LCMSIX(tmpResIp,nGig,2) + deallocate(oia,mrg) + end do + else !le repertoire n'existait pas encore => affectation directe + !call LCMEQU(tmpOtherIp,tmpResIp) + call LCMLEN(tmpOtherIp,'t ',lg,ty) + allocate(ora(lg),oia(lg)) + !t + call LCMGET(tmpOtherIp,'t ',oia) + call LCMPUT(tmpResIp,'t ',lg,1,oia) + !x + call LCMGET(tmpOtherIp,'x ',ora) + call LCMPUT(tmpResIp,'x ',lg,2,ora) + !y + call LCMGET(tmpOtherIp,'y ',ora) + call LCMPUT(tmpResIp,'y ',lg,2,ora) + deallocate(ora,oia) + do i = 1,lg + nGig = i2s(i) + call LCMSIX(tmpOtherIp,nGig,1) + call LCMLEN(tmpOtherIp,'gig ',lgO,ty) + allocate(oia(lgO),mrg(lgO)) + call LCMGET(tmpOtherIp,'gig ',oia) + call LCMGET(tmpOtherIp,'mrg ',mrg) + call LCMSIX(tmpOtherIp,nGig,2) + call LCMSIX(tmpResIp,nGig,1) + call LCMPUT(tmpResIp,'gig ',lgO,1,oia) + call LCMPUT(tmpResIp,'mrg ',lgO,1,mrg) + call LCMSIX(tmpResIp,nGig,2) + deallocate(oia,mrg) + end do + end if + end subroutine mergeDic + + subroutine compileResutats(geoIp,lCentreIp,lSideIp,boundDataIp) + type(c_ptr),intent(in) :: geoIp,lCentreIp,lSideIp,boundDataIp + + type(c_ptr) :: geo,centres + integer :: i,lg,ty,lgG + character*12 :: carreName,saveCarreName,gigName,num,mrgName + real,dimension(2) :: sideXY + real,dimension(4) :: minmaxXY + real,dimension(:),allocatable :: xx,yy + integer,dimension(:),allocatable :: tt,gigI,mrgI + + geo = geoIp ; centres = lCentreIp + call LCMSIX(geo,'NEW-DATA ',1) + carreName = ' ' + call LCMNXT(centres,carreName) + saveCarreName = carreName + do + call LCMSIX(geo,carreName,1) + !SIDEXY + call LCMGET(lSideIp,'x'//carreName,sideXY(1)) + call LCMGET(lSideIp,'y'//carreName,sideXY(2)) + if (geomTyp==RecTyp) then + call LCMPUT(geo,'SIDEXY ',2,2,sideXY) + else + call LCMPUT(geo,'SIDEXY ',2,2,sideXY(1)) + end if + !COORDX , COORDY , TURN , POSi + call LCMSIX(centres,carreName,1) + call LCMSIX(centres,'root ',1) + call LCMLEN(centres,'t ',lg,ty) + allocate(xx(lg),yy(lg),tt(lg)) + call LCMGET(centres,'x ',xx) + call LCMPUT(geo,'COORDX ',lg,2,xx) + call LCMGET(centres,'y ',yy) + call LCMPUT(geo,'COORDY ',lg,2,yy) + call LCMGET(centres,'t ',tt) + call LCMPUT(geo,'TURN ',lg,1,tt) + do i = 1,lg + num = i2s(i) + gigName = num + call LCMSIX(centres,gigName,1) + call LCMLEN(centres,'gig ',lgG,ty) + allocate(gigI(lgG),mrgI(lgG)) + call LCMGET(centres,'gig ',gigI) + call LCMGET(centres,'mrg ',mrgI) + call LCMSIX(centres,gigName,2) + gigName = 'POS' // num(:9) + call LCMPUT(geo,gigName,lgG,1,gigI) + mrgName = 'MRG' // num(:9) + call LCMPUT(geo,mrgName,lgG,1,mrgI) + deallocate(gigI,mrgI) + end do + deallocate(xx,yy,tt) + call LCMSIX(centres,'root ',2) + call LCMSIX(centres,carreName,2) + !on cycle + call LCMSIX(geo,carreName,2) + call LCMNXT(centres,carreName) + if (carreName==saveCarreName) exit + end do + !donnes sur les CL + call LCMSIX(geo,'BOUND-DATA ',1) + !SIDEXY + call LCMGET(lSideIp,'xroot ',sideXY(1)) + call LCMGET(lSideIp,'yroot ',sideXY(2)) + call LCMPUT(geo,'SIDEXY ',2,2,sideXY) + !MINMAXXY + call LCMGET(boundDataIp,'minX ',minmaxXY(1)) + call LCMGET(boundDataIp,'minY ',minmaxXY(2)) + call LCMGET(boundDataIp,'maxX ',minmaxXY(3)) + call LCMGET(boundDataIp,'maxY ',minmaxXY(4)) + call LCMPUT(geo,'MINMAXXY ',4,2,minmaxXY) + end subroutine compileResutats + + subroutine creerEtChargerNewDataHex(geoIp,szB,szP) + type(c_ptr),intent(in) :: geoIp + integer,intent(inout) :: szB,szP + + real,dimension(:),allocatable :: txx,tyy + integer,dimension(:),allocatable :: turn,mix,merg + character*12,dimension(:),allocatable :: cells + real :: side + double precision :: sdX,sdY + integer,dimension(40) :: sv + type(c_ptr) :: ip + integer :: nbH,lg,typ,iHex,i,j,k,l + integer :: a,b,aa,bb,nbC,ind,lgSide + integer,dimension(6) :: da,db + real,dimension(4) :: fooData + + ip = geoIp + !ajout de donnees pour le cas d'une seule cellule + call LCMSIX(ip,'NEW-DATA ',1) + call LCMSIX(ip,'BOUND-DATA ',1) + fooData = 0. + call LCMPUT(ip,'MINMAXXY ',4,2,fooData) + call LCMSIX(ip,'BOUND-DATA ',2) + call LCMSIX(ip,'NEW-DATA ',2) + !placement des cellules + call LCMGET(ip,'STATE-VECTOR',sv) + nbH=sv(3) + nbC=sv(9) + if ((sv(8)==0).and.(nbH==1)) then !pas de sous-cellules + call giveSide(ip,szB,side) + fooData(1)=2.0*side + fooData(2)=2.0*sqrt(0.75)*side + call LCMPUT(ip,'SIDEXY ',4,2,fooData) + call chargerNewData(geoIp,szB,szP) + return + end if + allocate(txx(nbH),tyy(nbH),turn(nbH),mix(nbH),merg(nbH)) + call LCMGET(ip,'MIX ',mix) + call LCMLEN(ip,'TURN ',lg,typ) + if (lg/=0) then + call LCMGET(ip,'TURN ',turn) + else + turn = 1 + end if + call LCMLEN(ip,'MERGE ',lg,typ) + if (lg/=0) then + call LCMGET(ip,'MERGE ',merg) + else + merg(:nbH) = (/(i,i=1,nbH)/) + end if + call LCMGET(ip,'IHEX ',iHex) + lgSide = 0 + select case(iHex) + case(H_S30) + a = -1 ; l = 0 ; i = 0 + do + i = i + 1 + do j = 0,1 + a = a + 1 ; b = j + do k = 1,i + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = 2*b ; b = b + 2 + if (l==nbH) goto 10 + end do + end do + end do + case(H_SA60) + l = 0 ; i = -1 + do + i = i + 1 ; a = i ; b = -i + do j = 0,i + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = 2*b ; b = b + 2 + if (l==nbH) goto 10 + end do + end do + case(H_SB60) + a = -1 ; l = 0 ; i = 0 + do + i = i + 1 + do j = 0,1 + a = a + 1 ; b = j + do k = 0,i + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = 2*b ; b = b + 2 + if (l==nbH) goto 10 + end do + aa = a ; bb = b + do k = 0,i-1 + l = l + 1 ; aa = aa - 1 ; bb = bb + 1 + txx(l) = aa ; tyy(l) = bb + if (l==nbH) goto 10 + end do + end do + end do + case(H_S90) + da = (/0,-1,0,0,0,0/) ; db = (/2,1,0,0,0,0/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = i ; b = -i + do j = 1,2 + do k = 1,i + a = a + da(j) ; b = b + db(j) + if (b>=0) then + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = b + if (l==nbH) goto 10 + end if + end do + end do + end do + case(H_R120) + da = (/0,-1,0,0,0,0/) ; db = (/2,1,0,0,0,0/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = i ; b = -i + do j = 1,2 + do k = 1,i + a = a + da(j) ; b = b + db(j) + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = b + if (l==nbH) goto 10 + end do + end do + end do + case(H_R180) + da = (/1,0,-1,0,0,0/) ; db = (/1,2,1,0,0,0/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = 0 ; b = -2*i + do j = 1,3 + do k = 1,i + a = a + da(j) ; b = b + db(j) + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = b + if (l==nbH) goto 10 + end do + end do + end do + case(H_SA180) + da = (/1,0,-1,0,0,0/) ; db = (/1,2,1,0,0,0/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = 0 ; b = -2*i + l = l + 1 ; txx(l) = a ; tyy(l) = b + if (l==nbH) goto 10 + do j = 1,3 + do k = 1,i + a = a + da(j) ; b = b + db(j) + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = b + if (l==nbH) goto 10 + end do + end do + end do + case(H_SB180) + da = (/0,-1,-1,0,0,0/) ; db = (/2,1,-1,-2,0,0/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = i ; b = -i + do j = 1,4 + do k = 1,i + if (b>=0) then + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = -2*a + if (l==nbH) goto 10 + end if + a = a + da(j) ; b = b + db(j) + end do + end do + end do + case(H_Complete) + da = (/-1,-1,0,1,1,0/) ; db = (/1,-1,-2,-1,1,2/) ; l = 1 ; i = 0 + txx(l) = 0 ; tyy(l) = 0 + if (l==nbH) goto 10 + do + i = i + 1 ; a = i ; b = i + do j = 1,6 + do k = 1,i + l = l + 1 ; txx(l) = a ; tyy(l) = b ; lgSide = 2*a + if (l==nbH) goto 10 + a = a + da(j) ; b = b + db(j) + end do + end do + end do + end select +10 continue !sortie de la boucle de remplissage des positions (bhaaa, un goto) + call giveSide(ip,szB,side) + sdX = side*1.5d0 ; sdY = side*0.5d0*sqrt(3.d0) + bCData%sidexy(1) = side ; bCData%sidexy(2) = lgSide * sdY + if (sv(8)==0) then !no sub-geometries + do j = 1,nbH + szP = szP + 1 + tabCellulePlaced(szP)%indice = 1 + tabCellulePlaced(szP)%xcenter = txx(j)*sdX + tabCellulePlaced(szP)%ycenter = tyy(j)*sdY + tabCellulePlaced(szP)%turn = turn(j) + allocate(tabCellulePlaced(szP)%gig(1)) + tabCellulePlaced(szP)%gig(1) = j + allocate(tabCellulePlaced(szP)%mrg(1)) + tabCellulePlaced(szP)%mrg(1) = merg(j) + enddo + else + allocate(cells(nbC)) + call LCMGTC(ip,'CELL ',12,nbC,cells) + do i = 1,szB + ind = 0 + do j = 1,szB + if (cells(i)==tabCelluleBase(j)%name) then + ind = j + exit + end if + end do + if (ind==0) call XABORT("G2S: internal error in function& + & creerEtChargerNewDataHex") + do j = 1,nbH + if (mix(j)==-i) then + !on va creer une cellule placee d'indice ind en position j + szP = szP + 1 + tabCellulePlaced(szP)%indice = ind + tabCellulePlaced(szP)%xcenter = txx(j)*sdX + tabCellulePlaced(szP)%ycenter = tyy(j)*sdY + tabCellulePlaced(szP)%turn = turn(j) + allocate(tabCellulePlaced(szP)%gig(1)) + tabCellulePlaced(szP)%gig(1) = j + allocate(tabCellulePlaced(szP)%mrg(1)) + tabCellulePlaced(szP)%mrg(1) = merg(j) + end if + end do + end do + deallocate(cells) + endif + deallocate(txx,tyy,turn,mix,merg) + end subroutine creerEtChargerNewDataHex + + subroutine giveSide(ip,szB,side) + type(c_ptr),intent(in) :: ip + integer,intent(in) :: szB + real,intent(inout) :: side + + integer :: i,lg,typ + logical :: gotIt + double precision :: dside + + gotIt = .false. + call LCMLEN(ip,'SIDE ',lg,typ) + if (lg/=0) then + call LCMGET(ip,'SIDE ',side) + gotIt = .true. + dside = side + end if + do i = 1,szB + if (tabCelluleBase(i)%ok(n_side)) then + if (gotIt) then + if (.not. isEqual(tabCelluleBase(i)%side,dside)) & + & call XABORT("G2S: Error in the value of argument SIDE& + & in the geometry") + else + gotIt = .true. + dside = tabCelluleBase(i)%side + end if + else if (gotIt) then + tabCelluleBase(i)%side = dside + end if + end do + if (.not. gotIt) then + call XABORT("G2S: Error in the value of argument SIDE in the geometry") + end if + do i = 1,szB + tabCelluleBase(i)%ok(n_side) = .true. + tabCelluleBase(i)%side = dside + end do + side = real(dside) + end subroutine giveSide + + subroutine creerEtChargerNewDataTri(geoIp,szB,szP) + type(c_ptr),intent(in):: geoIp + integer,intent(inout) :: szB,szP + + type(c_ptr) :: ip,lcIp,dicoIp,sidesIp,lCentreIp,centreCalculesIp + integer :: szLc + + ip = geoIp + !creation de la liste des triangles + call LCMOP(lcIp,'lc ',0,1,0) + call creeListeCellules(lcIp,ip,szLc) + + !creation d'une structure pour le calcul des longueurs des cotes + call LCMOP(dicoIp,'dico ',0,1,0) + call creeStructure(dicoIp,lcIp,szLc) + + !creation de la liste des longueurs des cotes et resolution du systeme + !et destruction de la structure + call LCMOP(sidesIp,'sides ',0,1,0) + call resoudStructure(dicoIp,sidesIp) + call LCMCL(dicoIp,2) + + !cree et prepare la liste des centres des cellules de base, + !avec leur gigogne d'origine + call LCMOP(lCentreIp,'centres ',0,1,0) + call prepareCentresTri(lCentreIp,lcIp,szLc,sidesIp) + + !resolution du systeme donnant les coordonnees des centres des cellules + !de base et stockage dans une nouvelle structure + call LCMOP(centreCalculesIp,'resCentres ',0,1,0) + call calculeCentre(centreCalculesIp,lCentreIp,0) + + !compilation des resultats et ajout dans la geometrie des nouvelles donnees + call compileResutatsTri(ip,lcIp,szLc,centreCalculesIp,sidesIp) + call LCMCL(lcIp,2) + call LCMCL(centreCalculesIp,2) + call LCMCL(sidesIp,2) + !utilisation des donnees (on accroche le traitement classique developpe + !pour python) + call chargerNewData(geoIp,szB,szP) + end subroutine creerEtChargerNewDataTri + + subroutine creeStructure(dicoIp,lcIp,szLc) + type(c_ptr),intent(in) :: dicoIp,lcIp + integer,intent(in) :: szLc + + character*12 :: triName + type(c_ptr) :: triIp,dIp + integer :: i,j,k,lg,ty,lenMix,lignNbr + real :: side,value + integer,dimension(40) :: sv + integer,dimension(:),allocatable :: mix + character*12,dimension(:),allocatable :: cell + + !initialisation + dIp = dicoIp + triName = ' ' + do i = 1,szLc + call LCMNXT(lcIp,triName) + triIp=LCMGID(lcIp,triName) + side = -1. + call LCMLEN(triIp,'SIDE ',lg,ty) + if (lg/=0) call LCMGET(triIp,'SIDE ',side) + call LCMSIX(dIp,triName,1) + call LCMPUT(dIp,'value',1,2,side) + call LCMSIX(dIp,triName,2) + end do + !remplissage + dIp = dicoIp + triName = ' ' + lignNbr = 0 + do k = 1,szLc + !tavail sur la kieme cellule + call LCMNXT(lcIp,triName) + triIp=LCMGID(lcIp,triName) + call LCMGET(triIp,'STATE-VECTOR',sv) + call LCMLEN(triIp,'MIX ',lenMix,ty) + if ((lenMix==1).or.(sv(9)==0)) cycle !pas d'info a retirer + allocate(mix(lenMix),cell(sv(9))) + call LCMGET(triIp,'MIX ',mix) + !on teste si tous les milieux sont bien negatifs + do i = 1,lenMix + if (mix(i)<0) cycle + call XABORT("G2S: error, meltig MIX and CELL not supported") + end do + call LCMGTC(triIp,'CELL ',12,sv(9),cell) + !on traite les donnees + !dans le triangle considere + call LCMSIX(dIp,triName,1) + value = 1.*sv(3) + do i = 1,lenMix + call LCMPUT(dIp,cell(-mix(i)),1,2,value) + end do + call LCMSIX(dIp,triName,2) + !dans les autres + value = 1./sv(3) + do i = 1,sv(9) + call LCMSIX(dIp,cell(i),1) + do j = 1,sv(9) + if (j==i) cycle + call LCMPUT(dIp,cell(j),1,2,1.) + end do + call LCMPUT(dIp,triName,1,2,value) + call LCMSIX(dIp,cell(i),2) + end do + deallocate(mix,cell) + end do + end subroutine creeStructure + + subroutine resoudStructure(dicoIp,resIp) + type(c_ptr),intent(in) :: dicoIp,resIp + + type(c_ptr) :: dIp + logical :: fini,flag1,flag2 + integer :: long,typ + character*12 :: triName,saveTriName + character*131:: hsmg + real :: val + + fini = .false. + flag1 = .true. + do + if (fini) exit + fini = .true. + flag2 = flag1 !garde la valeur de flag1 a l'essai precedent + flag1 = .false. + triName = ' ' + call LCMNXT(dicoIp,triName) + saveTriName = triName + do !boucle sur toutes les donnees du dico + fini = fini .and. evaluateLine(dicoIp,triName) + !test si fini a ete vrai au moins une fois + if (fini) flag1 = .true. + call LCMNXT(dicoIp,triName) + if (triName == saveTriName) exit + end do + if (.not.(flag1 .or. flag2)) then + !il ne s'est rien passe au cours des 2 tentatives + !precedentes => on abandonne + call XABORT("G2S: not enought data in the geometry(2)") + end if + end do + !recopie des resultats + dIp = dicoIp + triName = ' ' + call LCMNXT(dIp,triName) + saveTriName = triName + do + call LCMSIX(dIp,triName,1) + call LCMLEN(dIp,'value',long,typ) + if(long == 0) then + write(hsmg,'(36hG2S: value key missing for triName=,a,5h (3).)') trim(triName) + call LCMLIB(dIp) + call LCMLIB(dicoIp) + call XABORT(hsmg) + endif + call LCMGET(dIp,'value',val) + call LCMPUT(resIp,triName,1,2,val) + call LCMSIX(dIp,triName,2) + call LCMNXT(dIp,triName) + if (triName == saveTriName) exit + end do + end subroutine resoudStructure + + function evaluateLine(dicoIp,triName) + type(c_ptr),intent(in) :: dicoIp + character*12,intent(in) :: triName + logical :: evaluateLine + + type(c_ptr) :: triIp + integer :: long,typ + character*12 :: eqName,saveEqName + character*131:: hsmg + real :: factor,val,res + + evaluateLine = .true. + triIp = dicoIp + call LCMSIX(triIp,triName,1) + eqName = ' ' + call LCMNXT(triIp,eqName) + saveEqName = eqName + do !boucle sur toutes les equations de triName + if (eqName=='value') then + call LCMNXT(triIp,eqName) + if (eqName==saveEqName) exit + end if + call getValueOf(eqName,dicoIp,val) + if (val<0.) then + evaluateLine = .false. + else + call LCMGET(triIp,eqName,factor) + res = val * factor + call LCMLEN(triIp,'value',long,typ) + if(long == 0) then + write(hsmg,'(36hG2S: value key missing for eqName=,a,5h (4).)') trim(eqName) + call LCMLIB(triIp) + call LCMLIB(dicoIp) + call XABORT(hsmg) + endif + call LCMGET(triIp,'value',val) + if (val<0.) then + call LCMPUT(triIp,'value',1,2,res) + else if (abs(val-res)>geomPrec) then + call XABORT("G2S: error, incoherent geometry(3)") + end if + end if + call LCMNXT(triIp,eqName) + if (eqName==saveEqName) exit + end do + end function evaluateLine + + subroutine prepareCentresTri(lCentreIp,lcIp,szLc,sidesIp) + type(c_ptr),intent(in) :: lCentreIp,lcIp,sidesIp + integer,intent(in) :: szLc + + character*12 :: triName + type(c_ptr) :: triIp + integer :: i,j,k,lenMix,ind,iTri,lg,ty + real :: dsidex,dsidey,xi,yj + integer,dimension(40) :: sv + real,dimension(:),allocatable :: xx,yy + integer,dimension(:),allocatable :: mix,turn,merg + character*12,dimension(:),allocatable :: cell + + triName = ' ' + do k = 1,szLc + !tavail sur la kieme cellule + call LCMNXT(lcIp,triName) + triIp=LCMGID(lcIp,triName) + call LCMGET(triIp,'STATE-VECTOR',sv) + call LCMLEN(triIp,'MIX ',lenMix,ty) + if ((lenMix==1).or.(sv(9)==0)) cycle !pas d'info a retirer + allocate(mix(lenMix),turn(lenMix),merg(lenMix), & + xx(lenMix),yy(lenMix),cell(sv(9))) + call LCMGET(triIp,'MIX ',mix) + call LCMGET(triIp,'TURN ',turn) + call LCMLEN(triIp,'MERGE ',lg,ty) + if (lg/=0) then + call LCMGET(triIp,'MERGE ',merg) + else + merg = (/(i,i=1,lenMix)/) + end if + call LCMGTC(triIp,'CELL ',12,sv(9),cell) + call LCMGET(sidesIp,cell(1),dsidex) + dsidex = 0.5*dsidex + dsidey = sqrt_3f*dsidex + iTri = T_ST60 + if (triName=='root') call LCMGET(triIp,'ITRI ',iTri) + select case(iTri) + case(T_ST60) !triangle normal ou 'root en ST60' + ind = 0 + yj = (1-sv(4))*0.5*dsidey + do j = 1,sv(4) + xi = (j-sv(3))*dsidex + do i = 1,2*(sv(3)-j)+1 + ind = ind + 1 + xx(ind) = xi + yy(ind) = yj + xi = xi + dsidex + end do + yj = yj + dsidey + end do + case(T_S30) + ind = 0 + yj = 0.5*dsidey + do j = 1,sv(4) + xi = (3*j-2)*dsidex + do i = 1,sv(3)-3*(j-1) + ind = ind + 1 + xx(ind) = xi + yy(ind) = yj + xi = xi + dsidex + end do + yj = yj + dsidey + end do + case(T_SA60) + ind = 0 + yj = -(2*sv(4)-1)*0.5*dsidey + do j = 1,sv(4) + xi = (3*(sv(4)-j)+1)*dsidex + do i = 1,3*j-2+mod(sv(3)-1,3) + ind = ind + 1 + xx(ind) = xi + yy(ind) = yj + xi = xi + dsidex + end do + yj = yj + dsidey + end do + yj = 0.5*dsidey + do j = 1,sv(4) + xi = (3*j-2)*dsidex + do i = 1,sv(3)-3*(j-1) + ind = ind + 1 + xx(ind) = xi + yy(ind) = yj + xi = xi + dsidex + end do + yj = yj + dsidey + end do + case(T_Complete) + ind = 0 + yj = (1-sv(4))*0.5*dsidey + do j = 1,sv(4) + xi = (j-sv(3))*dsidex + do i = 1,2*sv(3) + ind = ind + 1 + xx(ind) = xi + yy(ind) = yj + xi = xi + dsidex + end do + yj = yj + dsidey + end do + end select + do i = 1,lenMix + call addInRes(lCentreIp,cell(-mix(i)),triName,xx(i), & + yy(i),turn(i),i,merg(i)) + end do + deallocate(mix,turn,merg,xx,yy,cell) + end do + end subroutine prepareCentresTri + + subroutine compileResutatsTri(geoIp,lcIp,szLc,centreCalculesIp,sidesIp) + type(c_ptr),intent(in) :: geoIp,lcIp,centreCalculesIp,sidesIp + integer,intent(in) :: szLc + + type(c_ptr) :: lMinMaxIp,sidexyIp + integer :: i,iTri + character*12 :: cellName + real :: val + integer,dimension(40) :: sv + + call LCMOP(lMinMaxIp,'minAndMaxXY ',0,1,0) + call LCMPUT(lMinMaxIp,'minX ',1,2,0.) + call LCMPUT(lMinMaxIp,'maxX ',1,2,0.) + call LCMPUT(lMinMaxIp,'minY ',1,2,0.) + call LCMPUT(lMinMaxIp,'maxY ',1,2,0.) + + call LCMOP(sidexyIp,'sidesXY ',0,1,0) + cellName = ' ' + do i = 1,szLc + call LCMNXT(lcIp,cellName) + call LCMGET(sidesIp,cellName,val) + call LCMPUT(sidexyIp,'x'//cellName(1:11),1,2,val) + call LCMPUT(sidexyIp,'y'//cellName(1:11),1,2,val) + if (cellName=='root') then !modification de xroot et yroot si besoin + call LCMGET(geoIp,'STATE-VECTOR',sv) + call LCMGET(geoIp,'ITRI ',iTri) + select case(iTri) + case(T_Complete) + val=val*sv(4)/real(sv(3)) + call LCMPUT(sidexyIp,'y'//cellName(1:11),1,2,val) + case(T_S30,T_SA60) + val=val/real(sv(3)) + call LCMPUT(sidexyIp,'x'//cellName(1:11),1,2,val) + val=val*int((sv(3)+1)/2) + call LCMPUT(sidexyIp,'y'//cellName(1:11),1,2,val) + end select + end if + end do + + call compileResutats(geoIp,centreCalculesIp,sidexyIp,lMinMaxIp) + + call LCMCL(lMinMaxIp,2) + call LCMCL(sidexyIp,2) + end subroutine compileResutatsTri + +end module pretraitement |