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

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