Sources Esope | Cast3M

geoli3
C GEOLI3    SOURCE    GOUNAND   05/12/21    21:27:17     5281
      SUBROUTINE GEOLI3(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC,
     $     JDTJAC,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : GEOLI3
C PROJET      : Noyau linéaire NLIN
C DESCRIPTION : Calcul du jacobien dans le cas où la matrice
C               jacobienne A n'est pas carrée. A(n,m)
C               On calcule sqrt (det (transpose(A) * A))
C               Ceci est effectué pour chaque point de Gauss d'un
C               élément.
C
C               Pour calculer det (tA A), on utilise la formule suivante
C               déduite de celle de Binet-Cauchy :
C
C               det(tAA) = SUM (det(A_I))^2
C                           I
C
C               où I parcourt tous les sous-ensembles {i_1&lt;...&lt;i_m} dans
C               {1,...,N} et A_I est la matrice (m,m) obtenue de A en
C               retenant seulement ses lignes de numéro i_1,...,i_k.
C
C LANGAGE     : Fortran 77 (sauf E/S)
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : -
C APPELE PAR       : GEOLIN
C***********************************************************************
C ENTREES            : * IESREL (type entier) : dimension de l'espace
C                        réel (i.e. géométrique).
C                      * IESREF (type entier) : dimension de l'espace de
C                        référence.
C                      * NDNOEU (type entier) : nombre de ddl (par
C                        élément) de la transformation géométrique.
C                      * NDPOGO (type entier) : nombre de points
C                        d'intégration.
C                      * NDELEM (type entier) : nombre d'éléments du
C                        maillage élémentaire courant.
C                      * JMAJAC (type MCHEVA) : valeurs de la matrice
C                        jacobienne aux points de Gauss sur le maillage
C                        élémentaire courant.
C ENTREES/SORTIES    : * JDTJAC (type MCHEVA) : valeurs de
C                        sqrt(det(trans(J).J)) aux points de Gauss sur
C                        le maillage élémentaire courant.
C SORTIES            : -
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 12/08/99, version initiale
C HISTORIQUE : v1, 12/08/99, création
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
      INTEGER IESREF,IESREL,NDPOGO,NDELEM
      REAL*8 JMAJAC(IESREL,IESREF,NDPOGO,NDELEM)
      REAL*8 JDTJAC(NDPOGO,NDELEM)
*
      INTEGER IMPR,IRET
*
      REAL*8 ZERO
      PARAMETER (ZERO=0.D0)
*
      INTEGER IELEM,IPOGO,IREEL,JREEL
      REAL*8 DETCAR,DET
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli3'
      IF (IESREF.EQ.1) THEN
         DO 1 IELEM=1,NDELEM
            DO 12 IPOGO=1,NDPOGO
               DETCAR=ZERO
               DO 122 IREEL=1,IESREL
                  DETCAR=DETCAR+
     $                 (JMAJAC(IREEL,1,IPOGO,IELEM)
     $                 *JMAJAC(IREEL,1,IPOGO,IELEM))
 122           CONTINUE
               IF (DETCAR.GT.ZERO) THEN
                  JDTJAC(IPOGO,IELEM)=SQRT(DETCAR)
               ELSE
                  GOTO 9998
               ENDIF
 12         CONTINUE
 1       CONTINUE
      ELSEIF (IESREF.EQ.2) THEN
         DO 3 IELEM=1,NDELEM
            DO 32 IPOGO=1,NDPOGO
               DETCAR=ZERO
               DO 322 IREEL=1,IESREL-1
                  DO 3222 JREEL=IREEL+1,IESREL
                     DET=(JMAJAC(IREEL,1,IPOGO,IELEM)
     $                   *JMAJAC(JREEL,2,IPOGO,IELEM))
     $                  -(JMAJAC(JREEL,1,IPOGO,IELEM)
     $                   *JMAJAC(IREEL,2,IPOGO,IELEM))
                     DETCAR=DETCAR+(DET*DET)
 3222             CONTINUE
 322           CONTINUE
               IF (DETCAR.GT.ZERO) THEN
                  JDTJAC(IPOGO,IELEM)=SQRT(DETCAR)
               ELSE
                  GOTO 9998
               ENDIF
 32         CONTINUE
 3       CONTINUE
      ELSE
         WRITE(IOIMP,*) 'Je ne sais pas calculer la racine'
         WRITE(IOIMP,*) 'du det. de la matrice tAA de dimension '
         WRITE(IOIMP,*) 'IESREF=',IESREF
      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Format handling
*
*
* Error handling
*
 9998 CONTINUE
      WRITE(IOIMP,*) 'Déterminant de la matrice tAA nul ou négatif'
      WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
      GOTO 9999
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine geoli3'
      RETURN
*
* End of subroutine GEOLI3
*
      END