Télécharger geolik.eso

Retour à la liste

Numérotation des lignes :

  1. C GEOLIK SOURCE GOUNAND 14/05/28 21:15:07 8056
  2. SUBROUTINE GEOLIK(DFFPG,JCOOR,NBELEM,
  3. $ JMAJAC,JMIJAC,JDTJAC,LERJ,IMREG,
  4. $ IMPR,IRET)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM : GEOLIK
  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 geolik'
  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 (IMREG.EQ.1) THEN
  152. NDELEM=1
  153. NDPOGO=1
  154. ENDIF
  155. IF (IESREF.LT.IESREL) THEN
  156. LCARRE=.FALSE.
  157. ELSEIF (IESREF.EQ.IESREL) THEN
  158. LCARRE=.TRUE.
  159. ELSE
  160. WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???'
  161. WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL
  162. GOTO 9999
  163. ENDIF
  164. *
  165. * Initialisations...
  166. *
  167. NBLIG=1
  168. NBCOL=1
  169. N2LIG=IESREL
  170. N2COL=IESREF
  171. NBPOI=NDPOGO
  172. NBELM=NDELEM
  173. SEGINI JMAJAC
  174. *
  175. * On effectue le calcul de la matrice jacobienne aux points de Gauss
  176. *
  177. CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM,
  178. $ DFFPG.VELCHE,JCOOR.VELCHE,
  179. $ JMAJAC.VELCHE,
  180. $ IMPR,IRET)
  181. IF (IRET.NE.0) GOTO 9999
  182. IF (IMPR.GT.4) THEN
  183. WRITE(IOIMP,*) 'On a créé',
  184. $ ' JMAJAC(élément , poi.gauss ,',
  185. $ ' dksi(i) , dx(j) ,1,1)'
  186. CALL PRCHVA(JMAJAC,IMPR,IRET)
  187. IF (IRET.NE.0) GOTO 9999
  188. ENDIF
  189. *
  190. * Si la matrice est carrée, on tente de calculer
  191. * son inverse et son déterminant...
  192. *
  193. IF (LCARRE) THEN
  194. NBLIG=1
  195. NBCOL=1
  196. N2LIG=IESREF
  197. N2COL=IESREL
  198. NBPOI=NDPOGO
  199. NBELM=NDELEM
  200. SEGINI JMIJAC
  201. NBLIG=1
  202. NBCOL=1
  203. N2LIG=1
  204. N2COL=1
  205. NBPOI=NDPOGO
  206. NBELM=NDELEM
  207. SEGINI JDTJAC
  208. CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE,
  209. $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ,
  210. $ IMPR,IRET)
  211. IF (IRET.NE.0) THEN
  212. IF (LERJ) GOTO 9666
  213. GOTO 9999
  214. ENDIF
  215. SEGDES JDTJAC
  216. SEGDES JMIJAC
  217. IF (IMPR.GT.4) THEN
  218. WRITE(IOIMP,*) 'On a créé',
  219. $ ' JMIJAC(élément , poi.gauss ,',
  220. $ ' dx(i) , dksi(j) ,1,1)'
  221. CALL PRCHVA(JMIJAC,IMPR,IRET)
  222. IF (IRET.NE.0) GOTO 9999
  223. ENDIF
  224. ELSE
  225. *
  226. * Sinon on calcule : SQRT ( det (transpose(A) * A))
  227. * et la pseudo-inverse (JtJ)^-1 Jt
  228. *
  229. NBLIG=1
  230. NBCOL=1
  231. N2LIG=IESREF
  232. N2COL=IESREL
  233. NBPOI=NDPOGO
  234. NBELM=NDELEM
  235. SEGINI JMIJAC
  236. *
  237. NBLIG=1
  238. NBCOL=1
  239. N2LIG=1
  240. N2COL=1
  241. NBPOI=NDPOGO
  242. NBELM=NDELEM
  243. SEGINI JDTJAC
  244. * Objets temporaire
  245. NBLIG=1
  246. NBCOL=1
  247. N2LIG=IESREF
  248. N2COL=IESREF
  249. NBPOI=NDPOGO
  250. NBELM=NDELEM
  251. SEGINI JJTJ
  252. SEGINI JJTJM1
  253. CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE,
  254. $ JJTJ.VELCHE,JJTJM1.VELCHE,
  255. $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ,
  256. $ IMPR,IRET)
  257. IF (IRET.NE.0) GOTO 9999
  258. SEGSUP JJTJM1
  259. SEGSUP JJTJ
  260. SEGDES JMIJAC
  261. SEGDES JDTJAC
  262. ENDIF
  263. SEGDES JMAJAC
  264. SEGDES JCOOR
  265. SEGDES DFFPG
  266. IF (IMPR.GT.3) THEN
  267. * IF (IMREG.GT.0) THEN
  268. WRITE(IOIMP,*) 'On a créé',
  269. $ ' JDTJAC(élément , poi.gauss ,',
  270. $ '1,1,1,1)'
  271. CALL PRCHVA(JDTJAC,IMPR,IRET)
  272. * CALL PRCHVA(JDTJAC,6,IRET)
  273. IF (IRET.NE.0) GOTO 9999
  274. ENDIF
  275. * ENDIF
  276. *
  277. * Normal termination
  278. *
  279. IRET=0
  280. RETURN
  281. *
  282. * Format handling
  283. *
  284. *
  285. * Error handling
  286. *
  287. 9666 CONTINUE
  288. IRET=666
  289. RETURN
  290. 9999 CONTINUE
  291. IRET=1
  292. WRITE(IOIMP,*) 'An error was detected in subroutine geolik'
  293. RETURN
  294. *
  295. * End of subroutine GEOLIK
  296. *
  297. END
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  

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