Télécharger ccgahu.eso

Retour à la liste

Numérotation des lignes :

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

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