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

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