Télécharger ccgahu.eso

Retour à la liste

Numérotation des lignes :

ccgahu
  1. C CCGAHU SOURCE PV 22/04/22 23:00:16 11344
  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 LBID,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. LBID=.FALSE.
  206. CALL GEOLI2(IREF,1,1,H.XMAT,HM1.XMAT,DETH,LBID,IMPR,IRET)
  207. IF (IRET.NE.0) GOTO 9999
  208. *!!! XH=1.0D0
  209. XH=1.D0/(SQRT(DETH(1,1)))
  210. *!! XH=1.D0/(DETH(1,1))
  211. *
  212. LFIRST = .FALSE.
  213. DO ILFC=1,NLFC
  214. IF (NLJM.EQ.1) THEN
  215. ILJM=1
  216. ELSE
  217. ILJM=ILFC
  218. ENDIF
  219. IF (NLJA.EQ.1) THEN
  220. ILJA=1
  221. ELSE
  222. ILJA=ILFC
  223. ENDIF
  224. IF (NLJT.EQ.1) THEN
  225. ILJT=1
  226. ELSE
  227. ILJT=ILFC
  228. ENDIF
  229. IF (NLJG.EQ.1) THEN
  230. ILJG=1
  231. ELSE
  232. ILJG=ILFC
  233. ENDIF
  234. DO IPFC=1,NPFC
  235. IF (NPJM.EQ.1) THEN
  236. IPJM=1
  237. ELSE
  238. IPJM=IPFC
  239. ENDIF
  240. IF (NPJA.EQ.1) THEN
  241. IPJA=1
  242. ELSE
  243. IPJA=IPFC
  244. ENDIF
  245. IF (NPJT.EQ.1) THEN
  246. IPJT=1
  247. ELSE
  248. IPJT=IPFC
  249. ENDIF
  250. IF (NPJG.EQ.1) THEN
  251. IPJG=1
  252. ELSE
  253. IPJG=IPFC
  254. ENDIF
  255. *
  256. * Récupération de l'exposant multiplié par 2
  257. * On suppose que c'est un entier
  258. *
  259. XTHE=JTHE.WELCHE(1,1,1,1,IPJT,ILJT)
  260. XGAM=JGAM.WELCHE(1,1,1,1,IPJG,ILJG)
  261. XCO1=XTHE*0.5D0
  262. IXP1=NINT(IREF*XGAM)
  263. XCO2=(1.D0-XTHE)*((SQRT(DBLE(IREF)))**IXP1)
  264. *!!!
  265. *!!! XCO2=1.D0
  266. IXP2=NINT(1-XGAM)
  267. IF (LFIRST) THEN
  268. WRITE(IOIMP,*) 'XTHE=',XTHE
  269. WRITE(IOIMP,*) 'XGAM=',XGAM
  270. WRITE(IOIMP,*) 'XCO1=',XCO1
  271. WRITE(IOIMP,*) 'XCO2=',XCO2
  272. WRITE(IOIMP,*) 'IXP1=',IXP1,'/2'
  273. WRITE(IOIMP,*) 'IXP2=',IXP2
  274. LFIRST=.FALSE.
  275. ENDIF
  276. *
  277. * Copie des coefficients de la métrique
  278. *
  279. DO IIDIM=1,IDIM
  280. JMET=MET.COEF(IIDIM,IIDIM)
  281. M.XMAT(IIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  282. ENDDO
  283. DO IIDIM=1,IIDIM
  284. NJ=IDIM-IIDIM
  285. IF (NJ.GE.1) THEN
  286. DO JIDIM=IIDIM+1,IDIM
  287. JMET=MET.COEF(IIDIM,JIDIM)
  288. M.XMAT(IIDIM,JIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  289. M.XMAT(JIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM)
  290. ENDDO
  291. ENDIF
  292. ENDDO
  293. * SEGPRT,M
  294. * Calcul de l'inverse
  295. *
  296. LBID=.FALSE.
  297. CALL GEOLI2(IREL,1,1,M.XMAT,MM1.XMAT,DETM,LBID,IMPR,IRET)
  298. IF (IRET.NE.0) GOTO 9999
  299. XDM=1.D0/(SQRT(DETM(1,1)))
  300. *1 XDM=1.D0
  301. ****** XDM=SQRT(DETM(1,1))
  302. *
  303. * Copie du jacobien
  304. *
  305. CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF,
  306. $ 'COPIE ',
  307. $ JAC.XMAT,IREL,IREF,
  308. $ IMPR,IRET)
  309. IF (IRET.NE.0) GOTO 9999
  310. * SEGPRT,JAC
  311. *
  312. * Calcul de la métrique G
  313. *
  314. * Calcul de Jt
  315. CALL MAMA(JAC.XMAT,IREL,IREF,
  316. $ 'TRANSPOS',JT.XMAT,IREF,IREL,
  317. $ IMPR,IRET)
  318. IF (IRET.NE.0) GOTO 9999
  319. * Calcul de MJ
  320. CALL MAMAMA(M.XMAT,IREL,IREL,JAC.XMAT,IREL,IREF,
  321. $ 'FOIS ',MJ.XMAT,IREL,IREF,
  322. $ IMPR,IRET)
  323. IF (IRET.NE.0) GOTO 9999
  324. * Calcul de JTM
  325. CALL MAMA(MJ.XMAT,IREL,IREF,
  326. $ 'TRANSPOS',JTM.XMAT,IREF,IREL,
  327. $ IMPR,IRET)
  328. IF (IRET.NE.0) GOTO 9999
  329. * Calcul de G=JtMJ
  330. CALL MAMAMA(JT.XMAT,IREF,IREL,MJ.XMAT,IREL,IREF,
  331. $ 'FOIS ',G.XMAT,IREF,IREF,
  332. $ IMPR,IRET)
  333. IF (IRET.NE.0) GOTO 9999
  334. * Calcul de l'inverse, du déterminant et trace de l'inverse de g
  335. LBID=.FALSE.
  336. CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET)
  337. IF (IRET.NE.0) GOTO 9999
  338. XM=SQRT(DETG(1,1))
  339. *! WRITE(IOIMP,*) 'XM=',XM
  340. * Calcul de hg-1 et de sa trace
  341. CALL MAMAMA(H.XMAT,IREF,IREF,IG.XMAT,IREF,IREF,
  342. $ 'FOIS ',HIG.XMAT,IREF,IREF,IMPR,IRET)
  343. IF (IRET.NE.0) GOTO 9999
  344. CALL MARE(HIG.XMAT,IREF,IREF,'TRACE ',
  345. $ XL,IMPR,IRET)
  346. IF (IRET.NE.0) GOTO 9999
  347. * Calcul de A = (J+)t = MJg-1 (pseudo-inverse)
  348. CALL MAMAMA(MJ.XMAT,IREL,IREF,IG.XMAT,IREF,IREF,
  349. $ 'FOIS ',A.XMAT,IREL,IREF,IMPR,IRET)
  350. IF (IRET.NE.0) GOTO 9999
  351. * Calcul de B = (J+)t hg-1 = MJg-1hg-1
  352. CALL MAMAMA(A.XMAT,IREL,IREF,HIG.XMAT,IREF,IREF,
  353. $ 'FOIS ',B.XMAT,IREL,IREF,IMPR,IRET)
  354. IF (IRET.NE.0) GOTO 9999
  355. * Calcul de C = g-1hg-1
  356. CALL MAMAMA(IG.XMAT,IREF,IREF,HIG.XMAT,IREF,IREF,
  357. $ 'FOIS ',C.XMAT,IREF,IREF,IMPR,IRET)
  358. IF (IRET.NE.0) GOTO 9999
  359. * Calcul de D = MJ g-1hg-1 JtM = B JtM
  360. CALL MAMAMA(B.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
  361. $ 'FOIS ',D.XMAT,IREL,IREL,IMPR,IRET)
  362. IF (IRET.NE.0) GOTO 9999
  363. * Calcul de F = MJ g-1 JtM = A JtM
  364. CALL MAMAMA(A.XMAT,IREL,IREF,JTM.XMAT,IREF,IREL,
  365. $ 'FOIS ',F.XMAT,IREL,IREL,IMPR,IRET)
  366. IF (IRET.NE.0) GOTO 9999
  367. *
  368. * Calcul de la fonctionnelle
  369. *
  370. IF (NOMLOI.EQ.'AHUF ') THEN
  371. XL1=SQRT(XL)**IXP1
  372. FONC1=XL1*XM
  373. *
  374. XM1=XM**IXP2
  375. FONC2=XM1
  376. *
  377. CONTRI=(XCO1*FONC1)+(XCO2*FONC2*(XH**(IXP2-1)))
  378. *!!!!
  379. CONTRI=CONTRI*XDM
  380. *
  381. * Calcul du résidu
  382. *
  383. ELSEIF (NOMLOI(1:4).EQ.'AHUR') THEN
  384. CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  385. IF (IRET.NE.0) GOTO 9999
  386. CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  387. IF (IRET.NE.0) GOTO 9999
  388. K=IDIM1
  389. L=IDIM2
  390. XL1=SQRT(XL)**IXP1
  391. XL2=IXP1*(SQRT(XL)**(IXP1-2))
  392. AKL1=XL1*XM*A.XMAT(K,L)
  393. AKL2=-1.D0*XL2*XM*B.XMAT(K,L)
  394. AKL=AKL1+AKL2
  395. *
  396. XM2=IXP2*(XM**IXP2)
  397. CKL=XM2*A.XMAT(K,L)
  398. *
  399. CONTRI=(XCO1*AKL)+(XCO2*CKL*(XH**(IXP2-1)))
  400. *!!!!
  401. CONTRI=CONTRI*XDM
  402. *
  403. * Calcul du jacobien
  404. *
  405. ELSEIF (NOMLOI(1:4).EQ.'AHUJ') THEN
  406. CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET)
  407. IF (IRET.NE.0) GOTO 9999
  408. CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET)
  409. IF (IRET.NE.0) GOTO 9999
  410. CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET)
  411. IF (IRET.NE.0) GOTO 9999
  412. CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET)
  413. IF (IRET.NE.0) GOTO 9999
  414. I=IDIM1
  415. J=IDIM2
  416. K=IDIM3
  417. L=IDIM4
  418. XL1=SQRT(XL)**IXP1
  419. XL2=IXP1*(SQRT(XL)**(IXP1-2))
  420. XL3=IXP1*(IXP1-2)*(SQRT(XL)**(IXP1-4))
  421. B11=XM*XL1*A.XMAT(I,J)*A.XMAT(K,L)
  422. B12=-1.D0*XL2*XM*B.XMAT(I,J)*A.XMAT(K,L)
  423. B13=XM*XL1*M.XMAT(K,I)*IG.XMAT(J,L)
  424. B14=-1.D0*XM*XL1*F.XMAT(K,I)*IG.XMAT(J,L)
  425. B15=-1.D0*XM*XL1*A.XMAT(K,J)*A.XMAT(I,L)
  426. BIJKL1=B11+B12+B13+B14+B15
  427. B20=+1.D0*XL3*XM*B.XMAT(I,J)*B.XMAT(K,L)
  428. B21=-1.D0*XL2*XM*A.XMAT(I,J)*B.XMAT(K,L)
  429. B22=-1.D0*XL2*XM*M.XMAT(K,I)*C.XMAT(J,L)
  430. B23=+1.D0*XL2*XM*F.XMAT(K,I)*C.XMAT(J,L)
  431. B24=+1.D0*XL2*XM*A.XMAT(K,J)*B.XMAT(I,L)
  432. B25=+1.D0*XL2*XM*D.XMAT(K,I)*IG.XMAT(J,L)
  433. B26=+1.D0*XL2*XM*B.XMAT(K,J)*A.XMAT(I,L)
  434. BIJKL2=B20+B21+B22+B23+B24+B25+B26
  435. BIJKL=BIJKL1+BIJKL2
  436. *
  437. XM2=IXP2*(XM**IXP2)
  438. XM3=IXP2*IXP2*(XM**IXP2)
  439. D1=XM3*A.XMAT(I,J)*A.XMAT(K,L)
  440. D2=XM2*M.XMAT(K,I)*IG.XMAT(J,L)
  441. D3=-1.D0*XM2*F.XMAT(K,I)*IG.XMAT(J,L)
  442. D4=-1.D0*XM2*A.XMAT(K,J)*A.XMAT(I,L)
  443. DIJKL=D1+D2+D3+D4
  444. *
  445. CONTRI=(XCO1*BIJKL)+(XCO2*DIJKL*(XH**(IXP2-1)))
  446. *!!!!
  447. CONTRI=CONTRI*XDM
  448. *! WRITE(IOIMP,*) 'CONTRI=',CONTRI
  449. ELSE
  450. WRITE(IOIMP,*) 'Erreur grave'
  451. GOTO 9999
  452. ENDIF
  453. FC.WELCHE(1,1,1,1,IPFC,ILFC)=
  454. $ FC.WELCHE(1,1,1,1,IPFC,ILFC)+
  455. $ CONTRI
  456. ENDDO
  457. ENDDO
  458. SEGSUP,F
  459. SEGSUP,D
  460. SEGSUP,MM1
  461. SEGSUP,M
  462. SEGSUP,C
  463. SEGSUP,HIG
  464. SEGSUP,HM1
  465. SEGSUP,H
  466. SEGSUP,IG
  467. SEGSUP,G
  468. SEGSUP,JTM
  469. SEGSUP,JT
  470. SEGSUP,B
  471. SEGSUP,A
  472. SEGSUP,MJ
  473. SEGSUP,JAC
  474. SEGSUP,MET
  475. *
  476. * Normal termination
  477. *
  478. IRET=0
  479. RETURN
  480. *
  481. * Format handling
  482. *
  483. *
  484. * Error handling
  485. *
  486. 9999 CONTINUE
  487. IRET=1
  488. WRITE(IOIMP,*) 'An error was detected in subroutine ccgahu'
  489. RETURN
  490. *
  491. * End of subroutine CCGAHU
  492. *
  493. END
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  

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