Télécharger la plaquette Cast3M


fr en

Télécharger geolik.eso

Retour à la liste

Numérotation des lignes : 

geolik

  1. C GEOLIK    SOURCE    GOUNAND   26/01/09    21:15:26     12441          
  2.       SUBROUTINE GEOLIK(DFFPG,JCOOR,NBELEM,
  3.      $     JMAJAC,JMIJAC,JDTJAC,IMREG,
  4.      $     IMPR,IRET)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM         : GEOLIK
  9. C PROJET      : Noyau linéaire NLIN
  10. C DESCRIPTION : Calcul de la matrice jacobienne d'une transformation
  11. C               géométrique aux points de Gauss d'un élément de
  12. C               référence pour chaque élément réel.
  13. C               On a :
  14. C    Fonction f :      R^m      ->    R^n
  15. C                 a=(a1...am)  |->  f(a)=(f1(a)...fn(a)
  16. C      =>         matjac(i,j)(a)=Dj fi (a) = dfi / dxj (a)
  17. C
  18. C  Par exemple, pour une surface en 3D, la matrice jacobienne
  19. C  s'exprime comme :
  20. C
  21. C                [  <xn>  ]
  22. C         [J] =  [  <yn>  ]   .  [ { dNg/d(ksi) } { dNg/d(eta) } ]
  23. C                [  <zn>  ]
  24. C
  25. C        (3x2)   (3 x Nnoeuds)      (Nnoeuds x 2)
  26. C
  27. C    Ici, le nb de ddl est égal aux nbs de noeuds car l'interpolation
  28. C    pour la géométrie est de type nodale.
  29. C
  30. C               Si la matrice jacobienne J est carrée, on calcule
  31. C               également l'inverse de la matrice jacobienne j ainsi que
  32. C               son déterminant : det J (appelé Jacobien).
  33. C               Si la matrice jacobienne n'est pas carrée, le Jacobien
  34. C               est calculé comme :
  35. C                                 sqrt(det(transpose(J)*J))
  36. C               Le but de ce sous-programme est d'effectuer toutes les
  37. C               opérations de calcul dépendant uniquement de la
  38. C               géométrie (ie de la variable : coordonnées des points).
  39. C
  40. C LANGAGE     : ESOPE
  41. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  42. C               mél : gounand@semt2.smts.cea.fr
  43. C***********************************************************************
  44. C APPELES          : PRCHVA (impression d'un segment de type MCHEVA)
  45. C                    GEOLI1 (calcul de la matrice jacobienne (fortran 77
  46. C                            ))
  47. C                    GEOLI2 (calcul de l'inverse et du déterminant de la
  48. C                            matrice jacobienne (fortran 77))
  49. C                    GEOLI3 (calcul de sqrt(det(transpose(J)*J)), J
  50. C                            étant la matrice jacobienne (fortran 77))
  51. C APPELE PAR       : NLIN
  52. C***********************************************************************
  53. C ENTREES            : * DFFPG (type MCHEVA) : valeurs des dérivées
  54. C                        premières des fonctions d'interpolation de la
  55. C                        transformation géométrique aux points de gauss
  56. C                        sur l'élément de référence.
  57. C                        Structure (cf.include SMCHAEL) :
  58. C                        (1, nb. ddl, 1, dim.esp.réf, nb. poi. gauss, 1)
  59. C                      * JCOOR  (type MCHEVA) : valeurs des ddl de la
  60. C                        transformation géométrique sur le maillage
  61. C                        élémentaire courant.
  62. C                        Structure (cf.include SMCHAEL) :
  63. C                        (1, nb. ddl, 1, dim. esp. réel,
  64. C                         1, nb. éléments)
  65. C                      * NBELEM (type entier) : nombre d'éléments du
  66. C                        maillage élémentaire courant.
  67. C ENTREES/SORTIES    : -
  68. C SORTIES            : * JMAJAC (type MCHEVA) : valeurs de la matrice
  69. C                        jacobienne aux points de Gauss sur le maillage
  70. C                        élémentaire courant.
  71. C                        Structure (cf.include SMCHAEL) :
  72. C                        (1, 1, dim. esp. réel, dim. esp. référence,
  73. C                         nb. poi. gauss, nb. éléments)
  74. C                      * JMIJAC (type MCHEVA) : valeurs de l'inverse de
  75. C                        la matrice jacobienne aux points de Gauss sur
  76. C                        le maillage élémentaire courant.
  77. C                        Structure (cf.include SMCHAEL) :
  78. C                        (1, 1, dim. esp. référence, dim. esp. réel,
  79. C                         nb. poi. gauss, nb. éléments)
  80. C                        JMIJAC n'existe que si dim.esp.réf=dim.esp.réel
  81. C                      * JDTJAC (type MCHEVA) : valeurs du déterminant
  82. C                        de la matrice jacobienne aux points de Gauss
  83. C                        sur le maillage élémentaire courant.
  84. C                        Structure (cf.include SMCHAEL) :
  85. C                        (1, 1, 1, 1, nb. poi. gauss, nb. éléments)
  86. C                        Si la matrice jacobienne n'est pas carrée, on
  87. C                        a calculé sqrt(det(trans(J).J))
  88. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  89. C***********************************************************************
  90. C VERSION    : v1, 09/08/99, version initiale
  91. C HISTORIQUE : v1, 09/08/99, création
  92. C HISTORIQUE : 2025/12/24 : JMIJAC=-666 si det nul ou -667 si chgt
  93. C              signe sur un element
  94. C HISTORIQUE :
  95. C***********************************************************************
  96. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  97. C en cas de modification de ce sous-programme afin de faciliter
  98. C la maintenance !
  99. C***********************************************************************
  100. -INC PPARAM
  101. -INC CCOPTIO
  102. -INC TNLIN
  103. *-INC SMCHAEL
  104.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
  105.       POINTEUR JCOOR.MCHEVA
  106.       POINTEUR JMAJAC.MCHEVA
  107.       POINTEUR JMIJAC.MCHEVA
  108.       POINTEUR JDTJAC.MCHEVA
  109.       POINTEUR JJTJ.MCHEVA
  110.       POINTEUR JJTJM1.MCHEVA
  111. * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
  112.       POINTEUR DFFPG.MCHEVA
  113. *
  114.       INTEGER NBELEM
  115.       INTEGER IMPR,IRET
  116. *
  117.       LOGICAL LCARRE
  118.       INTEGER IESREF,IESREL,NDNOEU,NDCOL,NDPOGO,NDELEM
  119. *
  120. * Executable statements
  121. *
  122.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolik'
  123.       SEGACT DFFPG
  124.       NDNOEU=DFFPG.WELCHE(/2)
  125.       IESREF=DFFPG.WELCHE(/4)
  126.       NDPOGO=DFFPG.WELCHE(/5)
  127.       SEGACT JCOOR
  128.       NDCOL=JCOOR.WELCHE(/2)
  129.       IESREL=JCOOR.WELCHE(/4)
  130.       NDELEM=JCOOR.WELCHE(/6)
  131.       IF (NDCOL.NE.NDNOEU) THEN
  132.          WRITE(IOIMP,*) 'Incompatibilité fns d''interpolation-géométrie'
  133.          WRITE(IOIMP,*) 'NDNOEU=',NDNOEU,' NDCOL=',NDCOL
  134.          GOTO 9999
  135.       ENDIF
  136. * On pourrait envisager un cas où tous les éléments seraient
  137. * identiques (ex. carré) NDELEM=1
  138.       IF (NDELEM.NE.NBELEM) THEN
  139.          WRITE(IOIMP,*) 'Incompatibilité coordonnées-géométrie'
  140.          WRITE(IOIMP,*) 'NDELEM=',NDELEM,' NBELEM=',NBELEM
  141.          GOTO 9999
  142.       ENDIF
  143.       IF (IMREG.EQ.1) THEN
  144.          NDELEM=1
  145.          NDPOGO=1
  146.       ENDIF
  147.       IF (IESREF.LT.IESREL) THEN
  148.          LCARRE=.FALSE.
  149.       ELSEIF (IESREF.EQ.IESREL) THEN
  150.          LCARRE=.TRUE.
  151.       ELSE
  152.          WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  153.          WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  154.          GOTO 9999
  155.       ENDIF
  156. *
  157. * Initialisations...
  158. *
  159.       NBLIG=1
  160.       NBCOL=1
  161.       N2LIG=IESREL
  162.       N2COL=IESREF
  163.       NBPOI=NDPOGO
  164.       NBELM=NDELEM
  165.       SEGINI JMAJAC
  166. *
  167. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  168. *
  169.       CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM,
  170.      $     DFFPG.WELCHE,JCOOR.WELCHE,
  171.      $     JMAJAC.WELCHE,
  172.      $     IMPR,IRET)
  173.       IF (IRET.NE.0) GOTO 9999
  174.       IF (IMPR.GT.4) THEN
  175.          WRITE(IOIMP,*) 'On a créé',
  176.      $   ' JMAJAC(élément , poi.gauss ,',
  177.      $        ' dksi(i) , dx(j) ,1,1)'
  178.          CALL PRCHVA(JMAJAC,IMPR,IRET)
  179.          IF (IRET.NE.0) GOTO 9999
  180.       ENDIF
  181. *
  182. * Si la matrice est carrée, on tente de calculer
  183. * son inverse et son déterminant...
  184. *
  185.       IF (LCARRE) THEN
  186.          NBLIG=1
  187.          NBCOL=1
  188.          N2LIG=IESREF
  189.          N2COL=IESREL
  190.          NBPOI=NDPOGO
  191.          NBELM=NDELEM
  192.          SEGINI JMIJAC
  193.          NBLIG=1
  194.          NBCOL=1
  195.          N2LIG=1
  196.          N2COL=1
  197.          NBPOI=NDPOGO
  198.          NBELM=NDELEM
  199.          SEGINI JDTJAC
  200.          CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  201.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,
  202.      $        IMPR,IRET)
  203.          IF (IRET.LT.0) THEN
  204.             SEGSUP JMIJAC
  205.             JMIJAC=0
  206.          ELSEIF (IRET.GT.0) THEN
  207.             GOTO 9999
  208.          ENDIF
  209. *         SEGDES JDTJAC
  210. *         SEGDES JMIJAC
  211.          IF (IMPR.GT.4) THEN
  212.             IF (JMIJAC.GT.0) THEN
  213.                WRITE(IOIMP,*) 'On a créé',
  214.      $              ' JMIJAC(élément , poi.gauss ,',
  215.      $              ' dx(i) , dksi(j) ,1,1)'
  216.                CALL PRCHVA(JMIJAC,IMPR,IRET)
  217.                IF (IRET.NE.0) GOTO 9999
  218.             ENDIF
  219.          ENDIF
  220.       ELSE
  221. *
  222. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  223. * et la pseudo-inverse (JtJ)^-1 Jt
  224. *
  225.          NBLIG=1
  226.          NBCOL=1
  227.          N2LIG=IESREF
  228.          N2COL=IESREL
  229.          NBPOI=NDPOGO
  230.          NBELM=NDELEM
  231.          SEGINI JMIJAC
  232. *
  233.          NBLIG=1
  234.          NBCOL=1
  235.          N2LIG=1
  236.          N2COL=1
  237.          NBPOI=NDPOGO
  238.          NBELM=NDELEM
  239.          SEGINI JDTJAC
  240. * Objets temporaire
  241.          NBLIG=1
  242.          NBCOL=1
  243.          N2LIG=IESREF
  244.          N2COL=IESREF
  245.          NBPOI=NDPOGO
  246.          NBELM=NDELEM
  247.          SEGINI JJTJ
  248.          SEGINI JJTJM1
  249.          CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  250.      $        JJTJ.WELCHE,JJTJM1.WELCHE,
  251.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,
  252.      $        IMPR,IRET)
  253.          IF (IRET.LT.0) THEN
  254.             SEGSUP JMIJAC
  255.             JMIJAC=0
  256.          ELSEIF (IRET.GT.0) THEN
  257.             GOTO 9999
  258.          ENDIF
  259.          SEGSUP JJTJM1
  260.          SEGSUP JJTJ
  261. *         SEGDES JMIJAC
  262. *         SEGDES JDTJAC
  263.       ENDIF
  264. *      SEGDES JMAJAC
  265. *      SEGDES JCOOR
  266. *      SEGDES DFFPG
  267.       IF (IMPR.GT.3) THEN
  268. *      IF (IMREG.GT.0) THEN
  269.          WRITE(IOIMP,*) 'On a créé',
  270.      $        ' JDTJAC(élément , poi.gauss ,',
  271.      $        '1,1,1,1)'
  272.          CALL PRCHVA(JDTJAC,IMPR,IRET)
  273. *         CALL PRCHVA(JDTJAC,6,IRET)
  274.          IF (IRET.NE.0) GOTO 9999
  275.       ENDIF
  276. *      ENDIF
  277. *
  278. * Normal termination
  279. *
  280.       IRET=0
  281.       RETURN
  282. *
  283. * Format handling
  284. *
  285. *
  286. * Error handling
  287. *
  288.  9999 CONTINUE
  289.       IRET=1
  290.       WRITE(IOIMP,*) 'An error was detected in subroutine geolik'
  291.       RETURN
  292. *
  293. * End of subroutine GEOLIK
  294. *
  295.       END
  296.  
  297.  

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