Télécharger ccgnor.eso

Retour à la liste

Numérotation des lignes :

  1. C CCGNOR SOURCE GOUNAND 11/04/29 21:15:18 6947
  2. SUBROUTINE CCGNOR(LCOF,NOMLOI,
  3. $ FC,
  4. $ IMPR,IRET)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. IMPLICIT INTEGER (I-N)
  7. C***********************************************************************
  8. C NOM : CCGNOR
  9. C DESCRIPTION : Calcul des composantes d'un vecteur normal
  10. C
  11. C
  12. C LANGAGE : ESOPE
  13. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  14. C mél : gounand@semt2.smts.cea.fr
  15. C***********************************************************************
  16. C APPELES :
  17. C APPELE PAR :
  18. C***********************************************************************
  19. C ENTREES :
  20. C ENTREES/SORTIES :
  21. C SORTIES : -
  22. C TRAVAIL :
  23. C***********************************************************************
  24. C VERSION : v1, 09/03/07, version initiale
  25. C HISTORIQUE : v1, 09/03/07, création
  26. C HISTORIQUE :
  27. C HISTORIQUE :
  28. C***********************************************************************
  29. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  30. C en cas de modification de ce sous-programme afin de faciliter
  31. C la maintenance !
  32. C***********************************************************************
  33. -INC CCOPTIO
  34. -INC CCREEL
  35. CBEGININCLUDE SMCHAEL
  36. SEGMENT MCHAEL
  37. POINTEUR IMACHE(N1).MELEME
  38. POINTEUR ICHEVA(N1).MCHEVA
  39. ENDSEGMENT
  40. SEGMENT MCHEVA
  41. REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
  42. ENDSEGMENT
  43. SEGMENT LCHEVA
  44. POINTEUR LISCHE(NBCHE).MCHEVA
  45. ENDSEGMENT
  46. CENDINCLUDE SMCHAEL
  47. INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  48. POINTEUR FC.MCHEVA
  49. POINTEUR LCOF.LCHEVA
  50. POINTEUR JMAJAC.MCHEVA
  51. POINTEUR JMIJAC.MCHEVA
  52. POINTEUR JDTJAC.MCHEVA
  53. CHARACTER*8 NOMLOI
  54. INTEGER ICOF
  55. *
  56. -INC TMXMAT
  57. * Objets temporaires
  58. POINTEUR JAC.MXMAT,JT.MXMAT,JP.MXMAT
  59. POINTEUR G.MXMAT,IG.MXMAT
  60. *
  61. SEGMENT MIMAT2
  62. INTEGER IMAT2(2,2)
  63. ENDSEGMENT
  64. POINTEUR EIJ.MIMAT2
  65. SEGMENT MIMAT3
  66. INTEGER IMAT3(3,3,3)
  67. ENDSEGMENT
  68. POINTEUR EIJK.MIMAT3
  69. *
  70. LOGICAL LBID
  71. INTEGER LAXSP
  72. *
  73. INTEGER IMPR,IRET
  74. *
  75. * Executable statements
  76. *
  77. * IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgnor'
  78. * WRITE(IOIMP,*) 'Entrée dans ccgnor'
  79. * WRITE(IOIMP,*) 'NOMLOI=',NOMLOI
  80. C IF (.NOT.(IDIM.EQ.1)) THEN
  81. C WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
  82. C GOTO 9999
  83. C ENDIF
  84. NLFC=FC.VELCHE(/6)
  85. NPFC=FC.VELCHE(/5)
  86. ICOF=0
  87. *
  88. ICOF=ICOF+1
  89. JMAJAC=LCOF.LISCHE(ICOF)
  90. NLJA=JMAJAC.VELCHE(/6)
  91. NPJA=JMAJAC.VELCHE(/5)
  92. IREF=JMAJAC.VELCHE(/4)
  93. IREL=JMAJAC.VELCHE(/3)
  94. * SEGPRT,JMAJAC
  95. *
  96. * WRITE(IOIMP,*) 'IREL=',IREL
  97. * WRITE(IOIMP,*) 'IREF=',IREF
  98. IF (IREL.NE.IDIM) THEN
  99. WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  100. GOTO 9999
  101. ENDIF
  102. IF (IREL.NE.IREF+1) THEN
  103. WRITE(IOIMP,*) 'Le maillage donne nest pas une surface'
  104. GOTO 9999
  105. ENDIF
  106. IF ((IREL.NE.2).AND.(IREL.NE.3)) THEN
  107. WRITE(IOIMP,*) 'Ne marche quen dimension despace 2 ou 3'
  108. GOTO 9999
  109. ENDIF
  110. *
  111. * Objets temporaires et à préconditionner
  112. *
  113. LDIM1=IREL
  114. LDIM2=IREF
  115. SEGINI,JAC
  116. LDIM1=IREF
  117. LDIM2=IREL
  118. SEGINI,JT
  119. * SEGINI,JP
  120. LDIM1=IREF
  121. LDIM2=IREF
  122. SEGINI,G
  123. SEGINI,IG
  124. *
  125. * Initialisation des tenseurs de permutation
  126. *
  127. SEGINI,EIJ
  128. IMULT=1
  129. ICPT=0
  130. DO I=1,2
  131. DO J=1,2
  132. IF (I.NE.J) THEN
  133. ICPT=ICPT+1
  134. IF (ICPT.EQ.2) THEN
  135. ICPT=0
  136. IMULT=IMULT*(-1)
  137. ENDIF
  138. EIJ.IMAT2(I,J)=IMULT
  139. ENDIF
  140. ENDDO
  141. ENDDO
  142. SEGINI,EIJK
  143. IMULT=1
  144. ICPT=0
  145. DO I=1,3
  146. DO J=1,3
  147. IF (I.NE.J) THEN
  148. DO K=1,3
  149. IF ((K.NE.I).AND.(K.NE.J)) THEN
  150. ICPT=ICPT+1
  151. IF (ICPT.EQ.2) THEN
  152. ICPT=0
  153. IMULT=IMULT*(-1)
  154. ENDIF
  155. EIJK.IMAT3(I,J,K)=IMULT
  156. ENDIF
  157. ENDDO
  158. ENDIF
  159. ENDDO
  160. ENDDO
  161. * SEGPRT,EIJ
  162. * SEGPRT,EIJK
  163. *
  164. DO ILFC=1,NLFC
  165. IF (NLJA.EQ.1) THEN
  166. ILJA=1
  167. ELSE
  168. ILJA=ILFC
  169. ENDIF
  170. DO IPFC=1,NPFC
  171. IF (NPJA.EQ.1) THEN
  172. IPJA=1
  173. ELSE
  174. IPJA=IPFC
  175. ENDIF
  176. *
  177. * Copie du jacobien
  178. *
  179. CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  180. $ 'COPIE ',
  181. $ JAC.XMAT,IREL,IREF,
  182. $ IMPR,IRET)
  183. IF (IRET.NE.0) GOTO 9999
  184. * Calcul de Jt
  185. CALL MAMA(JAC.XMAT,IREL,IREF,
  186. $ 'TRANSPOS',JT.XMAT,IREF,IREL,
  187. $ IMPR,IRET)
  188. IF (IRET.NE.0) GOTO 9999
  189. * Calcul de G=JtJ
  190. CALL MAMAMA(JT.XMAT,IREF,IREL,JAC.XMAT,IREL,IREF,
  191. $ 'FOIS ',G.XMAT,IREF,IREF,
  192. $ IMPR,IRET)
  193. IF (IRET.NE.0) GOTO 9999
  194. * Calcul de l'inverse, du déterminant et trace de l'inverse de g
  195. LBID=.FALSE.
  196. CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET)
  197. IF (IRET.NE.0) GOTO 9999
  198. XM=SQRT(DETG)
  199. **
  200. ** Calcul de la pseudo-inverse J+ = g-1 Jt
  201. **
  202. * CALL MAMAMA(IG.XMAT,IREF,IREF,JT.XMAT,IREF,IREL,
  203. * $ 'FOIS ',JP.XMAT,IREF,IREL,
  204. * $ IMPR,IRET)
  205. * IF (IRET.NE.0) GOTO 9999
  206.  
  207. IF (NOMLOI(1:4).EQ.'VNOR') THEN
  208. CALL CH2INT(NOMLOI(5:5),I,IMPR,IRET)
  209. IF (IRET.NE.0) GOTO 9999
  210. **
  211. ** Calcul de la pseudo-inverse J+ = g-1 Jt
  212. **
  213. * CALL MAMAMA(IG.XMAT,IREF,IREF,JT.XMAT,IREF,IREL,
  214. * $ 'FOIS ',JP.XMAT,IREF,IREL,
  215. * $ IMPR,IRET)
  216. * IF (IRET.NE.0) GOTO 9999
  217.  
  218. *
  219. IF (IREL.EQ.2) THEN
  220. CONTRI=0.D0
  221. DO J=1,2
  222. CONTRI=CONTRI-(EIJ.IMAT2(I,J)*
  223. $ JMAJAC.VELCHE(1,1,J,1,IPJA,ILJA))
  224. ENDDO
  225. CONTRI=CONTRI/XM
  226. ELSEIF (IREL.EQ.3) THEN
  227. C XNUM=XM
  228. C XDENO=0.D0
  229. C* C'est louche parce que II ne varie pas de 1 à IREL
  230. C DO II=1,IREF
  231. C DO IJ=1,IREF
  232. C* DO IN=1,IREF
  233. C* DO IO=1,IREF
  234. C XDENO=XDENO+
  235. C $ (EIJ.IMAT2(II,IJ)*
  236. C $ JAC.XMAT(II,1)*JAC.XMAT(IJ,2))**2
  237. CC $ (EIJ.IMAT2(II,IJ)*EIJ.IMAT2(IN,IO)*
  238. CC $ JAC.XMAT(II,1)*JAC.XMAT(IJ,2)*
  239. CC $ JAC.XMAT(IN,1)*JAC.XMAT(IO,2))
  240. CC ENDDO
  241. CC ENDDO
  242. C ENDDO
  243. C ENDDO
  244. C XAL=XNUM/XDENO
  245. C CONTRI=0.D0
  246. C DO J=1,IREL
  247. C DO K=1,IREL
  248. C CONTRI=CONTRI+
  249. C $ (EIJK.IMAT3(I,J,K)*
  250. C $ JAC.XMAT(J,1)*JAC.XMAT(K,2))
  251. C ENDDO
  252. C ENDDO
  253. C CONTRI=CONTRI*XAL
  254. CONTRI=0.D0
  255. DO IA=1,IREF
  256. DO IB=1,IREF
  257. DO J=1,IREL
  258. DO K=1,IREL
  259. CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  260. $ EIJK.IMAT3(I,J,K)*
  261. $ JMAJAC.VELCHE(1,1,J,IA,IPJA,ILJA)*
  262. $ JMAJAC.VELCHE(1,1,K,IB,IPJA,ILJA))
  263. C $ JP.XMAT(IA,J)*JP.XMAT(IB,K))
  264. ENDDO
  265. ENDDO
  266. ENDDO
  267. ENDDO
  268. * SEGPRT,JAC
  269. * SEGPRT,JP
  270. * WRITE(IOIMP,*) 'XM=',XM
  271. CONTRI=CONTRI/(2.D0*XM)
  272. * WRITE(IOIMP,*) 'CONTRI=',CONTRI
  273. ELSE
  274. WRITE(IOIMP,*) 'Erreur grave IREL=',IREL
  275. GOTO 9999
  276. ENDIF
  277. ELSEIF (NOMLOI(1:4).EQ.'VNOJ') THEN
  278. CALL CH2INT(NOMLOI(5:5),I,IMPR,IRET)
  279. IF (IRET.NE.0) GOTO 9999
  280. CALL CH2INT(NOMLOI(6:6),J,IMPR,IRET)
  281. IF (IRET.NE.0) GOTO 9999
  282. CALL CH2INT(NOMLOI(7:7),L,IMPR,IRET)
  283. IF (IRET.NE.0) GOTO 9999
  284. IF (IREL.EQ.2) THEN
  285. CONTRI=0.D0
  286. CONTRI=CONTRI-(EIJ.IMAT2(I,J)*
  287. $ JMAJAC.VELCHE(1,1,L,1,IPJA,ILJA))
  288. CONTRI=CONTRI/XM
  289. ELSEIF (IREL.EQ.3) THEN
  290. CONTRI=0.D0
  291. DO IA=1,IREF
  292. DO IB=1,IREF
  293. * DO J=1,IREL
  294. DO K=1,IREL
  295. CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  296. $ EIJK.IMAT3(I,J,K)*
  297. $ JMAJAC.VELCHE(1,1,L,IA,IPJA,ILJA)*
  298. $ JMAJAC.VELCHE(1,1,K,IB,IPJA,ILJA))
  299. CONTRI=CONTRI+(EIJ.IMAT2(IA,IB)*
  300. $ EIJK.IMAT3(I,K,J)*
  301. $ JMAJAC.VELCHE(1,1,K,IA,IPJA,ILJA)*
  302. $ JMAJAC.VELCHE(1,1,L,IB,IPJA,ILJA))
  303. ENDDO
  304. * ENDDO
  305. ENDDO
  306. ENDDO
  307. * SEGPRT,JAC
  308. * SEGPRT,JP
  309. * WRITE(IOIMP,*) 'XM=',XM
  310. CONTRI=CONTRI/(2.D0*XM)
  311. * WRITE(IOIMP,*) 'CONTRI=',CONTRI
  312. ELSE
  313. WRITE(IOIMP,*) 'Erreur grave IREL=',IREL
  314. GOTO 9999
  315. ENDIF
  316. ELSE
  317. WRITE(IOIMP,*) 'Erreur grave NOMLOI=',NOMLOI
  318. GOTO 9999
  319. ENDIF
  320. FC.VELCHE(1,1,1,1,IPFC,ILFC)=
  321. $ FC.VELCHE(1,1,1,1,IPFC,ILFC)+
  322. $ CONTRI
  323. ENDDO
  324. ENDDO
  325. SEGSUP,EIJK
  326. SEGSUP,EIJ
  327. * SEGPRT,FC
  328. SEGSUP,JAC
  329. * SEGSUP,JP
  330. SEGSUP,JT
  331. SEGSUP,G
  332. SEGSUP,IG
  333. *
  334. * Normal termination
  335. *
  336. IRET=0
  337. RETURN
  338. *
  339. * Format handling
  340. *
  341. *
  342. * Error handling
  343. *
  344. 9999 CONTINUE
  345. IRET=1
  346. WRITE(IOIMP,*) 'An error was detected in subroutine ccgnor'
  347. RETURN
  348. *
  349. * End of subroutine CCGNOR
  350. *
  351. END
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  

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