Esope source files

   1 : C TRINO6    SOURCE    CHAT      05/01/13    03:47:35     5004           
   2 :       SUBROUTINE TRINO6(ICM,IDIM,XE,XEC,XX,SHPC,SGA,QSI,IG)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)  
   5 :       DIMENSION A(2,6),AJ(2,2),X(2),XE(3,*),XEC(3,2),XX(3,2),
   6 :      & BJ(2,2),F(2),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(2),SIB(2)
   7 :       DATA EPSI,FR/1.D-7,4.D0/
   8 : C --- IG: NUMERO DU POINT DE GAUSS
   9 :       DO 10 I=1,IDIM
  10 :       A(I,1)=XE(I,1)-XX(I,IG)
  11 :       A(I,2)=-3.D0*XE(I,1)+4.D0*XE(I,2)-XE(I,3)
  12 :       A(I,3)=-3.D0*XE(I,1)-XE(I,5)+4.D0*XE(I,6)
  13 :       A(I,4)= 2.D0*XE(I,1)-4.D0*XE(I,2)+2.D0*XE(I,3)
  14 :       A(I,5)= 4.D0*(XE(I,1)-XE(I,2)+XE(I,4)-XE(I,6))
  15 :   10  A(I,6)= 2.D0*XE(I,1)+2.D0*XE(I,5)-4.D0*XE(I,6)
  16 : C --- Initialisation  
  17 :       DO 15 I=1,IDIM
  18 :   15  X(I)=0.D0
  19 : C
  20 :       ITER=0
  21 :       ICRI=0
  22 :   16  CONTINUE
  23 :  100  CONTINUE
  24 : C --- LES DERIVEES DES FONCTIONS DE TRANSFORMATION
  25 :       DO 20 I=1,IDIM
  26 :       AJ(I,1)=A(I,2)+2.0*A(I,4)*X(1)+A(I,5)*X(2)
  27 :   20  AJ(I,2)=A(I,3)+A(I,5)*X(1)+2.0*A(I,6)*X(2)
  28 : C --- Calcul de l'inverse de la matrice AJ et son determinant :BJ=inverse(AJ)
  28 :       CALL INVER2(AJ,BJ,DET,ICRI)
  29 :       DO 40 I=1,IDIM
  30 :   40  F(I)=A(I,1)+A(I,2)*X(1)+A(I,3)*X(2)+A(I,4)*X(1)*X(1)
  31 :      &    +A(I,5)*X(1)*X(2)+A(I,6)*X(2)*X(2)
  32 :       DD=SQRT(F(1)*F(1)+F(2)*F(2))
  33 :       IF(DD.GT.EPSI) THEN
  34 :       ITER=ITER+1
  35 :       DO 50 I=1,IDIM
  36 :   50  X(I)=X(I)-BJ(I,1)*F(1)-BJ(I,2)*F(2)
  37 :       GOTO 100
  38 :       END IF
  39 :       DO 60 I=1,IDIM
  40 :   60  QSI(I,IG)=X(I)
  41 :       C=1.D0-X(1)-X(2)
  42 : C --- Calcul des fonctions d'interpolation leur derivees
  43 : C --- Fonctions d'interpolation Ni
  44 :       SHPC(1,1,IG)=-C*(1.D0-2.D0*C)
  45 :       SHPC(1,2,IG)= FR*X(1)*C
  46 :       SHPC(1,3,IG)=-X(1)*(1.D0-2.D0*X(1))
  47 :       SHPC(1,4,IG)= FR*X(1)*X(2)
  48 :       SHPC(1,5,IG)=-X(2)*(1.D0-2.D0*X(2))
  49 :       SHPC(1,6,IG)= FR*X(2)*C
  50 :       IF(ICM.EQ.0) GOTO 1000
  51 : C --- Derivee de Ni par rapport a QSI
  52 :       SHPC(2,1,IG)= 1.D0-FR*C
  53 :       SHPC(2,2,IG)= FR*(C-X(1))
  54 :       SHPC(2,3,IG)=-1.D0+FR*X(1)
  55 :       SHPC(2,4,IG)= FR*X(2)
  56 :       SHPC(2,5,IG)= 0.D0
  57 :       SHPC(2,6,IG)=-FR*X(2)
  58 : C --- Derivee de Ni par rapport a ETA
  59 :       SHPC(3,1,IG)= SHPC(2,1,IG)
  60 :       SHPC(3,2,IG)=-FR*X(1)
  61 :       SHPC(3,3,IG)= 0.D0
  62 :       SHPC(3,4,IG)= FR*X(1)
  63 :       SHPC(3,5,IG)=-1.D0+FR*X(2)
  64 :       SHPC(3,6,IG)= FR*(C-X(2))
  65 : C
  66 :       DJ=0.0
  67 :       DO 70 I=1,IDIM
  68 :       DT(I)=0.5D0*(XEC(I,2)-XEC(I,1))
  69 :       SIB(I)=DT(I)
  70 :       DJ=DJ+DT(I)*DT(I)
  71 :       DO 70 J=1,IDIM
  72 :   70  BJ(I,J)=0.0
  73 :       DO 90 I=1,IDIM
  74 :       DO 90 J=1,IDIM
  75 :       DO 80 K=1,6
  76 :   80  BJ(I,J)=BJ(I,J)+XE(I,K)*SHPC(J+1,K,IG)
  77 :   90  CONTINUE
  78 : 
  79 : C --- LES DRIVEES
  79 :       CALL LINEQ1(BJ,DT,X)
  80 :       DJ=SQRT(DJ)
  81 :       DO 120 I=1,IDIM
  82 :  120  SIB(I)=SIB(I)/DJ
  83 : 1000  CONTINUE
  84 :       DO 2000 IS=1,IDIM
  85 : C --- Coordonnees du point d'integration dans le systeme reel (x,y,z)
  86 :       SGA(IS,1,IG)=XX(IS,IG)
  87 : C --- Coordonnees du point d'integration dans le systeme de reference 
  88 : C     (qsi,eta,dze)
  89 :       SGA(IS,2,IG)=QSI(IS,IG)
  90 : C --- Derivee de (qsi,eta) par rapport a miu
  91 :       SGA(IS,3,IG)=X(IS)
  92 : C --- Cosinus du cable avec les axes
  93 :       SGA(IS,4,IG)=SIB(IS)
  94 : C --- Jacobien
  95 :       SGA(IS,5,IG)=DJ
  96 :  2000 CONTINUE
  97 :       END
  98 :       
  99 :       
 100 : 
 101 :  
 102 :