Télécharger la plaquette Cast3M


fr en

Télécharger geoli2.eso

Retour à la liste

Numérotation des lignes : 

geoli2

  1. C GEOLI2    SOURCE    GOUNAND   26/01/09    21:15:23     12441          
  2.       SUBROUTINE GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC,
  3.      $     JMIJAC,JDTJAC,
  4.      $     IMPR,IRET)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM         : GEOLI2
  9. C PROJET      : Noyau linéaire NLIN
  10. C DESCRIPTION : Calcul de l'inverse et du déterminant d'une matrice
  11. C               jacobienne dans le cas où celle-ci est carrée.
  12. C               Ceci est effectué pour chaque point de Gauss d'un
  13. C               élément.
  14. C               On teste également si le déterminant s'annule ou change
  15. C               de signe sur l'élément auxquels cas on génère une
  16. C               erreur.
  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            : * IESREF (type entier) : dimension de l'espace de
  26. C                        référence.
  27. C                      * NDPOGO (type entier) : nombre de points
  28. C                        d'intégration.
  29. C                      * NDELEM (type entier) : nombre d'éléments du
  30. C                        maillage élémentaire courant.
  31. C                      * JMAJAC (type MCHEVA) : valeurs de la matrice
  32. C                        jacobienne aux points de Gauss sur le maillage
  33. C                        élémentaire courant.
  34. C ENTREES/SORTIES    : * JMIJAC (type MCHEVA) : valeurs de l'inverse de
  35. C                        la matrice jacobienne aux points de Gauss sur
  36. C                        le maillage élémentaire courant.
  37. C                        JMIJAC n'existe que si dim.esp.réf=dim.esp.réel
  38. C                      * JDTJAC (type MCHEVA) : valeurs du déterminant
  39. C                        de la matrice jacobienne aux points de Gauss
  40. C                        sur le maillage élémentaire courant.
  41. C                        Si la matrice jacobienne n'est pas carrée, on
  42. C                        a calculé sqrt(det(trans(J).J))
  43. C SORTIES            : -
  44. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  45. C***********************************************************************
  46. C VERSION    : v1, 10/08/99, version initiale
  47. C HISTORIQUE : v1, 10/08/99, création
  48. C HISTORIQUE :
  49. C HISTORIQUE :
  50. C***********************************************************************
  51. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  52. C en cas de modification de ce sous-programme afin de faciliter
  53. C la maintenance !
  54. C***********************************************************************
  55. -INC CCREEL
  56.  
  57. -INC PPARAM
  58. -INC CCOPTIO
  59.       INTEGER IESREF,NDPOGO,NDELEM
  60.       REAL*8 JMAJAC(IESREF,IESREF,NDPOGO,NDELEM)
  61.       REAL*8 JMIJAC(IESREF,IESREF,NDPOGO,NDELEM)
  62.       REAL*8 JDTJAC(NDPOGO,NDELEM)
  63. *
  64.       INTEGER IMPR,IRET
  65. *
  66.       REAL*8 ZERO,UN
  67.       PARAMETER (ZERO=0.D0,UN=1.D0)
  68. *
  69.       REAL*8 DET,INVDET
  70.       INTEGER IELEM,IPOGO,IREFER,JREFER
  71.       LOGICAL LSIGN,LDETZ,LCHSIG
  72. *
  73. * Executable statements
  74. *
  75. *     La ligne suivante est inutile mais le compilateur du jour bugge si
  76. *     on ne la mat pas
  77.       LSIGN=.FALSE.
  78.       LDETZ=.FALSE.
  79.       LCHSIG=.FALSE.
  80.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli2'
  81.       IF (IESREF.EQ.1) THEN
  82.          DO 1 IELEM=1,NDELEM
  83.             DO 12 IPOGO=1,NDPOGO
  84.                DET=JMAJAC(1,1,IPOGO,IELEM)
  85.                JDTJAC(IPOGO,IELEM)=DET
  86. *               IF (DET.NE.ZERO) THEN
  87.                IF (ABS(DET).GT.(SQRT(XPETIT))) THEN
  88.                   IF (IPOGO.EQ.1) THEN
  89.                      LSIGN=(DET.GT.ZERO)
  90.                   ELSE
  91.                      IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE.
  92.                   ENDIF
  93.                   INVDET=UN/DET
  94.                   JMIJAC(1,1,IPOGO,IELEM)=INVDET
  95.                ELSE
  96.                   LDETZ=.TRUE.
  97.                ENDIF
  98.  12         CONTINUE
  99.  1       CONTINUE
  100.       ELSEIF (IESREF.EQ.2) THEN
  101.          DO 3 IELEM=1,NDELEM
  102.             DO 32 IPOGO=1,NDPOGO
  103.                DET=(JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  104.      $           - (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,2,IPOGO,IELEM))
  105.                JDTJAC(IPOGO,IELEM)=DET
  106. *               IF (DET.NE.ZERO) THEN
  107.                IF (ABS(DET).GT.(SQRT(XPETIT))) THEN
  108.                   IF (IPOGO.EQ.1) THEN
  109.                      LSIGN=(DET.GT.ZERO)
  110.                   ELSE
  111.                      IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE.
  112.                   ENDIF
  113.                   INVDET=UN/DET
  114.                   JMIJAC(1,1,IPOGO,IELEM)=
  115.      $                 (+JMAJAC(2,2,IPOGO,IELEM))*INVDET
  116.                   JMIJAC(2,1,IPOGO,IELEM)=
  117.      $                 (-JMAJAC(2,1,IPOGO,IELEM))*INVDET
  118.                   JMIJAC(1,2,IPOGO,IELEM)=
  119.      $                 (-JMAJAC(1,2,IPOGO,IELEM))*INVDET
  120.                   JMIJAC(2,2,IPOGO,IELEM)=
  121.      $                 (+JMAJAC(1,1,IPOGO,IELEM))*INVDET
  122.                ELSE
  123.                   LDETZ=.TRUE.
  124.                ENDIF
  125.  32         CONTINUE
  126.  3       CONTINUE
  127.       ELSEIF (IESREF.EQ.3) THEN
  128.          DO 5 IELEM=1,NDELEM
  129.             DO 52 IPOGO=1,NDPOGO
  130.                JMIJAC(1,1,IPOGO,IELEM)=
  131.      $              (JMAJAC(2,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  132.      $             -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  133.                JMIJAC(2,1,IPOGO,IELEM)=
  134.      $              (JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  135.      $             -(JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  136.                JMIJAC(3,1,IPOGO,IELEM)=
  137.      $              (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM))
  138.      $             -(JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  139.                JMIJAC(1,2,IPOGO,IELEM)=
  140.      $              (JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM))
  141.      $             -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  142.                JMIJAC(2,2,IPOGO,IELEM)=
  143.      $              (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  144.      $             -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM))
  145.                JMIJAC(3,2,IPOGO,IELEM)=
  146.      $              (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM))
  147.      $             -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM))
  148.                JMIJAC(1,3,IPOGO,IELEM)=
  149.      $              (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  150.      $             -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  151.                JMIJAC(2,3,IPOGO,IELEM)=
  152.      $              (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,3,IPOGO,IELEM))
  153.      $             -(JMAJAC(2,3,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM))
  154.                JMIJAC(3,3,IPOGO,IELEM)=
  155.      $              (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  156.      $             -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,1,IPOGO,IELEM))
  157.                DET=(JMAJAC(1,1,IPOGO,IELEM)*JMIJAC(1,1,IPOGO,IELEM))
  158.      $            +(JMAJAC(1,2,IPOGO,IELEM)*JMIJAC(2,1,IPOGO,IELEM))
  159.      $            +(JMAJAC(1,3,IPOGO,IELEM)*JMIJAC(3,1,IPOGO,IELEM))
  160.                JDTJAC(IPOGO,IELEM)=DET
  161. *               IF (DET.NE.ZERO) THEN
  162.                IF (ABS(DET).GT.(SQRT(XPETIT))) THEN
  163.                   IF (IPOGO.EQ.1) THEN
  164.                      LSIGN=(DET.GT.ZERO)
  165.                   ELSE
  166.                      IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE.
  167.                   ENDIF
  168.                   INVDET=UN/DET
  169.                   DO 522 JREFER=1,IESREF
  170.                      DO 5222 IREFER=1,IESREF
  171.                         JMIJAC(IREFER,JREFER,IPOGO,IELEM)=
  172.      $                       JMIJAC(IREFER,JREFER,IPOGO,IELEM)
  173.      $                       * INVDET
  174.  5222                CONTINUE
  175.  522              CONTINUE
  176.                ELSE
  177.                   LDETZ=.TRUE.
  178.                ENDIF
  179.  52         CONTINUE
  180.  5       CONTINUE
  181.       ELSE
  182.          WRITE(IOIMP,*) 'Je ne sais pas calculer l''inverse d''une'
  183.          WRITE(IOIMP,*) 'matrice jacobienne de dimension '
  184.          WRITE(IOIMP,*) 'IESREF=',IESREF
  185.       ENDIF
  186.       IF (LDETZ) GOTO 9998
  187.       IF (LCHSIG) GOTO 9997
  188. *
  189. * Normal termination
  190. *
  191.       IRET=0
  192.       RETURN
  193. *
  194. * Format handling
  195. *
  196. *
  197. * Error handling
  198. *
  199.  9997 CONTINUE
  200. *      IF (LERJ) THEN
  201.          IRET=-666
  202. *         CALL ECRENT(IRET)
  203.          RETURN
  204. *      ENDIF
  205. *      WRITE(IOIMP,*) 'Le déterminant de la matrice jacobienne'
  206. *      WRITE(IOIMP,*) 'change de signe sur l''élément.'
  207. *      WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
  208. *      GOTO 9999
  209.  9998 CONTINUE
  210. *      IF (LERJ) THEN
  211.          IRET=-667
  212. *         CALL ECRENT(IRET)
  213.          RETURN
  214. *      ENDIF
  215. *      WRITE(IOIMP,*) 'Déterminant de la matrice jacobienne nul'
  216. *      WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
  217. *      GOTO 9999
  218.  9999 CONTINUE
  219.       IRET=1
  220.       WRITE(IOIMP,*) 'An error was detected in subroutine geoli2'
  221.       RETURN
  222. *
  223. * End of subroutine GEOLI2
  224. *
  225.       END
  226.  
  227.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales