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.  
  102. -INC PPARAM
  103. -INC CCOPTIO
  104. CBEGININCLUDE SFACTIV
  105. SEGMENT FACTIV
  106. POINTEUR IFACTI(NBSOUV).SFACTI
  107. ENDSEGMENT
  108. SEGMENT SFACTI
  109. POINTEUR ISFACT(NBSOFV).SSFACT
  110. ENDSEGMENT
  111. SEGMENT SSFACT
  112. LOGICAL LFACTI(NBELFV,NBELEV)
  113. ENDSEGMENT
  114. CENDINCLUDE SFACTIV
  115. CBEGININCLUDE SMCHAEL
  116. SEGMENT MCHAEL
  117. POINTEUR IMACHE(N1).MELEME
  118. POINTEUR ICHEVA(N1).MCHEVA
  119. ENDSEGMENT
  120. SEGMENT MCHEVA
  121. REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
  122. ENDSEGMENT
  123. SEGMENT LCHEVA
  124. POINTEUR LISCHE(NBCHE).MCHEVA
  125. ENDSEGMENT
  126. CENDINCLUDE SMCHAEL
  127. INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM
  128. POINTEUR JCOOR.MCHEVA
  129. POINTEUR JMAJAC.MCHEVA
  130. POINTEUR JMIJAC.MCHEVA
  131. POINTEUR JDTJAC.MCHEVA
  132. * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss
  133. POINTEUR DFFPG.MCHEVA
  134. *
  135. INTEGER NBELEV,NBELEF,NBELFV
  136. INTEGER IMPR,IRET
  137. *
  138. LOGICAL LCARRE,LERJ
  139. INTEGER IESREF,IESREL,NBPOGO
  140. *
  141. * Executable statements
  142. *
  143. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolif'
  144. SEGACT SSFACT
  145. NBELFV=SSFACT.LFACTI(/1)
  146. NBELEV=SSFACT.LFACTI(/2)
  147. SEGACT JCOOR
  148. NBLIG=JCOOR.VELCHE(/1)
  149. NBCOL=JCOOR.VELCHE(/2)
  150. N2LIG=JCOOR.VELCHE(/3)
  151. N2COL=JCOOR.VELCHE(/4)
  152. NBPOI=JCOOR.VELCHE(/5)
  153. NBELM=JCOOR.VELCHE(/6)
  154. * On pourrait envisager un cas où tous les éléments seraient
  155. * identiques (ex. carré) NBELM=1
  156. IF (NBLIG.NE.1.OR.N2LIG.NE.1
  157. $ .OR.NBPOI.NE.1.OR.NBELM.NE.NBELEV) THEN
  158. WRITE(IOIMP,*) 'Erreur dims JCOOR'
  159. GOTO 9999
  160. ENDIF
  161. NDDL=NBCOL
  162. IESREL=N2COL
  163. SEGACT DFFPG
  164. NBLIG=DFFPG.VELCHE(/1)
  165. NBCOL=DFFPG.VELCHE(/2)
  166. N2LIG=DFFPG.VELCHE(/3)
  167. N2COL=DFFPG.VELCHE(/4)
  168. NBPOI=DFFPG.VELCHE(/5)
  169. NBELM=DFFPG.VELCHE(/6)
  170. IF (NBLIG.NE.1.OR.NBCOL.NE.NDDL.OR.N2LIG.NE.1
  171. $ .OR.(NBELM.NE.NBELFV.AND.NBELM.NE.1)) THEN
  172. WRITE(IOIMP,*) 'Erreur dims DFFPG'
  173. GOTO 9999
  174. ENDIF
  175. IESREF=N2COL
  176. NBPOGO=NBPOI
  177. NLFVDF=NBELM
  178. IF (IESREF.LT.IESREL) THEN
  179. LCARRE=.FALSE.
  180. ELSEIF (IESREF.EQ.IESREL) THEN
  181. LCARRE=.TRUE.
  182. ELSE
  183. WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  184. WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  185. GOTO 9999
  186. ENDIF
  187. *
  188. * Initialisations...
  189. *
  190. NBLIG=1
  191. NBCOL=1
  192. N2LIG=IESREL
  193. N2COL=IESREF
  194. NBPOI=NBPOGO
  195. NBELM=NBELEF
  196. SEGINI JMAJAC
  197. * SEGPRT,DFFPG
  198. *
  199. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  200. *
  201. CALL GEOLF1(IESREL,IESREF,NDDL,NBPOGO,NBELEV,NBELFV,NBELEF,
  202. $ NLFVDF,
  203. $ DFFPG.VELCHE,JCOOR.VELCHE,SSFACT.LFACTI,
  204. $ JMAJAC.VELCHE,
  205. $ IMPR,IRET)
  206. IF (IRET.NE.0) GOTO 9999
  207. * SEGPRT,JMAJAC
  208. *
  209. * Si la matrice est carrée, on tente de calculer
  210. * son inverse et son déterminant...
  211. *
  212. IF (LCARRE) THEN
  213. NBLIG=1
  214. NBCOL=1
  215. N2LIG=IESREF
  216. N2COL=IESREL
  217. NBPOI=NBPOGO
  218. NBELM=NBELEF
  219. SEGINI JMIJAC
  220. NBLIG=1
  221. NBCOL=1
  222. N2LIG=1
  223. N2COL=1
  224. NBPOI=NBPOGO
  225. NBELM=NBELEF
  226. SEGINI JDTJAC
  227. CALL GEOLI2(IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE,
  228. $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ,
  229. $ IMPR,IRET)
  230. * SEGPRT,JDTJAC
  231. IF (IRET.NE.0) THEN
  232. IF (LERJ) GOTO 9666
  233. GOTO 9999
  234. ENDIF
  235. SEGDES JDTJAC
  236. SEGDES JMIJAC
  237. ELSE
  238. *
  239. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  240. *
  241. * Pas de matrice inverse
  242. JMIJAC=0
  243. *
  244. NBLIG=1
  245. NBCOL=1
  246. N2LIG=1
  247. N2COL=1
  248. NBPOI=NBPOGO
  249. NBELM=NBELEF
  250. SEGINI JDTJAC
  251. CALL GEOLI3(IESREL,IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE,
  252. $ JDTJAC.VELCHE,
  253. $ IMPR,IRET)
  254. IF (IRET.NE.0) GOTO 9999
  255. SEGDES JDTJAC
  256. ENDIF
  257. SEGDES JMAJAC
  258. SEGDES SSFACT
  259. SEGDES JCOOR
  260. SEGDES DFFPG
  261. *
  262. * Normal termination
  263. *
  264. IRET=0
  265. RETURN
  266. *
  267. * Format handling
  268. *
  269. *
  270. * Error handling
  271. *
  272. 9666 CONTINUE
  273. IRET=666
  274. RETURN
  275. 9999 CONTINUE
  276. IRET=1
  277. WRITE(IOIMP,*) 'An error was detected in subroutine geolif'
  278. RETURN
  279. *
  280. * End of subroutine GEOLIF
  281. *
  282. END
  283.  
  284.  
  285.  

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