Esope source files

   1 : C GEOLI4    SOURCE    GOUNAND   06/11/24    21:15:01     5574           
   2 :       SUBROUTINE GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC,
   3 :      $     JJTJ,JJTJM1,
   4 :      $     JMIJAC,JDTJAC,
   5 :      $     IMPR,IRET)
   6 :       IMPLICIT REAL*8 (A-H,O-Z)
   7 :       IMPLICIT INTEGER (I-N)
   8 : C***********************************************************************
   9 : C NOM         : GEOLI4
  10 : C PROJET      : Noyau linéaire NLIN
  11 : C DESCRIPTION : Calcul du jacobien dans le cas où la matrice 
  12 : C               jacobienne A n'est pas carrée. A(n,m)
  13 : C               On calcule sqrt (det (transpose(A) * A))
  14 : C               Ceci est effectué pour chaque point de Gauss d'un
  15 : C               élément.
  16 : C               On calcule aussi le pseudo-inverse de J
  17 : C             
  18 : C LANGAGE     : Fortran 77 (sauf E/S)
  19 : C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  20 : C               mél : gounand@semt2.smts.cea.fr
  21 : C***********************************************************************
  22 : C APPELES          : -
  23 : C APPELE PAR       : GEOLIN
  24 : C***********************************************************************
  25 : C ENTREES            : * IESREL (type entier) : dimension de l'espace
  26 : C                        réel (i.e. géométrique).
  27 : C                      * IESREF (type entier) : dimension de l'espace de
  28 : C                        référence.
  29 : C                      * NDNOEU (type entier) : nombre de ddl (par
  30 : C                        élément) de la transformation géométrique.
  31 : C                      * NDPOGO (type entier) : nombre de points
  32 : C                        d'intégration.
  33 : C                      * NDELEM (type entier) : nombre d'éléments du
  34 : C                        maillage élémentaire courant.
  35 : C                      * JMAJAC (type MCHEVA) : valeurs de la matrice
  36 : C                        jacobienne aux points de Gauss sur le maillage
  37 : C                        élémentaire courant. 
  38 : C ENTREES/SORTIES    : * JDTJAC (type MCHEVA) : valeurs de
  39 : C                        sqrt(det(trans(J).J)) aux points de Gauss sur
  40 : C                        le maillage élémentaire courant.
  41 : C SORTIES            : -
  42 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  43 : C***********************************************************************
  44 : C VERSION    : v1, 26/10/06, version initiale
  45 : C HISTORIQUE : v1, 26/10/06, création
  46 : C HISTORIQUE :
  47 : C HISTORIQUE :
  48 : C***********************************************************************
  49 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
  50 : C en cas de modification de ce sous-programme afin de faciliter
  51 : C la maintenance !
  52 : C***********************************************************************
  53 : -INC CCOPTIO
  54 :       INTEGER IESREF,IESREL,NDPOGO,NDELEM
  55 :       REAL*8 JMAJAC(IESREL,IESREF,NDPOGO,NDELEM)
  56 :       REAL*8 JMIJAC(IESREF,IESREL,NDPOGO,NDELEM)
  57 :       REAL*8 JJTJ(IESREF,IESREF,NDPOGO,NDELEM)
  58 :       REAL*8 JJTJM1(IESREF,IESREF,NDPOGO,NDELEM)
  59 :       REAL*8 JDTJAC(NDPOGO,NDELEM)
  60 : *
  61 :       INTEGER IMPR,IRET
  62 : *
  63 :       INTEGER IELEM,IPOGO,IKF,IIF,IJL,IKL,IJF
  64 :       LOGICAL LERJ
  65 : *     
  66 : * Executable statements
  67 : *
  68 :       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli4'
  69 : *
  70 : * Calcul de JtJ
  71 : *
  72 :       DO IELEM=1,NDELEM
  73 :          DO IPOGO=1,NDPOGO
  74 :             DO IKF=1,IESREF
  75 :                DO IIF=1,IESREF
  76 :                   DO IJL=1,IESREL
  77 :                      JJTJ(IIF,IKF,IPOGO,IELEM)=
  78 :      $                    JJTJ(IIF,IKF,IPOGO,IELEM)+
  79 :      $                    (JMAJAC(IJL,IIF,IPOGO,IELEM)
  80 :      $                    *JMAJAC(IJL,IKF,IPOGO,IELEM))
  81 :                   ENDDO
  82 :                ENDDO
  83 :             ENDDO
  84 :          ENDDO
  85 :       ENDDO
  86 : *
  87 : * Inverse et déterminant
  88 : *
  88 :       CALL GEOLI2(IESREF,NDPOGO,NDELEM,JJTJ,
  89 :      $     JJTJM1,JDTJAC,LERJ,
  90 :      $     IMPR,IRET)
  91 :       IF (IRET.NE.0) THEN
  92 :          IF (LERJ) THEN
  93 :             WRITE(IOIMP,*) 'Erreur totalement anormale'
  94 :          ENDIF
  95 :          GOTO 9999
  96 :       ENDIF
  97 : *
  98 : * Pseudo-inverse
  99 : *      
 100 :       DO IELEM=1,NDELEM
 101 :          DO IPOGO=1,NDPOGO
 102 :             DO IKL=1,IESREL
 103 :                DO IIF=1,IESREF
 104 :                   DO IJF=1,IESREF
 105 :                      JMIJAC(IIF,IKL,IPOGO,IELEM)=
 106 :      $                    JMIJAC(IIF,IKL,IPOGO,IELEM)+
 107 :      $                    (JJTJM1(IIF,IJF,IPOGO,IELEM)
 108 :      $                    *JMAJAC(IKL,IJF,IPOGO,IELEM))
 109 :                   ENDDO
 110 :                ENDDO
 111 :             ENDDO
 112 :          ENDDO
 113 :       ENDDO
 114 : *
 115 : * Racine du déterminant
 116 : *
 117 :       DO IELEM=1,NDELEM
 118 :          DO IPOGO=1,NDPOGO
 119 :             JDTJAC(IPOGO,IELEM)=SQRT(JDTJAC(IPOGO,IELEM))
 120 :          ENDDO
 121 :       ENDDO
 122 : *
 123 : * Normal termination
 124 : *
 125 :       IRET=0
 126 :       RETURN
 127 : *
 128 : * Format handling
 129 : *
 130 : *
 131 : * Error handling
 132 : *
 133 :  9998 CONTINUE
 134 :       WRITE(IOIMP,*) 'Déterminant de la matrice tAA nul ou négatif'
 135 :       WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
 136 :       GOTO 9999
 137 :  9999 CONTINUE
 138 :       IRET=1
 139 :       WRITE(IOIMP,*) 'An error was detected in subroutine geoli4'
 140 :       RETURN
 141 : *
 142 : * End of subroutine GEOLI4
 143 : *
 144 :       END
 145 : 
 146 :  
 147 :