diff options
Diffstat (limited to 'Dragon/src/SAL_GEOMETRY_MOD.f90')
| -rw-r--r-- | Dragon/src/SAL_GEOMETRY_MOD.f90 | 3756 |
1 files changed, 3756 insertions, 0 deletions
diff --git a/Dragon/src/SAL_GEOMETRY_MOD.f90 b/Dragon/src/SAL_GEOMETRY_MOD.f90 new file mode 100644 index 0000000..373c95d --- /dev/null +++ b/Dragon/src/SAL_GEOMETRY_MOD.f90 @@ -0,0 +1,3756 @@ +! +!--------------------------------------------------------------------- +! +!Purpose: +! Support module used to create a 2D geometry. +! +!Copyright: +! Copyright (C) 2001 Ecole Polytechnique de Montreal. +! +!Author(s): +! X. Warin +! +!--------------------------------------------------------------------- +! +MODULE SAL_GEOMETRY_MOD + + USE SAL_GEOMETRY_TYPES + USE PRECISION_AND_KINDS, ONLY : PDB, PI,TWOPI,HALFPI + USE SAL_NUMERIC_MOD, ONLY : SAL141 + USE SALGET_FUNS_MOD + +CONTAINS + + SUBROUTINE SAL100(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! perform input data and first allocation for geometry OBJECT T_G_BASIC + ! + !Parameters: input/output + ! GG geometry basic information. + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,TYPGEO,NBFOLD,EPS,ANGGEO, & + INDEX,KNDEX,PREC,ISPEC,LGSPEC + USE SAL_TRACKING_TYPES, ONLY : PRTIND + IMPLICIT NONE + TYPE(T_G_BASIC) :: GG + !***** + ! local variable + ! ************** + INTEGER, PARAMETER :: N_DATAIN=25, N_DATARE=20 + INTEGER, DIMENSION (N_DATAIN) :: DATAIN + REAL, DIMENSION (N_DATARE) :: DATARE + INTEGER :: OK + INTEGER, PARAMETER :: FOUT =6 + !***** + ! TYPGEO = type of geometry: + ! NBFOLD = n in angle definition of rotation or symmetry geometry + ! NB_NODE = number of nodes + ! NBELEM = number of elements + ! + CALL SALGET(DATAIN,6,F_GEO,FOUT0,'dimensions for geometry') + TYPGEO=DATAIN(1) + NBFOLD=DATAIN(2) + GG%NB_NODE=DATAIN(3) + GG%NB_ELEM=DATAIN(4) + GG%NB_MACRO=DATAIN(5) + GG%NB_FLUX=DATAIN(6) + ANGGEO=0.0D0 + SELECT CASE(TYPGEO) + CASE(1) + IF(NBFOLD.EQ.3) THEN + ANGGEO=TWOPI/8.0D0 ! half diagonal symmetry case + ELSE + ANGGEO=TWOPI/NBFOLD + ENDIF + CASE(2) + ANGGEO=TWOPI/NBFOLD + CASE(5:6) + ANGGEO=HALFPI + CASE(7) + ANGGEO=PI*0.25 + CASE(8:11) + ANGGEO=PI/3. + CASE(12) + ANGGEO=PI/6. + END SELECT + IF((PRTIND >= 1).AND.(TYPGEO.EQ.0)) THEN + WRITE(FOUT,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD + ELSE IF(PRTIND >= 1) THEN + WRITE(FOUT,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD, & + & ' ANGGEO=',ANGGEO,' radians' + ENDIF + LGSPEC=(TYPGEO/=0).AND.(NBFOLD==0) + IF(LGSPEC) THEN + IF(ISPEC==0) THEN + WRITE(*,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD + CALL XABORT('SAL100: TISO option is incompatible with the surfacic file') + ENDIF + ELSE + IF(ISPEC==1) THEN + WRITE(*,*) 'SAL100: TYPGEO=',TYPGEO,' NBFOLD=',NBFOLD + CALL XABORT('SAL100: TSPC option is incompatible with the surfacic file') + ENDIF + ENDIF + ! + !* read printing indexes for general domain data and topological deformations + ! INDEX = to print general domain data + ! KNDEX = to print motions of topological adjustment + ! PREC = if 0 then read RPAR & BCDATA with e20.12) + CALL SALGET(DATAIN,3,F_GEO,FOUT0,'index kndex prec') + INDEX=DATAIN(1) + KNDEX=DATAIN(2) + PREC=DATAIN(3) + ! + !* read epsilons for topological deformations + ! EPS = if the distance of two ends of elements < eps, + ! they will be united to one point + CALL SALGET(DATARE,1,F_GEO,FOUT0,'eps') + EPS=DATARE(1) + ! + IF(PRTIND >= 1) THEN + WRITE(FOUT,'(//,5X,''domain checkout:'',/, & + & 5X,''elements are in contact for distance < '',1P,E12.4,//)') EPS + ENDIF + ! + ALLOCATE(GG%NUM_MERGE(GG%NB_NODE),STAT =OK) + IF(OK /= 0) CALL XABORT('SAL100: failure to allocate NB_NODE') + CALL SALGET(GG%NUM_MERGE,GG%NB_NODE,F_GEO,FOUT0,'FLUX INDEX PER NODE') + IF(MAXVAL(GG%NUM_MERGE) /= GG%NB_FLUX) CALL XABORT('SAL100: inconsistent NBFLUX') + ALLOCATE(GG%NAME_MACRO(GG%NB_MACRO),STAT =OK) + IF(OK /= 0) CALL XABORT('SAL100: failure to allocate NB_MACRO') + CALL SALGET(GG%NAME_MACRO,GG%NB_MACRO,F_GEO,FOUT0,'NAMES OF MACROS') + ALLOCATE(GG%NUM_MACRO(GG%NB_FLUX)) + CALL SALGET(GG%NUM_MACRO,GG%NB_FLUX,F_GEO,FOUT0,'macro order number per flux region') + !* do the work (SAL100_2 is called here!): + CALL SAL110(GG) + ! + END SUBROUTINE SAL100 + ! + SUBROUTINE SAL110(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! constructs geometrical domain + ! + !Parameters: input/output + ! GG geometry basic information. + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,NBFOLD,NIPAR,NRPAR,ALLSUR,NANIS,IC,ISPEC,ANGGEO + USE SAL_TRACKING_TYPES, ONLY : PRTIND + !**** + IMPLICIT NONE + ! in variable + ! ************ + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + + ! local variable + ! *************** + INTEGER :: ELEM,OK,I,TYPE + CHARACTER(LEN=4) :: HSYM + REAL(PDB),PARAMETER :: CONV=PI/180._PDB + INTEGER, PARAMETER :: FOUT =6 + !**** + ! allocate node arrays + HSYM=' ' + IF((TYPGEO==0).OR.((TYPGEO==6).AND.(NBFOLD==0))) THEN + CONTINUE + ELSE IF((TYPGEO==1).AND.(NBFOLD==3)) THEN + HSYM='QUAR' + ELSE IF(((TYPGEO==1).AND.(NBFOLD==8)).OR.((TYPGEO==7).AND.(NBFOLD==0))) THEN + HSYM='EIGH' + ELSE IF(((TYPGEO==1).AND.(NBFOLD==6)).OR.((TYPGEO==8).AND.(NBFOLD==0))) THEN + HSYM='SIXT' + ELSE IF(((TYPGEO==1).AND.(NBFOLD==12)).OR.((TYPGEO==12).AND.(NBFOLD==0))) THEN + HSYM='S30' + ELSE IF(((TYPGEO==1).AND.(NBFOLD==4)).OR.((TYPGEO==3).AND.(NBFOLD==0))) THEN + HSYM='SYME' + ELSE IF((TYPGEO==1).AND.(NBFOLD==2)) THEN + HSYM='SY1D' + ELSE IF(((TYPGEO==2).AND.(NBFOLD==6)).OR.((TYPGEO==10).AND.(NBFOLD==0))) THEN + HSYM='RA60' + ELSE IF((TYPGEO==5).AND.(NBFOLD==0)) THEN + HSYM='TRAN' + ELSE IF((TYPGEO==9).AND.(NBFOLD==0)) THEN + HSYM='TRAN' + ELSE IF(((TYPGEO==2).AND.(NBFOLD==3)).OR.(TYPGEO==11).AND.(NBFOLD==0)) THEN + HSYM='R120' + ELSE + WRITE(*,*) "TYPGEO=",TYPGEO," NBFOLD=",NBFOLD + CALL XABORT('SAL110: non supported type of symmetry') + ENDIF + ALLOCATE (GG%IPAR(NIPAR,GG%NB_ELEM), GG%RPAR(NRPAR,GG%NB_ELEM), STAT=OK) + IF(OK/=0) CALL XABORT('SAL110: not enough memory I,R') + + !* read surfacic file + CALL SALINP(GG) + ! + !* unite the ends of elements, redefine elements + CALL SAL128(GG%RPAR,GG%IPAR,GG%NB_ELEM) + ! + IF((ISPEC==0).AND.(IC.EQ.4)) THEN + !* unfold domain + IF(HSYM=='QUAR') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: DIAG unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + ANGGEO=2.0*ANGGEO + CALL SALFOLD_1('DIAG',GG) + ELSE IF(HSYM=='EIGH') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: DIAG + SYME unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + ANGGEO=2.0*ANGGEO + CALL SALFOLD_1('DIAG',GG) + CALL SALFOLD_1('SYMX',GG) + CALL SALFOLD_1('SYMY',GG) + ELSE IF(HSYM=='SYME') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: SYME unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_1('SYMX',GG) + CALL SALFOLD_1('SYMY',GG) + ELSE IF(HSYM=='SY1D') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: SY1D unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_1('SYMX',GG) + ELSE IF(HSYM=='SIXT') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: SA60 unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_1('SA60',GG) + CALL SALFOLD_1('SYMH',GG) + ELSE IF(HSYM=='S30') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: S30 unfold" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_1('S30 ',GG) + CALL SALFOLD_1('SB60',GG) + CALL SALFOLD_1('SYMX',GG) + ELSE IF(HSYM=='RA60') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: SR60 unfold with rotation" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_2('SR60',GG) + CALL SALFOLD_2('R180',GG) + ELSE IF(HSYM=='R120') THEN + IF(PRTIND>0) WRITE(*,*) "SAL110: SR120 unfold with rotation" + IF(NANIS>1) CALL XABORT('SAL110: unfold unsupported with NANIS>1') + CALL SALFOLD_2('R120',GG) + ELSE + IF(PRTIND>0) WRITE(*,*) "SAL110: no unfold" + ENDIF + IF((TYPGEO/=0).AND.(NBFOLD==0)) THEN + TYPGEO=6 + ELSE + TYPGEO=0; NBFOLD=0 + ENDIF + ENDIF + IF(PRTIND>0) WRITE(FOUT,*) 'SAL110: after unfolding -- NB_ELEM=',GG%NB_ELEM, & + & ' NB_PERIM=',GG%NBBCDA + + IF(PRTIND>5) THEN + !* print surfacic file + WRITE(FOUT,'(5H--cut,75(1H-))') + WRITE(FOUT,'(5HBEGIN)') + WRITE(FOUT,'(42H* typgeo nbfold nbnode nbelem nbmacr nbreg)') + WRITE(FOUT,'(6I7)') TYPGEO,NBFOLD,GG%NB_NODE,GG%NB_ELEM,GG%NB_MACRO,GG%NB_NODE + WRITE(FOUT,'(20H* index kndex prec)') + WRITE(FOUT,'(4I7)') 0,0,1 + WRITE(FOUT,'(18H* eps eps0)') + WRITE(FOUT,'(1P,2E18.9)') 1.0E-03,1.0E-05 + WRITE(FOUT,'(20H* num_of_region/mesh)') + WRITE(FOUT,'(10I7)') (GG%NUM_MERGE(I),I=1,GG%NB_NODE) + WRITE(FOUT,'(13H* macro names)') + WRITE(FOUT,'(4(3x,a10,2x))') (GG%NAME_MACRO(I),I=1,GG%NB_MACRO) + WRITE(FOUT,'(35H* macro_order_index_per_flux_region)') + WRITE(FOUT,'(10I7)') (GG%NUM_MACRO(I),I=1,GG%NB_FLUX) + DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + WRITE(FOUT,'(7h elem =,I6)') ELEM + WRITE(FOUT,'(22H*type node- node+)') + WRITE(FOUT,'(3I6)') (GG%IPAR(I,ELEM),I=1,3) + WRITE(FOUT,'(63H*cx cy ex_or_R ey_or_theta1 theta2)') + IF(TYPE<=2) THEN + WRITE(FOUT,'(1P,5E18.9)') (GG%RPAR(I,ELEM),I=1,5) + ELSE IF(TYPE==3) THEN + WRITE(FOUT,'(1P,5E18.9)') (GG%RPAR(I,ELEM),I=1,3),GG%RPAR(4,ELEM)/CONV, & + (GG%RPAR(5,ELEM)-GG%RPAR(4,ELEM))/CONV + ENDIF + ENDDO + WRITE(FOUT,'(40H*defaul nbbcda allsur divsur ndivsur)') + WRITE(FOUT,'(1P,5I8)') GG%DEFAUL,GG%NBBCDA,ALLSUR,0,0 + WRITE(FOUT,'(17H*albedo deltasur)') + WRITE(FOUT,'(1P,2E18.9)') GG%ALBEDO,0.0 + DO ELEM=1,GG%NBBCDA + WRITE(FOUT,'(37H particular boundary condition number,i12)') ELEM + WRITE(FOUT,'(13H*type nber)') + WRITE(FOUT,'(1P,2I8)') GG%BCDATAREAD(ELEM)%SALTYPE,GG%BCDATAREAD(ELEM)%NBER + WRITE(FOUT,'(14H*elems(1,nber))') + WRITE(FOUT,'(1P,10I8)') (GG%BCDATAREAD(ELEM)%ELEMNB(I),I=1,GG%BCDATAREAD(ELEM)%NBER) + IF(GG%BCDATAREAD(ELEM)%SALTYPE==0) THEN + WRITE(FOUT,'(7H*albedo)') + WRITE(FOUT,'(1P,E18.9)') GG%BCDATAREAD(ELEM)%BCDATA(6) + ELSE + WRITE(FOUT,'(22H*cx cy angle)') + WRITE(FOUT,'(1P,3E18.9)') GG%BCDATAREAD(ELEM)%BCDATA(1:2),GG%BCDATAREAD(ELEM)%BCDATA(5)*180._PDB/PI + ENDIF + ENDDO + ENDIF + + ! allocate media and element arrays + ALLOCATE (GG%VOL_NODE(GG%NB_NODE),GG%PPERIM_NODE(GG%NB_NODE+1),GG%IBC2_ELEM(GG%NB_ELEM), & + GG%ISURF2_ELEM(GG%NB_ELEM),GG%MED(GG%NB_NODE), STAT=OK) + IF(OK/=0) CALL XABORT('SAL110: not enough memory VOL') + GG%ISURF2_ELEM(:GG%NB_ELEM)=0 + + !* 2D boundary conditions and macro contacts: + ! - defines NB_BC2, NBSUR2 + ! - defines surface strctures for each 2D macro + ! - defines perimeter structure for each 2D macro + ! - read 2D boundary conditions + CALL SAL130(GG) + ! + !* topological check + CALL SAL140(GG%NB_NODE,GG%RPAR,GG%IPAR,GG%PPERIM_NODE,GG%PERIM_NODE) + ! + !* volumes, surfaces, put local nbers in node, and read media: + CALL SAL160(GG) + IF(PRTIND>5) THEN + WRITE(FOUT,'(12H* mil(nbreg))') + WRITE(FOUT,'(10I7)') (GG%MED(I),I=1,GG%NB_NODE) + WRITE(FOUT,'(3HEND)') + WRITE(FOUT,'(5H--cut,75(1H-))') + ENDIF + ! + !* printout basic domain + CALL SAL170(GG) + ! + END SUBROUTINE SAL110 + ! + SUBROUTINE SALINP(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! read surfacic file + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : ALLSUR,PREC + !**** + IMPLICIT NONE + ! in variable + ! ************ + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + !**** + INTEGER, PARAMETER :: N_DATAIN=25 + INTEGER, DIMENSION (N_DATAIN) :: DATAIN + INTEGER :: ELEM,I,TYPE,NBER + INTEGER, PARAMETER, DIMENSION(0:4) :: READ_BC_LEN=(/1,1,3,3,3/) + INTEGER, PARAMETER :: FOUT =6 + REAL(PDB) :: ANGLE,BCDATA_TDT(3) + ! + !* read element data + DO ELEM=1,GG%NB_ELEM + CALL SAL126(GG%RPAR(:,ELEM),GG%IPAR(:,ELEM)) + ENDDO + !* read + ! DEFAUL = default bc condition + ! 0 = all external surfaces are isotropically reflected + ! with the same ALBEDO = ALBED0 + ! 1 = all external surfaces are reflected + ! NBBCDA = number of groups of bc conditions + ! + CALL SALGET(DATAIN,3,F_GEO,FOUT0,'general bc data') + GG%DEFAUL=DATAIN(1) + GG%NBBCDA=DATAIN(2) + ALLSUR=DATAIN(3) + ! we can define only two default bc's + IF(GG%DEFAUL>4.OR.GG%DEFAUL<0) THEN + WRITE(FOUT,'(8H defaul=,I5)') GG%DEFAUL + CALL XABORT('SALINP: wrong default bc type') + ENDIF + !* read albedo : defaul bcdata + CALL SALGET(GG%ALBEDO,F_GEO,FOUT0,PREC,'GENERAL ALBEDO') + ! + !* read detailed bcdata if required (motions) + LBCDIAG=.FALSE. + IF(GG%NBBCDA>0)THEN + ALLOCATE(GG%BCDATAREAD(GG%NBBCDA)) + DO I=1,GG%NBBCDA + GG%BCDATAREAD(I)%BCDATA(:6)=0.0 + CALL SALGET(DATAIN,2,F_GEO,FOUT0,'SPECIFIC BC: TYPE NBER') + TYPE=DATAIN(1) + NBER=DATAIN(2) + GG%BCDATAREAD(I)%SALTYPE=TYPE + GG%BCDATAREAD(I)%NBER=NBER + IF(TYPE>5.OR.TYPE<0) CALL XABORT('SALINP: wrong bc type') + IF(NBER>GG%NB_ELEM) CALL XABORT('SALINP: bc def exceeds nber of elements') + ! + ! read order numbers of elements affected + ALLOCATE(GG%BCDATAREAD(I)%ELEMNB(NBER)) + CALL SALGET(GG%BCDATAREAD(I)%ELEMNB,NBER,F_GEO,FOUT0,'BC ELEMENTS') + ! read bc motion + CALL SALGET(BCDATA_TDT,READ_BC_LEN(TYPE),F_GEO,FOUT0,PREC,'data for specific bc condition') + IF(READ_BC_LEN(TYPE).EQ.1) THEN + GG%BCDATAREAD(I)%BCDATA(1:5)=0._PDB + GG%BCDATAREAD(I)%BCDATA(6)=BCDATA_TDT(1) + ELSE + GG%BCDATAREAD(I)%BCDATA(1:2)=BCDATA_TDT(1:2) + ANGLE=BCDATA_TDT(3)*PI/180._PDB + GG%BCDATAREAD(I)%BCDATA(3)=COS(ANGLE) + GG%BCDATAREAD(I)%BCDATA(4)=SIN(ANGLE) + GG%BCDATAREAD(I)%BCDATA(5)=ANGLE + GG%BCDATAREAD(I)%BCDATA(6)=1._PDB + LBCDIAG=LBCDIAG.OR.((GG%BCDATAREAD(I)%BCDATA(1)==0._PDB).AND.(GG%BCDATAREAD(I)%BCDATA(2)==0._PDB) & + .AND.(GG%BCDATAREAD(I)%BCDATA(5)==PI/4._PDB)) + ENDIF + ENDDO + ENDIF + END SUBROUTINE SALINP + ! + SUBROUTINE SAL126(RPAR,IPAR) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! reads one element data + ! + !Parameters: input/output + ! RPAR floating point geometry descriptors + ! IPAR integer geometry descriptors + ! + !--------------------------------------------------------------------- + ! + !**** + USE SAL_GEOMETRY_TYPES, ONLY : PREC,G_ELE_TYPE + !**** + IMPLICIT NONE + ! in variable + INTEGER, INTENT(OUT), DIMENSION(:) :: IPAR + REAL(PDB), INTENT(INOUT), DIMENSION(:) :: RPAR + !**** + ! local variable + INTEGER :: TYPE,NBER,IAUX + REAL(PDB) :: PHI1,PHI2,ANGMAX,DELPHI + !**** + ! IPAR(1) = 1 (segment), 2 (circle), 3(arc of circle) + ! IPAR(2) = order nber of node in - side + ! IPAR(3) = order nber of node in + side + ! + !* read integer element data + CALL SALGET(IPAR,3,F_GEO,FOUT0,'integer descriptors') + TYPE=IPAR(1) + ! + ! ** segment: + ! R = C+T*F with T in (0,1) + ! RPAR(1~5) = CX CY FX FY F + ! + ! ** arc of circle: + ! R = C+R*F(THETA) with + ! THETA in (THETA1 < THETA2) + ! with THETA1 in (0,360) + ! RPAR(1~5) = CX CY R THETA1 THETA2 (in degrees) + ! + !* read real element data + + IF(TYPE == G_ELE_TYPE(1)) THEN + NBER=4 + ELSE IF(TYPE == G_ELE_TYPE(2)) THEN + NBER=3 + ELSE IF(TYPE == G_ELE_TYPE(3)) THEN + NBER=5 + ANGMAX=360._PDB + ELSE + WRITE(FOUT0,'(1X,''==> SAL126: unknown type '',I3)') TYPE + CALL XABORT('SAL126: unknown element type') + ENDIF + + CALL SALGET(RPAR,NBER,F_GEO,FOUT0,PREC,'real descriptors') + + IF(TYPE == G_ELE_TYPE(1)) THEN + ! segment: compute length + RPAR(5)=SQRT(RPAR(3)*RPAR(3)+RPAR(4)*RPAR(4)) + RPAR(6)=0._PDB + ELSE IF(TYPE == G_ELE_TYPE(2)) THEN + ! full circle: set angles + RPAR(4)=0._PDB + RPAR(5)=TWOPI + RPAR(6)=0._PDB + ELSE IF((TYPE == G_ELE_TYPE(3)).OR.(TYPE == G_ELE_TYPE(4))) THEN + ! check angles + PHI1=RPAR(NBER-1) + DELPHI=RPAR(NBER) + ! order angles in increasing values: + IF(DELPHI>0._PDB)THEN + IF(DELPHI>ANGMAX)THEN + WRITE(FOUT0,'(1X,''==> SAL126: DELPHI = '',1P,E12.4, & + &'' > '',1P,E12.4,'' FOR TYPE'',I3)')DELPHI,ANGMAX,TYPE + CALL XABORT('SAL126: invalid value of delphi') + ENDIF + PHI2=PHI1+DELPHI + ELSE + IF(DELPHI<-ANGMAX)THEN + WRITE(FOUT0,'(1X,''==> SAL126: DELPHI = '',1P,E12.4, & + &'' < '',1P,E12.4,'' FOR TYPE'',I3)')DELPHI,-ANGMAX,TYPE + CALL XABORT('SAL126: invalid value of delphi') + ENDIF + PHI2=PHI1 + PHI1=PHI1+DELPHI + ENDIF + IF(TYPE==G_ELE_TYPE(3))THEN + ! arc of circle: put phi1 within 0 and 360. + IF(PHI1>360._PDB)THEN + IAUX=INT(PHI1/360._PDB) + DELPHI=360._PDB*IAUX + PHI2=PHI2-DELPHI + PHI1=PHI1-DELPHI + ELSEIF(PHI1<0._PDB)THEN + IAUX=INT((-PHI1+1.D-7)/360._PDB)+1 + DELPHI=360._PDB*IAUX + PHI2=PHI2+DELPHI + PHI1=PHI1+DELPHI + ENDIF + RPAR(6)=0._PDB + ELSE + CALL XABORT('SAL126: unsupported option') + ENDIF + ! convert to radians + RPAR(NBER-1)=PHI1*(TWOPI/360._PDB) + RPAR(NBER)=PHI2*(TWOPI/360._PDB) + ENDIF + ! + END SUBROUTINE SAL126 + ! + SUBROUTINE SAL128(RPAR,IPAR,NB_ELEM) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! create library of ends of elements and redefine element discriptors + ! + !Parameters: input/output + ! RPAR floating point geometry descriptors + ! IPAR integer geometry descriptors + ! NB_ELEM number of surfacic elements + ! + !--------------------------------------------------------------------- + ! + !**** + USE SAL_GEOMETRY_TYPES, ONLY : G_ELE_TYPE,EPS,KNDEX,TYPGEO + ! in variable + !************ + IMPLICIT NONE + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: RPAR + INTEGER , INTENT(IN) :: NB_ELEM + ! local variable + ! *************** + REAL(PDB), DIMENSION (:,:), ALLOCATABLE :: POINT ! coordinates of the ends of elements + INTEGER, DIMENSION (:,:), ALLOCATABLE :: ELEM_NEW ! point nbers of two ends of elements + !**** + ! NN = total number of points + ! NBP = count of number of ends having the same coordinates + INTEGER :: TYPE,NN,ELEM,I,P,POS,I1,I2,OK + REAL(PDB) :: X(2),Y(2),D,DX,DY + REAL :: EPS2 + LOGICAL :: LGONE + INTEGER, DIMENSION(NB_ELEM) :: NBP + INTEGER, PARAMETER :: FOUT =6 + !**** + ! - POINT(2,NB_ELEM) = coordinates of the ends of elements + ! - ELEM_NEW(2,NB_ELEM) = point nbers of two ends of elements + ! + ALLOCATE(POINT(2,NB_ELEM),ELEM_NEW(2,NB_ELEM),STAT=OK) + IF(OK/=0) CALL XABORT('SAL100_3_1: not enough memory I,CH') + + EPS2=EPS*EPS + IF(KNDEX/=0)WRITE(FOUT,'(//,& + &" reunite the ends of elements of distance less than ",E13.4)') EPS + ! + !* nn counts the nber of points + NBP(1:NB_ELEM)=0 + NN=0 + DO ELEM=1,NB_ELEM + TYPE=IPAR(1,ELEM) + IF(TYPE/=G_ELE_TYPE(2)) THEN + DO I=1,2 + ! get coordinates of end + CALL SAL141(TYPE,RPAR(:,ELEM),X(I),Y(I),I) + ! find if there is a defined neighbouring point + LGONE=.FALSE. + DO P=1,NN + DX=POINT(1,P)-X(I) + DY=POINT(2,P)-Y(I) + D=DX*DX+DY*DY + IF(D<EPS2) THEN + LGONE=.TRUE. + POS=P + EXIT + ENDIF + ENDDO + IF(LGONE) THEN + ! choose average value + X(I)=(X(I)+POINT(1,POS)*NBP(POS))/REAL(NBP(POS)+1) + Y(I)=(Y(I)+POINT(2,POS)*NBP(POS))/REAL(NBP(POS)+1) + ELSE + ! add points + NN=NN+1 + IF(NN > NB_ELEM) CALL XABORT('SAL128: point overflow') + POS=NN + ENDIF + POINT(1,POS)=X(I) + POINT(2,POS)=Y(I) + NBP(POS)=NBP(POS)+1 + ! define the end of element + ELEM_NEW(I,ELEM)=POS + ENDDO + ELSE + ELEM_NEW(1:2,ELEM)=0 + ENDIF + ENDDO + !* adjust points on the axes + SELECT CASE(TYPGEO) + CASE(5:6) + CALL SAL128_3(NN,POINT,EPS,EPS2) + CASE(7) + CALL SAL128_4(NN,POINT,EPS,EPS2) + CASE(8,12) + CALL SAL128_5(NN,POINT,EPS,EPS2) + CASE(9) + CALL SAL128_6(NN,POINT,EPS,EPS2) + CASE(10) + CALL SAL128_7(NN,POINT,EPS,EPS2) + CASE(11) + CALL SAL128_8(NN,POINT,EPS,EPS2) + END SELECT + + !* redefine the elements in rpar + IF(KNDEX/=0)WRITE(FOUT,'(/," redefine elements : ",/)') + DO ELEM=1,NB_ELEM + TYPE=IPAR(1,ELEM) + IF(TYPE/=G_ELE_TYPE(2)) THEN + I1=ELEM_NEW(1,ELEM) + I2=ELEM_NEW(2,ELEM) + CALL SAL129(POINT(1,I1),POINT(2,I1),POINT(1,I2),POINT(2,I2), & + TYPE,RPAR(:,ELEM),ELEM,KNDEX,FOUT) + ENDIF + ENDDO + ! + DEALLOCATE(POINT,ELEM_NEW) + ! + END SUBROUTINE SAL128 + ! + SUBROUTINE SAL128_3(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process retangular geometry with translation or symmetry boundary condition + ! (TYPGEO=5,6): adjusts points on the axes, computes the lengths of the + ! retangle sides; in case of with translation,element lengths on the opposite + ! axes will be adjusted to be the same + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + ! DIMENSION POINT(2,NN) + !**** + INTEGER :: NP,NAXIS,I,J,K,M,P,IPOINT(4,NN),BP(4),IP(4) + LOGICAL :: LGADJUST + REAL(PDB) :: D,AUX,EX(4),EY(4),DIS(4),AX(4),AY(4),D_AXIS(4,NN) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! flag to adjust the element lengths on the opposite axes + LGADJUST=TYPGEO==5 + ! compute sides of the rectangle + LX=0.; LY=0. + DO P=1,NN + IF(ABS(POINT(1,P))<EPS) THEN + IF(LY<POINT(2,P)) LY=POINT(2,P) + ENDIF + IF(ABS(POINT(2,P))<EPS) THEN + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDIF + ENDDO + !* axis directions and their distance to (0,0) + NAXIS=4 + EX(1)=1.; EY(1)=0.; DIS(1)=0. + EX(2)=1.; EY(2)=0.; DIS(2)=-LY + EX(3)=0.; EY(3)=1.; DIS(3)=0. + EX(4)=0.; EY(4)=1.; DIS(4)=LX + NP=4 + AX(1)=0.; AY(1)=0. + AX(2)=0.; AY(2)=LY + AX(3)=LX; AY(3)=0. + AX(4)=LX; AY(4)=LY + ! beginning point of the axes + IF(LGADJUST) THEN + ! axis numbering: + ! M=2 + ! ************ + ! * * + ! M=3 * * M=4 + ! ************ + ! M=1 + ! + BP(1)=1; BP(2)=2; BP(3)=1; BP(4)=3 + IP=0 + ENDIF + ITER0:DO P=1,NN + DO I=1,NP + ! if it is the 4 corners of the rectangle + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the corners + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + IF(LGADJUST) THEN + ! put distance to the axis origin + SELECT CASE(I) + CASE(2) + ! end of axis 3,keep its distance and number + IP(3)=IP(3)+1 + D_AXIS(3,IP(3))=AY(I) + IPOINT(3,IP(3))=P + CASE(3) + ! end of axis 1,keep its distance and number + IP(1)=IP(1)+1 + D_AXIS(1,IP(1))=AX(I) + IPOINT(1,IP(1))=P + CASE(4) + ! end of axes 2&4,keep its distance and number + IP(2)=IP(2)+1 + D_AXIS(2,IP(2))=AX(I) + IPOINT(2,IP(2))=P + IP(4)=IP(4)+1 + D_AXIS(4,IP(4))=AY(I) + IPOINT(4,IP(4))=P + END SELECT + ENDIF + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + IF(LGADJUST) THEN + ! compute distance to the axis origin + D=SQRT((POINT(1,P)-AX(BP(I)))**2+(POINT(2,P)-AY(BP(I)))**2) + IP(I)=IP(I)+1 + D_AXIS(I,IP(I))=D + IPOINT(I,IP(I))=P + ENDIF + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + IF(LGADJUST) THEN + DO I=1,NAXIS,2 + IF(IP(I)/=IP(I+1)) & + CALL XABORT('SAL128_3: axial points nber not the same,axis') + ! sort the 'd_axis' table + DO J=I,I+1 + DO P=1,IP(J) + DO K=P+1,IP(J) + IF(D_AXIS(J,P)>D_AXIS(J,K)) THEN + D=D_AXIS(J,P) + D_AXIS(J,P)=D_AXIS(J,K) + D_AXIS(J,K)=D + M=IPOINT(J,P) + IPOINT(J,P)=IPOINT(J,K) + IPOINT(J,K)=M + ENDIF + ENDDO + ENDDO + ENDDO + DO P=1,IP(I) + IF(ABS(D_AXIS(I,P)-D_AXIS(I+1,P))>EPS) THEN + IF(KNDEX/=0) & + WRITE(FOUT,'(" warning: too great axial length difference",& + & 2(/,2X,"axis = ",I3," point = ",I3," d = ",E13.6))')& + I,P,D_AXIS(I,P),I+1,P,D_AXIS(I+1,P) + ENDIF + D=(D_AXIS(I,P)+D_AXIS(I+1,P))*0.5 + DO J=I,I+1 + K=IPOINT(J,P) + POINT(1,K)=D*EX(J)+AX(BP(J)) + POINT(2,K)=D*EY(J)+AY(BP(J)) + ENDDO + ENDDO + ENDDO + ENDIF + ! + END SUBROUTINE SAL128_3 + ! + SUBROUTINE SAL128_4(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process a 1/8 square geometry (typgeo=7): adjusts points on the axes, + ! computes the square side length + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + !**** + INTEGER :: NP,NAXIS,I,P + REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! compute the square sides + LX=0.; LY=0. + DO P=1,NN + IF(LY<POINT(2,P)) LY=POINT(2,P) + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDDO + IF(KNDEX>0.AND.ABS(LY-LX)>EPS) & + WRITE(FOUT,'(" warning: the square sides are not the same!", & + & " lx = ",E12.4," ly = ",E12.4)')LX,LY + LX=(LX+LY)*0.5 + LY=LX + !* adjust points on the axes + NAXIS=3 + ! axis directions + EX(1)=1.; EY(1)=0.; DIS(1)=0. + EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0. + EX(3)=0.; EY(3)=1.; DIS(3)=LX + NP=3 + ! axis ends + AX(1)=0.; AY(1)=0. + AX(2)=LX; AY(2)=0. + AX(3)=LX; AY(3)=LY + ITER0:DO P=1,NN + DO I=1,NP + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the axis origin + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + END SUBROUTINE SAL128_4 + ! + SUBROUTINE SAL128_5(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process a trianglar geometry (TYPGEO=8 or 12) with specular reflection: + ! adjusts points on the axes, computes the triangular side length + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + !**** + INTEGER :: NP,NAXIS,I,P + REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! compute the triangular side + LX=0.; LY=0. + DO P=1,NN + IF(ABS(POINT(2,P))<EPS) THEN + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDIF + IF(LY<POINT(2,P)) LY=POINT(2,P) + ENDDO + IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) & + WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", & + & " a = ",E12.4," h = ",E12.4)')LX,LY + ! adjust h according a + LY=LX*SIN(ANGGEO) + !* adjust points on the axes + NAXIS=3 + ! axis directions + EX(1)=1.; EY(1)=0.; DIS(1)=0. + EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0. + EX(3)=COS(ANGGEO*2.); EY(3)=SIN(ANGGEO*2.); DIS(3)=LY + NP=3 + ! axis ends + AX(1)=0.; AY(1)=0. + AX(2)=LX; AY(2)=0. + AX(3)=LX*0.5; AY(3)=LY + ITER0:DO P=1,NN + DO I=1,NP + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the axis origin + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + END SUBROUTINE SAL128_5 + ! + SUBROUTINE SAL128_6(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process a hexagonal geometry with translations on all sides (typgeo=9): + ! adjusts points on the axes, computes the lengths of hexagon sides; + ! element lengths on the opposite axes will be adjusted to be the same + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,KNDEX,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + !**** + INTEGER :: NP,NAXIS,I,J,K,M,P,IPOINT(6,NN),BP(6),IP(6) + REAL(PDB) :: D,AUX,EX(6),EY(6),DIS(6),AX(6),AY(6),D_AXIS(6,NN) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! compute sides of the hexagon: LX=SIDE LENGTH,LY=LX*SIN(A) + LX=0.; LY=0. + DO P=1,NN + IF(ABS(POINT(2,P))<EPS) THEN + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDIF + IF(LY<POINT(2,P)) LY=POINT(2,P) + ENDDO + IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) & + WRITE(FOUT,'(" warning: h is different from A*SIN(ANGGEO)!", & + & " a = ",E12.4," h = ",E12.4)')LX,LY + !* axis directions and their distance to (0,0) + NAXIS=6 + EX(1)=1.; EY(1)=0.; DIS(1)=LY + EX(2)=1.; EY(2)=0.; DIS(2)=-LY + EX(3)=COS(ANGGEO); EY(3)=SIN(ANGGEO); DIS(3)=LY + EX(4)=EX(3); EY(4)=EY(3); DIS(4)=-LY + EX(5)=COS(2.*ANGGEO); EY(5)=SIN(2.*ANGGEO); DIS(5)=LY + EX(6)=EX(5); EY(6)=EY(5); DIS(6)=-LY + NP=6 + AX(1)=-LX*0.5; AY(1)=-LY + AX(2)=LX*0.5; AY(2)=-LY + AX(3)=LX; AY(3)=0. + AX(4)=LX*0.5; AY(4)=LY + AX(5)=-LX*0.5; AY(5)=LY + AX(6)=-LX; AY(6)=0. + ! axis numbering: + ! M=2 + ! ***** + ! M=4 * * M=5 + ! * * * + ! M=6 * * M=3 + ! ***** + ! M=1 + ! + BP(1)=1; BP(2)=5; BP(3)=2; BP(4)=6; BP(5)=3; BP(6)=1 + IP=0 + ITER0:DO P=1,NN + DO I=1,NP + ! if it is the 6 corners of the hexagon + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the corner + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + ! put distance to the axis origin + J=0 + SELECT CASE(I) + CASE(1) + CYCLE ITER0 + CASE(2) + ! END OF AXIS 1 + J=1 + CASE(3) + ! END OF AXIS 3 + J=3 + CASE(4) + ! END OF AXES 2&5 + J=2 + CASE(5) + ! END OF AXIS 4 + J=4 + CASE(6) + ! END OF AXIS 6 + J=6 + END SELECT + ! keep its distance and number +10 IP(J)=IP(J)+1 + D_AXIS(J,IP(J))=LX + IPOINT(J,IP(J))=P + IF(J==2) THEN + J=5 + GOTO 10 + ENDIF + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=(POINT(1,P)-AX(BP(I)))*EX(I)+(POINT(2,P)-AY(BP(I)))*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + ! compute distance to the axis origin + D=SQRT((POINT(1,P)-AX(BP(I)))**2+(POINT(2,P)-AY(BP(I)))**2) + IP(I)=IP(I)+1 + D_AXIS(I,IP(I))=D + IPOINT(I,IP(I))=P + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + DO I=1,NAXIS,2 + IF(IP(I)/=IP(I+1)) & + CALL XABORT('SAL128_6: axial points nber not the same,axis') + ! sort the 'd_axis' table + DO J=I,I+1 + DO P=1,IP(J) + DO K=P+1,IP(J) + IF(D_AXIS(J,P)>D_AXIS(J,K)) THEN + D=D_AXIS(J,P) + D_AXIS(J,P)=D_AXIS(J,K) + D_AXIS(J,K)=D + M=IPOINT(J,P) + IPOINT(J,P)=IPOINT(J,K) + IPOINT(J,K)=M + ENDIF + ENDDO + ENDDO + ENDDO + DO P=1,IP(I) + IF(ABS(D_AXIS(I,P)-D_AXIS(I+1,P))>EPS) THEN + IF(KNDEX/=0) & + WRITE(FOUT,'(" warning: too great axial length difference",& + & 2(/,2X,"axis = ",I3," point = ",I3," d = ",E13.6))')& + I,P,D_AXIS(I,P),I+1,P,D_AXIS(I+1,P) + ENDIF + D=(D_AXIS(I,P)+D_AXIS(I+1,P))*0.5 + DO J=I,I+1 + K=IPOINT(J,P) + POINT(1,K)=D*EX(J)+AX(BP(J)) + POINT(2,K)=D*EY(J)+AY(BP(J)) + ENDDO + ENDDO + ENDDO + ! + END SUBROUTINE SAL128_6 + ! + SUBROUTINE SAL128_7(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process an isocel trianglar geometry (TYPGEO=10) with RA60 rotation: + ! adjusts points on the axes, computes the triangular side length + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + !**** + INTEGER :: NP,NAXIS,I,P + REAL(PDB) :: D,AUX,EX(3),EY(3),DIS(3),AX(3),AY(3) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! compute the triangular side + LX=0.; LY=0. + DO P=1,NN + IF(ABS(POINT(2,P))<EPS) THEN + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDIF + IF(LY<POINT(2,P)) LY=POINT(2,P) + ENDDO + IF(KNDEX>0.AND.ABS(LY-LX*SIN(ANGGEO))>EPS) & + WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", & + & " a = ",E12.4," h = ",E12.4)')LX,LY + ! adjust h according a + LY=LX*SIN(ANGGEO) + !* adjust points on the axes + NAXIS=3 + ! axis directions + EX(1)=1.; EY(1)=0.; DIS(1)=0. + EX(2)=COS(ANGGEO); EY(2)=SIN(ANGGEO); DIS(2)=0. + EX(3)=COS(ANGGEO*2.); EY(3)=SIN(ANGGEO*2.); DIS(3)=LY + NP=3 + ! axis ends + AX(1)=0.; AY(1)=0. + AX(2)=LX; AY(2)=0. + AX(3)=LX*0.5; AY(3)=LY + ITER0:DO P=1,NN + DO I=1,NP + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the axis origin + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + END SUBROUTINE SAL128_7 + ! + SUBROUTINE SAL128_8(NN,POINT,EPS,EPS2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! process R120 lozenge geometry (TYPGEO=11): adjusts points on the axes, + ! computes the lozenge side length + ! + !Parameters: input + ! NN total number of points + ! EPS when distance of points to axes less than EPS, displace + ! the points onto the axes + ! EPS2 variable set to EPS*EPS + ! + !Parameters: input/output + ! POINT coordinates of points + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,LX=>LENGTHX,LY=>LENGTHY + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NN + REAL, INTENT(IN) :: EPS,EPS2 + REAL(PDB), INTENT(INOUT), DIMENSION(:,:) :: POINT + !**** + INTEGER :: NP,NAXIS,I,P + REAL(PDB) :: D,AUX,EX(4),EY(4),DIS(4),AX(4),AY(4) + INTEGER, PARAMETER :: FOUT =6 + !**** + ! compute the triangular side + LX=0.; LY=0. + DO P=1,NN + IF(ABS(POINT(2,P))<EPS) THEN + IF(LX<POINT(1,P)) LX=POINT(1,P) + ENDIF + IF(LY<POINT(2,P)) LY=POINT(2,P) + ENDDO + IF(ABS(LY-LX*SIN(ANGGEO))>EPS) & + WRITE(FOUT,'(" warning: h is different from a*sin(anggeo)!", & + & " a = ",E12.4," h = ",E12.4)') LX,LY + ! adjust h according to base + LY=LX*SIN(ANGGEO) + !* adjust points on the axes + NAXIS=4 + ! axis directions + EX(1)=1.; EY(1)=0.; DIS(1)=0. + EX(2)=1.; EY(2)=0.; DIS(2)=-LY + EX(3)=COS(ANGGEO); EY(3)=SIN(ANGGEO); DIS(3)=0. + EX(4)=COS(ANGGEO); EY(4)=SIN(ANGGEO); DIS(4)=LY + NP=4 + ! axis ends + AX(1)=0. ; AY(1)=0. + AX(2)=LX*0.5; AY(2)=LY + AX(3)=LX ; AY(3)=0. + AX(4)=LX*1.5; AY(4)=LY + ITER0:DO P=1,NN + DO I=1,NP + D=(POINT(1,P)-AX(I))**2+(POINT(2,P)-AY(I))**2 + IF(D.LT.EPS2) THEN + ! move point to the axis origin + POINT(1,P)=AX(I) ; POINT(2,P)=AY(I) + CYCLE ITER0 + ENDIF + ENDDO + DO I=1,NAXIS + ! distance to the axis + D=POINT(1,P)*EY(I)-POINT(2,P)*EX(I)-DIS(I) + ! when aux>0,the point is near the half axis + AUX=POINT(1,P)*EX(I)+POINT(2,P)*EY(I) + IF(ABS(D).LT.EPS.AND.AUX>=0) THEN + ! move point to the axis + POINT(1,P)=POINT(1,P)-D*EY(I) + POINT(2,P)=POINT(2,P)+D*EX(I) + CYCLE ITER0 + ENDIF + ENDDO + ENDDO ITER0 + ! + END SUBROUTINE SAL128_8 + ! + SUBROUTINE SAL129(X1,Y1,X2,Y2,TYPE,RPAR,II,KNDEX,FOUT) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! recompute element data by giving two ends of the element + ! + !Parameters: input + ! X1 new position in X for the end 1 + ! Y1 new position in Y for the end 1 + ! X2 new position in X for the end 2 + ! Y2 new position in Y for the end 2 + ! TYPE type of element 1 (segment) 3 (arc of circle) + ! II element to be changed + ! KNDEX if not 0 print out data when motion>EPS + ! FOUT printout file + ! + !Parameters: input/output + ! RPAR descriptors of the element + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : NRPAR,EPS,G_ELE_TYPE,G_ELE_LEN + USE SAL_NUMERIC_MOD, ONLY : SALACO + !**** + IMPLICIT NONE + ! in variable + ! *********** + INTEGER, INTENT(IN) :: TYPE,KNDEX,FOUT,II + REAL(PDB), INTENT(IN) :: X1,Y1,X2,Y2 + REAL(PDB), INTENT(INOUT), DIMENSION(:) :: RPAR + !**** + + ! local variable + ! ************** + REAL(PDB) :: CX,CY,DELX,DELY,R,THETA,SX,SY,VX,VY,V2,DIST,MOV1,MOV2 + INTEGER :: I + REAL(PDB) :: RPAR_OLD(NRPAR) + !**** + ! keep old parameters: + RPAR_OLD(1:NRPAR)=RPAR(1:NRPAR) + ! + MOV1=0._PDB; MOV2=0._PDB; + IF(TYPE == G_ELE_TYPE(1)) THEN + !* segment : + ! compute the motions of two ends + DELX=X1-RPAR(1) + DELY=Y1-RPAR(2) + MOV1=SQRT(DELX*DELX+DELY*DELY) + DELX=X2-(RPAR(1)+RPAR(3)) + DELY=Y2-(RPAR(2)+RPAR(4)) + MOV2=SQRT(DELX*DELX+DELY*DELY) + ! set new ends, compute direction + RPAR(1)=X1 + RPAR(2)=Y1 + RPAR(3)=X2-X1 + RPAR(4)=Y2-Y1 + RPAR(5)=SQRT(RPAR(3)*RPAR(3)+RPAR(4)*RPAR(4)) + ELSE IF(TYPE == G_ELE_TYPE(3)) THEN + !* arc of circle : + ! get old radius + CX=RPAR(1); CY=RPAR(2) + ! compute the motion of two ends + R=RPAR(3) + THETA=RPAR(4) + DELX=CX+R*COS(THETA)-X1 + DELY=CY+R*SIN(THETA)-Y1 + MOV1=SQRT(DELX*DELX+DELY*DELY) + THETA=RPAR(5) + DELX=CX+R*COS(THETA)-X2 + DELY=CY+R*SIN(THETA)-Y2 + MOV2=SQRT(DELX*DELX+DELY*DELY) + ! compute new radius, new center and new angles + SX=X1-CX + SY=Y1-CY + ! get demi-vector from (x1,y1) to (x2,y2) + VX=(X2-X1)/2.0 + VY=(Y2-Y1)/2.0 + V2=VX*VX+VY*VY + ! get center closest to old center, change radius + ! compute new r + R=SQRT(V2+(SX*VY-SY*VX)**2/V2) + RPAR(3)=R + ! compute new center + DIST=1.0_PDB+(VX*SX+VY*SY)/V2 + ! get new center and compute angles + CX=CX+VX*DIST + CY=CY+VY*DIST + RPAR(1)=CX + RPAR(2)=CY + THETA=SALACO((X1-CX)/R,Y1-CY) + RPAR(4)=THETA + THETA=SALACO((X2-CX)/R,Y2-CY) + RPAR(5)=THETA + ! if phi1 > phi2 put phi2 equal to phi2+twopi + IF(RPAR(4)>RPAR(5)) RPAR(5)=RPAR(5)+TWOPI + ENDIF + IF(KNDEX/=0) THEN + IF(MOV1>EPS.OR.MOV2>EPS) THEN + WRITE(FOUT,'(/," old element ",I4," : ",6(1P,E13.4))') & + II,(RPAR_OLD(I),I=1,G_ELE_LEN(TYPE)) + WRITE(FOUT,'(" new element ",I4," : ",6(1P,E13.4))') & + II,(RPAR(I),I=1,G_ELE_LEN(TYPE)) + ENDIF + ENDIF + ! + END SUBROUTINE SAL129 + ! + SUBROUTINE SAL130(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! perform domain definition and set boundary conditions + ! - computes node perimeters per 2d macro + ! - reads boundary conditions + ! - computes external perimeters + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE,TYPGEO + !**** + IMPLICIT NONE + ! in variable + ! ************ + TYPE(T_G_BASIC), INTENT(INOUT):: GG + ! local variable + !**** + INTEGER :: NN,OK + INTEGER, DIMENSION(GG%NB_ELEM*2) :: AUX_ARR + !**** + !* compute node perimeters for the macro + CALL SAL130_2(GG%NB_ELEM,GG%NB_NODE,GG%IPAR,GG%PPERIM_NODE, & + GG%PERIM_NODE,AUX_ARR) + ! + !* - compute number of bc's per 2D macro, + ! NB_BC2 counts total nber of 2D bc's + ! - compute IBC2_ELEM, keep relative 2D bc nber to elements + ! - allocation : 2D bc structures + ! 2D perimeter structure for a macro + ! - get list of elements in 2d macro perimeter + + CALL SAL130_4(GG%NB_ELEM,NN,GG%IPAR,GG%IBC2_ELEM,AUX_ARR) + + GG%NB_BC2=NN + + ALLOCATE (GG%TYPE_BC2(NN),GG%IDATA_BC2(NN),GG%PERIM_MAC2(NN), & + GG%PPERIM_MAC2(7),STAT=OK) + IF(OK/=0) CALL XABORT('SAL130: not enough memory I,R') + GG%PPERIM_MAC2(:7)=0 + GG%PERIM_MAC2(1:NN)=AUX_ARR(1:NN) + GG%NPERIM_MAC2=NN + ! + !* read boundary conditions: + CALL SAL131(GG) + ! + ! - define IELEM_SURF2 + IF(GG%NB_SURF2>0) THEN + ! allocate + ALLOCATE(GG%SURF2(GG%NB_SURF2),STAT = OK) + IF(OK /= 0) CALL XABORT('SAL130: not enough memory I,R') + CALL SAL130_6(GG%NB_SURF2,GG%IBC2_SURF2,GG%PERIM_MAC2,GG%IELEM_SURF2) + ENDIF + ! + !* construct perimeter structures for rotative or symmetrical geometry + SELECT CASE(TYPGEO) + CASE(1:2) + CALL SAL130_8(GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%PPERIM_MAC2,GG%DIST_AXIS, & + GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2,GG%BCDATA,GG%IPAR,GG%RPAR) + CASE(5:12) + CALL SAL130_10(GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%PPERIM_MAC2,GG%DIST_AXIS, & + GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2,GG%BCDATA,GG%IPAR,GG%RPAR) + END SELECT + ! + END SUBROUTINE SAL130 + ! + SUBROUTINE SAL130_2(NB_ELEM,NB_NODE,IPAR,PPERIM,PERIM,LIST) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! compute node perimeters for one 2D macro + ! + !Parameters: input + ! NB_ELEM number of elements + ! NB_NODE number of nodes + ! IPAR integer descriptors for elements + ! + !Parameters: output + ! PPERIM array pointer to list of elements in perimeter per node + ! PERIM elements in perimeters per node + ! LIST temporary array + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + INTEGER, INTENT(IN) :: NB_ELEM,NB_NODE + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + INTEGER, INTENT(OUT), DIMENSION(:) :: LIST + INTEGER, INTENT(OUT), DIMENSION(:) :: PPERIM + INTEGER, POINTER, DIMENSION(:) :: PERIM + !**** + INTEGER :: NT,NN,NODE,ELEM,OK + CHARACTER(LEN=131) :: HSMG + !**** + NT=0 + PPERIM(1)=1 + DO NODE=1,NB_NODE + DO ELEM=1, NB_ELEM + IF(IPAR(2,ELEM)==NODE.OR.IPAR(3,ELEM)==NODE) THEN + NT=NT+1 + LIST(NT)=ELEM + ENDIF + ENDDO + NN=NT+1-PPERIM(NODE) + IF(NN>0) THEN + PPERIM(NODE+1)=NT+1 + ELSE + WRITE(HSMG,'(15HSAL130_2: node=,i5,19H without perimeter.)') NODE + CALL XABORT(HSMG) + ENDIF + ENDDO + IF(NT>0) THEN + ALLOCATE (PERIM(NT), STAT=OK) + IF(OK/=0) CALL XABORT('SAL130_2: not enough memory I') + PERIM(1:NT)=LIST(1:NT) + ENDIF + ! + END SUBROUTINE SAL130_2 + ! + SUBROUTINE SAL130_4(NB_ELEM,NB_BC,IPAR,IBC2_ELEM,LIST_BC) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! count number of bc's in a 2D macro + ! + !Parameters: input + ! NB_ELEM number of elements + ! + !Parameters: input/output + ! IPAR integer descriptors for elements + ! + !Parameters: output + ! NB_BC nber of bc's + ! IBC2_ELEM relative 2D bc index per element + ! LIST_BC list of elements in boundary + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + INTEGER, INTENT(IN) :: NB_ELEM + INTEGER, INTENT(OUT) :: NB_BC + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + INTEGER, INTENT(OUT), DIMENSION(:) :: IBC2_ELEM,LIST_BC + !**** + INTEGER :: ELEM + !**** + ! initiation + IBC2_ELEM(:)=0 + NB_BC=0 + DO ELEM=1, NB_ELEM + IF(IPAR(2,ELEM)<=0.AND.IPAR(3,ELEM)<=0) THEN + CALL XABORT('SAL130_4: two boundaries for element') + ELSE IF(IPAR(2,ELEM)<=0.OR.IPAR(3,ELEM)<=0) THEN + NB_BC=NB_BC+1 + LIST_BC(NB_BC)=ELEM + IF(IBC2_ELEM(ELEM)/=0) CALL XABORT('SAL130_4: two surfaces to element') + IBC2_ELEM(ELEM)=NB_BC + ! set mark for the macro connection surface : + ENDIF + ENDDO + ! + END SUBROUTINE SAL130_4 + ! + SUBROUTINE SAL130_6(NSURF,IBC2_SURF2,IELEM_BC2,IELEM_SURF2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! get element order indices per surface + ! + !Parameters: input + ! NSURF number of surfaces + ! IBC2_SURF2 2D bc order index per surface + ! IELEM_BC2 element order index per bc + ! + !Parameters: output + ! IELEM_SURF2 element order index per surface + ! + !--------------------------------------------------------------------- + ! + !**** + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NSURF + INTEGER, INTENT(IN), DIMENSION(:) :: IBC2_SURF2,IELEM_BC2 + INTEGER, INTENT(OUT), DIMENSION(:) :: IELEM_SURF2 + !**** + INTEGER :: SURF,IBC,ELEM + !**** + IF(NSURF > 0) THEN + DO SURF=1, NSURF + ! get relative bc order number + IBC=IBC2_SURF2(SURF) + ! get relative element + ELEM=IELEM_BC2(IBC) + ! define ielem_surf2 + IELEM_SURF2(SURF)=ELEM + ENDDO + ENDIF + ! + END SUBROUTINE SAL130_6 + ! + SUBROUTINE SAL130_8(NPERIM,PERIM,PPERIM,DIST_AXIS,IBC2_ELEM,TYPE_BC2,IDATA_BC2, & + BCDATA,IPAR,RPAR) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! calculate PERIM_MAC2,DIST_AXIS for rotative or symmetrical geometry + ! (TYPGEO=1 & 2) + ! + !Parameters: input + ! NPERIM number of elements in perimeter + ! IBC2_ELEM 2D bc order number per element + ! TYPE_BC2 type of boundary conditions per 2D bc + ! IDATA_BC2 position in the 'bcdata' table for each 2D bc + ! BCDATA table of boundary conditions data + ! IPAR integer element descriptor table + ! RPAR floating point element descriptor table + ! + !Parameters: input/output + ! PERIM list of elements in perimeter,in return it will be in following + ! order: (elems on axis 1)+(elems on axis 2)+(other elems) + ! + !Parameters: output + ! PPERIM pointers to the table 'perim': + ! (1): first elem on axis 1 (2): first elem on axis 2 + ! (3): first of other elems (4): NPERIM + 1 + ! DIST_AXIS distance of points on axis to the center (0,0),in order of: + ! (distances of points on axis 1)+(distances of points on axis 2) + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : ANGGEO,EPS,G_BC_TYPE + USE SAL_NUMERIC_MOD, ONLY : SAL141 + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NPERIM + INTEGER, INTENT(IN), DIMENSION(:) :: IBC2_ELEM,TYPE_BC2,IDATA_BC2 + INTEGER, INTENT(INOUT), DIMENSION(:) :: PERIM + INTEGER, INTENT(OUT), DIMENSION(:) :: PPERIM + REAL(PDB), POINTER, DIMENSION(:) :: DIST_AXIS + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + REAL(PDB), INTENT(INOUT),DIMENSION(:,:) :: RPAR + REAL(PDB), INTENT(IN),DIMENSION(:,:) :: BCDATA + !**** + ! LIST_ELEMS = table of elements,elements on axis are in increasing order + ! of distance to (0,0): + ! 1=elements on axis 1; 2=elements on axis 2; 3=other elements + ! AUX_DIST = max distance of element ends to the beginnings of the axes: + ! 1=distances on axis 1; 2=distances on axis 2 + INTEGER, DIMENSION(NPERIM,3) :: LIST_ELEMS + REAL(PDB), DIMENSION(NPERIM,2) :: AUX_DIST + INTEGER :: I,J,K,M,ELEM,TYPBC,IBC,IDATA,NBE(3),OK,NAXES + REAL(PDB) :: ANGLE,X,Y,D + !**** + NAXES=2 + !* calculate number of elements on axis 1 & 2 + ! and their distances to the beginning of the axes + NBE=0 + DO I=1,NPERIM + ELEM=PERIM(I) + IBC=IBC2_ELEM(ELEM) + TYPBC=TYPE_BC2(IBC) + IDATA=IDATA_BC2(IBC) + ! default is the elements not on axis + M=NAXES+1 + IF(TYPBC==G_BC_TYPE(3)) THEN + !* rotation: + ANGLE=BCDATA(5,IDATA) + IF(ABS(ANGLE-ANGGEO)<EPS) THEN + ! element is on the axis 1 (angle=anggeo) + M=1 + ELSEIF(ABS(ANGLE+ANGGEO)<EPS) THEN + ! element is on the axis 2 (angle=-anggeo) + M=2 + ENDIF + ELSEIF(TYPBC==G_BC_TYPE(4)) THEN + !* symmetry: + ANGLE=BCDATA(5,IDATA) + IF(ABS(ANGLE)<EPS) THEN + ! element is on the axis 1 (angle=0) + M=1 + ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN + ! element is on the axis 2 (angle=anggeo) + M=2 + ENDIF + ENDIF + NBE(M)=NBE(M)+1 + LIST_ELEMS(NBE(M),M)=ELEM + IF(M/=NAXES+1) THEN + ! sort the element list according their distance to (0,0) + D=0. + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT(X*X+Y*Y)) + ENDDO + AUX_DIST(NBE(M),M)=D + IF(NBE(M)>1) THEN + DO K=1,NBE(M)-1 + IF(AUX_DIST(NBE(M),M)<AUX_DIST(K,M)) THEN + ! insert the last element to position k + D=AUX_DIST(NBE(M),M) + ELEM=LIST_ELEMS(NBE(M),M) + DO J=NBE(M),K+1,-1 + AUX_DIST(J,M)=AUX_DIST(J-1,M) + LIST_ELEMS(J,M)=LIST_ELEMS(J-1,M) + ENDDO + AUX_DIST(K,M)=D + LIST_ELEMS(K,M)=ELEM + EXIT + ENDIF + ENDDO + ENDIF + ENDIF + ENDDO + ! nber of blocks in 'perim' + M=NAXES+1 + !* organize elements as: + ! (elems on axis 1)+(elems on axis 2)+(other elems) + PPERIM(:)=0 + PPERIM(1)=1 + DO I=1,M + PPERIM(I+1)=PPERIM(I)+NBE(I) + ENDDO + DO I=1,M + PERIM(PPERIM(I):PPERIM(I+1)-1)=LIST_ELEMS(1:NBE(I),I) + ENDDO + !* allocate and set dist_axis + IF(PPERIM(NAXES+1)-1>0) THEN + ALLOCATE(DIST_AXIS(PPERIM(NAXES+1)-1),STAT=OK) + IF(OK.NE.0) CALL XABORT('SAL130_8: not enough memory R') + DO I=1,NAXES + DIST_AXIS(PPERIM(I):PPERIM(I+1)-1)=AUX_DIST(1:NBE(I),I) + ENDDO + ELSE + NULLIFY(DIST_AXIS) + ENDIF + ! + end subroutine sal130_8 + ! + SUBROUTINE SAL130_10(NPERIM,PERIM,PPERIM,DIST_AXIS,IBC2_ELEM,TYPE_BC2,IDATA_BC2, & + BCDATA,IPAR,RPAR) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! calculate PERIM_MAC2 and DIST_AXIS for the cyclical geometries: + ! type 5&6: retangle with translations or symmetries on all sides + ! type 7 : 1/8 assembly with symmetries on all sides + ! type 8 : equilateral triangle with symmetries on all sides + ! type 9 : hexagon with translations on all sides + ! type 10 : hexagon with RA60 rotation + ! type 11 : hexagon with R120 rotation + ! type 12 : rectangular S30 triangle geometry + ! + !Parameters: input + ! NPERIM number of elements in perimeter + ! IBC2_ELEM 2D bc order number per element + ! TYPE_BC2 type of boundary conditions per 2D bc + ! IDATA_BC2 position in the 'bcdata' table for each 2D bc + ! BCDATA table of boundary conditions data + ! IPAR integer element descriptor table + ! RPAR floating point element descriptor table + ! + !Parameters: input/output + ! PERIM list of elements in perimeter,in return it will be in following + ! order: (elems on axis 1)+(elems on axis 2)+(other elems) + ! + !Parameters: output + ! PPERIM pointers to the table 'perim': + ! DIST_AXIS distance of points on axis to the axis origin + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE + USE SAL_NUMERIC_MOD, ONLY : SAL141 + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NPERIM + INTEGER, INTENT(IN), DIMENSION(:) :: IBC2_ELEM,TYPE_BC2,IDATA_BC2 + INTEGER, INTENT(INOUT), DIMENSION(:) :: PERIM + INTEGER, INTENT(OUT), DIMENSION(:) :: PPERIM + REAL(PDB), POINTER, DIMENSION(:) :: DIST_AXIS + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + REAL(PDB), INTENT(INOUT),DIMENSION(:,:) :: RPAR + REAL(PDB), INTENT(IN),DIMENSION(:,:) :: BCDATA + !**** + ! LIST_ELEMS = table of elements,elements on axis are in increasing order + ! of distance to axis origin: + ! 1=elements on axis 1; 2=elements on axis 2; ... + ! AUX_DIST = max distance of element ends to axis origin: + ! 1=distances on axis 1; 2=distances on axis 2; ... + INTEGER, DIMENSION(NPERIM,6) :: LIST_ELEMS + REAL(PDB), DIMENSION(NPERIM,6) :: AUX_DIST + INTEGER :: I,J,K,M,ELEM,TYPBC,IBC,IDATA,NBE(6),OK,NAXES + REAL(PDB) :: ANGLE,X,Y,D + !**** + NAXES=0 + IF((TYPGEO==5).OR.(TYPGEO==6).OR.(TYPGEO==11)) THEN + NAXES=4 + ELSEIF((TYPGEO==7).OR.(TYPGEO==8).OR.(TYPGEO==10).OR.(TYPGEO==12)) THEN + NAXES=3 + ELSEIF(TYPGEO==9) THEN + NAXES=6 + ENDIF + !* calculate number of elements on axes, + ! and their distances to the origin of the axis + NBE=0 + + DO I=1,NPERIM + ELEM=PERIM(I) + IBC=IBC2_ELEM(ELEM) + TYPBC=TYPE_BC2(IBC) + IDATA=IDATA_BC2(IBC) + ENDDO + DO I=1,NPERIM + ELEM=PERIM(I) + IBC=IBC2_ELEM(ELEM) + TYPBC=TYPE_BC2(IBC) + IDATA=IDATA_BC2(IBC) + ANGLE=BCDATA(5,IDATA) + ! default is the elements not on axis + M=NAXES+1 + SELECT CASE(TYPGEO) + CASE(5) + !* rectangle with translations: + ! axes definition: + ! M=3 + ! ************ + ! * * + ! M=2 * * M=4 + ! ************ + ! M=1 + IF(TYPBC==G_BC_TYPE(2)) THEN + IF(BCDATA(2,IDATA)>0) THEN + M=1 + ELSEIF(BCDATA(2,IDATA)<0) THEN + M=3 + ELSEIF(BCDATA(1,IDATA)>0) THEN + M=2 + ELSEIF(BCDATA(1,IDATA)<0) THEN + M=4 + ENDIF + ELSE + CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=5.') + ENDIF + CASE(6) + !* rectangle with symmetries: + ! axes definition: + ! M=3 + ! ************ + ! * * + ! M=2 * * M=4 + ! ************ + ! M=1 + IF(TYPBC==G_BC_TYPE(4)) THEN + IF(ABS(ANGLE)<EPS) THEN + IF(BCDATA(2,IDATA)>0) THEN + ! cy>0: element is on the axis 3 (angle=0) + M=3 + ELSE + ! cy=0: element is on the axis 1 (angle=0) + M=1 + ENDIF + ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN + IF(BCDATA(1,IDATA)>0) THEN + ! cx>0: element is on the axis 4 (angle=anggeo) + M=4 + ELSE + ! cx=0: element is on the axis 2 (angle=anggeo) + M=2 + ENDIF + ENDIF + ELSE + CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=6.') + ENDIF + CASE(7,8,12) + !* triangles with symmetries: + ! typgeo=7 axes definition: + ! * + ! M=2 * * M=3 + ! * * + ! ******* + ! M=1 + ! typgeo=8 axes definition: + ! * + ! M=2 * * M=3 + ! * * + ! ************* + ! M=1 + IF(TYPBC==G_BC_TYPE(4)) THEN + IF(ABS(ANGLE)<EPS) THEN + ! element is on the axis 1 (angle=0) + M=1 + ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN + ! element is on the axis 2 (angle=anggeo) + M=2 + ELSEIF(ABS(ANGLE-HALFPI)<EPS) THEN + ! typgeo=7:element is on the axis 3 (angle=pi/2) + M=3 + ELSEIF(ABS(ANGLE-2.*ANGGEO)<EPS) THEN + ! typgeo=8:element is on the axis 3 (angle=2*anggeo) + M=3 + ELSEIF(ABS(ANGLE-2.*PI/3.)<EPS) THEN + ! typgeo=12:element is on the axis 3 (angle=2*pi/3) + M=3 + ELSE + CALL XABORT('SAL130_10: boundary condition data error in element for TYPGEO=7,8,12.') + ENDIF + ELSE + CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=7,8.') + ENDIF + CASE(9) + !* hexagon with translations: + ! axes definition: + ! M=4 (0,-2ly) + ! ***** + ! M=3 (3/2lx,-ly) * * M=5 (-3/2lx,-ly) + ! * * + ! M=2 (3/2lx, ly) * * M=6 (-3/2lx, ly) + ! ***** + ! M=1 (0, 2ly) + ! origins of axes: + ! axis 1&2: (-lx/2,-ly) + ! axis 3: (-lx , 0) + ! axis 4: (-lx/2, ly) + ! axis 5: ( lx , 0) + ! axis 6: ( lx/2,-ly) + IF(TYPBC==G_BC_TYPE(2)) THEN + IF(ABS(BCDATA(1,IDATA))<EPS) THEN + IF(BCDATA(2,IDATA)>0) THEN + M=1 + ELSE + M=4 + ENDIF + ELSEIF(BCDATA(1,IDATA)>0.AND.BCDATA(2,IDATA)>0) THEN + M=2 + ELSEIF(BCDATA(1,IDATA)>0.AND.BCDATA(2,IDATA)<0) THEN + M=3 + ELSEIF(BCDATA(1,IDATA)<0.AND.BCDATA(2,IDATA)>0) THEN + M=6 + ELSEIF(BCDATA(1,IDATA)<0.AND.BCDATA(2,IDATA)<0) THEN + M=5 + ENDIF + ELSE + CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=9.') + ENDIF + CASE(10) + !* triangles with rotations and translations: + ! typgeo=10 axes definition: + ! * + ! M=2 * * M=3 + ! * * + ! ************* + ! M=1 + IF((TYPBC==G_BC_TYPE(2)).OR.(TYPBC==G_BC_TYPE(3))) THEN + IF(ABS(ANGLE)<EPS) THEN + ! element is on the axis 1 (angle=0) + M=1 + ELSEIF(ABS(ANGLE-ANGGEO)<EPS) THEN + ! element is on the axis 2 (angle=anggeo) + M=2 + ELSEIF(ABS(ANGLE-2.*ANGGEO)<EPS) THEN + ! element is on the axis 3 (angle=2*anggeo) + M=3 + ELSE + CALL XABORT('SAL130_10: boundary condition data error in element for TYPGEO=10.') + ENDIF + ELSE + CALL XABORT('SAL130_10: wrong boundary condition for TYPGEO=10.') + ENDIF + CASE(11) + !* lozenge with rotations and translations: + ! axes definition: + ! M=3 + ! ************ + ! * * + ! M=2 * * M=4 + ! ************ + ! M=1 + IF((TYPBC==G_BC_TYPE(2)).OR.(TYPBC==G_BC_TYPE(3))) THEN + IF(ABS(ANGLE)<EPS) THEN + IF(BCDATA(2,IDATA)>0) THEN + ! cy>0: element is on the axis 3 (angle=0) + M=3 + ELSE + ! cy=0: element is on the axis 1 (angle=0) + M=1 + ENDIF + ELSE + IF(BCDATA(1,IDATA)>0) THEN + ! cx>0: element is on the axis 4 (angle>0) + M=4 + ELSE + ! cx=0: element is on the axis 2 (angle>0) + M=2 + ENDIF + ENDIF + ELSE + CALL XABORT('SAL130_10: TYPGEO=11:wrong boundary condition type.') + ENDIF + CASE DEFAULT + CALL XABORT('SAL130_10: boundary condition not implemented.') + END SELECT + IF(M>=NAXES+1) CALL XABORT('SAL130_10: element not on axes') + NBE(M)=NBE(M)+1 + LIST_ELEMS(NBE(M),M)=ELEM + ! sort the element list according their distance to + ! the origins of the axes + D=0. + SELECT CASE(TYPGEO) + CASE(5:6) + SELECT CASE(M) + CASE(1:2) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT(X*X+Y*Y)) + ENDDO + CASE(3) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,X) + ENDDO + CASE(4) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,Y) + ENDDO + END SELECT + CASE(7) + SELECT CASE(M) + CASE(1:2) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT(X*X+Y*Y)) + ENDDO + CASE(3) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,Y) + ENDDO + END SELECT + CASE(8,10,12) + SELECT CASE(M) + CASE(1:2) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT(X*X+Y*Y)) + ENDDO + CASE(3) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + ! axis origin is (lx,0) + D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y)) + ENDDO + END SELECT + CASE(9) + ! origins of axes: + ! axis 1&2: (-lx/2,-ly) + ! axis 3: (-lx , 0) + ! axis 4: (-lx/2, ly) + ! axis 5: ( lx , 0) + ! axis 6: ( lx/2,-ly) + SELECT CASE(M) + CASE(1:2) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X+LX*0.5)*(X+LX*0.5)+(Y+LY)*(Y+LY))) + ENDDO + CASE(3) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X+LX)*(X+LX)+Y*Y)) + ENDDO + CASE(4) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X+LX*0.5)*(X+LX*0.5)+(Y-LY)*(Y-LY))) + ENDDO + CASE(5) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y)) + ENDDO + CASE(6) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X-LX*0.5)*(X-LX*0.5)+(Y+LY)*(Y+LY))) + ENDDO + END SELECT + CASE(11) + SELECT CASE(M) + CASE(1:2) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT(X*X+Y*Y)) + ENDDO + CASE(3) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X-LX*0.5)*(X-LX*0.5)+(Y-LY)*(Y-LY))) + ENDDO + CASE(4) + DO K=1,2 + CALL SAL141(IPAR(1,ELEM),RPAR(:,ELEM),X,Y,K) + D=MAX(D,SQRT((X-LX)*(X-LX)+Y*Y)) + ENDDO + END SELECT + END SELECT + AUX_DIST(NBE(M),M)=D + IF(NBE(M)>1) THEN + DO K=1,NBE(M)-1 + IF(AUX_DIST(NBE(M),M)<AUX_DIST(K,M)) THEN + ! insert the last element to position k + D=AUX_DIST(NBE(M),M) + ELEM=LIST_ELEMS(NBE(M),M) + DO J=NBE(M),K+1,-1 + AUX_DIST(J,M)=AUX_DIST(J-1,M) + LIST_ELEMS(J,M)=LIST_ELEMS(J-1,M) + ENDDO + AUX_DIST(K,M)=D + LIST_ELEMS(K,M)=ELEM + EXIT + ENDIF + ENDDO + ENDIF + ENDDO + ! nber of blocks in 'perim' + M=NAXES + !* organize elements as: (elems on axis 1)+(elems on axis 2)+... + PPERIM(1)=1 + DO I=1,M + PPERIM(I+1)=PPERIM(I)+NBE(I) + ENDDO + DO I=1,M + PERIM(PPERIM(I):PPERIM(I+1)-1)=LIST_ELEMS(1:NBE(I),I) + ENDDO + !* allocate and set dist_axis + IF(PPERIM(NAXES+1)-1>0) THEN + ALLOCATE(DIST_AXIS(PPERIM(NAXES+1)-1),STAT=OK) + IF(OK.NE.0) CALL XABORT('SAL130_10: not enough memory R') + DO I=1,NAXES + DIST_AXIS(PPERIM(I):PPERIM(I+1)-1)=AUX_DIST(1:NBE(I),I) + ENDDO + ELSE + NULLIFY(DIST_AXIS) + ENDIF + ! + END SUBROUTINE SAL130_10 + ! + SUBROUTINE SAL131(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! modifies GG%IPAR and constructs bcdata (modifies GG%IPAR and + ! constructs BCDATA) + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_MAX_LEN,G_BC_LEN,G_BC_TYPE,NIPAR,EPS, & + ANGGEO,TYPGEO,LENGTHX,LENGTHY,ALLSUR + !**** + IMPLICIT NONE + ! in variable + ! ************ + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + !**** + LOGICAL :: LGBC,LGALLS + INTEGER :: TYPE,II,NN,NBER,I,J,ITBC,IELEM + REAL(PDB) :: ANGLE + INTEGER :: ELEM,IBC,OK,IDATA(2) + INTEGER, POINTER, DIMENSION(:) :: TYPE_BC2,ISURF2_ELEM,IELEM_BC2 + !***** + !** TMP_BCDATA = description of motions at the boundary + !** IPAR = pointer to geometry descriptors + INTEGER, DIMENSION(GG%NB_ELEM) :: AUX_ARR + REAL(PDB), DIMENSION(G_BC_MAX_LEN,GG%NB_ELEM+5) :: TMP_BCDATA + !***** + !* initialization + IF(TYPGEO==1.OR.TYPGEO==2) IDATA(:)=0 + !* BCDATA for surfaces of type G_BC_TYPE(-1) + ! + LGALLS=ALLSUR/=0 + ITBC=0 ! the first bc data + ! + !* treat approximate boundary condictions + IF(LGALLS)THEN + ITBC=1 + ! ALL BC'S PRODUIT SURFACES + IF(GG%DEFAUL==1)THEN + ! SPECULAR REFLEXION -> ISOTROPIC REFLEXION WITH ALBEDO=1 + GG%DEFAUL=0 ; TMP_BCDATA(6,ITBC)=1._PDB + ENDIF + ENDIF + ! + !* put default value in all bc elements: + CALL SAL131_2(GG%NB_ELEM,GG%DEFAUL,GG%IPAR,GG%IBC2_ELEM,GG%TYPE_BC2,GG%IDATA_BC2) + ! + IF(GG%NBBCDA>0)THEN + ANGLE=0._PDB + DO I=1,GG%NBBCDA + ITBC=ITBC+1 + IF(ITBC.GT.GG%NB_ELEM+5) CALL XABORT('SAL131: BCDATA overflow') + TMP_BCDATA(:,ITBC)=GG%BCDATAREAD(I)%BCDATA(:) + TYPE=GG%BCDATAREAD(I)%SALTYPE + SELECT CASE(TYPE) + CASE(1) + IF(LGALLS)THEN + ! specular reflexion -> isotropic reflexion with albedo=1 + TYPE=0 ; TMP_BCDATA(6,ITBC)=1._PDB + ENDIF + CASE(2) + ANGLE=TMP_BCDATA(5,ITBC) + SELECT CASE(TYPGEO) + CASE(5) + ! adjust translation data according to the sides of rectangle + IF(TMP_BCDATA(1,ITBC)/=0) TMP_BCDATA(1,ITBC)=SIGN(LENGTHX,TMP_BCDATA(1,ITBC)) + IF(TMP_BCDATA(2,ITBC)/=0) TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC)) + CASE(9) + ! adjust translation data according to the hexagonal sides + IF(ABS(TMP_BCDATA(1,ITBC))<EPS) THEN + ! axes 1 & 4 + TMP_BCDATA(2,ITBC)=SIGN(2.*LENGTHY,TMP_BCDATA(2,ITBC)) + ELSE + TMP_BCDATA(1,ITBC)=SIGN(LENGTHX*1.5,TMP_BCDATA(1,ITBC)) + TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC)) + ENDIF + TMP_BCDATA(3,ITBC)=COS(ANGLE) + TMP_BCDATA(4,ITBC)=SIN(ANGLE) + TMP_BCDATA(5,ITBC)=ANGLE + CASE(11) + ! adjust translation data according to the hexagonal sides + IF(ABS(ANGLE)<EPS) THEN + IF(TMP_BCDATA(2,ITBC)>0.) THEN + ! cy>0 and angle=0 + TMP_BCDATA(1,ITBC)=LENGTHX/2.0 + TMP_BCDATA(2,ITBC)=LENGTHY + ELSE + ! cy=0 and angle=0 + TMP_BCDATA(1,ITBC)=0.0 + TMP_BCDATA(2,ITBC)=0.0 + ENDIF + ELSE + IF(TMP_BCDATA(1,ITBC)>0.) THEN + ! cx>0 and angle=pi/3 rad + TMP_BCDATA(1,ITBC)=LENGTHX + TMP_BCDATA(2,ITBC)=0.0 + ELSE + ! cx=0 and angle=pi/3 rad + TMP_BCDATA(1,ITBC)=0.0 + TMP_BCDATA(2,ITBC)=0.0 + ENDIF + ENDIF + TMP_BCDATA(3,ITBC)=COS(ANGLE) + TMP_BCDATA(4,ITBC)=SIN(ANGLE) + TMP_BCDATA(5,ITBC)=ANGLE + END SELECT + IF(LGALLS)THEN + ! translation -> isotropic translation with albedo=1 + TYPE=TYPE+4 + ! perform array shift + TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1) + ! albedo=1 + TMP_BCDATA(6,ITBC)=1._PDB + ENDIF + CASE(3) + ! cases of rotation: + ! read angle, compute cos and sin + ANGLE=TMP_BCDATA(5,ITBC) + IF(TYPGEO==1.OR.TYPGEO==2) THEN + IF(ABS(ANGLE-ANGGEO)<EPS) THEN + ! for the axis 1:keep anggeo,cos(anggeo),sin(anggeo) + ANGLE=ANGGEO + IF(IDATA(1)==0) IDATA(1)=ITBC + ELSEIF(ABS(ANGLE+ANGGEO)<EPS) THEN + ! for the axis 2:keep -anggeo,cos(-anggeo),sin(-anggeo) + ANGLE=-ANGGEO + IF(IDATA(2)==0) IDATA(2)=ITBC + ELSE + CALL XABORT('SAL131: error in angle of rotative axis') + ENDIF + ENDIF + TMP_BCDATA(3,ITBC)=COS(ANGLE) + TMP_BCDATA(4,ITBC)=SIN(ANGLE) + TMP_BCDATA(5,ITBC)=ANGLE + IF(LGALLS)THEN + ! rotation -> isotropic rotation with albedo=1 + ! axial symmetry -> isotropic axial symmetry with albedo=1 + TYPE=TYPE+4 + ! perform array shift + TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1) + ! albedo=1 + TMP_BCDATA(6,ITBC)=1._PDB + ENDIF + CASE(4) + ! cases of specular symmetry: + ! read center and angle, compute cos and sin + ANGLE=TMP_BCDATA(5,ITBC) + SELECT CASE(TYPGEO) + CASE(1,2) + IF(ABS(ANGLE-0.)<EPS) THEN + IF(IDATA(1)==0) IDATA(1)=ITBC + ELSE + IF(IDATA(2)==0) IDATA(2)=ITBC + ENDIF + CASE(6:8) + ! adjust translation data according to the rectangle/triangle sides + IF(TMP_BCDATA(1,ITBC)/=0) TMP_BCDATA(1,ITBC)=SIGN(LENGTHX,TMP_BCDATA(1,ITBC)) + IF(TMP_BCDATA(2,ITBC)/=0) TMP_BCDATA(2,ITBC)=SIGN(LENGTHY,TMP_BCDATA(2,ITBC)) + END SELECT + TMP_BCDATA(3,ITBC)=COS(ANGLE) + TMP_BCDATA(4,ITBC)=SIN(ANGLE) + TMP_BCDATA(5,ITBC)=ANGLE + IF(LGALLS)THEN + ! rotation -> isotropic rotation with albedo=1 + ! axial symmetry -> isotropic axial symmetry with albedo=1 + TYPE=TYPE+4 + ! perform array shift + TMP_BCDATA(:,ITBC)=CSHIFT(TMP_BCDATA(:,ITBC),1) + ! albedo=1 + TMP_BCDATA(6,ITBC)=1._PDB + ENDIF + END SELECT + ! + ! modify notation for boundary conditions + NBER=GG%BCDATAREAD(I)%NBER + DO J=1,NBER + ELEM=GG%BCDATAREAD(I)%ELEMNB(J) + IF(ELEM>GG%NB_ELEM.OR.ELEM<=0) CALL XABORT('SAL131: unknown bc element') + ! get local surface nber + IBC=GG%IBC2_ELEM(ELEM) + LGBC=GG%IPAR(2,ELEM)<=0 + II=0 + IF(LGBC)THEN + II=2 + ELSE + LGBC=GG%IPAR(3,ELEM)<=0 + IF(LGBC) II=3 + ENDIF + IF(.NOT.LGBC) THEN + WRITE(*,*) 'elem :',ELEM + WRITE(*,*) 'GG%IPAR(:,ELEM) :',GG%IPAR(:,ELEM) + CALL XABORT('SAL131: wrong bc element') + ENDIF + ! put bc type + GG%IPAR(II,ELEM)=G_BC_TYPE(TYPE) + GG%TYPE_BC2(IBC)=G_BC_TYPE(TYPE) + ! put bc data position : + GG%IDATA_BC2(IBC)=ITBC + ENDDO + ENDDO + ENDIF + ! + !* - set BCDATA position for surfaces of type G_BC_TYPE(-1) + ! - compute the nber of surfaces (type -1,0,-12,-13,-14,-15) : nbsur2 + ! - allocate structures for the surfaces + ! - compute surf_mac2 + TYPE_BC2=>GG%TYPE_BC2 + ISURF2_ELEM=>GG%ISURF2_ELEM + IELEM_BC2=>GG%PERIM_MAC2 + NN=0 + DO IBC=1,GG%NB_BC2 + IF(TYPE_BC2(IBC)==G_BC_TYPE(-1)) THEN + ! macro contact surfaces : set bcdata position to 1 + GG%IDATA_BC2(IBC)=1 + ENDIF + ! relative element nber + IELEM=IELEM_BC2(IBC) + ! count 2D surfaces number + IF(TYPE_BC2(IBC)==G_BC_TYPE(-1) .OR. TYPE_BC2(IBC)==G_BC_TYPE(0) .OR. & + TYPE_BC2(IBC)==G_BC_TYPE(1)) THEN + NN=NN+1 + AUX_ARR(NN)=IBC + ISURF2_ELEM(IELEM)=NN + ELSE + ISURF2_ELEM(IELEM)=0 + ENDIF + ENDDO + GG%NB_SURF2=NN + IF(NN>0) THEN + ALLOCATE (GG%IBC2_SURF2(NN),GG%IELEM_SURF2(NN),STAT=OK) + IF(OK/=0) CALL XABORT('SAL131: NOT ENOUGH MEMORY I,R') + GG%IBC2_SURF2(1:NN)=AUX_ARR(1:NN) + ELSE + NULLIFY(GG%IBC2_SURF2,GG%IELEM_SURF2) + ENDIF + ! + !* allocate idata_axis + IF(TYPGEO==1.OR.TYPGEO==2) THEN + DO I=1,2 + IF(IDATA(I)==0) THEN + ! there is no elements on this axis,add a bcdata for this axis + ITBC=ITBC+1 + ANGLE=0._PDB + IF(TYPGEO==1) THEN + ! symmetry + IF(I==1) THEN + ANGLE=0._PDB + ELSE + ANGLE=ANGGEO + ENDIF + ELSEIF(TYPGEO==2) THEN + ! rotation + IF(I==1) THEN + ANGLE=ANGGEO + ELSE + ANGLE=-ANGGEO + ENDIF + ENDIF + TMP_BCDATA(1,ITBC)=0._PDB + TMP_BCDATA(2,ITBC)=0._PDB + TMP_BCDATA(3,ITBC)=COS(ANGLE) + TMP_BCDATA(4,ITBC)=SIN(ANGLE) + TMP_BCDATA(5,ITBC)=ANGLE + IDATA(I)=ITBC + ENDIF + ENDDO + ALLOCATE(GG%IDATA_AXIS(2), STAT=OK) + IF(OK.NE.0) CALL XABORT('SAL131: not enough memory I') + GG%IDATA_AXIS=IDATA + ENDIF + !* allocate bcdata + ALLOCATE (GG%BCDATA(G_BC_MAX_LEN,ITBC), STAT=OK) + IF(OK/=0) CALL XABORT('SAL131: not enough memory R') + GG%BCDATA(:,1:ITBC)=TMP_BCDATA(:,1:ITBC) + GG%NALBG=ITBC + ! + END SUBROUTINE SAL131 + ! + SUBROUTINE SAL131_2(NB_ELEM,DEFAUL,IPAR,IBC2_ELEM,TYPE_BC2,IDATA_BC2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! put default bc type and BCDATA to a 2D macro + ! + !Parameters: input + ! NB_ELEM nber of elements + ! DEFAUL default bc type + ! IBC2_ELEM relative 2D bc number per element + ! + !Parameters: input/output + ! IPAR integer descriptors for elements + ! + !Parameters: input + ! TYPE_BC2 2D boundary condiction type + ! IDATA_BC2 position in bcdata par 2D bc + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NB_ELEM,DEFAUL + INTEGER, INTENT(INOUT), DIMENSION(:,:) :: IPAR + INTEGER, INTENT(IN), DIMENSION(:) :: IBC2_ELEM + INTEGER, INTENT(OUT), DIMENSION(:) :: TYPE_BC2 + INTEGER, INTENT(OUT), DIMENSION(:) :: IDATA_BC2 + !**** + INTEGER :: ELEM,II,IBC + LOGICAL :: LGBC + !**** + ! initiation + TYPE_BC2=G_BC_TYPE(-1) + ! + DO ELEM=1,NB_ELEM + LGBC=IPAR(2,ELEM)==0 + II=0 + IF(LGBC) THEN + II=2 + IF(IPAR(3,ELEM)<=0) CALL XABORT('SAL131_2: element with 2 bc''s') + ELSE + LGBC=IPAR(3,ELEM)==0 + IF(LGBC) II=3 + ENDIF + IF(LGBC) THEN + ! put bc type value to ipar + IPAR(II,ELEM)=G_BC_TYPE(DEFAUL) + ! put bc type value to bc type structure + IBC=IBC2_ELEM(ELEM) + IF(IBC==0) CALL XABORT('SAL131_2: surf-element relation error') + IF(TYPE_BC2(IBC)/=G_BC_TYPE(-1)) CALL XABORT('SAL131_2: two elements to a surface') + TYPE_BC2(IBC)=G_BC_TYPE(DEFAUL) + ! put position of "defaul" in bcdata table (use default albedo) + IDATA_BC2(IBC)=0 + ENDIF + ENDDO + ! + END SUBROUTINE SAL131_2 + ! + SUBROUTINE SAL140(NB_NODE,RPAR,IPAR,PPERIM,PERIM) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! checks domain topology in a 2D macro + ! + !Parameters: input + ! NB_NODE number of nodes in macro + ! + !Parameters: input/output + ! RPAR floating point geometry descriptors + ! IPAR integer descriptors for elements + ! PPERIM pointer to list of elements in perimeter + ! PERIM list of perimeters + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : NRPAR,NIPAR + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NB_NODE + INTEGER, INTENT(IN), DIMENSION (:) :: PPERIM + INTEGER, INTENT(INOUT), DIMENSION (:,:) :: IPAR + INTEGER, INTENT(INOUT), DIMENSION (:) :: PERIM + REAL(PDB), INTENT(INOUT), DIMENSION (:,:) :: RPAR + !**** + INTEGER :: ELEM,NODE,ID,L,FIRST,LAST,NEXT,NLOOP, & + NEWEND,KEEP,IEND,I1,I2 + REAL(PDB) :: X,Y,XNEW,YNEW,DIST + LOGICAL :: LGOPEN,LGON,LGERR + INTEGER, PARAMETER :: MXKEEP = 20 + REAL, DIMENSION(MXKEEP) :: KEEPD + INTEGER, PARAMETER :: FOUT =6 + REAL(PDB),PARAMETER :: EPS2=1.0E-7_PDB + !***** + NEXT=0 + !* checks topology of each node + DO NODE=1,NB_NODE + I1=PPERIM(NODE) + ! ID counts elements in the perimeter already treated + ID=I1-1 + I2=PPERIM(NODE+1)-1 + IF(I2<I1) CALL XABORT('SAL140: node without perimeter') + ! mark elements in perimeter of node as untreated + ! (except full circle: ipar=2) + DO L=I1,I2 + ELEM=PERIM(L) + IF(IPAR(1,ELEM)==2)THEN + ID=ID+1 + IF(ID<L)THEN + ! move circles at beginning + PERIM(L)=PERIM(ID) + PERIM(ID)=ELEM + ENDIF + ENDIF + ENDDO + ID=ID+1 + ! nloop counts nber of loops for node + NLOOP=ID-1 + ! treat all elements in perimeter until they have been all done: + DO WHILE (ID<I2) + ! an element in the perimeter that has not been done + FIRST=PERIM(ID) + NLOOP=NLOOP+1 + ! get first end for this element: + IEND=1 + ! determine the loop defined by element last: + LGOPEN=.TRUE. + LAST=FIRST + DO WHILE(LGOPEN) + ! get coordinates x and y of end of element last + CALL SAL141(IPAR(1,LAST),RPAR(:,LAST),X,Y,IEND) + ! get next element in the perimeter in touch with + ! last element at its end iend + LGON=.TRUE. + L=ID+1 + KEEP=0 + DO WHILE(LGON.AND.L<=I2) + NEXT=PERIM(L) + ! check whether next is really 'next' to last + CALL SAL142(X,Y,XNEW,YNEW,IPAR(1,NEXT),RPAR(:,NEXT),NEWEND,EPS2,DIST) + IF(KEEP<MXKEEP)KEEP=KEEP+1 + KEEPD(KEEP)=REAL(DIST) + LGON=NEWEND<0 + IF(LGON)L=L+1 + ENDDO + LGOPEN=.NOT.LGON + IF(LGOPEN)THEN + ! replace last element in order in the perimeter + ID=ID+1 + IF(ID<L)THEN + PERIM(L)=PERIM(ID) + PERIM(ID)=NEXT + ENDIF + ELSE + LGERR=.FALSE. + ! try to close with first element of loop: + CALL SAL142(X,Y,XNEW,YNEW,IPAR(1,FIRST),RPAR(:,FIRST),NEWEND,EPS2,DIST) + LGOPEN=NEWEND<0 + IF(LGOPEN)THEN + ! fatal error: cannot close loop for this node + WRITE(FOUT,'(//,''==> cannot close node '',i5, & + &'' AT ELEMENT'',I5,/)')NODE,LAST + LGERR=.TRUE. + ELSEIF(FIRST==LAST)THEN + WRITE(FOUT,'(//,''==> node '',i5,'' with '', & + &''isolated element'',i5,/)')NODE,LAST + LGERR=.TRUE. + ENDIF + IF(LGERR) CALL XABORT('SAL140: node not closed') + NEXT=FIRST + ENDIF + ! define last = next element and proceed + LAST=NEXT + IEND=3-NEWEND + ENDDO + ID=ID+1 + ENDDO + ENDDO + ! + END SUBROUTINE SAL140 + ! + SUBROUTINE SAL142(X,Y,XNEW,YNEW,TYPE,RPAR,IEND,EPS2,DIST) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! checks whether element is very close to point (X,Y) + ! + !Parameters: input + ! X abscissa coordinate + ! Y ordinate coordinate + ! TYPE type of element 1 (segment), 3 (arc of circle) + ! EPS2 criterium for closeness + ! + !Parameters: input/output + ! XNEW abscissa coordinate of end of element close to point + ! YNEW ordinate coordinate of end of element close to point + ! RPAR floating point geometry descriptors + ! IEND = 1 (beginning is close to point dist < EPS2) + ! 2 (end is close to point DIST < EPS2) + ! -1 (beginning is close to point dist > EPS2) + ! -2 (end is close to point DIST > EPS2) + ! DIST distance from end of element to point (X,Y) + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_ELE_TYPE + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: TYPE + INTEGER, INTENT(OUT) :: IEND + REAL(PDB), INTENT(IN) :: X,Y,EPS2 + REAL(PDB), INTENT(OUT) :: XNEW,YNEW,DIST + REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR + ! DIMENSION RPAR(NRPAR) + !**** + REAL(PDB) :: CX,CY,THETA,R,DIST2 + !**** + !* function giving distance between two points: + REAL(PDB) :: FUNC,X1,Y1,X2,Y2 + FUNC(X1,Y1,X2,Y2)=(X1-X2)**2+(Y1-Y2)**2 + !**** + CX=RPAR(1) + CY=RPAR(2) + DIST2=0._PDB + IF(TYPE==G_ELE_TYPE(1))THEN + !* segment + XNEW=CX + YNEW=CY + DIST=FUNC(XNEW,YNEW,X,Y) + IF(DIST<=EPS2)THEN + IEND=1 + RETURN + ELSE + XNEW=CX+RPAR(3) + YNEW=CY+RPAR(4) + DIST2=FUNC(XNEW,YNEW,X,Y) + ENDIF + ELSEIF(TYPE<=G_ELE_TYPE(3))THEN + !* arc of circle + R=RPAR(3) + THETA=RPAR(4) + XNEW=CX+R*COS(THETA) + YNEW=CY+R*SIN(THETA) + DIST=FUNC(XNEW,YNEW,X,Y) + IF(DIST<=EPS2)THEN + IEND=1 + RETURN + ELSE + THETA=RPAR(5) + XNEW=CX+R*COS(THETA) + YNEW=CY+R*SIN(THETA) + DIST2=FUNC(XNEW,YNEW,X,Y) + ENDIF + ELSE + CALL XABORT('SAL142: not implemented') + ENDIF + ! + IF(DIST2<=EPS2)THEN + DIST=DIST2 + IEND=2 + RETURN + ELSE + IF(DIST<=DIST2)THEN + IEND=-1 + ELSE + DIST=DIST2 + IEND=-2 + ENDIF + ENDIF + ! + END SUBROUTINE SAL142 + ! + SUBROUTINE SAL160(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! analyse domain definition: 2D volumes, surfaces + ! - compute node volumes + ! - compute areas of 2d surfaces + ! - read medium data + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE,NBMED + !**** + IMPLICIT NONE + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + !**** + INTEGER :: INB,IMED + INTEGER, DIMENSION(GG%NB_NODE) :: DATAIN + !**** + ! SUBROUTINE SAL160_2(NB_ELEM,IPAR,RPAR,VOL2,ISURF2_ELEM,NB_SURF2,SURF2) + IF(GG%NB_SURF2==0) THEN + CALL SAL160_2(GG%NB_ELEM,GG%IPAR,GG%RPAR,GG%VOL_NODE,GG%ISURF2_ELEM, & + GG%NB_SURF2) + ELSE + CALL SAL160_2(GG%NB_ELEM,GG%IPAR,GG%RPAR,GG%VOL_NODE,GG%ISURF2_ELEM, & + GG%NB_SURF2,GG%SURF2) + ENDIF + ! + !* read medium per region + CALL SALGET(DATAIN,GG%NB_NODE,F_GEO,FOUT0,'media per node') + + ! number of media fixed to maximum of datain + NBMED = MAXVAL(DATAIN(1:GG%NB_NODE)) + ! + !* check and define med for code regions + DO INB=1,GG%NB_NODE + IMED=DATAIN(INB) + IF(IMED>NBMED.OR.IMED<0)THEN + WRITE(*,*) 'medium : ',IMED + WRITE(*,*) 'inb, nbmed, nbnode : ',INB,NBMED,GG%NB_NODE + CALL XABORT('SAL160: wrong medium in a region') + ENDIF + GG%MED(INB)=IMED + ENDDO + ! + END SUBROUTINE SAL160 + ! + SUBROUTINE SAL160_2(NB_ELEM,IPAR,RPAR,VOL2,ISURF2_ELEM,NB_SURF2,SURF2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! compute 2D volumes and surfaces + ! - compute node volumes + ! - compute 2D surface areas + ! + !Parameters: input + ! NB_ELEM number of elements + ! NB_NODE number of nodes + ! IPAR integer descriptors for elements + ! RPAR floating point descriptors for elements + ! ISURF2_ELEM 2D surface nber per elem + ! NB_SURF2 number of 2D surface + ! + !Parameters: output + ! VOL2 2D volumes of node + ! SURF2 2D areas of node + ! + !--------------------------------------------------------------------- + ! + USE PRECISION_AND_KINDS, ONLY : PDB + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE + !**** + IMPLICIT NONE + INTEGER, INTENT(IN) :: NB_ELEM + INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR + REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR + REAL(PDB), INTENT(OUT), DIMENSION(:) :: VOL2 + INTEGER, INTENT(IN), DIMENSION(:) :: ISURF2_ELEM + INTEGER, INTENT(IN) :: NB_SURF2 + REAL(PDB), INTENT(OUT), OPTIONAL, DIMENSION(:) :: SURF2 + !**** + INTEGER :: ELEM,NODEBC,NODE,ISURF + LOGICAL :: LGBC + REAL(PDB) :: AUX + !**** + ! initiation + VOL2=0._PDB + IF(NB_SURF2 > 0) SURF2=0._PDB + DO ELEM=1,NB_ELEM + !* compute volume of node and add to volume of region (local) + CALL SAL161(IPAR(1,ELEM),RPAR(:,ELEM),AUX) + NODEBC=IPAR(2,ELEM) + LGBC=NODEBC<=0 + IF(.NOT.LGBC) VOL2(NODEBC)=VOL2(NODEBC)+AUX + NODE=IPAR(3,ELEM) + IF(NODE>0)THEN + VOL2(NODE)=VOL2(NODE)-AUX + ELSE + IF(LGBC) CALL XABORT('SAL160_2: isolated element') + LGBC=.TRUE. + NODEBC=NODE + ENDIF + IF(LGBC) THEN + IF(NODEBC==G_BC_TYPE(-1).OR.NODEBC==G_BC_TYPE(0).OR.NODEBC==G_BC_TYPE(1))THEN + ! compute external surface (loaded in total order number) + ISURF=ISURF2_ELEM(ELEM) + IF(ISURF==0) CALL XABORT('SAL160_2: wrong relation of element - surf ') + CALL SAL162(IPAR(1,ELEM),RPAR(:,ELEM),SURF2(ISURF)) + ENDIF + ENDIF + ENDDO + ! + END SUBROUTINE SAL160_2 + ! + SUBROUTINE SAL161(TYPE,RPAR,VOL2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! computes 'volume' between an element and the x axis + ! + !Parameters: input + ! TYPE type of element + ! RPAR floating point descriptors for elements + ! + !Parameters: output + ! VOL2 '2D area' between the element and the horizontal axis + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + INTEGER, INTENT(IN) :: TYPE + REAL(PDB), INTENT(IN), DIMENSION (:) :: RPAR + REAL(PDB), INTENT(OUT) :: VOL2 + ! DIMENSION RPAR(*) + !**** + REAL(PDB) :: YC,EX,EY,R,PHI1,PHI2,COS1,COS2 + !**** + ! volume is added to node- and substracted from node+ + YC=RPAR(2) + IF(TYPE==1)THEN + ! segment: + EX=RPAR(3) + EY=RPAR(4) + VOL2=EX*(YC+EY/2.) + ELSEIF(TYPE==2)THEN + ! whole circle + VOL2=PI*RPAR(3)*RPAR(3) + ELSEIF(TYPE==3)THEN + ! arc of circle: + R=RPAR(3) + PHI1=RPAR(4) + PHI2=RPAR(5) + COS1=COS(PHI1) + COS2=COS(PHI2) + VOL2=R*(YC*(COS1-COS2)+(R/2.0)*(PHI2-PHI1+ & + & (SIN(PHI1)*COS1-SIN(PHI2)*COS2))) + ELSE + CALL XABORT('SAL161: not implemented') + ENDIF + ! + END SUBROUTINE SAL161 + ! + SUBROUTINE SAL162(TYPE,RPAR,SURF2) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! computes surface of a boundary element + ! + !Parameters: input + ! TYPE type of element + ! RPAR floating point descriptors for elements + ! + !Parameters: output + ! SURF2 '2D length' of element + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + INTEGER, INTENT(IN) :: TYPE + REAL(PDB), INTENT(IN), DIMENSION (:) :: RPAR + REAL(PDB), INTENT(OUT) :: SURF2 + !**** + IF(TYPE==1)THEN + ! segment + SURF2=RPAR(5) + ELSEIF(TYPE<=3)THEN + ! arc of circle + SURF2=(RPAR(5)-RPAR(4))*RPAR(3) + ELSE + CALL XABORT('SAL162: not implemented') + ENDIF + ! + END SUBROUTINE SAL162 + ! + SUBROUTINE SAL170(GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! prints out volume, surface and medium info for 2D macros + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE + USE SAL_TRACKING_TYPES, ONLY : PRTIND + !**** + IMPLICIT NONE + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + INTEGER, PARAMETER :: FOUT =6 + + !**** + IF(PRTIND == 1) THEN + WRITE(FOUT,'(//,10x,''2D geometry'',/,10X,11(''=''),//, & + &I8,'' regions'',/, & + &I8,'' external surfaces'',//)') GG%NB_NODE,GG%NB_SURF2 + ENDIF + ! + END SUBROUTINE SAL170 + ! + SUBROUTINE SALSYM(X1,Y1,X2,Y2,X0,Y0,X4,Y4) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! set the symmetric point (x4,y4) of (x0,y0) relative to symmetry axis + ! (x1,y1)->(x2,y2) + ! + !Parameters: input + ! X1 abscissa coordinate of a point on the symmetry axis + ! Y1 ordinate coordinate of a point on the symmetry axis + ! X2 abscissa coordinate of a point on the symmetry axis + ! Y2 ordinate coordinate of a point on the symmetry axis + ! X0 abscissa coordinate of the point to mirror + ! Y0 ordinate coordinate of the point to mirror + ! + !Parameters: output + ! X4 abscissa coordinate of the symmetric point + ! Y4 ordinate coordinate of the symmetric point + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + REAL(PDB),INTENT(IN) :: X1,Y1,X2,Y2,X0,Y0 + REAL(PDB),INTENT(OUT) :: X4,Y4 + + REAL(PDB) :: A,B,C,DEN,X3,Y3 + A=Y2-Y1; B=X1-X2; C=X2*Y1-X1*Y2; + DEN=A*A+B*B; + IF(DEN==0._PDB) CALL XABORT('SALSYM: division by zero') + X3=(B*(B*X0-A*Y0)-A*C)/DEN; + Y3=(A*(-B*X0+A*Y0)-B*C)/DEN; + X4=X0+2._PDB*(X3-X0); Y4=Y0+2._PDB*(Y3-Y0); + END SUBROUTINE SALSYM + ! + SUBROUTINE SALROT(X0,Y0,THETA,X3,Y3,X4,Y4) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! set the symmetric point (x4,y4) of (x3,y3) relative to the rotation + ! center (x1,y1) with angle theta. + ! + !Parameters: input + ! X0 abscissa coordinate of the rotation center + ! Y0 ordinate coordinate of the rotation center + ! THETA rotation angle + ! X3 abscissa coordinate of the point to mirror + ! Y3 ordinate coordinate of the point to mirror + ! + !Parameters: output + ! X4 abscissa coordinate of the symmetric point + ! Y4 ordinate coordinate of the symmetric point + ! + !--------------------------------------------------------------------- + ! + IMPLICIT NONE + REAL(PDB),INTENT(IN) :: X0,Y0,THETA,X3,Y3 + REAL(PDB),INTENT(OUT) :: X4,Y4 + + REAL(PDB) :: A,B + IF(THETA==0.0_PDB) THEN + X4=X3; Y4=Y3 + ELSE + A=X3-X0; B=Y3-Y0 + X4=X0+A*COS(THETA)-B*SIN(THETA); Y4=Y0+A*SIN(THETA)+B*COS(THETA) + ENDIF + END SUBROUTINE SALROT + ! + SUBROUTINE SALFOLD_0(GG,IPASS,IB,NBBCDA,ALIGN,LFOLD,IFOLD) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! unfold the domain with reflection relative to axis AXIS_XY + ! + !Parameters: input + ! IB actual unfolding axis + ! NBBCDA number of perimeters before unfolding + ! ALIGN unfolding axes + ! LFOLD identification flag to all unfolding axes + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !Parameters: output + ! IFOLD folded element indices corresponding to unfolded ones + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY + USE SAL_NUMERIC_MOD, ONLY : FINDLC,DET_ROSETTA + IMPLICIT NONE + INTEGER, INTENT(IN) :: IPASS,IB,NBBCDA + LOGICAL, DIMENSION(NBBCDA), INTENT(IN) :: LFOLD + REAL(PDB), DIMENSION(3,3,NBBCDA), INTENT(INOUT) :: ALIGN + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + INTEGER, DIMENSION(:), INTENT(OUT) :: IFOLD + ! + ! AXIS_XY values of AXIS_X1, AXIS_Y1, AXIS_X2 and AXIS_Y2 for the + ! reflecting axis + INTEGER :: ELEM,TYPE,OK,TMP_NB_ELEM,TMP_NBBCDA,I,J,IBC,INDBC,IAUX + INTEGER, DIMENSION(3) :: IPAR_TMP + REAL(PDB), DIMENSION(4) :: AXIS_XY + REAL(PDB), DIMENSION(6) :: RPAR_TMP + REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB + REAL(PDB) :: X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, & + Y4B,XMIN,YMIN,XMAX,YMAX,PHI1,PHI2,DELPHI,DET1,DET2 + ! + ! allocatable arrays + INTEGER, ALLOCATABLE, DIMENSION(:) :: PERIM_ELEM + LOGICAL, ALLOCATABLE, DIMENSION(:) :: ISPERIM + INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2 + REAL(PDB), ALLOCATABLE, DIMENSION(:) :: ANGLE,ALBEDO + REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR + REAL(PDB), ALLOCATABLE, DIMENSION(:,:,:) :: ALIGN2 + TYPE(T_SALBCDATA), POINTER, DIMENSION(:) :: TMP_BCDATAREAD + ! + ! compute size of the unfold geometry + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + IF(TYPE==1) THEN + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2); + ENDIF + ENDDO + LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN; + ! + ! allocate new surfacic element containers + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), & + I2(2,GG%NB_ELEM), STAT=OK) + IF(OK/=0) CALL XABORT('SALFOLD_0: not enough memory') + TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB; + ! + ! loop over old elements + TMP_NB_ELEM=0 + THETA1=0._PDB; THETA2=0._PDB; + I2(:2,:GG%NB_ELEM)=0 + AXIS_XY(1)=ALIGN(1,1,IB) ; AXIS_XY(2)=ALIGN(1,2,IB) + AXIS_XY(3)=ALIGN(2,1,IB) ; AXIS_XY(4)=ALIGN(2,2,IB) + OUT1: DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM) + IF(TYPE==1) THEN + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM) + ! Cycle if this element is sitting on an unfolding axe + DO IBC=1,NBBCDA + IF(.NOT.LFOLD(IBC)) CYCLE + IF((IPASS.EQ.1).AND.(IBC.NE.IB)) CYCLE + ALIGN(3,1,IBC)=X1; ALIGN(3,2,IBC)=Y1 + DET1 = DET_ROSETTA(ALIGN(1,1,IBC),3) + ALIGN(3,1,IBC)=X2; ALIGN(3,2,IBC)=Y2; + DET2 = DET_ROSETTA(ALIGN(1,1,IBC),3) + IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) CYCLE OUT1 + ENDDO + ! + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4) + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X2,Y2,DX4,DY4) + ELSE IF(TYPE==2) THEN + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4) + THETA1=0._PDB; THETA2=0._PDB; + ELSE IF(TYPE==3) THEN + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1,Y1,X4,Y4) + X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM)); + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1B,Y1B,X4B,Y4B) + IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN + THETA1=PI/2._PDB + ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN + THETA1=3._PDB*PI/2._PDB + ELSE IF(X4B-X4 > 0._PDB) THEN + THETA1=ATAN((Y4B-Y4)/(X4B-X4)) + ELSE + THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI + ENDIF + X1B=X1+RAD*COS(GG%RPAR(5,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(5,ELEM)); + CALL SALSYM(AXIS_XY(1),AXIS_XY(2),AXIS_XY(3),AXIS_XY(4),X1B,Y1B,X4B,Y4B) + IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN + THETA2=PI/2._PDB + ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN + THETA2=3._PDB*PI/2._PDB + ELSE IF(X4B-X4 > 0._PDB) THEN + THETA2=ATAN((Y4B-Y4)/(X4B-X4)) + ELSE + THETA2=ATAN((Y4B-Y4)/(X4B-X4))+PI + ENDIF + ELSE + WRITE(*,*) " elem=",ELEM," type=",TYPE + CALL XABORT('SALFOLD_0: invalid type of surfacic element') + ENDIF + IPAR_TMP(:3)=0 + RPAR_TMP(:6)=0_PDB + RPAR_TMP(1)=X4; RPAR_TMP(2)=Y4; + IPAR_TMP(1)=TYPE; + IF(TYPE==1) THEN + RPAR_TMP(3)=DX4-X4; RPAR_TMP(4)=DY4-Y4; + RPAR_TMP(5)=SQRT(RPAR_TMP(3)**2+RPAR_TMP(4)**2) + XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4); + XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4); + IPAR_TMP(2)=GG%IPAR(3,ELEM); IPAR_TMP(3)=GG%IPAR(2,ELEM); + ELSE IF((TYPE==2).OR.(TYPE==3)) THEN + RPAR_TMP(3)=GG%RPAR(3,ELEM) ! RADIUS + IF(THETA2>THETA1) THETA1=THETA1+2._PDB*PI + PHI1=THETA2; DELPHI=THETA1-THETA2; + IF(DELPHI>0._PDB)THEN + PHI2=PHI1+DELPHI + ELSE + PHI2=PHI1 + PHI1=PHI1+DELPHI + ENDIF + IF(TYPE==3)THEN + ! arc of circle: put phi1 within 0 and 2*pi + IF(PHI1>2._PDB*PI)THEN + IAUX=INT(PHI1/(2._PDB*PI)) + DELPHI=(2._PDB*PI)*IAUX + PHI1=PHI1-DELPHI ; PHI2=PHI2-DELPHI + ELSEIF(PHI1<0._PDB)THEN + IAUX=INT((-PHI1+1.D-7)/(2._PDB*PI))+1 + DELPHI=(2._PDB*PI)*IAUX + PHI1=PHI1+DELPHI ; PHI2=PHI2+DELPHI + ENDIF + ENDIF + RPAR_TMP(4)=PHI1; RPAR_TMP(5)=PHI2; ! ANGLES + IPAR_TMP(2)=GG%IPAR(2,ELEM); IPAR_TMP(3)=GG%IPAR(3,ELEM) + ENDIF + RPAR_TMP(6)=0._PDB + IF(IPASS==2) THEN + ! remove identical elements at pass 2 + DO I=1,TMP_NB_ELEM + IF((ABS(TMP_RPAR(1,I)-RPAR_TMP(1))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(2,I)-RPAR_TMP(2))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(3,I)-RPAR_TMP(3))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(4,I)-RPAR_TMP(4))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(5,I)-RPAR_TMP(5))<=10.0*EPS)) THEN + CYCLE OUT1 + ENDIF + ENDDO + ENDIF + TMP_NB_ELEM=TMP_NB_ELEM+1 + IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_0: tmp_nb_elem overflow(1)') + TMP_IPAR(:3,TMP_NB_ELEM)=IPAR_TMP(:3) + TMP_RPAR(:6,TMP_NB_ELEM)=RPAR_TMP(:6) + I2(2,ELEM)=TMP_NB_ELEM + IF(IPASS==2) THEN + ! remove identical elements at pass 2 + DO I=1,TMP_NB_ELEM + IF((ABS(TMP_RPAR(1,I)-GG%RPAR(1,ELEM))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(2,I)-GG%RPAR(2,ELEM))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(3,I)-GG%RPAR(3,ELEM))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(4,I)-GG%RPAR(4,ELEM))<=10.0*EPS).AND. & + (ABS(TMP_RPAR(5,I)-GG%RPAR(5,ELEM))<=10.0*EPS)) THEN + CYCLE OUT1 + ENDIF + ENDDO + ENDIF + TMP_NB_ELEM=TMP_NB_ELEM+1 + IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_0: tmp_nb_elem overflow(2)') + TMP_IPAR(:,TMP_NB_ELEM)=GG%IPAR(:,ELEM) + TMP_RPAR(:,TMP_NB_ELEM)=GG%RPAR(:,ELEM) + I2(1,ELEM)=TMP_NB_ELEM + ENDDO OUT1 + DEALLOCATE(GG%IPAR,GG%RPAR) + DO ELEM=1,GG%NB_ELEM + IF(I2(1,ELEM).GT.2*GG%NB_ELEM) CALL XABORT('SALFOLD_0: IFOLD overflow') + IF(I2(1,ELEM).NE.0) IFOLD(I2(1,ELEM))=ELEM + IF(I2(2,ELEM).NE.0) IFOLD(I2(2,ELEM))=ELEM + ENDDO + GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR; + GG%NB_ELEM=TMP_NB_ELEM + ! + ! loop over boundary conditions + ALLOCATE(ISPERIM(GG%NB_ELEM),ALIGN2(3,3,GG%NB_ELEM),ANGLE(GG%NB_ELEM), & + & ALBEDO(GG%NB_ELEM),PERIM_ELEM(GG%NB_ELEM)) + ALIGN2(:3,3,:GG%NB_ELEM)=1.0_PDB + PERIM_ELEM(:GG%NB_ELEM)=0 + ISPERIM(:GG%NB_ELEM)=.FALSE. + TMP_NBBCDA=0 + DO IBC=1,GG%NBBCDA + DO I=1,GG%BCDATAREAD(IBC)%NBER + INDBC=GG%BCDATAREAD(IBC)%ELEMNB(I) + IF(INDBC==0) CYCLE + IF(I2(1,INDBC)/=0) ISPERIM(I2(1,INDBC))=.TRUE. + IF(I2(2,INDBC)/=0) ISPERIM(I2(2,INDBC))=.TRUE. + ENDDO + ENDDO + ITER0: DO ELEM=1,GG%NB_ELEM + IF(.NOT.ISPERIM(ELEM)) CYCLE + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + DO J=1,TMP_NBBCDA + ALIGN2(3,1,J)=X1; ALIGN2(3,2,J)=Y1; + DET1 = DET_ROSETTA(ALIGN2(1,1,J),3) + ALIGN2(3,1,J)=X2; ALIGN2(3,2,J)=Y2; + DET2 = DET_ROSETTA(ALIGN2(1,1,J),3) + IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) THEN + PERIM_ELEM(ELEM) = J + CYCLE ITER0 + ENDIF + ENDDO + TMP_NBBCDA=TMP_NBBCDA+1 + PERIM_ELEM(ELEM) = TMP_NBBCDA + ANGLE(TMP_NBBCDA)=ATAN((Y2-Y1)/(X2-X1)) + IF(ABS(ANGLE(TMP_NBBCDA)).LE.1.0E-5) ANGLE(TMP_NBBCDA)=0.0 + ALIGN2(1,1,TMP_NBBCDA)=X1; ALIGN2(1,2,TMP_NBBCDA)=Y1 + ALIGN2(2,1,TMP_NBBCDA)=X2; ALIGN2(2,2,TMP_NBBCDA)=Y2 + ! Recover albedo from folded geometry + ALBEDO(TMP_NBBCDA)=1.0 + DO IBC=1,GG%NBBCDA + J = FINDLC(GG%BCDATAREAD(IBC)%ELEMNB,ELEM) + IF(J.EQ.1) THEN + ALBEDO(TMP_NBBCDA)=GG%BCDATAREAD(IBC)%BCDATA(6) + EXIT + ENDIF + ENDDO + ENDDO ITER0 + ALLOCATE(TMP_BCDATAREAD(TMP_NBBCDA)) + DO IBC=1,TMP_NBBCDA + TMP_BCDATAREAD(IBC)%NBER = COUNT(PERIM_ELEM(:GG%NB_ELEM) == IBC) + ALLOCATE(TMP_BCDATAREAD(IBC)%ELEMNB(TMP_BCDATAREAD(IBC)%NBER)) + TMP_BCDATAREAD(IBC)%SALTYPE = 0 + J=0 + DO I=1,GG%NB_ELEM + IF(PERIM_ELEM(I) == IBC) THEN + J=J+1 + TMP_BCDATAREAD(IBC)%ELEMNB(J) = I + ENDIF + ENDDO + TMP_BCDATAREAD(IBC)%BCDATA(1) = ALIGN2(1,1,IBC) + TMP_BCDATAREAD(IBC)%BCDATA(2) = ALIGN2(1,2,IBC) + TMP_BCDATAREAD(IBC)%BCDATA(3) = COS(ANGLE(IBC)) + TMP_BCDATAREAD(IBC)%BCDATA(4) = SIN(ANGLE(IBC)) + TMP_BCDATAREAD(IBC)%BCDATA(5) = ANGLE(IBC) + TMP_BCDATAREAD(IBC)%BCDATA(6) = ALBEDO(IBC) + ENDDO + DEALLOCATE(I2,PERIM_ELEM,ALBEDO,ANGLE,ALIGN2,ISPERIM) + DEALLOCATE(GG%BCDATAREAD) + GG%BCDATAREAD=>TMP_BCDATAREAD + GG%NBBCDA=TMP_NBBCDA + GG%ALBEDO=1.D0 + END SUBROUTINE SALFOLD_0 + ! + SUBROUTINE SALFOLD_1(HSYM,GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! unfold the domain with reflection + ! + !Parameters: input + ! HSYM: type of symmetry: 'DIAG': diagonal symmetry; 'SYMX': symmetry + ! relative to X axis; 'SYMY': symmetry relative to Y axis; + ! 'SA60': SA60 symmetry; 'S30': S30 symmetry + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY + USE SAL_NUMERIC_MOD, ONLY : FINDLC,DET_ROSETTA + IMPLICIT NONE + CHARACTER(LEN=4),INTENT(IN) :: HSYM + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + ! + INTEGER :: ELEM,TYPE,OK,TMP_NB_ELEM,TMP_NBBCDA,I,J,IBC,INDBC,IAUX,ISYM,NSYM + REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB + REAL(PDB) :: AXIS_X1(2),AXIS_X2(2),AXIS_Y1(2),AXIS_Y2(2) + REAL(PDB) :: X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, & + Y4B,XMIN,YMIN,XMAX,YMAX,PHI1,PHI2,DELPHI,DET1,DET2 + LOGICAL :: NOCOPY(2) + ! + ! allocatable arrays + INTEGER, ALLOCATABLE, DIMENSION(:) :: PERIM_ELEM + LOGICAL, ALLOCATABLE, DIMENSION(:) :: ISPERIM + INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2 + REAL(PDB), ALLOCATABLE, DIMENSION(:) :: ANGLE,ALBEDO + REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR + REAL(PDB), ALLOCATABLE, DIMENSION(:,:,:) :: ALIGN + TYPE(T_SALBCDATA), POINTER, DIMENSION(:) :: TMP_BCDATAREAD + ! + ! compute size of the unfold geometry + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + IF(TYPE==1) THEN + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2); + ENDIF + ENDDO + LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN; + ! + NSYM=1 + IF(HSYM=='SYMX') THEN ! define symmetry axis + AXIS_X1(1)=0._PDB; AXIS_X2(1)=100._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=0._PDB; + ELSE IF(HSYM=='SYMY') THEN + AXIS_X1(1)=0._PDB; AXIS_X2(1)=0._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=100._PDB; + ELSE IF(HSYM=='DIAG') THEN + AXIS_X1(1)=0._PDB; AXIS_X2(1)=100._PDB; AXIS_Y1(1)=0._PDB; AXIS_Y2(1)=100._PDB; + ELSE IF(HSYM=='SA60') THEN + ! the hexagon side is on south + NSYM=2 + AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.5_PDB*LENGTHX + AXIS_Y2(1)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX + AXIS_X1(2)=XMIN+LENGTHX; AXIS_Y1(2)=YMIN; AXIS_X2(2)=XMIN+0.5_PDB*LENGTHX + AXIS_Y2(2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX + ELSE IF(HSYM=='SB60') THEN + ! the hexagon side is on north-east + NSYM=2 + AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.5_PDB*LENGTHX + AXIS_Y2(1)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX + AXIS_X1(2)=XMIN; AXIS_Y1(2)=YMIN; AXIS_X2(2)=XMIN + AXIS_Y2(2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX + ELSE IF(HSYM=='S30') THEN + AXIS_X1(1)=XMIN; AXIS_Y1(1)=YMIN; AXIS_X2(1)=XMIN+0.75_PDB*LENGTHX + AXIS_Y2(1)=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX + ELSE IF(HSYM=='SYMH') THEN + AXIS_X1(1)=XMIN ; AXIS_Y1(1)=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX + AXIS_X2(1)=XMIN+LENGTHX ; AXIS_Y2(1)=AXIS_Y1(1) + ELSE + CALL XABORT('SALFOLD_1: invalid type of symmetry axis') + ENDIF + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + ! + ! allocate new surfacic element containers + ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), & + I2(3,GG%NB_ELEM), STAT=OK) + IF(OK/=0) CALL XABORT('SALFOLD_1: not enough memory') + TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB; + ! + ! loop over old elements + TMP_NB_ELEM=0 + THETA1=0._PDB; THETA2=0._PDB; + DO ELEM=1,GG%NB_ELEM + I2(:,ELEM)=0 + TYPE=GG%IPAR(1,ELEM) + NOCOPY(:2)=.FALSE. + DO ISYM=1,NSYM + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM) + IF(TYPE==1) THEN + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4) + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X2,Y2,DX4,DY4) + NOCOPY(ISYM)=ABS(X1-X4)<10.0*EPS .AND. ABS(Y1-Y4)<10.0*EPS .AND. & + ABS(X2-DX4)<10.0*EPS .AND. ABS(Y2-DY4)<10.0*EPS + IF((HSYM=='SB60').AND.(ISYM==2)) THEN + NOCOPY(2)=(ABS(2._PDB*ABS(Y1-YMIN)-ABS(X4-X1)*SQRT(3._PDB))<10.0*EPS .AND. & + ABS(Y1-Y4)<10.0*EPS .AND. ABS(2._PDB*ABS(Y2-YMIN)-ABS(DX4-X2)*SQRT(3._PDB))<10.0*EPS & + .AND. ABS(Y2-DY4)<10.0*EPS) + ENDIF + IF(NOCOPY(ISYM)) CYCLE + ELSE IF(TYPE==2) THEN + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4) + THETA1=0._PDB; THETA2=0._PDB; + ELSE IF(TYPE==3) THEN + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1,Y1,X4,Y4) + X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM)); + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1B,Y1B,X4B,Y4B) + IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN + THETA1=PI/2._PDB + ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN + THETA1=3._PDB*PI/2._PDB + ELSE IF(X4B-X4 > 0._PDB) THEN + THETA1=ATAN((Y4B-Y4)/(X4B-X4)) + ELSE + THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI + ENDIF + X1B=X1+RAD*COS(GG%RPAR(5,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(5,ELEM)); + CALL SALSYM(AXIS_X1(ISYM),AXIS_Y1(ISYM),AXIS_X2(ISYM),AXIS_Y2(ISYM),X1B,Y1B,X4B,Y4B) + IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN + THETA2=PI/2._PDB + ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN + THETA2=3._PDB*PI/2._PDB + ELSE IF(X4B-X4 > 0._PDB) THEN + THETA2=ATAN((Y4B-Y4)/(X4B-X4)) + ELSE + THETA2=ATAN((Y4B-Y4)/(X4B-X4))+PI + ENDIF + ELSE + WRITE(*,*) " elem=",ELEM," type=",TYPE," isym=",ISYM + CALL XABORT('SALFOLD_1: invalid type of surfacic element') + ENDIF + IF((TYPE==1).AND.(HSYM=='SB60').AND.(ISYM==1).AND.(ABS(Y1)<10.0*EPS).AND.(ABS(Y2)<10.0*EPS)) CYCLE + TMP_NB_ELEM=TMP_NB_ELEM+1 + IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_1: tmp_nb_elem overflow(1)') + I2(ISYM+1,ELEM)=TMP_NB_ELEM + TMP_RPAR(1,TMP_NB_ELEM)=X4; TMP_RPAR(2,TMP_NB_ELEM)=Y4; + TMP_IPAR(1,TMP_NB_ELEM)=TYPE; + IF(TYPE==1) THEN + TMP_RPAR(3,TMP_NB_ELEM)=DX4-X4; TMP_RPAR(4,TMP_NB_ELEM)=DY4-Y4; + TMP_RPAR(5,TMP_NB_ELEM)=SQRT(TMP_RPAR(3,TMP_NB_ELEM)**2+TMP_RPAR(4,TMP_NB_ELEM)**2) + XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4); + XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4); + TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(3,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(2,ELEM); + ELSE IF((TYPE==2).OR.(TYPE==3)) THEN + TMP_RPAR(3,TMP_NB_ELEM)=GG%RPAR(3,ELEM) ! RADIUS + IF(THETA2>THETA1) THETA1=THETA1+2._PDB*PI + PHI1=THETA2; DELPHI=THETA1-THETA2; + IF(DELPHI>0._PDB)THEN + PHI2=PHI1+DELPHI + ELSE + PHI2=PHI1 + PHI1=PHI1+DELPHI + ENDIF + IF(TYPE==3)THEN + ! arc of circle: put phi1 within 0 and 2*pi + IF(PHI1>2._PDB*PI)THEN + IAUX=INT(PHI1/(2._PDB*PI)) + DELPHI=(2._PDB*PI)*IAUX + PHI1=PHI1-DELPHI ; PHI2=PHI2-DELPHI + ELSEIF(PHI1<0._PDB)THEN + IAUX=INT((-PHI1+1.D-7)/(2._PDB*PI))+1 + DELPHI=(2._PDB*PI)*IAUX + PHI1=PHI1+DELPHI ; PHI2=PHI2+DELPHI + ENDIF + ENDIF + TMP_RPAR(4,TMP_NB_ELEM)=PHI1; TMP_RPAR(5,TMP_NB_ELEM)=PHI2; ! ANGLES + TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM) + ENDIF + TMP_RPAR(6,TMP_NB_ELEM)=0._PDB + ENDDO + IF((.NOT.NOCOPY(1)).AND.(.NOT.NOCOPY(2))) THEN + TMP_NB_ELEM=TMP_NB_ELEM+1 + IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_1: tmp_nb_elem overflow(2)') + TMP_IPAR(:,TMP_NB_ELEM)=GG%IPAR(:,ELEM) + TMP_RPAR(:,TMP_NB_ELEM)=GG%RPAR(:,ELEM) + I2(1,ELEM)=TMP_NB_ELEM + ENDIF + ENDDO + DEALLOCATE(GG%IPAR,GG%RPAR) + GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR; + GG%NB_ELEM=TMP_NB_ELEM + ! + ! translate the domain + DO ELEM=1,GG%NB_ELEM + GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-XMIN + GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-YMIN + ENDDO + LENGTHX=XMAX-XMIN ; LENGTHY=YMAX-YMIN ; + ! + ! loop over boundary conditions + TMP_NBBCDA=0 + ALLOCATE(ISPERIM(GG%NB_ELEM),ALIGN(3,3,GG%NB_ELEM),ANGLE(GG%NB_ELEM), & + & ALBEDO(GG%NB_ELEM),PERIM_ELEM(GG%NB_ELEM)) + ALIGN(:3,3,:GG%NB_ELEM)=1.0_PDB + PERIM_ELEM(:GG%NB_ELEM)=0 + ISPERIM(:GG%NB_ELEM)=.FALSE. + DO IBC=1,GG%NBBCDA + DO I=1,GG%BCDATAREAD(IBC)%NBER + INDBC=GG%BCDATAREAD(IBC)%ELEMNB(I) + IF(INDBC==0) CYCLE + DO ISYM=1,NSYM+1 + IF(I2(ISYM,INDBC)/=0) ISPERIM(I2(ISYM,INDBC))=.TRUE. + ENDDO + ENDDO + ENDDO + ITER0: DO ELEM=1,GG%NB_ELEM + IF(.NOT.ISPERIM(ELEM)) CYCLE + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + DO J=1,TMP_NBBCDA + ALIGN(3,1,J)=X1; ALIGN(3,2,J)=Y1; + DET1 = DET_ROSETTA(ALIGN(1,1,J),3) + ALIGN(3,1,J)=X2; ALIGN(3,2,J)=Y2; + DET2 = DET_ROSETTA(ALIGN(1,1,J),3) + IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) THEN + PERIM_ELEM(ELEM) = J + CYCLE ITER0 + ENDIF + ENDDO + TMP_NBBCDA=TMP_NBBCDA+1 + PERIM_ELEM(ELEM) = TMP_NBBCDA + ANGLE(TMP_NBBCDA)=ATAN((Y2-Y1)/(X2-X1)) + IF(ABS(ANGLE(TMP_NBBCDA)).LE.1.0E-5) ANGLE(TMP_NBBCDA)=0.0 + ALIGN(1,1,TMP_NBBCDA)=X1; ALIGN(1,2,TMP_NBBCDA)=Y1 + ALIGN(2,1,TMP_NBBCDA)=X2; ALIGN(2,2,TMP_NBBCDA)=Y2 + ! Recover albedo from folded geometry + ALBEDO(TMP_NBBCDA)=1.0 + DO IBC=1,GG%NBBCDA + J = FINDLC(GG%BCDATAREAD(IBC)%ELEMNB,ELEM) + IF(J.EQ.1) THEN + ALBEDO(TMP_NBBCDA)=GG%BCDATAREAD(IBC)%BCDATA(6) + EXIT + ENDIF + ENDDO + ENDDO ITER0 + ALLOCATE(TMP_BCDATAREAD(TMP_NBBCDA)) + DO IBC=1,TMP_NBBCDA + TMP_BCDATAREAD(IBC)%NBER = COUNT(PERIM_ELEM(:GG%NB_ELEM) == IBC) + ALLOCATE(TMP_BCDATAREAD(IBC)%ELEMNB(TMP_BCDATAREAD(IBC)%NBER)) + TMP_BCDATAREAD(IBC)%SALTYPE = 0 + J=0 + DO I=1,GG%NB_ELEM + IF(PERIM_ELEM(I) == IBC) THEN + J=J+1 + TMP_BCDATAREAD(IBC)%ELEMNB(J) = I + ENDIF + ENDDO + TMP_BCDATAREAD(IBC)%BCDATA(1) = ALIGN(1,1,IBC) + TMP_BCDATAREAD(IBC)%BCDATA(2) = ALIGN(1,2,IBC) + TMP_BCDATAREAD(IBC)%BCDATA(3) = COS(ANGLE(IBC)) + TMP_BCDATAREAD(IBC)%BCDATA(4) = SIN(ANGLE(IBC)) + TMP_BCDATAREAD(IBC)%BCDATA(5) = ANGLE(IBC) + TMP_BCDATAREAD(IBC)%BCDATA(6) = ALBEDO(IBC) + ENDDO + DEALLOCATE(I2,PERIM_ELEM,ALBEDO,ANGLE,ALIGN,ISPERIM) + DEALLOCATE(GG%BCDATAREAD) + GG%BCDATAREAD=>TMP_BCDATAREAD + GG%NBBCDA=TMP_NBBCDA + GG%ALBEDO=1.D0 + END SUBROUTINE SALFOLD_1 + ! + SUBROUTINE SALFOLD_2(HSYM,GG) + ! + !--------------------------------------------------------------------- + ! + !Purpose: + ! unfold the domain with rotation + ! + !Parameters: input + ! HSYM: type of symmetry: SR60': dual 60-degree rotation; R180: 180-degree + ! rotation; R120: dual 120-degree rotation + ! + !Parameters: input/output + ! GG geometry descriptor + ! + !--------------------------------------------------------------------- + ! + USE SAL_GEOMETRY_TYPES, ONLY : NIPAR,NRPAR,LENGTHX,LENGTHY + IMPLICIT NONE + CHARACTER(LEN=4),INTENT(IN) :: HSYM + TYPE(T_G_BASIC), INTENT(INOUT) :: GG + ! + INTEGER :: ELEM,ELEM2,TYPE,OK,TMP_NB_ELEM,ISYM,NSYM,IAUX + REAL(PDB) :: THROT(3),X1,X2,X4,Y1,Y2,Y4,DX4,DY4,RAD,THETA1,THETA2,X1B,Y1B,X4B, & + Y4B,XMIN,YMIN,XMAX,YMAX,CENTER_X,CENTER_Y,DELPHI + REAL(PDB),PARAMETER :: EPS=1.0E-5_PDB + ! + ! allocatable arrays + LOGICAL, ALLOCATABLE, DIMENSION(:) :: ELEM_DUP + INTEGER, POINTER, DIMENSION(:,:) :: TMP_IPAR + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: I2 + REAL(PDB), POINTER, DIMENSION(:,:) :: TMP_RPAR + ! + ! compute size of the unfold geometry + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + IF(TYPE==1) THEN + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2); + ENDIF + ENDDO + LENGTHX=XMAX-XMIN; LENGTHY=YMAX-YMIN; + ! + NSYM=0 + IF(HSYM=='SR60') THEN ! define rotation center + NSYM=3 + CENTER_X=XMIN+0.5_PDB*LENGTHX; CENTER_Y=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX + THROT(1)=0.0_PDB; THROT(2)=-PI/3.0_PDB; THROT(3)=PI/3.0_PDB + ELSE IF(HSYM=='R180') THEN + NSYM=2 + CENTER_X=XMIN+0.5_PDB*LENGTHX; CENTER_Y=YMIN+0.25_PDB*SQRT(3._PDB)*LENGTHX + THROT(1)=0.0_PDB; THROT(2)=PI; + ELSE IF(HSYM=='R120') THEN + NSYM=3 + CENTER_X=XMIN+0.5_PDB*LENGTHX/1.5_PDB; CENTER_Y=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX/1.5_PDB + THROT(1)=0.0_PDB; THROT(2)=-2.0_PDB*PI/3.0_PDB; THROT(3)=2.0_PDB*PI/3.0_PDB + ELSE + CALL XABORT('SALFOLD_2: invalid type of symmetry axis') + ENDIF + XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB; + ! + ! allocate new surfacic element containers + ALLOCATE(I2(NSYM,GG%NB_ELEM)) + ALLOCATE(TMP_IPAR(NIPAR,3*GG%NB_ELEM), TMP_RPAR(NRPAR,3*GG%NB_ELEM), STAT=OK) + IF(OK/=0) CALL XABORT('SALFOLD_2: not enough memory') + I2(:NSYM,:GG%NB_ELEM)=0 + TMP_IPAR(:,:)=0; TMP_RPAR(:,:)=0._PDB; + ! + ! loop over old elements + TMP_NB_ELEM=0 + THETA1=0._PDB; THETA2=0._PDB; + DO ELEM=1,GG%NB_ELEM + TYPE=GG%IPAR(1,ELEM) + DO ISYM=1,NSYM + X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM); RAD=GG%RPAR(3,ELEM) + IF(TYPE==1) THEN + XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1); + X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM); + CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4) + CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X2,Y2,DX4,DY4) + ELSE IF(TYPE==2) THEN + CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4) + THETA1=0._PDB; THETA2=0._PDB; + ELSE IF(TYPE==3) THEN + DELPHI=GG%RPAR(5,ELEM)-GG%RPAR(4,ELEM) + X1B=X1+RAD*COS(GG%RPAR(4,ELEM)); Y1B=Y1+RAD*SIN(GG%RPAR(4,ELEM)); + CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1,Y1,X4,Y4) + CALL SALROT(CENTER_X,CENTER_Y,THROT(ISYM),X1B,Y1B,X4B,Y4B) + IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 > 0._PDB)) THEN + THETA1=PI/2._PDB + ELSE IF((ABS(X4B-X4)<EPS*ABS(RAD)).AND.(Y4B-Y4 < 0._PDB)) THEN + THETA1=3._PDB*PI/2._PDB + ELSE IF(X4B-X4 > 0._PDB) THEN + THETA1=ATAN((Y4B-Y4)/(X4B-X4)) + ELSE + THETA1=ATAN((Y4B-Y4)/(X4B-X4))+PI + ENDIF + THETA2=THETA1+DELPHI + ! put THETA1 within 0 and 2*pi + IF(THETA1>2._PDB*PI)THEN + IAUX=INT(THETA1/(2._PDB*PI)) + DELPHI=(2._PDB*PI)*IAUX + THETA1=THETA1-DELPHI ; THETA2=THETA2-DELPHI + ELSEIF(THETA1<0._PDB)THEN + IAUX=INT((-THETA1+1.D-7)/(2._PDB*PI))+1 + DELPHI=(2._PDB*PI)*IAUX + THETA1=THETA1+DELPHI ; THETA2=THETA2+DELPHI + ENDIF + ELSE + WRITE(*,*) " elem=",ELEM," type=",TYPE," isym=",ISYM + CALL XABORT('SALFOLD_2: invalid type of surfacic element') + ENDIF + TMP_NB_ELEM=TMP_NB_ELEM+1 + IF(TMP_NB_ELEM>3*GG%NB_ELEM) CALL XABORT('SALFOLD_2: TMP_NB_ELEM overflow') + I2(ISYM,ELEM)=TMP_NB_ELEM + TMP_RPAR(1,TMP_NB_ELEM)=X4; TMP_RPAR(2,TMP_NB_ELEM)=Y4; + TMP_IPAR(1,TMP_NB_ELEM)=TYPE; + IF(TYPE==1) THEN + TMP_RPAR(3,TMP_NB_ELEM)=DX4-X4; TMP_RPAR(4,TMP_NB_ELEM)=DY4-Y4; + TMP_RPAR(5,TMP_NB_ELEM)=SQRT(TMP_RPAR(3,TMP_NB_ELEM)**2+TMP_RPAR(4,TMP_NB_ELEM)**2) + XMIN=MIN(XMIN,X4); YMIN=MIN(YMIN,Y4); XMAX=MAX(XMAX,X4); YMAX=MAX(YMAX,Y4); + XMIN=MIN(XMIN,DX4); YMIN=MIN(YMIN,DY4); XMAX=MAX(XMAX,DX4); YMAX=MAX(YMAX,DY4); + TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM) + ELSE IF((TYPE==2).OR.(TYPE==3)) THEN + TMP_RPAR(3,TMP_NB_ELEM)=GG%RPAR(3,ELEM) ! RADIUS + TMP_RPAR(4,TMP_NB_ELEM)=THETA1; TMP_RPAR(5,TMP_NB_ELEM)=THETA2; ! ANGLES + TMP_IPAR(2,TMP_NB_ELEM)=GG%IPAR(2,ELEM); TMP_IPAR(3,TMP_NB_ELEM)=GG%IPAR(3,ELEM) + ENDIF + TMP_RPAR(6,TMP_NB_ELEM)=0._PDB + ENDDO + ENDDO + DEALLOCATE(GG%IPAR,GG%RPAR) + ! + ! eliminate duplicate elements + ALLOCATE(ELEM_DUP(TMP_NB_ELEM)) + ELEM_DUP(:TMP_NB_ELEM)=.FALSE. + DO ELEM=1,TMP_NB_ELEM + TYPE=TMP_IPAR(1,ELEM) + IF((TYPE==1).AND.((TMP_IPAR(2,ELEM)==0).OR.(TMP_IPAR(3,ELEM)==0))) THEN + X4=TMP_RPAR(1,ELEM); Y4=TMP_RPAR(2,ELEM) + DX4=X4+TMP_RPAR(3,ELEM); DY4=Y4+TMP_RPAR(4,ELEM) + DO ELEM2=ELEM+1,TMP_NB_ELEM + IF(TMP_IPAR(1,ELEM2)/=1) CYCLE + X1=TMP_RPAR(1,ELEM2); Y1=TMP_RPAR(2,ELEM2) + X2=X1+TMP_RPAR(3,ELEM2); Y2=Y1+TMP_RPAR(4,ELEM2) + IF(((ABS(X1-X4)<10.0*EPS .AND. ABS(Y1-Y4)<10.0*EPS .AND. & + ABS(X2-DX4)<10.0*EPS .AND. ABS(Y2-DY4)<10.0*EPS)).OR. & + ((ABS(X1-DX4)<10.0*EPS .AND. ABS(Y1-DY4)<10.0*EPS .AND. & + ABS(X2-X4)<10.0*EPS .AND. ABS(Y2-Y4)<10.0*EPS))) THEN + ELEM_DUP(ELEM2)=.TRUE. + IF(TMP_IPAR(2,ELEM)==0) TMP_IPAR(2,ELEM)=MAX(TMP_IPAR(2,ELEM2),TMP_IPAR(3,ELEM2)) + IF(TMP_IPAR(3,ELEM)==0) TMP_IPAR(3,ELEM)=MAX(TMP_IPAR(2,ELEM2),TMP_IPAR(3,ELEM2)) + IF(TMP_IPAR(2,ELEM)==TMP_IPAR(3,ELEM)) ELEM_DUP(ELEM)=.TRUE. + EXIT + ENDIF + ENDDO + ENDIF + ENDDO + ! + ELEM=1 + DO WHILE(ELEM<=TMP_NB_ELEM) + IF(ELEM_DUP(ELEM)) THEN + TMP_NB_ELEM=TMP_NB_ELEM-1 + DO ELEM2=ELEM,TMP_NB_ELEM + TMP_IPAR(:,ELEM2)=TMP_IPAR(:,ELEM2+1); TMP_RPAR(:,ELEM2)=TMP_RPAR(:,ELEM2+1) + ELEM_DUP(ELEM2)=ELEM_DUP(ELEM2+1) + ENDDO + ELSE + ELEM=ELEM+1 + ENDIF + ENDDO + GG%IPAR=>TMP_IPAR; GG%RPAR=>TMP_RPAR; + GG%NB_ELEM=TMP_NB_ELEM + ! + ! translate the domain + DO ELEM=1,GG%NB_ELEM + GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-XMIN + GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-YMIN + ENDDO + LENGTHX=XMAX-XMIN ; LENGTHY=YMAX-YMIN ; + GG%NBBCDA=0; GG%ALBEDO=1.D0 + ! + DEALLOCATE(ELEM_DUP,I2) + END SUBROUTINE SALFOLD_2 +END MODULE SAL_GEOMETRY_MOD |