Télécharger geoli2.eso

Retour à la liste

Numérotation des lignes :

geoli2
  1. C GEOLI2 SOURCE GOUNAND 21/06/02 21:16:08 11022
  2. SUBROUTINE GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC,
  3. $ JMIJAC,JDTJAC,LERJ,
  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,LERJ
  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. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli2'
  79. IF (IESREF.EQ.1) THEN
  80. DO 1 IELEM=1,NDELEM
  81. DO 12 IPOGO=1,NDPOGO
  82. DET=JMAJAC(1,1,IPOGO,IELEM)
  83. IF (DET.NE.ZERO) THEN
  84. INVDET=UN/DET
  85. JDTJAC(IPOGO,IELEM)=DET
  86. JMIJAC(1,1,IPOGO,IELEM)=INVDET
  87. ELSE
  88. GOTO 9998
  89. ENDIF
  90. 12 CONTINUE
  91. 1 CONTINUE
  92. ELSEIF (IESREF.EQ.2) THEN
  93. DO 3 IELEM=1,NDELEM
  94. DO 32 IPOGO=1,NDPOGO
  95. DET=(JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  96. $ - (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,2,IPOGO,IELEM))
  97. IF (DET.NE.ZERO) THEN
  98. IF (IPOGO.EQ.1) THEN
  99. LSIGN=(DET.GT.ZERO)
  100. ELSE
  101. IF (LSIGN.NEQV.(DET.GT.ZERO)) GOTO 9997
  102. ENDIF
  103. INVDET=UN/DET
  104. JDTJAC(IPOGO,IELEM)=DET
  105. JMIJAC(1,1,IPOGO,IELEM)=
  106. $ (+JMAJAC(2,2,IPOGO,IELEM))*INVDET
  107. JMIJAC(2,1,IPOGO,IELEM)=
  108. $ (-JMAJAC(2,1,IPOGO,IELEM))*INVDET
  109. JMIJAC(1,2,IPOGO,IELEM)=
  110. $ (-JMAJAC(1,2,IPOGO,IELEM))*INVDET
  111. JMIJAC(2,2,IPOGO,IELEM)=
  112. $ (+JMAJAC(1,1,IPOGO,IELEM))*INVDET
  113. ELSE
  114. GOTO 9998
  115. ENDIF
  116. 32 CONTINUE
  117. 3 CONTINUE
  118. ELSEIF (IESREF.EQ.3) THEN
  119. DO 5 IELEM=1,NDELEM
  120. DO 52 IPOGO=1,NDPOGO
  121. JMIJAC(1,1,IPOGO,IELEM)=
  122. $ (JMAJAC(2,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  123. $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  124. JMIJAC(2,1,IPOGO,IELEM)=
  125. $ (JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  126. $ -(JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  127. JMIJAC(3,1,IPOGO,IELEM)=
  128. $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM))
  129. $ -(JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  130. JMIJAC(1,2,IPOGO,IELEM)=
  131. $ (JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM))
  132. $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  133. JMIJAC(2,2,IPOGO,IELEM)=
  134. $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM))
  135. $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM))
  136. JMIJAC(3,2,IPOGO,IELEM)=
  137. $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM))
  138. $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM))
  139. JMIJAC(1,3,IPOGO,IELEM)=
  140. $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM))
  141. $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  142. JMIJAC(2,3,IPOGO,IELEM)=
  143. $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,3,IPOGO,IELEM))
  144. $ -(JMAJAC(2,3,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM))
  145. JMIJAC(3,3,IPOGO,IELEM)=
  146. $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM))
  147. $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,1,IPOGO,IELEM))
  148. DET=(JMAJAC(1,1,IPOGO,IELEM)*JMIJAC(1,1,IPOGO,IELEM))
  149. $ +(JMAJAC(1,2,IPOGO,IELEM)*JMIJAC(2,1,IPOGO,IELEM))
  150. $ +(JMAJAC(1,3,IPOGO,IELEM)*JMIJAC(3,1,IPOGO,IELEM))
  151. * IF (DET.NE.ZERO) THEN
  152. IF (ABS(DET).GT.(SQRT(XPETIT))) THEN
  153. IF (IPOGO.EQ.1) THEN
  154. LSIGN=(DET.GT.ZERO)
  155. ELSE
  156. IF (LSIGN.NEQV.(DET.GT.ZERO)) GOTO 9997
  157. ENDIF
  158. INVDET=UN/DET
  159. JDTJAC(IPOGO,IELEM)=DET
  160. DO 522 JREFER=1,IESREF
  161. DO 5222 IREFER=1,IESREF
  162. JMIJAC(IREFER,JREFER,IPOGO,IELEM)=
  163. $ JMIJAC(IREFER,JREFER,IPOGO,IELEM)
  164. $ * INVDET
  165. 5222 CONTINUE
  166. 522 CONTINUE
  167. ELSE
  168. GOTO 9998
  169. ENDIF
  170. 52 CONTINUE
  171. 5 CONTINUE
  172. ELSE
  173. WRITE(IOIMP,*) 'Je ne sais pas calculer l''inverse d''une'
  174. WRITE(IOIMP,*) 'matrice jacobienne de dimension '
  175. WRITE(IOIMP,*) 'IESREF=',IESREF
  176. ENDIF
  177. *
  178. * Normal termination
  179. *
  180. IRET=0
  181. RETURN
  182. *
  183. * Format handling
  184. *
  185. *
  186. * Error handling
  187. *
  188. 9997 CONTINUE
  189. IF (LERJ) THEN
  190. IRET=666
  191. CALL ECRENT(IRET)
  192. RETURN
  193. ENDIF
  194. WRITE(IOIMP,*) 'Le déterminant de la matrice jacobienne'
  195. WRITE(IOIMP,*) 'change de signe sur l''élément.'
  196. WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
  197. GOTO 9999
  198. 9998 CONTINUE
  199. IF (LERJ) THEN
  200. IRET=667
  201. CALL ECRENT(IRET)
  202. RETURN
  203. ENDIF
  204. WRITE(IOIMP,*) 'Déterminant de la matrice jacobienne nul'
  205. WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO
  206. GOTO 9999
  207. 9999 CONTINUE
  208. IRET=1
  209. WRITE(IOIMP,*) 'An error was detected in subroutine geoli2'
  210. RETURN
  211. *
  212. * End of subroutine GEOLI2
  213. *
  214. END
  215.  
  216.  
  217.  
  218.  

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