diff options
Diffstat (limited to 'Trivac/src/TRITCO.f')
| -rwxr-xr-x | Trivac/src/TRITCO.f | 252 |
1 files changed, 252 insertions, 0 deletions
diff --git a/Trivac/src/TRITCO.f b/Trivac/src/TRITCO.f new file mode 100755 index 0000000..87ec38d --- /dev/null +++ b/Trivac/src/TRITCO.f @@ -0,0 +1,252 @@ +*DECK TRITCO + SUBROUTINE TRITCO (NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,KK4,KK5, + 1 VOL0,MAT,MATN,DIF,DDF,SIDE,ZZ,QFR,IPR,A) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Mesh centered finite difference coefficients in hexagonal geometry +* with triangular sub meshing. +* +*Copyright: +* Copyright (C) 2002 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. Benaboud +* +*Parameters: input +* NEL total number of finite elements. +* LL4 order of the system matrices. +* ISPLH number of triangles (equal to 6*(ISPLH-1)**2). +* IR first dimension for matrices DIF and DDF. +* IQF index in array QFR. +* K index of finite element. +* KK1 first neighbour of the triangular finite element. +* KK2 second neighbour of the triangular finite element. +* KK3 third neighbour of the triangular finite element. +* KK4 fourth neighbour of the triangular finite element. +* KK5 fifth neighbour of the triangular finite element. +* VOL0 volume of the finite element. +* MAT mixture index assigned to each hexagon. +* MATN mixture index assigned to each triangle. +* DIF directional diffusion coefficients. +* DDF variation of directional diffusion coefficients. +* SIDE side of an hexagon. +* ZZ Z-directed mesh spacings. +* QFR element-ordered boundary conditions. +* IPR type of matrix assembly: +* =0: compute the system matrices; +* =1: compute the derivative of system matrices; +* =2 or =3: compute the variation of system matrices. +* +*Parameters: output +* A mesh centered finite difference coefficients. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,KK4,KK5,MAT(NEL), + 1 MATN(LL4),IPR + REAL VOL0,DIF(IR,3),DDF(IR,3),SIDE,ZZ(NEL),QFR(8) + DOUBLE PRECISION A(5) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION SHARM,DHARM,VHARM +* FORMULA WIHOUT VARIATION OR DERIVATIVE. + SHARM(X1,X2,DIF1,DIF2)=2.0D0*DIF1*DIF2/(X1*DIF2+X2*DIF1) +* FORMULA WITH DERIVATIVE. + DHARM(X1,X2,DIF1,DIF2,DDF1,DDF2)=2.0D0*(X1*DIF2*DIF2*DDF1+ + 1 X2*DIF1*DIF1*DDF2)/(X1*DIF2+X2*DIF1)**2 +* FORMULA WITH VARIATION. + VHARM(X1,X2,DIF1,DIF2,DDF1,DDF2)=2.0D0*((DIF1+DDF1)*(DIF2+DDF2) + 1 /(X1*(DIF2+DDF2)+X2*(DIF1+DDF1))-DIF1*DIF2/(X1*DIF2+X2*DIF1)) +* + L=MAT(K) + DZ=ZZ(K) + DS=SIDE/(SQRT(3.0)*(ISPLH-1)) + DT=SIDE/(ISPLH-1) + IF(IPR.EQ.0) THEN +* FORMULE DIRECTE. + IF(KK1.GT.0) THEN + A(1)=SHARM(DS,DS,DIF(L,1),DIF(MATN(KK1),1))*DT*DZ + ELSE IF(KK1.EQ.-1) THEN + A(1)=SHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0)*DT*DZ + ELSE IF(KK1.EQ.-2) THEN + A(1)=0.0D0 + ELSE IF(KK1.EQ.-3) THEN + A(1)=2.0D0*SHARM(DS,DS,DIF(L,1),DIF(L,1))*DT*DZ + ENDIF +* + IF(KK2.GT.0) THEN + A(2)=SHARM(DS,DS,DIF(L,1),DIF(MATN(KK2),1))*DT*DZ + ELSE IF(KK2.EQ.-1) THEN + A(2)=SHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0)*DT*DZ + ELSE IF(KK2.EQ.-2) THEN + A(2)=0.0D0 + ELSE IF(KK2.EQ.-3) THEN + A(2)=2.0D0*SHARM(DS,DS,DIF(L,1),DIF(L,1))*DT*DZ + ENDIF +* + IF(KK3.GT.0) THEN + A(3)=SHARM(DS,DS,DIF(L,1),DIF(MATN(KK3),1))*DT*DZ + ELSE IF(KK3.EQ.-1) THEN + A(3)=SHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0)*DT*DZ + ELSE IF(KK3.EQ.-2) THEN + A(3)=0.0D0 + ELSE IF(KK3.EQ.-3) THEN + A(3)=2.0D0*SHARM(DS,DS,DIF(L,1),DIF(L,1))*DT*DZ + ENDIF +* + IF(KK4.GT.0) THEN + A(4)=SHARM(DZ,ZZ(KK4),DIF(L,1),DIF(MAT(KK4),1))*VOL0/DZ + ELSE IF(KK4.EQ.-1) THEN + A(4)=SHARM(DZ,DZ,DIF(L,1),DZ*QFR(7)/2.0)*VOL0/DZ + ELSE IF(KK4.EQ.-2) THEN + A(4)=0.0D0 + ELSE IF(KK4.EQ.-3) THEN + A(4)=2.0D0*SHARM(DZ,DZ,DIF(L,1),DIF(L,1))*VOL0/DZ + ENDIF +* + IF(KK5.GT.0) THEN + A(5)=SHARM(DZ,ZZ(KK5),DIF(L,1),DIF(MAT(KK5),1))*VOL0/DZ + ELSE IF(KK5.EQ.-1) THEN + A(5)=SHARM(DZ,DZ,DIF(L,1),DZ*QFR(8)/2.0)*VOL0/DZ + ELSE IF(KK5.EQ.-2) THEN + A(5)=0.0D0 + ELSE IF(KK5.EQ.-3) THEN + A(5)=2.0D0*SHARM(DZ,DZ,DIF(L,1),DIF(L,1))*VOL0/DZ + ENDIF +* + ELSE IF(IPR.EQ.1) THEN +* FORMULE DE DERIVEE. + IF(KK1.GT.0) THEN + A(1)=DHARM(DS,DS,DIF(L,1),DIF(MATN(KK1),1),DDF(L,1), + 1 DDF(MATN(KK1),1))*DZ*DT + ELSE IF(KK1.EQ.-1) THEN + A(1)=DHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK1.EQ.-2) THEN + A(1)=0.0D0 + ELSE IF(KK1.EQ.-3) THEN + A(1)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK2.GT.0) THEN + A(2)=DHARM(DS,DS,DIF(L,1),DIF(MATN(KK2),1),DDF(L,1), + 1 DDF(MATN(KK2),1))*DZ*DT + ELSE IF(KK2.EQ.-1) THEN + A(2)=DHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK2.EQ.-2) THEN + A(2)=0.0D0 + ELSE IF(KK2.EQ.-3) THEN + A(2)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK3.GT.0) THEN + A(3)=DHARM(DS,DS,DIF(L,1),DIF(MATN(KK3),1),DDF(L,1), + 1 DDF(MATN(KK3),1))*DZ*DT + ELSE IF(KK3.EQ.-1) THEN + A(3)=DHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK3.EQ.-2) THEN + A(3)=0.0D0 + ELSE IF(KK3.EQ.-3) THEN + A(3)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK4.GT.0) THEN + A(4)=DHARM(DZ,ZZ(KK4),DIF(L,3),DIF(MAT(KK4),3),DDF(L,3), + 1 DDF(MAT(KK4),3))*VOL0/DZ + ELSE IF(KK4.EQ.-1) THEN + A(4)=DHARM(DZ,DZ,DIF(L,3),DZ*QFR(7)/2.0,DDF(L,3),0.0) + 1 *VOL0/DZ + ELSE IF(KK4.EQ.-2) THEN + A(4)=0.0D0 + ELSE IF(KK4.EQ.-3) THEN + A(4)=2.0D0*DDF(L,3)*VOL0/(DZ*DZ) + ENDIF +* + IF(KK5.GT.0) THEN + A(5)=DHARM(DZ,ZZ(KK5),DIF(L,3),DIF(MAT(KK5),3),DDF(L,3), + 1 DDF(MAT(KK5),3))*VOL0/DZ + ELSE IF(KK5.EQ.-1) THEN + A(5)=DHARM(DZ,DZ,DIF(L,3),DZ*QFR(8)/2.0,DDF(L,3),0.0) + 1 *VOL0/DZ + ELSE IF(KK5.EQ.-2) THEN + A(5)=0.0D0 + ELSE IF(KK5.EQ.-3) THEN + A(5)=2.0D0*DDF(L,3)*VOL0/(DZ*DZ) + ENDIF +* + ELSE IF(IPR.GE.2) THEN +* FORMULE DE VARIATION. + IF(KK1.GT.0) THEN + A(1)=VHARM(DS,DS,DIF(L,1),DIF(MATN(KK1),1),DDF(L,1), + 1 DDF(MATN(KK1),1))*DZ*DT + ELSE IF(KK1.EQ.-1) THEN + A(1)=VHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK1.EQ.-2) THEN + A(1)=0.0D0 + ELSE IF(KK1.EQ.-3) THEN + A(1)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK2.GT.0) THEN + A(2)=VHARM(DS,DS,DIF(L,1),DIF(MATN(KK2),1),DDF(L,1), + 1 DDF(MATN(KK2),1))*DZ*DT + ELSE IF(KK2.EQ.-1) THEN + A(2)=VHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK2.EQ.-2) THEN + A(2)=0.0D0 + ELSE IF(KK2.EQ.-3) THEN + A(2)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK3.GT.0) THEN + A(3)=VHARM(DS,DS,DIF(L,1),DIF(MATN(KK3),1),DDF(L,1), + 1 DDF(MATN(KK3),1))*DZ*DT + ELSE IF(KK3.EQ.-1) THEN + A(3)=VHARM(DS,DS,DIF(L,1),DS*QFR(IQF)/2.0,DDF(L,1),0.0) + 1 *DZ*DT + ELSE IF(KK3.EQ.-2) THEN + A(3)=0.0D0 + ELSE IF(KK3.EQ.-3) THEN + A(3)=2.0D0*DDF(L,1)*DZ*DT/DS + ENDIF +* + IF(KK4.GT.0) THEN + A(4)=VHARM(DZ,ZZ(KK4),DIF(L,3),DIF(MAT(KK4),3),DDF(L,3), + 1 DDF(MAT(KK4),3))*VOL0/DZ + ELSE IF(KK4.EQ.-1) THEN + A(4)=VHARM(DZ,DZ,DIF(L,3),DZ*QFR(7)/2.0,DDF(L,3),0.0) + 1 *VOL0/DZ + ELSE IF(KK4.EQ.-2) THEN + A(4)=0.0D0 + ELSE IF(KK4.EQ.-3) THEN + A(4)=2.0D0*DDF(L,3)*VOL0/(DZ*DZ) + ENDIF +* + IF(KK5.GT.0) THEN + A(5)=VHARM(DZ,ZZ(KK5),DIF(L,3),DIF(MAT(KK5),3),DDF(L,3), + 1 DDF(MAT(KK5),3))*VOL0/DZ + ELSE IF(KK5.EQ.-1) THEN + A(5)=VHARM(DZ,DZ,DIF(L,3),DZ*QFR(8)/2.0,DDF(L,3),0.0) + 1 *VOL0/DZ + ELSE IF(KK5.EQ.-2) THEN + A(5)=0.0D0 + ELSE IF(KK5.EQ.-3) THEN + A(5)=2.0D0*DDF(L,3)*VOL0/(DZ*DZ) + ENDIF +* + ENDIF + RETURN + END |