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

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