! !--------------------------------------------------------------------- ! !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 REAL(PDB) :: X1,Y1,X2,Y2,XMIN,XMAX,YMIN,YMAX 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) ! ! translate the domain for cyclic cases 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 IF(ABS(XMIN).LT.1.D-10) XMIN=0.D0 IF(ABS(YMIN).LT.1.D-10) YMIN=0.D0 SELECT CASE(TYPGEO) CASE(4,9) DO ELEM=1,GG%NB_ELEM GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-(XMIN+XMAX)/2.D0 GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-(XMIN+XMAX)/2.D0 ENDDO DO I=1,GG%NBBCDA GG%BCDATAREAD(I)%BCDATA(1)=GG%BCDATAREAD(I)%BCDATA(1)-(XMIN+XMAX)/2.D0 GG%BCDATAREAD(I)%BCDATA(2)=GG%BCDATAREAD(I)%BCDATA(2)-(XMIN+XMAX)/2.D0 IF(ABS(GG%BCDATAREAD(I)%BCDATA(1)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(1)=0.D0 IF(ABS(GG%BCDATAREAD(I)%BCDATA(2)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(2)=0.D0 ENDDO CASE(5,6,7,8,10,11,12) 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 DO I=1,GG%NBBCDA GG%BCDATAREAD(I)%BCDATA(1)=GG%BCDATAREAD(I)%BCDATA(1)-XMIN GG%BCDATAREAD(I)%BCDATA(2)=GG%BCDATAREAD(I)%BCDATA(2)-YMIN IF(ABS(GG%BCDATAREAD(I)%BCDATA(1)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(1)=0.D0 IF(ABS(GG%BCDATAREAD(I)%BCDATA(2)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(2)=0.D0 ENDDO END SELECT ! !* 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 : TYPGEO,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 CHARACTER(LEN=131) :: HSMG 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)) IF((TYPGEO.EQ.3).OR.(TYPGEO.EQ.6).OR.(TYPGEO.EQ.6)) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/2._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ELSE IF(TYPGEO.EQ.7) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/2._PDB).AND.(ANGLE/=PI/4._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ELSE IF((TYPGEO.EQ.8).OR.(TYPGEO.EQ.10)) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/3._PDB).AND.(ANGLE/=2._PDB*PI/3._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ELSE IF(TYPGEO.EQ.9) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/6._PDB).AND.(ANGLE/=2._PDB*PI/3._PDB).AND. & & (ANGLE/=PI).AND.(ANGLE/=7._PDB*PI/6._PDB).AND.(ANGLE/=5._PDB*PI/3._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ELSE IF(TYPGEO.EQ.11) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/3._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ELSE IF(TYPGEO.EQ.12) THEN IF((ANGLE/=0.0).AND.(ANGLE/=PI/6._PDB).AND.(ANGLE/=2._PDB*PI/3._PDB)) THEN WRITE(HSMG,100) ANGLE,I,TYPGEO CALL XABORT(HSMG) ENDIF ENDIF ENDIF ENDDO ENDIF 100 FORMAT(34HSALINP: FORBIDDEN PERIMETER ANGLE=,1P,E13.4,18H RADIANS FOR SIDE=,I3, & & 12H AND TYPGEO=,I3,1H.) 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 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))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(LY0.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))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))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))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) & 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)1) THEN DO K=1,NBE(M)-1 IF(AUX_DIST(NBE(M),M)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 CHARACTER(LEN=131) :: HSMG !**** 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) 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)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)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)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)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) THEN WRITE(HSMG,'(42HSAL130_10: element not on axes for typgeo=,i3,1h.)') TYPGEO CALL XABORT(HSMG) ENDIF 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)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))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) 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.) 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 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) 0._PDB)) THEN THETA1=PI/2._PDB ELSE IF((ABS(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) 0._PDB)) THEN THETA2=PI/2._PDB ELSE IF((ABS(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,ALIGN3(3,3) 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(NOCOPY(ISYM)) CYCLE IF(HSYM=='SB60') THEN IF(ISYM==1) THEN IF((ABS(Y1)<10.0*EPS).AND.(ABS(Y2)<10.0*EPS)) CYCLE ELSE IF((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)) CYCLE ALIGN3(1,2)=YMIN ; ALIGN3(2,2)=YMIN+0.5_PDB*SQRT(3._PDB)*LENGTHX ALIGN3(1,1)=-XMIN ALIGN3(2,1)=-XMIN-0.5_PDB*LENGTHX ALIGN3(:3,3)=1.0_PDB ALIGN3(3,1)=X4; ALIGN3(3,2)=Y4; DET1 = DET_ROSETTA(ALIGN3(1,1),3) ALIGN3(3,1)=DX4; ALIGN3(3,2)=DY4; DET2 = DET_ROSETTA(ALIGN3(1,1),3) IF((ABS(DET1).LE.1.0E-4).AND.(ABS(DET2).LE.1.0E-4)) CYCLE ENDIF ENDIF 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) 0._PDB)) THEN THETA1=PI/2._PDB ELSE IF((ABS(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) 0._PDB)) THEN THETA2=PI/2._PDB ELSE IF((ABS(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 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) 0._PDB)) THEN THETA1=PI/2._PDB ELSE IF((ABS(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