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   21/06/02    21:16:10     11022          
  2.       SUBROUTINE GEOLIK(DFFPG,JCOOR,NBELEM,
  3.      $     JMAJAC,JMIJAC,JDTJAC,LERJ,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 :
  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. -INC PPARAM
  100. -INC CCOPTIO
  101. -INC TNLIN      
  102. *-INC SMCHAEL
  103.       INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
  104.       POINTEUR JCOOR.MCHEVA
  105.       POINTEUR JMAJAC.MCHEVA
  106.       POINTEUR JMIJAC.MCHEVA
  107.       POINTEUR JDTJAC.MCHEVA
  108.       POINTEUR JJTJ.MCHEVA
  109.       POINTEUR JJTJM1.MCHEVA
  110. * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
  111.       POINTEUR DFFPG.MCHEVA
  112. *
  113.       INTEGER NBELEM
  114.       INTEGER IMPR,IRET
  115. *
  116.       LOGICAL LCARRE,LERJ
  117.       INTEGER IESREF,IESREL,NDNOEU,NDCOL,NDPOGO,NDELEM
  118. *
  119. * Executable statements
  120. *
  121.       IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolik'
  122.       SEGACT DFFPG
  123.       NDNOEU=DFFPG.WELCHE(/2)
  124.       IESREF=DFFPG.WELCHE(/4)
  125.       NDPOGO=DFFPG.WELCHE(/5)
  126.       SEGACT JCOOR
  127.       NDCOL=JCOOR.WELCHE(/2)
  128.       IESREL=JCOOR.WELCHE(/4)
  129.       NDELEM=JCOOR.WELCHE(/6)
  130.       IF (NDCOL.NE.NDNOEU) THEN
  131.          WRITE(IOIMP,*) 'Incompatibilité fns d''interpolation-géométrie'
  132.          WRITE(IOIMP,*) 'NDNOEU=',NDNOEU,' NDCOL=',NDCOL
  133.          GOTO 9999
  134.       ENDIF
  135. * On pourrait envisager un cas où tous les éléments seraient
  136. * identiques (ex. carré) NDELEM=1
  137.       IF (NDELEM.NE.NBELEM) THEN
  138.          WRITE(IOIMP,*) 'Incompatibilité coordonnées-géométrie'
  139.          WRITE(IOIMP,*) 'NDELEM=',NDELEM,' NBELEM=',NBELEM
  140.          GOTO 9999
  141.       ENDIF
  142.       IF (IMREG.EQ.1) THEN
  143.          NDELEM=1
  144.          NDPOGO=1
  145.       ENDIF
  146.       IF (IESREF.LT.IESREL) THEN
  147.          LCARRE=.FALSE.
  148.       ELSEIF (IESREF.EQ.IESREL) THEN
  149.          LCARRE=.TRUE.
  150.       ELSE
  151.          WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  152.          WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  153.          GOTO 9999
  154.       ENDIF
  155. *
  156. * Initialisations...
  157. *
  158.       NBLIG=1
  159.       NBCOL=1
  160.       N2LIG=IESREL
  161.       N2COL=IESREF
  162.       NBPOI=NDPOGO
  163.       NBELM=NDELEM
  164.       SEGINI JMAJAC
  165. *
  166. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  167. *
  168.       CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM,
  169.      $     DFFPG.WELCHE,JCOOR.WELCHE,
  170.      $     JMAJAC.WELCHE,
  171.      $     IMPR,IRET)
  172.       IF (IRET.NE.0) GOTO 9999
  173.       IF (IMPR.GT.4) THEN
  174.          WRITE(IOIMP,*) 'On a créé',
  175.      $   ' JMAJAC(élément , poi.gauss ,',
  176.      $        ' dksi(i) , dx(j) ,1,1)'
  177.          CALL PRCHVA(JMAJAC,IMPR,IRET)
  178.          IF (IRET.NE.0) GOTO 9999
  179.       ENDIF
  180. *
  181. * Si la matrice est carrée, on tente de calculer
  182. * son inverse et son déterminant...
  183. *
  184.       IF (LCARRE) THEN
  185.          NBLIG=1
  186.          NBCOL=1
  187.          N2LIG=IESREF
  188.          N2COL=IESREL
  189.          NBPOI=NDPOGO
  190.          NBELM=NDELEM
  191.          SEGINI JMIJAC
  192.          NBLIG=1
  193.          NBCOL=1
  194.          N2LIG=1
  195.          N2COL=1
  196.          NBPOI=NDPOGO
  197.          NBELM=NDELEM
  198.          SEGINI JDTJAC
  199.          CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  200.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,LERJ,
  201.      $        IMPR,IRET)
  202.          IF (IRET.NE.0) THEN
  203.             IF (LERJ) GOTO 9666
  204.             GOTO 9999
  205.          ENDIF
  206.          SEGDES JDTJAC
  207.          SEGDES JMIJAC
  208.          IF (IMPR.GT.4) 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.       ELSE
  216. *
  217. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  218. * et la pseudo-inverse (JtJ)^-1 Jt
  219. *
  220.          NBLIG=1
  221.          NBCOL=1
  222.          N2LIG=IESREF
  223.          N2COL=IESREL
  224.          NBPOI=NDPOGO
  225.          NBELM=NDELEM
  226.          SEGINI JMIJAC
  227. *
  228.          NBLIG=1
  229.          NBCOL=1
  230.          N2LIG=1
  231.          N2COL=1
  232.          NBPOI=NDPOGO
  233.          NBELM=NDELEM
  234.          SEGINI JDTJAC
  235. * Objets temporaire
  236.          NBLIG=1
  237.          NBCOL=1
  238.          N2LIG=IESREF
  239.          N2COL=IESREF
  240.          NBPOI=NDPOGO
  241.          NBELM=NDELEM
  242.          SEGINI JJTJ
  243.          SEGINI JJTJM1
  244.          CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.WELCHE,
  245.      $        JJTJ.WELCHE,JJTJM1.WELCHE,
  246.      $        JMIJAC.WELCHE,JDTJAC.WELCHE,LERJ,
  247.      $        IMPR,IRET)
  248.          IF (IRET.NE.0) GOTO 9999
  249.          SEGSUP JJTJM1
  250.          SEGSUP JJTJ
  251.          SEGDES JMIJAC
  252.          SEGDES JDTJAC
  253.       ENDIF
  254.       SEGDES JMAJAC
  255.       SEGDES JCOOR
  256.       SEGDES DFFPG
  257.       IF (IMPR.GT.3) THEN
  258. *      IF (IMREG.GT.0) 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. *         CALL PRCHVA(JDTJAC,6,IRET)
  264.          IF (IRET.NE.0) GOTO 9999
  265.       ENDIF
  266. *      ENDIF
  267. *
  268. * Normal termination
  269. *
  270.       IRET=0
  271.       RETURN
  272. *
  273. * Format handling
  274. *
  275. *
  276. * Error handling
  277. *
  278.  9666 CONTINUE
  279.       IRET=666
  280.       RETURN
  281.  9999 CONTINUE
  282.       IRET=1
  283.       WRITE(IOIMP,*) 'An error was detected in subroutine geolik'
  284.       RETURN
  285. *
  286. * End of subroutine GEOLIK
  287. *
  288.       END
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  

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