Télécharger ccgmus.eso

Retour à la liste

Numérotation des lignes :

  1. C CCGMUS SOURCE GOUNAND 15/06/26 21:15:00 8557
  2. SUBROUTINE CCGMUS(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 : CCGMUS
  9. C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss :
  10. C Taille d'un élément suivant une direction donnée
  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, 30/08/06, version initiale
  25. C HISTORIQUE : v1, 30/08/06, création
  26. C HISTORIQUE : 09/06/2015 : SG ajout de tests de positivite pour eviter
  27. C des pbs dans les SQRTs
  28. C Ajout du cas ou l'espace de reference est de dimension
  29. C strictement inferieure a l'espace reel
  30. C HISTORIQUE :
  31. C***********************************************************************
  32. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  33. C en cas de modification de ce sous-programme afin de faciliter
  34. C la maintenance !
  35. C***********************************************************************
  36. -INC CCOPTIO
  37. -INC CCREEL
  38. CBEGININCLUDE SMCHAEL
  39. SEGMENT MCHAEL
  40. POINTEUR IMACHE(N1).MELEME
  41. POINTEUR ICHEVA(N1).MCHEVA
  42. ENDSEGMENT
  43. SEGMENT MCHEVA
  44. REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
  45. ENDSEGMENT
  46. SEGMENT LCHEVA
  47. POINTEUR LISCHE(NBCHE).MCHEVA
  48. ENDSEGMENT
  49. CENDINCLUDE SMCHAEL
  50. INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
  51. POINTEUR FC.MCHEVA
  52. POINTEUR LCOF.LCHEVA
  53. POINTEUR JMAJAC.MCHEVA
  54. POINTEUR JMIJAC.MCHEVA
  55. POINTEUR JDTJAC.MCHEVA
  56. POINTEUR JMAREG.MCHEVA
  57. POINTEUR JDIAMA.MCHEVA
  58. POINTEUR JPC.MCHEVA
  59. POINTEUR JVIT.MCHEVA
  60. POINTEUR JRHO.MCHEVA
  61. POINTEUR JMU.MCHEVA
  62. POINTEUR JPEC.MCHEVA
  63. CHARACTER*8 NOMLOI
  64. INTEGER ICOF
  65. *
  66. -INC TMXMAT
  67. POINTEUR A.MXMAT
  68. POINTEUR AP.MXMAT
  69. POINTEUR JMA.MXMAT
  70. POINTEUR KJMA.MXMAT
  71. POINTEUR JM1.MXMAT,J.MXMAT
  72. POINTEUR K.MXMAT
  73. *
  74. SEGMENT MVIT
  75. POINTEUR MVCOMP(IDIM).MCHEVA
  76. ENDSEGMENT
  77. LOGICAL LCRIT,LFIRST
  78. LOGICAL LRRHO,LRMU,LRPEC
  79. *
  80. INTEGER IMPR,IRET
  81. *
  82. * Executable statements
  83. *
  84. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgmus'
  85. XPETI=SQRT(XPETIT)
  86. CALL CH2INT(NOMLOI(7:7),IMETH,IMPR,IRET)
  87. IF (IRET.NE.0) GOTO 9999
  88. CALL CH2INT(NOMLOI(8:8),IDIMI,IMPR,IRET)
  89. IF (IRET.NE.0) GOTO 9999
  90. NLFC=FC.VELCHE(/6)
  91. NPFC=FC.VELCHE(/5)
  92. ICOF=0
  93. *
  94. * Récupération de RHO et MU
  95. *
  96. ICOF=ICOF+1
  97. JRHO=LCOF.LISCHE(ICOF)
  98. NLRO=JRHO.VELCHE(/6)
  99. NPRO=JRHO.VELCHE(/5)
  100. ICOF=ICOF+1
  101. JMU=LCOF.LISCHE(ICOF)
  102. NLMU=JMU.VELCHE(/6)
  103. NPMU=JMU.VELCHE(/5)
  104. *
  105. * Récupération des coefficients de la metrique
  106. *
  107. ICOF1=ICOF+1
  108. SEGINI MVIT
  109. DO IIDIM=1,IDIM
  110. ICOF=ICOF+1
  111. JVIT=LCOF.LISCHE(ICOF)
  112. IF (ICOF.EQ.ICOF1) THEN
  113. NLJV=JVIT.VELCHE(/6)
  114. NPJV=JVIT.VELCHE(/5)
  115. ELSE
  116. NLJV2=JVIT.VELCHE(/6)
  117. NPJV2=JVIT.VELCHE(/5)
  118. IF (NLJV2.NE.NLJV.OR.NPJV2.NE.NPJV) THEN
  119. WRITE(IOIMP,*) 'Erreur grave dims JVIT'
  120. GOTO 9999
  121. ENDIF
  122. ENDIF
  123. MVIT.MVCOMP(IIDIM)=JVIT
  124. ENDDO
  125. *
  126. * Récupération du Péclet critique
  127. *
  128. ICOF=ICOF+1
  129. JPEC=LCOF.LISCHE(ICOF)
  130. NLPEC=JPEC.VELCHE(/6)
  131. NPPEC=JPEC.VELCHE(/5)
  132. *
  133. * sg : 09/06/2015
  134. * Verification que rho, mu et pec sont positifs.
  135. * Meme si les champs sont positifs, il peut y avoir un souci
  136. * avec les quadratiques aux points de Gauss
  137. * 220 0
  138. *Attention: une des composantes du champ est négative ou nulle
  139. * 220 0
  140. *sa puissance est mise à zéro
  141. LRRHO=.FALSE.
  142. DO ILRO=1,NLRO
  143. DO IPRO=1,NPRO
  144. XRHO=JRHO.VELCHE(1,1,1,1,IPRO,ILRO)
  145. IF (XRHO.LT.XZERO) THEN
  146. IF (.NOT.LRRHO) THEN
  147. LRRHO=.TRUE.
  148. WRITE(IOIMP,*) 'ccgmus.eso : RHO'
  149. CALL ERREUR(220)
  150. ENDIF
  151. ENDIF
  152. ENDDO
  153. ENDDO
  154. LRMU=.FALSE.
  155. DO ILMU=1,NLMU
  156. DO IPMU=1,NPMU
  157. XMU=JMU.VELCHE(1,1,1,1,IPMU,ILMU)
  158. IF (XMU.LT.XZERO) THEN
  159. IF (.NOT.LRMU) THEN
  160. LRMU=.TRUE.
  161. WRITE(IOIMP,*) 'ccgmus.eso : MU'
  162. CALL ERREUR(220)
  163. ENDIF
  164. ENDIF
  165. ENDDO
  166. ENDDO
  167. LRPEC=.FALSE.
  168. DO ILPEC=1,NLPEC
  169. DO IPPEC=1,NPPEC
  170. XPEC=JPEC.VELCHE(1,1,1,1,IPPEC,ILPEC)
  171. IF (XPEC.LT.XZERO) THEN
  172. IF (.NOT.LRPEC) THEN
  173. LRPEC=.TRUE.
  174. WRITE(IOIMP,*) 'ccgmus.eso : PEC'
  175. CALL ERREUR(220)
  176. ENDIF
  177. ENDIF
  178. ENDDO
  179. ENDDO
  180. *dbg IF (CONTRI.LT.0.D0) THEN
  181. *dbg WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU
  182. *dbg $ ,' XV=',XV,' XH=',XH,' XPEC=',XPEC
  183. *dbg WRITE(IOIMP,*) 'NPRO=',NPRO,'ILRO=',ILRO
  184. *dbg do ii=1,npro
  185. *dbg XRHO=JRHO.VELCHE(1,1,1,1,Ii,ILRO)
  186. *dbg WRITE(IOIMP,*) 'ii XRHO=',ii,XRHO
  187. *dbg enddo
  188. *dbg GOTO 9999
  189. *dbg ENDIF
  190. *
  191. ICOF=ICOF+1
  192. JMAJAC=LCOF.LISCHE(ICOF)
  193. C NLJA=JMAJAC.VELCHE(/6)
  194. C NPJA=JMAJAC.VELCHE(/5)
  195. IREF=JMAJAC.VELCHE(/4)
  196. IREL=JMAJAC.VELCHE(/3)
  197. *dbg WRITE(IOIMP,*) 'IREF=',IREF,' IREL=',IREL
  198. *dbg IF (IREL.NE.IDIM.OR.IREF.NE.IDIM) THEN
  199. IF (IREL.NE.IDIM) THEN
  200. WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  201. GOTO 9999
  202. ENDIF
  203. *
  204. ICOF=ICOF+1
  205. JMIJAC=LCOF.LISCHE(ICOF)
  206. IF (JMIJAC.EQ.0) THEN
  207. WRITE(IOIMP,*) 'Erreur JMIJAC=0'
  208. GOTO 9999
  209. ENDIF
  210. NLJI=JMIJAC.VELCHE(/6)
  211. NPJI=JMIJAC.VELCHE(/5)
  212. IREL2=JMIJAC.VELCHE(/4)
  213. IREF2=JMIJAC.VELCHE(/3)
  214. *
  215. IF (IREL2.NE.IREL.OR.IREF2.NE.IREF) THEN
  216. WRITE(IOIMP,*) 'Erreur dims JMIJAC'
  217. GOTO 9999
  218. ENDIF
  219. *
  220. ICOF=ICOF+1
  221. ICOF=ICOF+1
  222. JMAREG=LCOF.LISCHE(ICOF)
  223. NLJR=JMAREG.VELCHE(/6)
  224. NPJR=JMAREG.VELCHE(/5)
  225. I1 =JMAREG.VELCHE(/4)
  226. I2 =JMAREG.VELCHE(/3)
  227. IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
  228. $ THEN
  229. WRITE(IOIMP,*) 'Erreur dims JMAREG'
  230. GOTO 9999
  231. ENDIF
  232. ICOF=ICOF+1
  233. JDIAMA=LCOF.LISCHE(ICOF)
  234. NLJD=JDIAMA.VELCHE(/6)
  235. NPJD=JDIAMA.VELCHE(/5)
  236. I1 =JDIAMA.VELCHE(/4)
  237. I2 =JDIAMA.VELCHE(/3)
  238. IF ((NLJD.NE.1).OR.(NPJD.NE.1).OR.(I1.NE.1).OR.(I2.NE.1))
  239. $ THEN
  240. WRITE(IOIMP,*) 'Erreur dims JDIAMA'
  241. GOTO 9999
  242. ENDIF
  243. XDIAMA=JDIAMA.VELCHE(1,1,1,1,1,1)
  244. *dbg
  245. *dbg WRITE(IOIMP,*) 'XDIAMA=',XDIAMA
  246. *
  247. * Objets temporaires
  248. *
  249. LDIM1=IREL
  250. LDIM2=1
  251. SEGINI,A
  252. SEGINI,AP
  253. LDIM1=IREF
  254. LDIM2=1
  255. SEGINI,JMA
  256. LDIM1=IREF
  257. LDIM2=1
  258. SEGINI,KJMA
  259. LDIM1=IREL
  260. LDIM2=IREF
  261. SEGINI,J
  262. LDIM1=IREF
  263. LDIM2=IREL
  264. SEGINI,JM1
  265. LDIM1=IREF
  266. LDIM2=IREF
  267. SEGINI,K
  268. *
  269. * Calcul de la métrique des éléments réguliers
  270. *
  271. CALL MAMA(JMAREG.VELCHE,IREF,IREF,
  272. $ 'COPIE ',K.XMAT,IREF,IREF,
  273. $ IMPR,IRET)
  274. IF (IRET.NE.0) GOTO 9999
  275. * SEGPRT,H
  276. *
  277. DO ILFC=1,NLFC
  278. IF (NLRO.EQ.1) THEN
  279. ILRO=1
  280. ELSE
  281. ILRO=ILFC
  282. ENDIF
  283. IF (NLMU.EQ.1) THEN
  284. ILMU=1
  285. ELSE
  286. ILMU=ILFC
  287. ENDIF
  288. IF (NLJV.EQ.1) THEN
  289. ILJV=1
  290. ELSE
  291. ILJV=ILFC
  292. ENDIF
  293. IF (NLPEC.EQ.1) THEN
  294. ILPEC=1
  295. ELSE
  296. ILPEC=ILFC
  297. ENDIF
  298. C IF (NLJA.EQ.1) THEN
  299. C ILJA=1
  300. C ELSE
  301. C ILJA=ILFC
  302. C ENDIF
  303. IF (NLJI.EQ.1) THEN
  304. ILJI=1
  305. ELSE
  306. ILJI=ILFC
  307. ENDIF
  308. *
  309. DO IPFC=1,NPFC
  310. IF (NPRO.EQ.1) THEN
  311. IPRO=1
  312. ELSE
  313. IPRO=IPFC
  314. ENDIF
  315. IF (NPMU.EQ.1) THEN
  316. IPMU=1
  317. ELSE
  318. IPMU=IPFC
  319. ENDIF
  320. IF (NPJV.EQ.1) THEN
  321. IPJV=1
  322. ELSE
  323. IPJV=IPFC
  324. ENDIF
  325. IF (NPPEC.EQ.1) THEN
  326. IPPEC=1
  327. ELSE
  328. IPPEC=IPFC
  329. ENDIF
  330. C IF (NPJA.EQ.1) THEN
  331. C IPJA=1
  332. C ELSE
  333. C IPJA=IPFC
  334. C ENDIF
  335. IF (NPJI.EQ.1) THEN
  336. IPJI=1
  337. ELSE
  338. IPJI=IPFC
  339. ENDIF
  340. XRHO=JRHO.VELCHE(1,1,1,1,IPRO,ILRO)
  341. XMU=JMU.VELCHE(1,1,1,1,IPMU,ILMU)
  342. XPEC=JPEC.VELCHE(1,1,1,1,IPPEC,ILPEC)
  343. IF (LRRHO) XRHO=MAX(XRHO,XZERO)
  344. IF (LRMU) XMU=MAX(XMU,XZERO)
  345. IF (LRPEC) XPEC=MAX(XPEC,XZERO)
  346.  
  347. *dbg LFIRST=(ILFC.EQ.1).AND.(IPFC.EQ.1)
  348. LFIRST=.FALSE.
  349. IF (LFIRST) WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU,' XPEC='
  350. $ ,XPEC
  351. *
  352. * Copie de la vitesse
  353. *
  354. DO IIDIM=1,IDIM
  355. JVIT=MVIT.MVCOMP(IIDIM)
  356. A.XMAT(IIDIM,1)=JVIT.VELCHE(1,1,1,1,IPJV,ILJV)
  357. ENDDO
  358. IF (LFIRST) WRITE(IOIMP,*) 'A ',(A.XMAT(ii,1),ii=1,IDIM)
  359. * SEGPRT,A
  360. *
  361. * Copie de l'inverse du jacobien ref->reel
  362. *
  363. CALL MAMA(JMIJAC.VELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL,
  364. $ 'COPIE ',
  365. $ JM1.XMAT,IREF,IREL,
  366. $ IMPR,IRET)
  367. IF (IRET.NE.0) GOTO 9999
  368. C SEGPRT,JM1
  369. *
  370. * Copie du jacobien ref->reel
  371. *
  372. CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJI,ILJI),IREL,IREF,
  373. $ 'COPIE ',
  374. $ J.XMAT,IREL,IREF,
  375. $ IMPR,IRET)
  376. IF (IRET.NE.0) GOTO 9999
  377. C SEGPRT,J
  378. *
  379. * Vecteur A dans le repère de référence
  380. *
  381. CALL MAMAMA(JM1.XMAT,IREF,IREL,A.XMAT,IREL,1,
  382. $ 'FOIS ',JMA.XMAT,IREF,1,
  383. $ IMPR,IRET)
  384. IF (IRET.NE.0) GOTO 9999
  385. *
  386. * Vecteur A' dans le repère réel (projection sur le domaine considéré)
  387. *
  388. CALL MAMAMA(J.XMAT,IREL,IREF,JMA.XMAT,IREF,1,
  389. $ 'FOIS ',AP.XMAT,IREL,1,
  390. $ IMPR,IRET)
  391. IF (IRET.NE.0) GOTO 9999
  392. IF (LFIRST) WRITE(IOIMP,*) 'AP ',(AP.XMAT(ii,1),ii=1,IDIM)
  393. *
  394. * Vecteur A dans le repère de l'élément régulier
  395. *
  396. CALL MAMAMA(K.XMAT,IREF,IREF,JMA.XMAT,IREF,1,
  397. $ 'FOIS ',KJMA.XMAT,IREF,1,
  398. $ IMPR,IRET)
  399. IF (IRET.NE.0) GOTO 9999
  400. *
  401. * Normes de A
  402. *
  403. XAREG2=XZERO
  404. DO IIREF=1,IREF
  405. XAREG2=XAREG2+(KJMA.XMAT(IIREF,1)**2)
  406. ENDDO
  407. XA2=XZERO
  408. DO IIDIM=1,IDIM
  409. XA2=XA2+(A.XMAT(IIDIM,1)**2)
  410. ENDDO
  411. XAP2=XZERO
  412. DO IIDIM=1,IDIM
  413. XAP2=XAP2+(AP.XMAT(IIDIM,1)**2)
  414. ENDDO
  415. IF (LFIRST) WRITE(IOIMP,*) 'XAREG2=',XAREG2,' XA2=',XA2,
  416. $ ' XAP2=',XAP2
  417. *
  418. *old IF (XA2.LT.XPETI) THEN
  419. IF (XAP2.LT.XPETI) THEN
  420. CONTRI=XZERO
  421. ELSE
  422. *old XV=SQRT(XA2)
  423. *old* Pondération par un diamètre externe !
  424. *old* XH=SQRT(XA2/XAREG2)
  425. *old XH=SQRT(XA2/XAREG2)*XDIAMA
  426. * Pour que le décentrement marche sur une surface
  427. XV=SQRT(XAP2)
  428. XH=SQRT(XAP2/XAREG2)*XDIAMA
  429. *
  430. IF (IMETH.EQ.1) THEN
  431. * Diffusion artificielle UPWIND
  432. CONTRI=((XRHO*XV*XH)/XPEC)
  433. CONTRI=SQRT(CONTRI)
  434. IF (IDIMI.GT.0) THEN
  435. CONTRI=CONTRI*
  436. *old $ (A.XMAT(IDIMI,1)/XV)
  437. $ (AP.XMAT(IDIMI,1)/XV)
  438. ENDIF
  439. IF (LFIRST) WRITE(IOIMP,*) 'CONTRI=',CONTRI
  440. ELSEIF (IMETH.EQ.2) THEN
  441. * Diffusion artificielle comme Alberto (equiv critical approx)
  442. * Doit-on décentrer ?
  443. LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU))
  444. IF (LCRIT) THEN
  445. CONTRI=((XRHO*XV*XH)/XPEC)-XMU
  446. CONTRI=SQRT(CONTRI)
  447. IF (IDIMI.GT.0) THEN
  448. CONTRI=CONTRI*
  449. *old $ (A.XMAT(IDIMI,1)/XV)
  450. $ (AP.XMAT(IDIMI,1)/XV)
  451. ENDIF
  452. ELSE
  453. CONTRI=XZERO
  454. ENDIF
  455. ELSEIF (IMETH.EQ.3) THEN
  456. * Diffusion artificielle (equiv doubly asymptotic)
  457. LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU*3.D0))
  458. IF (LCRIT) THEN
  459. CONTRI=((XRHO*XV*XH)/XPEC)
  460. ELSE
  461. CONTRI=(((XRHO*XV*XH)/XPEC)**2)/(3.D0*XMU)
  462. ENDIF
  463. CONTRI=SQRT(CONTRI)
  464. IF (IDIMI.GT.0) THEN
  465. CONTRI=CONTRI*
  466. *old $ (A.XMAT(IDIMI,1)/XV)
  467. $ (AP.XMAT(IDIMI,1)/XV)
  468. ENDIF
  469. ELSE
  470. WRITE(IOIMP,*) 'IMETH=',IMETH,' unknown'
  471. GOTO 9999
  472. ENDIF
  473. ENDIF
  474. FC.VELCHE(1,1,1,1,IPFC,ILFC)=
  475. $ FC.VELCHE(1,1,1,1,IPFC,ILFC)+
  476. $ CONTRI
  477. ENDDO
  478. ENDDO
  479. SEGSUP,K
  480. SEGSUP,J
  481. SEGSUP,JM1
  482. SEGSUP,KJMA
  483. SEGSUP,JMA
  484. SEGSUP,AP
  485. SEGSUP,A
  486. SEGSUP,MVIT
  487. *
  488. * Normal termination
  489. *
  490. IRET=0
  491. RETURN
  492. *
  493. * Format handling
  494. *
  495. *
  496. * Error handling
  497. *
  498. 9999 CONTINUE
  499. IRET=1
  500. WRITE(IOIMP,*) 'An error was detected in subroutine ccgmus'
  501. RETURN
  502. *
  503. * End of subroutine CCGMUS
  504. *
  505. END
  506.  
  507.  
  508.  
  509.  
  510.  
  511.  

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