Télécharger ccgahu.eso

Retour à la liste

Numérotation des lignes :

  1. C CCGAHU SOURCE GOUNAND 07/07/05 21:15:08 5784
  2. SUBROUTINE CCGAHU(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 : CCGAHU
  9. C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss :
  10. C DEDU ADAP avec metrique
  11. C On implémente ici la première partie de la fonctionnelle
  12. C de Huang
  13. C
  14. C LANGAGE : ESOPE
  15. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  16. C mél : gounand@semt2.smts.cea.fr
  17. C***********************************************************************
  18. C APPELES :
  19. C APPELE PAR :
  20. C***********************************************************************
  21. C ENTREES :
  22. C ENTREES/SORTIES :
  23. C SORTIES : -
  24. C TRAVAIL :
  25. C***********************************************************************
  26. C VERSION : v1, 30/04/07, version initiale
  27. C HISTORIQUE : v1, 30/04/07, création
  28. C HISTORIQUE :
  29. C***********************************************************************
  30. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  31. C en cas de modification de ce sous-programme afin de faciliter
  32. C la maintenance !
  33. C***********************************************************************
  34.  
  35. -INC PPARAM
  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 JMET.MCHEVA
  58. POINTEUR JTHE.MCHEVA
  59. POINTEUR JGAM.MCHEVA
  60. CHARACTER*8 NOMLOI
  61. INTEGER ICOF
  62. *
  63. -INC TMXMAT
  64. * Objets temporaires
  65. POINTEUR JAC.MXMAT,JT.MXMAT
  66. POINTEUR G.MXMAT,IG.MXMAT,H.MXMAT,HIG.MXMAT,HM1.MXMAT
  67. POINTEUR MM1.MXMAT,JTM.MXMAT,MJ.MXMAT
  68. * Objets à préconditionner (éventuellement)
  69. POINTEUR A.MXMAT,B.MXMAT,C.MXMAT,D.MXMAT,F.MXMAT,M.MXMAT
  70. *
  71. SEGMENT MCOF
  72. POINTEUR COEF(IDIM,IDIM).MCHEVA
  73. ENDSEGMENT
  74. POINTEUR MET.MCOF
  75. *
  76. LOGICAL LBID,LFIRST
  77. INTEGER LAXSP
  78. REAL*8 DEUPI,XR
  79. REAL*8 XL,XM,XH
  80. *
  81. INTEGER IMPR,IRET
  82. *
  83. * Executable statements
  84. *
  85. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgahu'
  86. C IF (.NOT.(IDIM.EQ.1)) THEN
  87. C WRITE(IOIMP,*) 'IDIM=',IDIM,' ?'
  88. C GOTO 9999
  89. C ENDIF
  90. NLFC=FC.VELCHE(/6)
  91. NPFC=FC.VELCHE(/5)
  92. ICOF=0
  93. *
  94. * Récupération des coefficients de la metrique
  95. *
  96. SEGINI MET
  97. DO IIDIM=1,IDIM
  98. ICOF=ICOF+1
  99. JMET=LCOF.LISCHE(ICOF)
  100. IF (ICOF.EQ.1) THEN
  101. NLJM=JMET.VELCHE(/6)
  102. NPJM=JMET.VELCHE(/5)
  103. ELSE
  104. NLJM2=JMET.VELCHE(/6)
  105. NPJM2=JMET.VELCHE(/5)
  106. IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN
  107. WRITE(IOIMP,*) 'Erreur grave dims JMET'
  108. GOTO 9999
  109. ENDIF
  110. ENDIF
  111. MET.COEF(IIDIM,IIDIM)=JMET
  112. ENDDO
  113. DO IIDIM=1,IDIM
  114. NJ=IDIM-IIDIM
  115. IF (NJ.GE.1) THEN
  116. DO JIDIM=IIDIM+1,IDIM
  117. ICOF=ICOF+1
  118. JMET=LCOF.LISCHE(ICOF)
  119. NLJM2=JMET.VELCHE(/6)
  120. NPJM2=JMET.VELCHE(/5)
  121. IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN
  122. WRITE(IOIMP,*) 'Erreur grave dims JMET2'
  123. GOTO 9999
  124. ENDIF
  125. MET.COEF(IIDIM,JIDIM)=JMET
  126. ENDDO
  127. ENDIF
  128. ENDDO
  129. *
  130. ICOF=ICOF+1
  131. JTHE=LCOF.LISCHE(ICOF)
  132. NLJT=JTHE.VELCHE(/6)
  133. NPJT=JTHE.VELCHE(/5)
  134. I1 =JTHE.VELCHE(/4)
  135. I2 =JTHE.VELCHE(/3)
  136. IF ((I1.NE.1).OR.(I2.NE.1))
  137. $ THEN
  138. WRITE(IOIMP,*) 'Erreur dims JTHE'
  139. GOTO 9999
  140. ENDIF
  141. *
  142. ICOF=ICOF+1
  143. JGAM=LCOF.LISCHE(ICOF)
  144. NLJG=JTHE.VELCHE(/6)
  145. NPJG=JTHE.VELCHE(/5)
  146. I1 =JTHE.VELCHE(/4)
  147. I2 =JTHE.VELCHE(/3)
  148. IF ((I1.NE.1).OR.(I2.NE.1))
  149. $ THEN
  150. WRITE(IOIMP,*) 'Erreur dims JGAM'
  151. GOTO 9999
  152. ENDIF
  153. *
  154. ICOF=ICOF+1
  155. JMAJAC=LCOF.LISCHE(ICOF)
  156. NLJA=JMAJAC.VELCHE(/6)
  157. NPJA=JMAJAC.VELCHE(/5)
  158. IREF=JMAJAC.VELCHE(/4)
  159. IREL=JMAJAC.VELCHE(/3)
  160. *
  161. IF (IREL.NE.IDIM) THEN
  162. WRITE(IOIMP,*) 'Erreur dims JMAJAC'
  163. GOTO 9999
  164. ENDIF
  165. *
  166. ICOF=ICOF+1
  167. ICOF=ICOF+1
  168. ICOF=ICOF+1
  169. JMAREG=LCOF.LISCHE(ICOF)
  170. NLJR=JMAREG.VELCHE(/6)
  171. NPJR=JMAREG.VELCHE(/5)
  172. I1 =JMAREG.VELCHE(/4)
  173. I2 =JMAREG.VELCHE(/3)
  174. IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
  175. $ THEN
  176. WRITE(IOIMP,*) 'Erreur dims JMAREG'
  177. GOTO 9999
  178. ENDIF
  179. *
  180. * Objets temporaires et à préconditionner
  181. *
  182. LDIM1=IREL
  183. LDIM2=IREF
  184. SEGINI,JAC
  185. SEGINI,MJ
  186. SEGINI,A
  187. SEGINI,B
  188. LDIM1=IREF
  189. LDIM2=IREL
  190. SEGINI,JT
  191. SEGINI,JTM
  192. LDIM1=IREF
  193. LDIM2=IREF
  194. SEGINI,G
  195. SEGINI,IG
  196. SEGINI,H
  197. SEGINI,HM1
  198. SEGINI,HIG
  199. SEGINI,C
  200. LDIM1=IREL
  201. LDIM2=IREL
  202. SEGINI,M
  203. SEGINI,MM1
  204. SEGINI,D
  205. SEGINI,F
  206. *
  207. * Calcul de la métrique des éléments réguliers
  208. *
  209. CALL MAMA(JMAREG.VELCHE,IREF,IREF,
  210. $ 'JTJ ',H.XMAT,IREF,IREF,
  211. $ IMPR,IRET)
  212. IF (IRET.NE.0) GOTO 9999
  213. * SEGPRT,H
  214. LBID=.FALSE.
  215. CALL GEOLI2(IREF,1,1,H.XMAT,HM1.XMAT,DETH,LBID,IMPR,IRET)
  216. IF (IRET.NE.0) GOTO 9999
  217. *!!! XH=1.0D0
  218. XH=1.D0/(SQRT(DETH))
  219. *!! XH=1.D0/(DETH)
  220. *
  221. LFIRST = .FALSE.
  222. DO ILFC=1,NLFC
  223. IF (NLJM.EQ.1) THEN
  224. ILJM=1
  225. ELSE
  226. ILJM=ILFC
  227. ENDIF
  228. IF (NLJA.EQ.1) THEN
  229. ILJA=1
  230. ELSE
  231. ILJA=ILFC
  232. ENDIF
  233. IF (NLJT.EQ.1) THEN
  234. ILJT=1
  235. ELSE
  236. ILJT=ILFC
  237. ENDIF
  238. IF (NLJG.EQ.1) THEN
  239. ILJG=1
  240. ELSE
  241. ILJG=ILFC
  242. ENDIF
  243. DO IPFC=1,NPFC
  244. IF (NPJM.EQ.1) THEN
  245. IPJM=1
  246. ELSE
  247. IPJM=IPFC
  248. ENDIF
  249. IF (NPJA.EQ.1) THEN
  250. IPJA=1
  251. ELSE
  252. IPJA=IPFC
  253. ENDIF
  254. IF (NPJT.EQ.1) THEN
  255. IPJT=1
  256. ELSE
  257. IPJT=IPFC
  258. ENDIF
  259. IF (NPJG.EQ.1) THEN
  260. IPJG=1
  261. ELSE
  262. IPJG=IPFC
  263. ENDIF
  264. *
  265. * Récupération de l'exposant multiplié par 2
  266. * On suppose que c'est un entier
  267. *
  268. XTHE=JTHE.VELCHE(1,1,1,1,IPJT,ILJT)
  269. XGAM=JGAM.VELCHE(1,1,1,1,IPJG,ILJG)
  270. XCO1=XTHE*0.5D0
  271. IXP1=NINT(IREF*XGAM)
  272. XCO2=(1.D0-XTHE)*((SQRT(DBLE(IREF)))**IXP1)
  273. *!!!
  274. *!!! XCO2=1.D0
  275. IXP2=NINT(1-XGAM)
  276. IF (LFIRST) THEN
  277. WRITE(IOIMP,*) 'XTHE=',XTHE
  278. WRITE(IOIMP,*) 'XGAM=',XGAM
  279. WRITE(IOIMP,*) 'XCO1=',XCO1
  280. WRITE(IOIMP,*) 'XCO2=',XCO2
  281. WRITE(IOIMP,*) 'IXP1=',IXP1,'/2'
  282. WRITE(IOIMP,*) 'IXP2=',IXP2
  283. LFIRST=.FALSE.
  284. ENDIF
  285. *
  286. * Copie des coefficients de la métrique
  287. *
  288. DO IIDIM=1,IDIM
  289. JMET=MET.COEF(IIDIM,IIDIM)
  290. M.XMAT(IIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
  291. ENDDO
  292. DO IIDIM=1,IIDIM
  293. NJ=IDIM-IIDIM
  294. IF (NJ.GE.1) THEN
  295. DO JIDIM=IIDIM+1,IDIM
  296. JMET=MET.COEF(IIDIM,JIDIM)
  297. M.XMAT(IIDIM,JIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
  298. M.XMAT(JIDIM,IIDIM)=JMET.VELCHE(1,1,1,1,IPJM,ILJM)
  299. ENDDO
  300. ENDIF
  301. ENDDO
  302. * SEGPRT,M
  303. * Calcul de l'inverse
  304. *
  305. LBID=.FALSE.
  306. CALL GEOLI2(IREL,1,1,M.XMAT,MM1.XMAT,DETM,LBID,IMPR,IRET)
  307. IF (IRET.NE.0) GOTO 9999
  308. XDM=1.D0/(SQRT(DETM))
  309. *1 XDM=1.D0
  310. ****** XDM=SQRT(DETM)
  311. *
  312. * Copie du jacobien
  313. *
  314. CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  315. $ 'COPIE ',
  316. $ JAC.XMAT,IREL,IREF,
  317. $ IMPR,IRET)
  318. IF (IRET.NE.0) GOTO 9999
  319. * SEGPRT,JAC
  320. *
  321. * Calcul de la métrique G
  322. *
  323. * Calcul de Jt
  324. CALL MAMA(JAC.XMAT,IREL,IREF,
  325. $ 'TRANSPOS',JT.XMAT,IREF,IREL,
  326. $ IMPR,IRET)
  327. IF (IRET.NE.0) GOTO 9999
  328. * Calcul de MJ
  329. CALL MAMAMA(M.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
  330. $ 'FOIS ',MJ.XMAT,IREL,IREF,
  331. $ IMPR,IRET)
  332. IF (IRET.NE.0) GOTO 9999
  333. * Calcul de JTM
  334. CALL MAMA(MJ.XMAT,IREL,IREF,
  335. $ 'TRANSPOS',JTM.XMAT,IREF,IREL,
  336. $ IMPR,IRET)
  337. IF (IRET.NE.0) GOTO 9999
  338. * Calcul de G=JtMJ
  339. CALL MAMAMA(JT.XMAT,IREF,IREL,MJ.XMAT,IREL,IREF,
  340. $ 'FOIS ',G.XMAT,IREF,IREF,
  341. $ IMPR,IRET)
  342. IF (IRET.NE.0) GOTO 9999
  343. * Calcul de l'inverse, du déterminant et trace de l'inverse de g
  344. LBID=.FALSE.
  345. CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET)
  346. IF (IRET.NE.0) GOTO 9999
  347. XM=SQRT(DETG)
  348. *! WRITE(IOIMP,*) 'XM=',XM
  349. * Calcul de hg-1 et de sa trace
  350. CALL MAMAMA(H.XMAT,IREF,IREF,IG.XMAT,IREF,IREF,
  351. $ 'FOIS ',HIG.XMAT,IREF,IREF,IMPR,IRET)
  352. IF (IRET.NE.0) GOTO 9999
  353. CALL MARE(HIG.XMAT,IREF,IREF,'TRACE ',
  354. $ XL,IMPR,IRET)
  355. IF (IRET.NE.0) GOTO 9999
  356. * Calcul de A = (J+)t = MJg-1 (pseudo-inverse)
  357. CALL MAMAMA(MJ.XMAT,IREL,IREF,IG.XMAT,IREF,IREF,
  358. $ 'FOIS ',A.XMAT,IREL,IREF,IMPR,IRET)
  359. IF (IRET.NE.0) GOTO 9999
  360. * Calcul de B = (J+)t hg-1 = MJg-1hg-1
  361. CALL MAMAMA(A.XMAT,IREL,IREF,HIG.XMAT,IREF,IREF,
  362. $ 'FOIS ',B.XMAT,IREL,IREF,IMPR,IRET)
  363. IF (IRET.NE.0) GOTO 9999
  364. * Calcul de C = g-1hg-1
  365. CALL MAMAMA(IG.XMAT,IREF,IREF,HIG.XMAT,IREF,IREF,
  366. $ 'FOIS ',C.XMAT,IREF,IREF,IMPR,IRET)
  367. IF (IRET.NE.0) GOTO 9999
  368. * Calcul de D = MJ g-1hg-1 JtM = B JtM
  369. CALL MAMAMA(B.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
  370. $ 'FOIS ',D.XMAT,IREL,IREL,IMPR,IRET)
  371. IF (IRET.NE.0) GOTO 9999
  372. * Calcul de F = MJ g-1 JtM = A JtM
  373. CALL MAMAMA(A.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
  374. $ 'FOIS ',F.XMAT,IREL,IREL,IMPR,IRET)
  375. IF (IRET.NE.0) GOTO 9999
  376. *
  377. * Calcul de la fonctionnelle
  378. *
  379. IF (NOMLOI.EQ.'AHUF ') THEN
  380. XL1=SQRT(XL)**IXP1
  381. FONC1=XL1*XM
  382. *
  383. XM1=XM**IXP2
  384. FONC2=XM1
  385. *
  386. CONTRI=(XCO1*FONC1)+(XCO2*FONC2*(XH**(IXP2-1)))
  387. *!!!!
  388. CONTRI=CONTRI*XDM
  389. *
  390. * Calcul du résidu
  391. *
  392. ELSEIF (NOMLOI(1:4).EQ.'AHUR') THEN
  393. CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  394. IF (IRET.NE.0) GOTO 9999
  395. CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  396. IF (IRET.NE.0) GOTO 9999
  397. K=IDIM1
  398. L=IDIM2
  399. XL1=SQRT(XL)**IXP1
  400. XL2=IXP1*(SQRT(XL)**(IXP1-2))
  401. AKL1=XL1*XM*A.XMAT(K,L)
  402. AKL2=-1.D0*XL2*XM*B.XMAT(K,L)
  403. AKL=AKL1+AKL2
  404. *
  405. XM2=IXP2*(XM**IXP2)
  406. CKL=XM2*A.XMAT(K,L)
  407. *
  408. CONTRI=(XCO1*AKL)+(XCO2*CKL*(XH**(IXP2-1)))
  409. *!!!!
  410. CONTRI=CONTRI*XDM
  411. *
  412. * Calcul du jacobien
  413. *
  414. ELSEIF (NOMLOI(1:4).EQ.'AHUJ') THEN
  415. CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  416. IF (IRET.NE.0) GOTO 9999
  417. CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  418. IF (IRET.NE.0) GOTO 9999
  419. CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET)
  420. IF (IRET.NE.0) GOTO 9999
  421. CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET)
  422. IF (IRET.NE.0) GOTO 9999
  423. I=IDIM1
  424. J=IDIM2
  425. K=IDIM3
  426. L=IDIM4
  427. XL1=SQRT(XL)**IXP1
  428. XL2=IXP1*(SQRT(XL)**(IXP1-2))
  429. XL3=IXP1*(IXP1-2)*(SQRT(XL)**(IXP1-4))
  430. B11=XM*XL1*A.XMAT(I,J)*A.XMAT(K,L)
  431. B12=-1.D0*XL2*XM*B.XMAT(I,J)*A.XMAT(K,L)
  432. B13=XM*XL1*M.XMAT(K,I)*IG.XMAT(J,L)
  433. B14=-1.D0*XM*XL1*F.XMAT(K,I)*IG.XMAT(J,L)
  434. B15=-1.D0*XM*XL1*A.XMAT(K,J)*A.XMAT(I,L)
  435. BIJKL1=B11+B12+B13+B14+B15
  436. B20=+1.D0*XL3*XM*B.XMAT(I,J)*B.XMAT(K,L)
  437. B21=-1.D0*XL2*XM*A.XMAT(I,J)*B.XMAT(K,L)
  438. B22=-1.D0*XL2*XM*M.XMAT(K,I)*C.XMAT(J,L)
  439. B23=+1.D0*XL2*XM*F.XMAT(K,I)*C.XMAT(J,L)
  440. B24=+1.D0*XL2*XM*A.XMAT(K,J)*B.XMAT(I,L)
  441. B25=+1.D0*XL2*XM*D.XMAT(K,I)*IG.XMAT(J,L)
  442. B26=+1.D0*XL2*XM*B.XMAT(K,J)*A.XMAT(I,L)
  443. BIJKL2=B20+B21+B22+B23+B24+B25+B26
  444. BIJKL=BIJKL1+BIJKL2
  445. *
  446. XM2=IXP2*(XM**IXP2)
  447. XM3=IXP2*IXP2*(XM**IXP2)
  448. D1=XM3*A.XMAT(I,J)*A.XMAT(K,L)
  449. D2=XM2*M.XMAT(K,I)*IG.XMAT(J,L)
  450. D3=-1.D0*XM2*F.XMAT(K,I)*IG.XMAT(J,L)
  451. D4=-1.D0*XM2*A.XMAT(K,J)*A.XMAT(I,L)
  452. DIJKL=D1+D2+D3+D4
  453. *
  454. CONTRI=(XCO1*BIJKL)+(XCO2*DIJKL*(XH**(IXP2-1)))
  455. *!!!!
  456. CONTRI=CONTRI*XDM
  457. *! WRITE(IOIMP,*) 'CONTRI=',CONTRI
  458. ELSE
  459. WRITE(IOIMP,*) 'Erreur grave'
  460. GOTO 9999
  461. ENDIF
  462. FC.VELCHE(1,1,1,1,IPFC,ILFC)=
  463. $ FC.VELCHE(1,1,1,1,IPFC,ILFC)+
  464. $ CONTRI
  465. ENDDO
  466. ENDDO
  467. SEGSUP,F
  468. SEGSUP,D
  469. SEGSUP,MM1
  470. SEGSUP,M
  471. SEGSUP,C
  472. SEGSUP,HIG
  473. SEGSUP,HM1
  474. SEGSUP,H
  475. SEGSUP,IG
  476. SEGSUP,G
  477. SEGSUP,JTM
  478. SEGSUP,JT
  479. SEGSUP,B
  480. SEGSUP,A
  481. SEGSUP,MJ
  482. SEGSUP,JAC
  483. SEGSUP,MET
  484. *
  485. * Normal termination
  486. *
  487. IRET=0
  488. RETURN
  489. *
  490. * Format handling
  491. *
  492. *
  493. * Error handling
  494. *
  495. 9999 CONTINUE
  496. IRET=1
  497. WRITE(IOIMP,*) 'An error was detected in subroutine ccgahu'
  498. RETURN
  499. *
  500. * End of subroutine CCGAHU
  501. *
  502. END
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  

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