Télécharger la plaquette Cast3M


fr en

Télécharger geolin.eso

Retour à la liste

Numérotation des lignes : 

geolin

  1. C GEOLIN    SOURCE    GOUNAND   26/01/09    21:15:28     12441          
  2.       SUBROUTINE GEOLIN(DFFPG,JCOOR,NBELEM,
  3.      $     JMAJAC,JMIJAC,JDTJAC,
  4.      $     IMPR,IRET)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM         : GEOLIN
  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 :
  93. C HISTORIQUE :
  94. C***********************************************************************
  95. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  96. C en cas de modification de ce sous-programme afin de faciliter
  97. C la maintenance !
  98. C***********************************************************************
  99.  
  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 geolin'
  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 (IESREF.LT.IESREL) THEN
  144.          LCARRE=.FALSE.
  145.       ELSEIF (IESREF.EQ.IESREL) THEN
  146.          LCARRE=.TRUE.
  147.       ELSE
  148.          WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  149.          WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  150.          GOTO 9999
  151.       ENDIF
  152. *
  153. * Initialisations...
  154. *
  155.       NBLIG=1
  156.       NBCOL=1
  157.       N2LIG=IESREL
  158.       N2COL=IESREF
  159.       NBPOI=NDPOGO
  160.       NBELM=NDELEM
  161.       SEGINI JMAJAC
  162. *
  163. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  164. *
  165.       CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM,
  166.      $     DFFPG.WELCHE,JCOOR.WELCHE,
  167.      $     JMAJAC.WELCHE,
  168.      $     IMPR,IRET)
  169.       IF (IRET.NE.0) GOTO 9999
  170.       IF (IMPR.GT.4) THEN
  171.          WRITE(IOIMP,*) 'On a créé',
  172.      $   ' JMAJAC(élément , poi.gauss ,',
  173.      $        ' dksi(i) , dx(j) ,1,1)'
  174.          CALL PRCHVA(JMAJAC,IMPR,IRET)
  175.          IF (IRET.NE.0) GOTO 9999
  176.       ENDIF
  177. *
  178. * Si la matrice est carrée, on tente de calculer
  179. * son inverse et son déterminant...
  180. *
  181.       IF (LCARRE) THEN
  182.          NBLIG=1
  183.          NBCOL=1
  184.          N2LIG=IESREF
  185.          N2COL=IESREL
  186.          NBPOI=NDPOGO
  187.          NBELM=NDELEM
  188.          SEGINI JMIJAC
  189.          NBLIG=1
  190.          NBCOL=1
  191.          N2LIG=1
  192.          N2COL=1
  193.          NBPOI=NDPOGO
  194.          NBELM=NDELEM
  195.          SEGINI JDTJAC
  196.          CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  197.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,
  198.      $        IMPR,IRET)
  199.          IF (IRET.LT.0) THEN
  200.             SEGSUP JMIJAC
  201.             JMIJAC=0
  202.          ELSEIF (IRET.GT.0) THEN
  203.             GOTO 9999
  204.          ENDIF
  205. *         SEGDES JDTJAC
  206. *         SEGDES JMIJAC
  207.          IF (IMPR.GT.4) THEN
  208.             IF (JMIJAC.GT.0) THEN
  209.                WRITE(IOIMP,*) 'On a créé',
  210.      $              ' JMIJAC(élément , poi.gauss ,',
  211.      $              ' dx(i) , dksi(j) ,1,1)'
  212.                CALL PRCHVA(JMIJAC,IMPR,IRET)
  213.                IF (IRET.NE.0) GOTO 9999
  214.             ENDIF
  215.          ENDIF
  216.       ELSE
  217. *
  218. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  219. * et la pseudo-inverse (JtJ)^-1 Jt
  220. *
  221.          NBLIG=1
  222.          NBCOL=1
  223.          N2LIG=IESREF
  224.          N2COL=IESREL
  225.          NBPOI=NDPOGO
  226.          NBELM=NDELEM
  227.          SEGINI JMIJAC
  228. *
  229.          NBLIG=1
  230.          NBCOL=1
  231.          N2LIG=1
  232.          N2COL=1
  233.          NBPOI=NDPOGO
  234.          NBELM=NDELEM
  235.          SEGINI JDTJAC
  236. * Objets temporaire
  237.          NBLIG=1
  238.          NBCOL=1
  239.          N2LIG=IESREF
  240.          N2COL=IESREF
  241.          NBPOI=NDPOGO
  242.          NBELM=NDELEM
  243.          SEGINI JJTJ
  244.          SEGINI JJTJM1
  245.          CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  246.      $        JJTJ.WELCHE,JJTJM1.WELCHE,
  247.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,
  248.      $        IMPR,IRET)
  249.          IF (IRET.NE.0) GOTO 9999
  250.          SEGSUP JJTJM1
  251.          SEGSUP JJTJ
  252. *         SEGDES JMIJAC
  253. *         SEGDES JDTJAC
  254.       ENDIF
  255. *      SEGDES JMAJAC
  256. *      SEGDES JCOOR
  257. *      SEGDES DFFPG
  258.       IF (IMPR.GT.3) THEN
  259.          WRITE(IOIMP,*) 'On a créé',
  260.      $        ' JDTJAC(élément , poi.gauss ,',
  261.      $        '1,1,1,1)'
  262.          CALL PRCHVA(JDTJAC,IMPR,IRET)
  263.          IF (IRET.NE.0) GOTO 9999
  264.       ENDIF
  265. *
  266. * Normal termination
  267. *
  268.       IRET=0
  269.       RETURN
  270. *
  271. * Format handling
  272. *
  273. *
  274. * Error handling
  275. *
  276.  9666 CONTINUE
  277.       IRET=666
  278.       RETURN
  279.  9999 CONTINUE
  280.       IRET=1
  281.       WRITE(IOIMP,*) 'An error was detected in subroutine geolin'
  282.       RETURN
  283. *
  284. * End of subroutine GEOLIN
  285. *
  286.       END
  287.  
  288.  

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