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

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