Sources Esope | Cast3M

geolif
C GEOLIF    SOURCE    GOUNAND   21/06/02    21:16:09     11022          
      SUBROUTINE GEOLIF(DFFPG,JCOOR,SSFACT,NBELEF,
     $     JMAJAC,JMIJAC,JDTJAC,LERJ,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : GEOLIF
C PROJET      : Noyau linéaire NLIN
C DESCRIPTION : Calcul de la matrice jacobienne d'une transformation
C               géométrique aux points de Gauss d'un élément de
C               référence pour chaque élément réel.
C               On a :
C    Fonction f :      R^m      ->    R^n
C                 a=(a1...am)  |->  f(a)=(f1(a)...fn(a)
C      =>         matjac(i,j)(a)=Dj fi (a) = dfi / dxj (a)
C
C  Par exemple, pour une surface en 3D, la matrice jacobienne
C  s'exprime comme :
C
C                [  &lt;xn>  ]
C         [J] =  [  &lt;yn>  ]   .  [ { dNg/d(ksi) } { dNg/d(eta) } ]
C                [  &lt;zn>  ]
C
C        (3x2)   (3 x Nnoeuds)      (Nnoeuds x 2)
C
C    Ici, le nb de ddl est égal aux nbs de noeuds car l'interpolation
C    pour la géométrie est de type nodale.
C
C               Si la matrice jacobienne J est carrée, on calcule
C               également l'inverse de la matrice jacobienne j ainsi que
C               son déterminant : det J (appelé Jacobien).
C               Si la matrice jacobienne n'est pas carrée, le Jacobien
C               est calculé comme :
C                                 sqrt(det(transpose(J)*J))
C               Le but de ce sous-programme est d'effectuer toutes les
C               opérations de calcul dépendant uniquement de la
C               géométrie (ie de la variable : coordonnées des points).
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : RSET8  (copie de tableaux de CHARACTER*8)
C                    PRCHVA (impression d'un segment de type MCHEVA)
C                    GEOLI1 (calcul de la matrice jacobienne (fortran 77
C                            ))
C                    GEOLI2 (calcul de l'inverse et du déterminant de la
C                            matrice jacobienne (fortran 77))
C                    GEOLI3 (calcul de sqrt(det(transpose(J)*J)), J
C                            étant la matrice jacobienne (fortran 77))
C APPELE PAR       : NLIN
C***********************************************************************
C ENTREES            : * DFFPG (type MCHEVA) : valeurs des dérivées
C                        premières des fonctions d'interpolation de la
C                        transformation géométrique aux points de gauss
C                        sur l'élément de référence.
C                        Structure (cf.include SMCHAEL) :
C                        (1, nb. ddl, 1, dim.esp.réf, nb. poi. gauss, 1)
C                      * JCOOR  (type MCHEVA) : valeurs des ddl de la
C                        transformation géométrique sur le maillage
C                        élémentaire courant.
C                        Structure (cf.include SMCHAEL) :
C                        (1, nb. ddl, 1, dim. esp. réel,
C                         1, nb. éléments)
C                      * NBELEV (type entier) : nombre d'éléments du
C                        maillage élémentaire courant.
C ENTREES/SORTIES    : -
C SORTIES            : * JMAJAC (type MCHEVA) : valeurs de la matrice
C                        jacobienne aux points de Gauss sur le maillage
C                        élémentaire courant.
C                        Structure (cf.include SMCHAEL) :
C                        (1, 1, dim. esp. réel, dim. esp. référence,
C                         nb. poi. gauss, nb. éléments)
C                      * JMIJAC (type MCHEVA) : valeurs de l'inverse de
C                        la matrice jacobienne aux points de Gauss sur
C                        le maillage élémentaire courant.
C                        Structure (cf.include SMCHAEL) :
C                        (1, 1, dim. esp. référence, dim. esp. réel,
C                         nb. poi. gauss, nb. éléments)
C                        JMIJAC n'existe que si dim.esp.réf=dim.esp.réel
C                      * JDTJAC (type MCHEVA) : valeurs du déterminant
C                        de la matrice jacobienne aux points de Gauss
C                        sur le maillage élémentaire courant.
C                        Structure (cf.include SMCHAEL) :
C                        (1, 1, 1, 1, nb. poi. gauss, nb. éléments)
C                        Si la matrice jacobienne n'est pas carrée, on
C                        a calculé sqrt(det(trans(J).J))
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v2, 03/10/03, refonte complète (modif SMTNLIN)
C VERSION    : v1, 09/08/99, version initiale
C HISTORIQUE : v1, 09/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
-INC TNLIN      
*-INC SFACTIV
*-INC SMCHAEL
      INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
      POINTEUR JCOOR.MCHEVA
      POINTEUR JMAJAC.MCHEVA
      POINTEUR JMIJAC.MCHEVA
      POINTEUR JDTJAC.MCHEVA
* Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
      POINTEUR DFFPG.MCHEVA
*
      INTEGER NBELEV,NBELEF,NBELFV
      INTEGER IMPR,IRET
*
      LOGICAL LCARRE,LERJ
      INTEGER IESREF,IESREL,NBPOGO
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolif'
      SEGACT SSFACT
      NBELFV=SSFACT.LFACTI(/1)
      NBELEV=SSFACT.LFACTI(/2)
      SEGACT JCOOR
      NBLIG=JCOOR.WELCHE(/1)
      NBCOL=JCOOR.WELCHE(/2)
      N2LIG=JCOOR.WELCHE(/3)
      N2COL=JCOOR.WELCHE(/4)
      NBPOI=JCOOR.WELCHE(/5)
      NBELM=JCOOR.WELCHE(/6)
* On pourrait envisager un cas où tous les éléments seraient
* identiques (ex. carré) NBELM=1
      IF (NBLIG.NE.1.OR.N2LIG.NE.1
     $     .OR.NBPOI.NE.1.OR.NBELM.NE.NBELEV) THEN
         WRITE(IOIMP,*) 'Erreur dims JCOOR'
         GOTO 9999
      ENDIF
      NDDL=NBCOL
      IESREL=N2COL
      SEGACT DFFPG
      NBLIG=DFFPG.WELCHE(/1)
      NBCOL=DFFPG.WELCHE(/2)
      N2LIG=DFFPG.WELCHE(/3)
      N2COL=DFFPG.WELCHE(/4)
      NBPOI=DFFPG.WELCHE(/5)
      NBELM=DFFPG.WELCHE(/6)
      IF (NBLIG.NE.1.OR.NBCOL.NE.NDDL.OR.N2LIG.NE.1
     $     .OR.(NBELM.NE.NBELFV.AND.NBELM.NE.1)) THEN
         WRITE(IOIMP,*) 'Erreur dims DFFPG'
         GOTO 9999
      ENDIF
      IESREF=N2COL
      NBPOGO=NBPOI
      NLFVDF=NBELM
      IF (IESREF.LT.IESREL) THEN
         LCARRE=.FALSE.
      ELSEIF (IESREF.EQ.IESREL) THEN
         LCARRE=.TRUE.
      ELSE
         WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
         WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
         GOTO 9999
      ENDIF
*
* Initialisations...
*
      NBLIG=1
      NBCOL=1
      N2LIG=IESREL
      N2COL=IESREF
      NBPOI=NBPOGO
      NBELM=NBELEF
      SEGINI JMAJAC
*      SEGPRT,DFFPG
*
* On effectue le calcul de la matrice jacobienne aux points de Gauss
*
      CALL GEOLF1(IESREL,IESREF,NDDL,NBPOGO,NBELEV,NBELFV,NBELEF,
     $     NLFVDF,
     $     DFFPG.WELCHE,JCOOR.WELCHE,SSFACT.LFACTI,
     $     JMAJAC.WELCHE,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*      SEGPRT,JMAJAC
*
* Si la matrice est carrée, on tente de calculer
* son inverse et son déterminant...
*
      IF (LCARRE) THEN
         NBLIG=1
         NBCOL=1
         N2LIG=IESREF
         N2COL=IESREL
         NBPOI=NBPOGO
         NBELM=NBELEF
         SEGINI JMIJAC
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=1
         NBPOI=NBPOGO
         NBELM=NBELEF
         SEGINI JDTJAC
         CALL GEOLI2(IESREF,NBPOGO,NBELEF,JMAJAC.WELCHE,
     $        JMIJAC.WELCHE,JDTJAC.WELCHE,LERJ,
     $        IMPR,IRET)
*         SEGPRT,JDTJAC
         IF (IRET.NE.0) THEN
            IF (LERJ) GOTO 9666
            GOTO 9999
         ENDIF
         SEGDES JDTJAC
         SEGDES JMIJAC
      ELSE
*
* Sinon on calcule : SQRT ( det (transpose(A) * A))
*
* Pas de matrice inverse
         JMIJAC=0
*
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=1
         NBPOI=NBPOGO
         NBELM=NBELEF
         SEGINI JDTJAC
         CALL GEOLI3(IESREL,IESREF,NBPOGO,NBELEF,JMAJAC.WELCHE,
     $        JDTJAC.WELCHE,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGDES JDTJAC
      ENDIF
      SEGDES JMAJAC
      SEGDES SSFACT
      SEGDES JCOOR
      SEGDES DFFPG
*
* Normal termination
*
      IRET=0
      RETURN
*
* Format handling
*
*
* Error handling
*
 9666 CONTINUE
      IRET=666
      RETURN
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine geolif'
      RETURN
*
* End of subroutine GEOLIF
*
      END