summaryrefslogtreecommitdiff
path: root/Dragon/src
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src')
-rw-r--r--Dragon/src/EVODRV.f12
-rw-r--r--Dragon/src/SAL_GEOMETRY_MOD.f9056
-rw-r--r--Dragon/src/SAL_TRAJECTORY_MOD.f9012
3 files changed, 65 insertions, 15 deletions
diff --git a/Dragon/src/EVODRV.f b/Dragon/src/EVODRV.f
index afe2f35..cb9b40c 100644
--- a/Dragon/src/EVODRV.f
+++ b/Dragon/src/EVODRV.f
@@ -552,7 +552,11 @@
ENDIF
IF(ITIM.EQ.0) THEN
IF(NTIM.GT.0) THEN
- IF(XT(IP).LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+ IF(XT(IP).LT.TIMES(NTIM)) THEN
+ WRITE(HSMG,'(19HEVODRV: INVALID XT(,I4,2H)=,1P,E12.4,2H <,
+ > E12.4,1H.)') IP,XT(IP),TIMES(NTIM)
+ CALL XABORT(HSMG)
+ ENDIF
ENDIF
NTIM=NTIM+1
TIMES(NTIM)=XT(IP)
@@ -603,7 +607,11 @@
ENDIF
IF(ITIM.EQ.0) THEN
IF(NTIM.GT.0) THEN
- IF(XTI.LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+ IF(XTI.LT.TIMES(NTIM)) THEN
+ WRITE(HSMG,'(20HEVODRV: INVALID XTI=,1P,E12.4,2H <,E12.4,
+ > 1H.)') XTI,TIMES(NTIM)
+ CALL XABORT(HSMG)
+ ENDIF
ENDIF
NTIM=NTIM+1
TIMES(NTIM)=XTI
diff --git a/Dragon/src/SAL_GEOMETRY_MOD.f90 b/Dragon/src/SAL_GEOMETRY_MOD.f90
index 4bce6e8..a56376f 100644
--- a/Dragon/src/SAL_GEOMETRY_MOD.f90
+++ b/Dragon/src/SAL_GEOMETRY_MOD.f90
@@ -154,6 +154,7 @@ CONTAINS
! 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
@@ -188,10 +189,48 @@ CONTAINS
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)
!
@@ -1818,7 +1857,7 @@ CONTAINS
!
!---------------------------------------------------------------------
!
- USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE
+ USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE
USE SAL_NUMERIC_MOD, ONLY : SAL141
!****
IMPLICIT NONE
@@ -1840,6 +1879,7 @@ CONTAINS
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
@@ -1852,13 +1892,6 @@ CONTAINS
!* calculate number of elements on axes,
! and their distances to the origin of the axis
NBE=0
-
- DO I=1,NPERIM
- ELEM=PERIM(I)
- IBC=IBC2_ELEM(ELEM)
- TYPBC=TYPE_BC2(IBC)
- IDATA=IDATA_BC2(IBC)
- ENDDO
DO I=1,NPERIM
ELEM=PERIM(I)
IBC=IBC2_ELEM(ELEM)
@@ -2048,7 +2081,10 @@ CONTAINS
CASE DEFAULT
CALL XABORT('SAL130_10: boundary condition not implemented.')
END SELECT
- IF(M>=NAXES+1) CALL XABORT('SAL130_10: element not on axes')
+ 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
diff --git a/Dragon/src/SAL_TRAJECTORY_MOD.f90 b/Dragon/src/SAL_TRAJECTORY_MOD.f90
index ac1714d..9408d1b 100644
--- a/Dragon/src/SAL_TRAJECTORY_MOD.f90
+++ b/Dragon/src/SAL_TRAJECTORY_MOD.f90
@@ -68,7 +68,7 @@ CONTAINS
!
!---------------------------------------------------------------------
!
- USE SAL_GEOMETRY_TYPES, ONLY : ISPEC
+ USE SAL_GEOMETRY_TYPES, ONLY : ISPEC,TYPGEO
USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,ANGTAB,ELMTAB,CNT,CNT0,NB_TOT,DNEW,DINIT, &
NNEW,LNEW,IERR,LGMORE,DD0,NTRACK,EPS1,EX0,EY0,LGOK,IPART,DELX, &
N_AXIS
@@ -110,7 +110,7 @@ CONTAINS
ELSE
WRITE(*,*) 'PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS) :',PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS)
WRITE(*,*) 'DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1) :',DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1)
- WRITE(*,*) 'DELX :',DELX,' RADIA=',RADIA
+ WRITE(*,*) 'DELX :',DELX,' RADIA=',RADIA,' TYPGEO=',TYPGEO
CALL XABORT('SALTRA: Cant find entry point')
ENDIF
ENDIF
@@ -164,6 +164,7 @@ CONTAINS
!
!---------------------------------------------------------------------
!
+ USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO
USE SAL_TRACKING_TYPES, ONLY : IPART,N_AXIS,DNEW,DELX,NNEW,LNEW,COSINE,AX, &
AY,HX,HY,BX,BY,EX,EY
INTEGER, INTENT(IN) :: NPERIM
@@ -171,6 +172,7 @@ CONTAINS
REAL(PDB), INTENT(IN), DIMENSION(:) :: DIST_AXIS
INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR
INTEGER :: I,J
+ CHARACTER(LEN=131) :: HSMG
!***
LNEW=0
!* compute crossed element
@@ -180,7 +182,11 @@ CONTAINS
EXIT
ENDIF
ENDDO
- IF(LNEW==0) CALL XABORT('SAL241_2: Error of distances on the axis')
+ IF(LNEW==0) THEN
+ WRITE(HSMG,'(52HSAL241_2: Error of distances on the axis for typgeo=, &
+ & i3,1h.)') TYPGEO
+ CALL XABORT(HSMG)
+ ENDIF
!* get entered node
NNEW=IPAR(2,LNEW)
IF(NNEW<0) NNEW=IPAR(3,LNEW)