Sources Esope | Cast3M

geolin
C GEOLIN    SOURCE    GOUNAND   21/06/02    21:16:11     11022          
      SUBROUTINE GEOLIN(DFFPG,JCOOR,NBELEM,
     $     JMAJAC,JMIJAC,JDTJAC,LERJ,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : GEOLIN
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          : 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                      * NBELEM (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    : 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 SMCHAEL
      INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
      POINTEUR JCOOR.MCHEVA
      POINTEUR JMAJAC.MCHEVA
      POINTEUR JMIJAC.MCHEVA
      POINTEUR JDTJAC.MCHEVA
      POINTEUR JJTJ.MCHEVA
      POINTEUR JJTJM1.MCHEVA
* Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
      POINTEUR DFFPG.MCHEVA
*
      INTEGER NBELEM
      INTEGER IMPR,IRET
*
      LOGICAL LCARRE,LERJ
      INTEGER IESREF,IESREL,NDNOEU,NDCOL,NDPOGO,NDELEM
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolin'
      SEGACT DFFPG
      NDNOEU=DFFPG.WELCHE(/2)
      IESREF=DFFPG.WELCHE(/4)
      NDPOGO=DFFPG.WELCHE(/5)
      SEGACT JCOOR
      NDCOL=JCOOR.WELCHE(/2)
      IESREL=JCOOR.WELCHE(/4)
      NDELEM=JCOOR.WELCHE(/6)
      IF (NDCOL.NE.NDNOEU) THEN
         WRITE(IOIMP,*) 'Incompatibilité fns d''interpolation-géométrie'
         WRITE(IOIMP,*) 'NDNOEU=',NDNOEU,' NDCOL=',NDCOL
         GOTO 9999
      ENDIF
* On pourrait envisager un cas où tous les éléments seraient
* identiques (ex. carré) NDELEM=1
      IF (NDELEM.NE.NBELEM) THEN
         WRITE(IOIMP,*) 'Incompatibilité coordonnées-géométrie'
         WRITE(IOIMP,*) 'NDELEM=',NDELEM,' NBELEM=',NBELEM
         GOTO 9999
      ENDIF
      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=NDPOGO
      NBELM=NDELEM
      SEGINI JMAJAC
*
* On effectue le calcul de la matrice jacobienne aux points de Gauss
*
      CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM,
     $     DFFPG.WELCHE,JCOOR.WELCHE,
     $     JMAJAC.WELCHE,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      IF (IMPR.GT.4) THEN
         WRITE(IOIMP,*) 'On a créé',
     $   ' JMAJAC(élément , poi.gauss ,',
     $        ' dksi(i) , dx(j) ,1,1)'
         CALL PRCHVA(JMAJAC,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
      ENDIF
*
* 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=NDPOGO
         NBELM=NDELEM
         SEGINI JMIJAC
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=1
         NBPOI=NDPOGO
         NBELM=NDELEM
         SEGINI JDTJAC
         CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
     $        JMIJAC.WELCHE,JDTJAC.WELCHE,LERJ,
     $        IMPR,IRET)
         IF (IRET.NE.0) THEN
            IF (LERJ) GOTO 9666
            GOTO 9999
         ENDIF
         SEGDES JDTJAC
         SEGDES JMIJAC
         IF (IMPR.GT.4) THEN
            WRITE(IOIMP,*) 'On a créé',
     $           ' JMIJAC(élément , poi.gauss ,',
     $           ' dx(i) , dksi(j) ,1,1)'
            CALL PRCHVA(JMIJAC,IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
         ENDIF
      ELSE
*
* Sinon on calcule : SQRT ( det (transpose(A) * A))
* et la pseudo-inverse (JtJ)^-1 Jt
*
         NBLIG=1
         NBCOL=1
         N2LIG=IESREF
         N2COL=IESREL
         NBPOI=NDPOGO
         NBELM=NDELEM
         SEGINI JMIJAC
*
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=1
         NBPOI=NDPOGO
         NBELM=NDELEM
         SEGINI JDTJAC
* Objets temporaire
         NBLIG=1
         NBCOL=1
         N2LIG=IESREF
         N2COL=IESREF
         NBPOI=NDPOGO
         NBELM=NDELEM
         SEGINI JJTJ
         SEGINI JJTJM1
         CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
     $        JJTJ.WELCHE,JJTJM1.WELCHE,
     $        JMIJAC.WELCHE,JDTJAC.WELCHE,LERJ,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGSUP JJTJM1
         SEGSUP JJTJ
         SEGDES JMIJAC
         SEGDES JDTJAC
      ENDIF
      SEGDES JMAJAC
      SEGDES JCOOR
      SEGDES DFFPG
      IF (IMPR.GT.3) THEN
         WRITE(IOIMP,*) 'On a créé',
     $        ' JDTJAC(élément , poi.gauss ,',
     $        '1,1,1,1)'
         CALL PRCHVA(JDTJAC,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
      ENDIF
*
* 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 geolin'
      RETURN
*
* End of subroutine GEOLIN
*
      END