diff options
Diffstat (limited to 'Dragon/src/g2s_convert.f90')
| -rw-r--r-- | Dragon/src/g2s_convert.f90 | 530 |
1 files changed, 530 insertions, 0 deletions
diff --git a/Dragon/src/g2s_convert.f90 b/Dragon/src/g2s_convert.f90 new file mode 100644 index 0000000..c8f1fe9 --- /dev/null +++ b/Dragon/src/g2s_convert.f90 @@ -0,0 +1,530 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! Convert an Alamos surfacic file towards Salomon format. +! +!Copyright: +! Copyright (C) 2022 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): +! A. Hebert +! +!Parameters: input +! ipAl Alamos ascii file index +! +!Parameters: output +! ipSal Salomon ascii file index +! +!----------------------------------------------------------------------- +! +subroutine g2s_convert(impx,ipAl,ipZa,ipSal) + use constUtiles + use segArc + ! typgeo type of geometry (0: TISO tracking; 5: rectangle TRAN; + ! 6: rectangle REFL; 7: eight;8: SA60; 9: hexagon; 10: RA60; + ! 11: R120) + integer,intent(in) :: impx,ipAl,ipZa,ipSal + + ! local variables + integer :: i, j, k, nbfold, defautCl, nangles, saltype, idummy + double precision :: aStart, bStart, rStart, anglStart, aEnd, bEnd, rEnd, anglEnd, & + & radius, angle, angle0, dx, dy, dx0, dy0, r_xmin, r_xmax, & + & r_ymin, r_ymax + character(len=4) :: strDCL + character(len=12) :: text12, name_geom + character(len=40) :: text40 + character(len=131) :: hsmg + integer, allocatable, dimension(:) :: nber,merg,milTab,milTab2,iper + integer, allocatable, dimension(:,:) :: bcElems,iboundary + double precision, allocatable, dimension(:) :: ccx,ccy,ttx,tty + double precision, allocatable, dimension(:,:) :: apnodes + character(len=24), allocatable, dimension(:) :: propertyNames,mixNames,mixNames2 + type(t_segArc) :: sa + double precision, allocatable, dimension(:) :: list_angles + ! + ! set angles + nangles=-99 + if(typgeo == 0) then + nangles=0 + else if((abs(typgeo) == 5).or.(abs(typgeo) == 6).or.(abs(typgeo) == 11)) then + nangles=4 + else if((abs(typgeo) == 7).or.(abs(typgeo) == 8).or.(abs(typgeo) == 10)) then + nangles=3 + else if(abs(typgeo) == 9) then + nangles=6 + else + write(hsmg,'(35hg2s_convert: invalid symmetry type=,i3,1h.)') typgeo + call XABORT(hsmg) + endif + allocate(list_angles(nangles)) + if((abs(typgeo) == 5).or.(abs(typgeo) == 6)) then + list_angles=(/ 0.0, 0.0, 90.0, 90.0 /) + else if(abs(typgeo) == 7) then + list_angles=(/ 0.0, 45.0, 90.0 /) + else if((abs(typgeo) == 8).or.(abs(typgeo) == 10)) then + list_angles=(/ 0.0, 60.0, 120.0 /) + else if(abs(typgeo) == 9) then + list_angles=(/ 0.0, 120.0, 60.0, 0.0, 120.0, 60.0 /) + else if(abs(typgeo) == 11) then + list_angles=(/ 0.0, 0.0, 60.0, 60.0 /) + endif + + ! ---------------------------------- + ! Read the Alamos surfacic file + ! ---------------------------------- + read(ipAl,'(a12,a24)',end=100) text12,name_geom + 10 read(ipAl,'(a12,a24)',end=100) text12,name_geom + if (text12 /= ' Geometrie:') go to 10 + 20 read(ipAl,'(a12)',end=100) text12 + if (text12 /= ' Milieux:') go to 20 + read(ipAl,"(i7)",advance='no') nbmil + allocate(mixNames(nbmil)) + read(ipAl,*) mixNames(:) + 30 read(ipAl,'(a12)',end=100) text12 + if (text12 /= ' Mailles:') go to 30 + read(ipAl,"(i7)") nbNode + allocate(milTab(nbNode),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT('g2s_convert: allocation problem.') + read(ipAl,'(10i7)',end=100) milTab(:) + nbFlux=nbNode + if(ipZa /= -1) then + ! ---------------------------------- + ! Read the PropertyMap file + ! ---------------------------------- + allocate(mixNames2(nbNode),milTab2(nbNode),propertyNames(nbNode),iper(nbNode),stat=alloc_ok) + if (alloc_ok /= 0) call XABORT('g2s_convert: allocation problem.') + 40 read(ipZa,'(a12,a24)',end=110) text12 + if (text12 /= 'PropertyMap:') go to 40 + read(ipZA,*) propertyNames(:) + nbmil=0 + loop1: do i=1,nbNode + do j=1,nbmil + if(propertyNames(i) == mixNames2(j)) then + milTab2(i)=j; cycle loop1 + endif + enddo + nbmil=nbmil+1 + mixNames2(nbmil)=propertyNames(i) + milTab2(i)=nbmil + iper(nbmil)=milTab(i) + enddo loop1 + write(6,'(1x,a)') name_geom + do i=1,nbmil + write(6,'(1x,i5,2x,a,4h--> ,a)') i,mixNames2(i),mixNames(iper(i)) + enddo + deallocate(mixNames); allocate(mixNames(nbmil)) + milTab(:nbNode)=milTab2(:nbNode); mixNames(:nbmil)=mixNames2(:nbmil) + deallocate(iper,propertyNames,milTab2,mixNames2) + else + write(6,'(1x,a)') name_geom + do i=1,nbmil + write(6,'(1x,i5,2x,a)') i,mixNames(i) + enddo + endif + if(impx > 1) then + write(6,'(5h--cut,75(1h-))') + i=1 + do j=1,1+(nbmil-1)/4 + write(text40,'(15h(10h INTEGER ,,i1,12h(a,1x),2h:=,,i1,8hi5,2h ;))') & + & min(i+3,nbmil)-i+1,min(i+3,nbmil)-i+1 + write(6,text40) (trim(mixNames(k)),k=i,min(i+3,nbmil)),(k,k=i,min(i+3,nbmil)) + i=i+4 + enddo + write(6,'(5h--cut,75(1h-))') + endif + deallocate(mixNames) + + ! read surfacic elements + rewind(ipAl) + 50 read(ipAl,'(a12)',end=100) text12 + if (text12 /= ' Noeuds:') go to 50 + read(ipAl,'(i7)') nbPoints + allocate(apnodes(2,nbPoints)) + r_xmin=999.0 + r_ymin=999.0 + do i=1,nbPoints + if((abs(typgeo) == 8).or.(abs(typgeo) == 10).or.(abs(typgeo) == 11)) then + read(ipAl,*,end=100) apnodes(2,i),apnodes(1,i) + apnodes(2,i)=-apnodes(2,i) ! rotate Alamos geometry by 90 degres + else + read(ipAl,*,end=100) apnodes(1,i),apnodes(2,i) + endif + r_xmin=min(r_xmin,apnodes(1,i)) + r_ymin=min(r_ymin,apnodes(2,i)) + enddo + do i=1,nbPoints + apnodes(1,i)=apnodes(1,i)-r_xmin + apnodes(2,i)=apnodes(2,i)-r_ymin + enddo + r_xmin=999.0 + r_xmax=-999.0 + r_ymin=999.0 + r_ymax=-999.0 + do i=1,nbPoints + r_xmin=min(r_xmin,apnodes(1,i)) + r_ymin=min(r_ymin,apnodes(2,i)) + r_xmax=max(r_xmax,apnodes(1,i)) + r_ymax=max(r_ymax,apnodes(2,i)) + enddo + basex=real(r_xmax-r_xmin) + basey=real(r_ymax-r_ymin) + + ! compute boundary elements characteristics + if(abs(typgeo) > 0) then + allocate(ttx(nangles),tty(nangles),ccx(nangles),ccy(nangles)) + if(abs(typgeo) == 9) basex=basex/2.0 + if(abs(typgeo) == 11) basex=2.0*basex/3.0 + do i = 1,nangles + ccx(i)=0.0 ; ccy(i)=0.0 ; ttx(i)=0.0 ; tty(i)=0.0 + if((abs(typgeo) == 5).or.(abs(typgeo) == 6)) then + select case(i) + case(2) + ccy(i)=basey + case(4) + ccx(i)=basex + end select + else if((abs(typgeo) == 7).or.(abs(typgeo) == 8).or.(abs(typgeo) == 10)) then + select case(i) + case(3) + ccx(i)=basex + end select + else if(abs(typgeo) == 9) then + select case(i) + case(1:2) + ccx(i)=basex/2.0 + case(3) + ccy(i)=basex*sqrt(0.75) + case(4) + ccx(i)=basex/2.0 + ccy(i)=2.0*basex*sqrt(0.75) + case(5) + ccx(i)=2.0*basex + ccy(i)=basex*sqrt(0.75) + case(6) + ccx(i)=1.5*basex + end select + else if(abs(typgeo) == 11) then + select case(i) + case(2) + ccx(i)=basex/2.0 + ccy(i)=basex*sqrt(0.75) + case(4) + ccx(i)=basex + end select + endif + if(typgeo == 5) then + ! pure translation in Cartesian geometry + select case(i) + case(1) + tty(i)=basey + case(2) + tty(i)=-basey + case(3) + ttx(i)=basex + case(4) + ttx(i)=-basex + end select + else if(typgeo == 9) then + ! pure translation in hexagonal geometry + select case(i) + case(1) + tty(i)=2.0*basex*sqrt(0.75) + case(2) + ttx(i)=2.0*basex + tty(i)=basex*sqrt(0.75) + case(3) + ttx(i)=2.0*basex + tty(i)=-basex*sqrt(0.75) + case(4) + tty(i)=-2.0*basex*sqrt(0.75) + case(5) + ttx(i)=-2.0*basex + tty(i)=-basex*sqrt(0.75) + case(6) + ttx(i)=-2.0*basex + tty(i)=basex*sqrt(0.75) + end select + else + ttx(i)=ccx(i) ; tty(i)=ccy(i) + endif + end do + endif + + ! read surfacic elements + read(ipAl,'(a12)') text12 + if (text12 /= ' Aretes:') call XABORT('g2s_convert: keyword Aretes: expected.') + read(ipAl,'(i7)',end=100) iszSA + allocate(tabSegArc(iszSA),iboundary(iszSA,2)) + ibd=0 + do i=1,iszSA + read(ipAl,*,end=100) itype + backspace(ipAl) + if(itype == 0) then + ! line segment + tabSegArc(i)%typ=1 + read(ipAl,*) itype,ipt1,ipt2,nmoins,nplus + tabSegArc(i)%x=real(apnodes(1,ipt1)) + tabSegArc(i)%y=real(apnodes(2,ipt1)) + tabSegArc(i)%dx=real(apnodes(1,ipt2)) + tabSegArc(i)%dy=real(apnodes(2,ipt2)) + else if(itype == 1) then + ! full circle + tabSegArc(i)%typ=2 + read(ipAl,*) itype,ipt3,radius,nplus,nmoins + tabSegArc(i)%x=real(apnodes(1,ipt3)) + tabSegArc(i)%y=real(apnodes(2,ipt3)) + tabSegArc(i)%r=real(radius) + tabSegArc(i)%b=0.0 ; tabSegArc(i)%a=0.0 + else if(itype == 2) then + ! circular arc + tabSegArc(i)%typ=3 + read(ipAl,*) itype,ipt1,ipt2,ipt3,nplus,nmoins + aStart=real(apnodes(1,ipt1)-apnodes(1,ipt3)) + bStart=real(apnodes(2,ipt1)-apnodes(2,ipt3)) + rStart=sqrt(aStart**2+bStart**2) + if(bStart >= 0.0) then + anglStart=acos(aStart/rStart) + else + anglStart=dpi_c-acos(aStart/rStart) + endif + aEnd=real(apnodes(1,ipt2)-apnodes(1,ipt3)) + bEnd=real(apnodes(2,ipt2)-apnodes(2,ipt3)) + rEnd=sqrt(aEnd**2+bEnd**2) + if(bEnd >= 0.0) then + anglEnd=acos(aEnd/rEnd) + else + anglEnd=dpi_c-acos(aEnd/rEnd) + endif + tabSegArc(i)%x=real(apnodes(1,ipt3)) + tabSegArc(i)%y=real(apnodes(2,ipt3)) + tabSegArc(i)%r=real(0.5d0*(rStart+rEnd)) + tabSegArc(i)%b=real(anglEnd) + tabSegArc(i)%a=real(anglStart) + else + write(hsmg,'(35hg2s_convert: invalid element type=,i3,1h.)') itype + call XABORT(hsmg) + endif + tabSegArc(i)%nodeg=nplus + tabSegArc(i)%noded=nmoins + if(abs(typgeo) > 0) then + if(((nplus==0).or.(nmoins==0)).and.(itype == 0)) then + ibd=ibd+1 + if(ibd > iszSA) call XABORT('g2s_convert: boundary overflow') + iboundary(ibd,2)=0 + dx=tabSegArc(i)%x-tabSegArc(i)%dx + dy=tabSegArc(i)%y-tabSegArc(i)%dy + angle=atan(dy/dx)*rad2deg + if(angle < 0.0) angle=angle+180. + do j=1,nangles + if(abs(angle-list_angles(j)) > 0.5) cycle + if((abs(tabSegArc(i)%x-ccx(j)) < epsilon).and.(abs(tabSegArc(i)%y-ccy(j)) < epsilon)) then + iboundary(ibd,2)=j + else if((abs(tabSegArc(i)%dx-ccx(j)) < epsilon).and.(abs(tabSegArc(i)%dy-ccy(j)) < epsilon)) then + iboundary(ibd,2)=j + else + dx0=tabSegArc(i)%x-ccx(j) + dy0=tabSegArc(i)%y-ccy(j) + angle0=atan(dy0/dx0)*rad2deg + if(angle0 < 0.0) angle0=angle0+180. + if(abs(angle-angle0) <= 1.0) iboundary(ibd,2)=j ! look for 1 degree agreement + endif + enddo + iboundary(ibd,1)=i + endif + endif + enddo + if(abs(typgeo) > 0) deallocate(ccy,ccx) + + ! list boundary elements + allocate(bcElems(iszSA,nangles),nber(nangles)) + if(abs(typgeo) > 0) then + nber(:)=0 + bcElems(:,:)=-1 + do ibound=1,ibd + i1=iboundary(ibound,1) + i2=iboundary(ibound,2) + if(i2 == 0) then + write(hsmg,'(26hg2s_convert: boundary side,i7,18h is not allocated., & + & 33h Try a different value of typgeo.)') i1 + call XABORT(hsmg) + endif + if(i2 > nangles) call XABORT('g2s_convert: boundary angle overflow.') + nber(i2)=nber(i2)+1 + if(nber(i2) > iszSA) call XABORT('g2s_convert: boundary element overflow.') + bcElems(nber(i2),i2)=i1 + enddo + endif + deallocate(iboundary) + + ! read Bords data in some specific Alamos files + read(ipAl,'(a12)') text12 + if (text12 == ' Bords:') then + read(ipAl,*) nbBords + do ib=1,nbBords + read(ipAl,*) itype,ib1,ib2 + enddo + read(ipAl,'(a12)') text12 + endif + + ! skip milTab data + if (text12 /= ' Mailles:') call XABORT('g2s_convert: keyword Mailles: expected.') + read(ipAl,'(i7)') nbNode + read(ipAl,'(10i7)',end=100) (idummy, i=1,nbNode) + read(ipAl,'(a12)') text12 + if (text12 == ' ') read(ipAl,'(a12)') text12 + print *,'read=',text12 + if (text12 /= ' Fin:') call XABORT('g2s_convert: keyword Fin: expected.') + ! + ! set nbfold + nbfold=0 + if(typgeo == -5) then + typgeo=1 ; nbfold=4 + else if(typgeo == -6) then + typgeo=1 ; nbfold=4 + else if(typgeo == -7) then + typgeo=1 ; nbfold=8 + else if(typgeo == -8) then + typgeo=1 ; nbfold=6 + else if(typgeo == -9) then + typgeo=1 ; nbfold=1 + else if(typgeo == -10) then + typgeo=2 ; nbfold=6 + else if(typgeo == -11) then + typgeo=2 ; nbfold=3 + endif + ! ---------------------------------- + ! Generate the Salomon surfacic file + ! ---------------------------------- + write(ipSal,'(a5)') 'BEGIN' + + write(ipSal,'(/a24)') 'DEFINE DOMAINE' + write(ipSal,'(a24/)') '==============' + + write(ipSal,'(a28/)') '1.main dimensions:' + write(ipSal,'(a32)') '*typgeo nbfold nbNode nbelem ' + write(ipSal,'(10'//formati//')') typgeo,nbfold,nbNode,iszSA,1,nbFlux + + write(ipSal,'(/a37)') '2.impression and precision:' + write(ipSal,'(/a21)') '*index kndex prec' + write(ipSal,'(10'//formati//')') 0,0,1 + + write(ipSal,'(/a40)') '3.precision of geometry data:' + write(ipSal,'(/a4)') '*eps' + write(ipSal,'('//formatr//')') gSALeps + + write(ipSal,'(/a58)') '4.flux region number per geometry region (mesh):' + write(ipSal,'(/a6)') '*merge' + allocate(merg(nbNode)) + do i=1,nbNode + merg(i)=i + enddo + write(ipSal,'(10'//formati//')') merg(:nbNode) + deallocate(merg) + + write(ipSal,'(/a29)') '5.name of geometry:' + write(ipSal,'(/a11)') '*macro_name' + write(ipSal,'(4(3x,a12))') name_geom + + write(ipSal,'(/a46)') '6.macro order number per flux region:' + write(ipSal,'(/a14)') '*macro_indices' + write(ipSal,'(10'//formati//')') (1,i=1,nbFlux) + + write(ipSal,'(/a57)') '7.read integer and real data for each elements:' + do i=1,iszSA + sa = tabSegArc(i) + cx=real(sa%x) + cy=real(sa%y) + if(sa%typ == 1) then + ! line segment + tx= real(sa%dx-sa%x) + ty= real(sa%dy-sa%y) + delta=0. + else if(sa%typ == 2) then + ! full circle + tx=real(sa%r) + ty=0. + delta=0. + else if(sa%typ == 3) then + ! circular arc + tx=real(sa%r) + ty=real(sa%a*rad2deg) + if (sa%b>sa%a) then + delta=real((sa%b-sa%a)*rad2deg) + else + delta=real((sa%b-sa%a)*rad2deg+360.d0) + endif + endif + write(ipSal,*) + write(ipSal,*) 'elem =',i + write(ipSal,'(a22)') '*type node- node+' + write(ipSal,'(10'//formati//')') sa%typ,sa%nodeg,sa%noded + write(ipSal,'(a63)') '*cx cy ex or R & + &ey or theta1 theta2' + write(ipSal,'(5'//formatr//')') cx,cy,tx,ty,delta + enddo + deallocate(tabSegArc) + + ! write boundary elements lists + write(ipSal,'(/a63)') '8.read integer and real data for boundary conditions:' + defautCl = 0 + albedo = 0.0 + strDCL = 'VOID' + write(ipSal,'(/a40)') '*defaul nbbcda allsur divsur ndivsur' + write(ipSal,'(10'//formati//')') defautCl,nangles,0,0,0 + write(ipSal,'(/a24)') 'DEFAULT = ' // strDCL + write(ipSal,'(a24)') '==============' + write(ipSal,'(/a17)') '*albedo deltasur' + write(ipSal,'(5'//formatr//')') albedo,0.0 + if(typgeo /= 0) then + do i = 1,nangles + write(ipSal,*) + write(ipSal,*) 'particular boundary condition number',i + write(ipSal,'(/a13)') '*type nber' + if(typgeo <= 2) then + saltype = 1 ! isotropic reflexion + else if((typgeo == 5).or.(typgeo == 9)) then + saltype = 2 ! translation + else if(typgeo < 9) then + saltype = 4 ! axial symmetry (specular reflexion) + else if(typgeo == 10) then ! RA60 + if(i == 1) then + saltype = 2 ! translation + else + saltype = 3 ! rotation + endif + else if(typgeo == 11) then ! R120 + if((i == 1).or.(i == 4)) then + saltype = 2 ! translation + else + saltype = 3 ! rotation + endif + else + call XABORT('g2s_convert: unknown type of geometry.') + endif + write(ipSal,'(10'//formati//')') saltype, nber(i) + write(ipSal,'(/a14)') '*elems(1,nber)' + write(ipSal,'(10'//formati//')') bcElems(:nber(i),i) + write(ipSal,'(/a22)') '*cx cy angle' + write(ipSal,'(5'//formatr//')') ttx(i),tty(i), list_angles(i) + end do + deallocate(tty,ttx) + endif + deallocate(nber,bcElems,list_angles) + + write(ipSal,'(/a28)') '9.medium per node:' + write(ipSal,'(/a11)') '*mil(nbreg)' + write(ipSal,'(10'//formati//')') milTab(:nbNode) + deallocate(milTab) + + write(ipSal,'(/a3)') 'END' + rewind(ipSal) + return + ! + 100 call XABORT('g2s_convert: end of Alamos surfacic file encountered.') + 110 call XABORT('g2s_convert: end of PropertyMap file encountered.') + end subroutine g2s_convert |