Télécharger geolif.eso

Retour à la liste

Numérotation des lignes :

  1. C GEOLIF SOURCE GOUNAND 06/08/04 21:15:57 5520
  2. SUBROUTINE GEOLIF(DFFPG,JCOOR,SSFACT,NBELEF,
  3. $ JMAJAC,JMIJAC,JDTJAC,LERJ,
  4. $ IMPR,IRET)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM : GEOLIF
  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 : RSET8 (copie de tableaux de CHARACTER*8)
  45. C PRCHVA (impression d'un segment de type MCHEVA)
  46. C GEOLI1 (calcul de la matrice jacobienne (fortran 77
  47. C ))
  48. C GEOLI2 (calcul de l'inverse et du déterminant de la
  49. C matrice jacobienne (fortran 77))
  50. C GEOLI3 (calcul de sqrt(det(transpose(J)*J)), J
  51. C étant la matrice jacobienne (fortran 77))
  52. C APPELE PAR : NLIN
  53. C***********************************************************************
  54. C ENTREES : * DFFPG (type MCHEVA) : valeurs des dérivées
  55. C premières des fonctions d'interpolation de la
  56. C transformation géométrique aux points de gauss
  57. C sur l'élément de référence.
  58. C Structure (cf.include SMCHAEL) :
  59. C (1, nb. ddl, 1, dim.esp.réf, nb. poi. gauss, 1)
  60. C * JCOOR (type MCHEVA) : valeurs des ddl de la
  61. C transformation géométrique sur le maillage
  62. C élémentaire courant.
  63. C Structure (cf.include SMCHAEL) :
  64. C (1, nb. ddl, 1, dim. esp. réel,
  65. C 1, nb. éléments)
  66. C * NBELEV (type entier) : nombre d'éléments du
  67. C maillage élémentaire courant.
  68. C ENTREES/SORTIES : -
  69. C SORTIES : * JMAJAC (type MCHEVA) : valeurs de la matrice
  70. C jacobienne aux points de Gauss sur le maillage
  71. C élémentaire courant.
  72. C Structure (cf.include SMCHAEL) :
  73. C (1, 1, dim. esp. réel, dim. esp. référence,
  74. C nb. poi. gauss, nb. éléments)
  75. C * JMIJAC (type MCHEVA) : valeurs de l'inverse de
  76. C la matrice jacobienne aux points de Gauss sur
  77. C le maillage élémentaire courant.
  78. C Structure (cf.include SMCHAEL) :
  79. C (1, 1, dim. esp. référence, dim. esp. réel,
  80. C nb. poi. gauss, nb. éléments)
  81. C JMIJAC n'existe que si dim.esp.réf=dim.esp.réel
  82. C * JDTJAC (type MCHEVA) : valeurs du déterminant
  83. C de la matrice jacobienne aux points de Gauss
  84. C sur le maillage élémentaire courant.
  85. C Structure (cf.include SMCHAEL) :
  86. C (1, 1, 1, 1, nb. poi. gauss, nb. éléments)
  87. C Si la matrice jacobienne n'est pas carrée, on
  88. C a calculé sqrt(det(trans(J).J))
  89. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  90. C***********************************************************************
  91. C VERSION : v2, 03/10/03, refonte complète (modif SMPOUET)
  92. C VERSION : v1, 09/08/99, version initiale
  93. C HISTORIQUE : v1, 09/08/99, création
  94. C HISTORIQUE :
  95. C HISTORIQUE :
  96. C***********************************************************************
  97. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  98. C en cas de modification de ce sous-programme afin de faciliter
  99. C la maintenance !
  100. C***********************************************************************
  101. -INC CCOPTIO
  102. CBEGININCLUDE SFACTIV
  103. SEGMENT FACTIV
  104. POINTEUR IFACTI(NBSOUV).SFACTI
  105. ENDSEGMENT
  106. SEGMENT SFACTI
  107. POINTEUR ISFACT(NBSOFV).SSFACT
  108. ENDSEGMENT
  109. SEGMENT SSFACT
  110. LOGICAL LFACTI(NBELFV,NBELEV)
  111. ENDSEGMENT
  112. CENDINCLUDE SFACTIV
  113. CBEGININCLUDE SMCHAEL
  114. SEGMENT MCHAEL
  115. POINTEUR IMACHE(N1).MELEME
  116. POINTEUR ICHEVA(N1).MCHEVA
  117. ENDSEGMENT
  118. SEGMENT MCHEVA
  119. REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
  120. ENDSEGMENT
  121. SEGMENT LCHEVA
  122. POINTEUR LISCHE(NBCHE).MCHEVA
  123. ENDSEGMENT
  124. CENDINCLUDE SMCHAEL
  125. INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
  126. POINTEUR JCOOR.MCHEVA
  127. POINTEUR JMAJAC.MCHEVA
  128. POINTEUR JMIJAC.MCHEVA
  129. POINTEUR JDTJAC.MCHEVA
  130. * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
  131. POINTEUR DFFPG.MCHEVA
  132. *
  133. INTEGER NBELEV,NBELEF,NBELFV
  134. INTEGER IMPR,IRET
  135. *
  136. LOGICAL LCARRE,LERJ
  137. INTEGER IESREF,IESREL,NBPOGO
  138. *
  139. * Executable statements
  140. *
  141. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolif'
  142. SEGACT SSFACT
  143. NBELFV=SSFACT.LFACTI(/1)
  144. NBELEV=SSFACT.LFACTI(/2)
  145. SEGACT JCOOR
  146. NBLIG=JCOOR.VELCHE(/1)
  147. NBCOL=JCOOR.VELCHE(/2)
  148. N2LIG=JCOOR.VELCHE(/3)
  149. N2COL=JCOOR.VELCHE(/4)
  150. NBPOI=JCOOR.VELCHE(/5)
  151. NBELM=JCOOR.VELCHE(/6)
  152. * On pourrait envisager un cas où tous les éléments seraient
  153. * identiques (ex. carré) NBELM=1
  154. IF (NBLIG.NE.1.OR.N2LIG.NE.1
  155. $ .OR.NBPOI.NE.1.OR.NBELM.NE.NBELEV) THEN
  156. WRITE(IOIMP,*) 'Erreur dims JCOOR'
  157. GOTO 9999
  158. ENDIF
  159. NDDL=NBCOL
  160. IESREL=N2COL
  161. SEGACT DFFPG
  162. NBLIG=DFFPG.VELCHE(/1)
  163. NBCOL=DFFPG.VELCHE(/2)
  164. N2LIG=DFFPG.VELCHE(/3)
  165. N2COL=DFFPG.VELCHE(/4)
  166. NBPOI=DFFPG.VELCHE(/5)
  167. NBELM=DFFPG.VELCHE(/6)
  168. IF (NBLIG.NE.1.OR.NBCOL.NE.NDDL.OR.N2LIG.NE.1
  169. $ .OR.(NBELM.NE.NBELFV.AND.NBELM.NE.1)) THEN
  170. WRITE(IOIMP,*) 'Erreur dims DFFPG'
  171. GOTO 9999
  172. ENDIF
  173. IESREF=N2COL
  174. NBPOGO=NBPOI
  175. NLFVDF=NBELM
  176. IF (IESREF.LT.IESREL) THEN
  177. LCARRE=.FALSE.
  178. ELSEIF (IESREF.EQ.IESREL) THEN
  179. LCARRE=.TRUE.
  180. ELSE
  181. WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  182. WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  183. GOTO 9999
  184. ENDIF
  185. *
  186. * Initialisations...
  187. *
  188. NBLIG=1
  189. NBCOL=1
  190. N2LIG=IESREL
  191. N2COL=IESREF
  192. NBPOI=NBPOGO
  193. NBELM=NBELEF
  194. SEGINI JMAJAC
  195. * SEGPRT,DFFPG
  196. *
  197. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  198. *
  199. CALL GEOLF1(IESREL,IESREF,NDDL,NBPOGO,NBELEV,NBELFV,NBELEF,
  200. $ NLFVDF,
  201. $ DFFPG.VELCHE,JCOOR.VELCHE,SSFACT.LFACTI,
  202. $ JMAJAC.VELCHE,
  203. $ IMPR,IRET)
  204. IF (IRET.NE.0) GOTO 9999
  205. * SEGPRT,JMAJAC
  206. *
  207. * Si la matrice est carrée, on tente de calculer
  208. * son inverse et son déterminant...
  209. *
  210. IF (LCARRE) THEN
  211. NBLIG=1
  212. NBCOL=1
  213. N2LIG=IESREF
  214. N2COL=IESREL
  215. NBPOI=NBPOGO
  216. NBELM=NBELEF
  217. SEGINI JMIJAC
  218. NBLIG=1
  219. NBCOL=1
  220. N2LIG=1
  221. N2COL=1
  222. NBPOI=NBPOGO
  223. NBELM=NBELEF
  224. SEGINI JDTJAC
  225. CALL GEOLI2(IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE,
  226. $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ,
  227. $ IMPR,IRET)
  228. * SEGPRT,JDTJAC
  229. IF (IRET.NE.0) THEN
  230. IF (LERJ) GOTO 9666
  231. GOTO 9999
  232. ENDIF
  233. SEGDES JDTJAC
  234. SEGDES JMIJAC
  235. ELSE
  236. *
  237. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  238. *
  239. * Pas de matrice inverse
  240. JMIJAC=0
  241. *
  242. NBLIG=1
  243. NBCOL=1
  244. N2LIG=1
  245. N2COL=1
  246. NBPOI=NBPOGO
  247. NBELM=NBELEF
  248. SEGINI JDTJAC
  249. CALL GEOLI3(IESREL,IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE,
  250. $ JDTJAC.VELCHE,
  251. $ IMPR,IRET)
  252. IF (IRET.NE.0) GOTO 9999
  253. SEGDES JDTJAC
  254. ENDIF
  255. SEGDES JMAJAC
  256. SEGDES SSFACT
  257. SEGDES JCOOR
  258. SEGDES DFFPG
  259. *
  260. * Normal termination
  261. *
  262. IRET=0
  263. RETURN
  264. *
  265. * Format handling
  266. *
  267. *
  268. * Error handling
  269. *
  270. 9666 CONTINUE
  271. IRET=666
  272. RETURN
  273. 9999 CONTINUE
  274. IRET=1
  275. WRITE(IOIMP,*) 'An error was detected in subroutine geolif'
  276. RETURN
  277. *
  278. * End of subroutine GEOLIF
  279. *
  280. END
  281.  
  282.  
  283.  

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