Télécharger geolik.eso

Retour à la liste

Numérotation des lignes :

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

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