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

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