Télécharger geolin.eso

Retour à la liste

Numérotation des lignes :

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

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