ccgahu
C CCGAHU SOURCE PV 22/04/22 23:00:16 11344 $ FC, $ IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : CCGAHU C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : C DEDU ADAP avec metrique C On implémente ici la première partie de la fonctionnelle C de Huang C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : C APPELE PAR : C*********************************************************************** C ENTREES : C ENTREES/SORTIES : C SORTIES : - C TRAVAIL : C*********************************************************************** C VERSION : v1, 30/04/07, version initiale C HISTORIQUE : v1, 30/04/07, création C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC CCREEL -INC TNLIN *-INC SMCHAEL INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1 POINTEUR FC.MCHEVA POINTEUR LCOF.LCHEVA POINTEUR JMAJAC.MCHEVA POINTEUR JMIJAC.MCHEVA POINTEUR JDTJAC.MCHEVA POINTEUR JMAREG.MCHEVA POINTEUR JMET.MCHEVA POINTEUR JTHE.MCHEVA POINTEUR JGAM.MCHEVA CHARACTER*8 NOMLOI INTEGER ICOF * -INC TMXMAT * Objets temporaires POINTEUR JAC.MXMAT,JT.MXMAT POINTEUR G.MXMAT,IG.MXMAT,H.MXMAT,HIG.MXMAT,HM1.MXMAT POINTEUR MM1.MXMAT,JTM.MXMAT,MJ.MXMAT * Objets à préconditionner (éventuellement) POINTEUR A.MXMAT,B.MXMAT,C.MXMAT,D.MXMAT,F.MXMAT,M.MXMAT * SEGMENT MCOF POINTEUR COEF(IDIM,IDIM).MCHEVA ENDSEGMENT POINTEUR MET.MCOF * LOGICAL LBID,LFIRST INTEGER LAXSP REAL*8 DEUPI,XR REAL*8 XL,XM,XH real*8 deth(1,1),detg(1,1),detm(1,1) * INTEGER IMPR,IRET * * Executable statements * IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgahu' C IF (.NOT.(IDIM.EQ.1)) THEN C WRITE(IOIMP,*) 'IDIM=',IDIM,' ?' C GOTO 9999 C ENDIF NLFC=FC.WELCHE(/6) NPFC=FC.WELCHE(/5) ICOF=0 * * Récupération des coefficients de la metrique * SEGINI MET DO IIDIM=1,IDIM ICOF=ICOF+1 JMET=LCOF.LISCHE(ICOF) IF (ICOF.EQ.1) THEN NLJM=JMET.WELCHE(/6) NPJM=JMET.WELCHE(/5) ELSE NLJM2=JMET.WELCHE(/6) NPJM2=JMET.WELCHE(/5) IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN WRITE(IOIMP,*) 'Erreur grave dims JMET' GOTO 9999 ENDIF ENDIF MET.COEF(IIDIM,IIDIM)=JMET ENDDO DO IIDIM=1,IDIM NJ=IDIM-IIDIM IF (NJ.GE.1) THEN DO JIDIM=IIDIM+1,IDIM ICOF=ICOF+1 JMET=LCOF.LISCHE(ICOF) NLJM2=JMET.WELCHE(/6) NPJM2=JMET.WELCHE(/5) IF (NLJM2.NE.NLJM.OR.NPJM2.NE.NPJM) THEN WRITE(IOIMP,*) 'Erreur grave dims JMET2' GOTO 9999 ENDIF MET.COEF(IIDIM,JIDIM)=JMET ENDDO ENDIF ENDDO * ICOF=ICOF+1 JTHE=LCOF.LISCHE(ICOF) NLJT=JTHE.WELCHE(/6) NPJT=JTHE.WELCHE(/5) I1 =JTHE.WELCHE(/4) $ THEN WRITE(IOIMP,*) 'Erreur dims JTHE' GOTO 9999 ENDIF * ICOF=ICOF+1 JGAM=LCOF.LISCHE(ICOF) NLJG=JTHE.WELCHE(/6) NPJG=JTHE.WELCHE(/5) I1 =JTHE.WELCHE(/4) $ THEN WRITE(IOIMP,*) 'Erreur dims JGAM' GOTO 9999 ENDIF * ICOF=ICOF+1 JMAJAC=LCOF.LISCHE(ICOF) NLJA=JMAJAC.WELCHE(/6) NPJA=JMAJAC.WELCHE(/5) IREF=JMAJAC.WELCHE(/4) IREL=JMAJAC.WELCHE(/3) * IF (IREL.NE.IDIM) THEN WRITE(IOIMP,*) 'Erreur dims JMAJAC' GOTO 9999 ENDIF * ICOF=ICOF+1 ICOF=ICOF+1 ICOF=ICOF+1 JMAREG=LCOF.LISCHE(ICOF) NLJR=JMAREG.WELCHE(/6) NPJR=JMAREG.WELCHE(/5) I1 =JMAREG.WELCHE(/4) $ THEN WRITE(IOIMP,*) 'Erreur dims JMAREG' GOTO 9999 ENDIF * * Objets temporaires et à préconditionner * LDIM1=IREL LDIM2=IREF SEGINI,JAC SEGINI,MJ SEGINI,A SEGINI,B LDIM1=IREF LDIM2=IREL SEGINI,JT SEGINI,JTM LDIM1=IREF LDIM2=IREF SEGINI,G SEGINI,IG SEGINI,H SEGINI,HM1 SEGINI,HIG SEGINI,C LDIM1=IREL LDIM2=IREL SEGINI,M SEGINI,MM1 SEGINI,D SEGINI,F * * Calcul de la métrique des éléments réguliers * $ 'JTJ ',H.XMAT,IREF,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * SEGPRT,H LBID=.FALSE. IF (IRET.NE.0) GOTO 9999 *!!! XH=1.0D0 XH=1.D0/(SQRT(DETH(1,1))) *!! XH=1.D0/(DETH(1,1)) * LFIRST = .FALSE. DO ILFC=1,NLFC IF (NLJM.EQ.1) THEN ILJM=1 ELSE ILJM=ILFC ENDIF IF (NLJA.EQ.1) THEN ILJA=1 ELSE ILJA=ILFC ENDIF IF (NLJT.EQ.1) THEN ILJT=1 ELSE ILJT=ILFC ENDIF IF (NLJG.EQ.1) THEN ILJG=1 ELSE ILJG=ILFC ENDIF DO IPFC=1,NPFC IF (NPJM.EQ.1) THEN IPJM=1 ELSE IPJM=IPFC ENDIF IF (NPJA.EQ.1) THEN IPJA=1 ELSE IPJA=IPFC ENDIF IF (NPJT.EQ.1) THEN IPJT=1 ELSE IPJT=IPFC ENDIF IF (NPJG.EQ.1) THEN IPJG=1 ELSE IPJG=IPFC ENDIF * * Récupération de l'exposant multiplié par 2 * On suppose que c'est un entier * XTHE=JTHE.WELCHE(1,1,1,1,IPJT,ILJT) XGAM=JGAM.WELCHE(1,1,1,1,IPJG,ILJG) XCO1=XTHE*0.5D0 IXP1=NINT(IREF*XGAM) XCO2=(1.D0-XTHE)*((SQRT(DBLE(IREF)))**IXP1) *!!! *!!! XCO2=1.D0 IXP2=NINT(1-XGAM) IF (LFIRST) THEN WRITE(IOIMP,*) 'XTHE=',XTHE WRITE(IOIMP,*) 'XGAM=',XGAM WRITE(IOIMP,*) 'XCO1=',XCO1 WRITE(IOIMP,*) 'XCO2=',XCO2 WRITE(IOIMP,*) 'IXP1=',IXP1,'/2' WRITE(IOIMP,*) 'IXP2=',IXP2 LFIRST=.FALSE. ENDIF * * Copie des coefficients de la métrique * DO IIDIM=1,IDIM JMET=MET.COEF(IIDIM,IIDIM) M.XMAT(IIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM) ENDDO DO IIDIM=1,IIDIM NJ=IDIM-IIDIM IF (NJ.GE.1) THEN DO JIDIM=IIDIM+1,IDIM JMET=MET.COEF(IIDIM,JIDIM) M.XMAT(IIDIM,JIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM) M.XMAT(JIDIM,IIDIM)=JMET.WELCHE(1,1,1,1,IPJM,ILJM) ENDDO ENDIF ENDDO * SEGPRT,M * Calcul de l'inverse * LBID=.FALSE. IF (IRET.NE.0) GOTO 9999 XDM=1.D0/(SQRT(DETM(1,1))) *1 XDM=1.D0 ****** XDM=SQRT(DETM(1,1)) * * Copie du jacobien * $ 'COPIE ', $ JAC.XMAT,IREL,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * SEGPRT,JAC * * Calcul de la métrique G * * Calcul de Jt $ 'TRANSPOS',JT.XMAT,IREF,IREL, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de MJ $ 'FOIS ',MJ.XMAT,IREL,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de JTM $ 'TRANSPOS',JTM.XMAT,IREF,IREL, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de G=JtMJ $ 'FOIS ',G.XMAT,IREF,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de l'inverse, du déterminant et trace de l'inverse de g LBID=.FALSE. IF (IRET.NE.0) GOTO 9999 XM=SQRT(DETG(1,1)) *! WRITE(IOIMP,*) 'XM=',XM * Calcul de hg-1 et de sa trace $ 'FOIS ',HIG.XMAT,IREF,IREF,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 $ XL,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de A = (J+)t = MJg-1 (pseudo-inverse) $ 'FOIS ',A.XMAT,IREL,IREF,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de B = (J+)t hg-1 = MJg-1hg-1 $ 'FOIS ',B.XMAT,IREL,IREF,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de C = g-1hg-1 $ 'FOIS ',C.XMAT,IREF,IREF,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de D = MJ g-1hg-1 JtM = B JtM $ 'FOIS ',D.XMAT,IREL,IREL,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de F = MJ g-1 JtM = A JtM $ 'FOIS ',F.XMAT,IREL,IREL,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * * Calcul de la fonctionnelle * IF (NOMLOI.EQ.'AHUF ') THEN XL1=SQRT(XL)**IXP1 FONC1=XL1*XM * XM1=XM**IXP2 FONC2=XM1 * CONTRI=(XCO1*FONC1)+(XCO2*FONC2*(XH**(IXP2-1))) *!!!! CONTRI=CONTRI*XDM * * Calcul du résidu * ELSEIF (NOMLOI(1:4).EQ.'AHUR') THEN IF (IRET.NE.0) GOTO 9999 IF (IRET.NE.0) GOTO 9999 K=IDIM1 L=IDIM2 XL1=SQRT(XL)**IXP1 XL2=IXP1*(SQRT(XL)**(IXP1-2)) AKL1=XL1*XM*A.XMAT(K,L) AKL2=-1.D0*XL2*XM*B.XMAT(K,L) AKL=AKL1+AKL2 * XM2=IXP2*(XM**IXP2) CKL=XM2*A.XMAT(K,L) * CONTRI=(XCO1*AKL)+(XCO2*CKL*(XH**(IXP2-1))) *!!!! CONTRI=CONTRI*XDM * * Calcul du jacobien * ELSEIF (NOMLOI(1:4).EQ.'AHUJ') THEN IF (IRET.NE.0) GOTO 9999 IF (IRET.NE.0) GOTO 9999 IF (IRET.NE.0) GOTO 9999 IF (IRET.NE.0) GOTO 9999 I=IDIM1 J=IDIM2 K=IDIM3 L=IDIM4 XL1=SQRT(XL)**IXP1 XL2=IXP1*(SQRT(XL)**(IXP1-2)) XL3=IXP1*(IXP1-2)*(SQRT(XL)**(IXP1-4)) B11=XM*XL1*A.XMAT(I,J)*A.XMAT(K,L) B12=-1.D0*XL2*XM*B.XMAT(I,J)*A.XMAT(K,L) B13=XM*XL1*M.XMAT(K,I)*IG.XMAT(J,L) B14=-1.D0*XM*XL1*F.XMAT(K,I)*IG.XMAT(J,L) B15=-1.D0*XM*XL1*A.XMAT(K,J)*A.XMAT(I,L) BIJKL1=B11+B12+B13+B14+B15 B20=+1.D0*XL3*XM*B.XMAT(I,J)*B.XMAT(K,L) B21=-1.D0*XL2*XM*A.XMAT(I,J)*B.XMAT(K,L) B22=-1.D0*XL2*XM*M.XMAT(K,I)*C.XMAT(J,L) B23=+1.D0*XL2*XM*F.XMAT(K,I)*C.XMAT(J,L) B24=+1.D0*XL2*XM*A.XMAT(K,J)*B.XMAT(I,L) B25=+1.D0*XL2*XM*D.XMAT(K,I)*IG.XMAT(J,L) B26=+1.D0*XL2*XM*B.XMAT(K,J)*A.XMAT(I,L) BIJKL2=B20+B21+B22+B23+B24+B25+B26 BIJKL=BIJKL1+BIJKL2 * XM2=IXP2*(XM**IXP2) XM3=IXP2*IXP2*(XM**IXP2) D1=XM3*A.XMAT(I,J)*A.XMAT(K,L) D2=XM2*M.XMAT(K,I)*IG.XMAT(J,L) D3=-1.D0*XM2*F.XMAT(K,I)*IG.XMAT(J,L) D4=-1.D0*XM2*A.XMAT(K,J)*A.XMAT(I,L) DIJKL=D1+D2+D3+D4 * CONTRI=(XCO1*BIJKL)+(XCO2*DIJKL*(XH**(IXP2-1))) *!!!! CONTRI=CONTRI*XDM *! WRITE(IOIMP,*) 'CONTRI=',CONTRI ELSE WRITE(IOIMP,*) 'Erreur grave' GOTO 9999 ENDIF FC.WELCHE(1,1,1,1,IPFC,ILFC)= $ FC.WELCHE(1,1,1,1,IPFC,ILFC)+ $ CONTRI ENDDO ENDDO SEGSUP,F SEGSUP,D SEGSUP,MM1 SEGSUP,M SEGSUP,C SEGSUP,HIG SEGSUP,HM1 SEGSUP,H SEGSUP,IG SEGSUP,G SEGSUP,JTM SEGSUP,JT SEGSUP,B SEGSUP,A SEGSUP,MJ SEGSUP,JAC SEGSUP,MET * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ccgahu' RETURN * * End of subroutine CCGAHU * END
© Cast3M 2003 - Tous droits réservés.
Mentions légales