C CCGTSU SOURCE GOUNAND 26/01/09 21:15:13 12441 SUBROUTINE CCGTSU(LCOF,NOMLOI, $ FC, $ IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : CCGTSU C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : C Tension de surface C 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, 17/01/07, version initiale C HISTORIQUE : v1, 17/01/07, création C HISTORIQUE : 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 JPC.MCHEVA POINTEUR JTS.MCHEVA CHARACTER*8 NOMLOI INTEGER ICOF * -INC TMXMAT POINTEUR JAC.MXMAT,JT.MXMAT POINTEUR G.MXMAT POINTEUR IG.MXMAT POINTEUR JPT.MXMAT POINTEUR JPTJT.MXMAT * INTEGER LAXSP REAL*8 DEUPI,QUATPI,XR REAL*8 DETG(1) PARAMETER (DEUPI=2.D0*XPI,QUATPI=4.D0*XPI) * INTEGER IMPR,IRET * * Executable statements * IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgtsu' 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 du coefficent de tension de surface * ICOF=ICOF+1 JTS =LCOF.LISCHE(ICOF) NLTS=JTS.WELCHE(/6) NPTS=JTS.WELCHE(/5) * 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 * *! WRITE(IOIMP,*) 'IFOMOD=',IFOMOD IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN *! WRITE(IOIMP,*) 'MODE AXISYMETRIQUE !!!!!!' LAXSP=1 ELSEIF (IFOMOD.EQ.5) THEN *! WRITE(IOIMP,*) 'MODE SPHERIQUE !!!!!!' LAXSP=2 ELSE *! WRITE(IOIMP,*) 'MODE PLAN !!!!!!' LAXSP=0 ENDIF *! WRITE(IOIMP,*) 'NOMLOI=',NOMLOI *! WRITE(IOIMP,*) 'LAXSP=',LAXSP * IF (LAXSP.GT.0) THEN ICOF=ICOF+1 ICOF=ICOF+1 JPC=LCOF.LISCHE(ICOF) NLPC=JPC.WELCHE(/6) NPPC=JPC.WELCHE(/5) IF (((NLPC.NE.1).AND.(NLPC.NE.NLFC)).OR. $ ((NPPC.NE.1).AND.(NPPC.NE.NPFC))) THEN WRITE(IOIMP,*) 'Erreur dims JPC' GOTO 9999 ENDIF ELSE JPC=0 NLPC=0 NPPC=0 ENDIF * * Objets temporaires * LDIM1=IREL LDIM2=IREF SEGINI,JAC SEGINI,JPT LDIM1=IREF LDIM2=IREL SEGINI,JT LDIM1=IREF LDIM2=IREF SEGINI,G SEGINI,IG LDIM1=IREL LDIM2=IREL SEGINI,JPTJT * DO ILFC=1,NLFC IF (NLTS.EQ.1) THEN ILTS=1 ELSE ILTS=ILFC ENDIF IF (NLJA.EQ.1) THEN ILJA=1 ELSE ILJA=ILFC ENDIF IF (NLPC.EQ.1) THEN ILPC=1 ELSE ILPC=ILFC ENDIF DO IPFC=1,NPFC IF (NPTS.EQ.1) THEN IPTS=1 ELSE IPTS=IPFC ENDIF IF (NPJA.EQ.1) THEN IPJA=1 ELSE IPJA=IPFC ENDIF IF (NPPC.EQ.1) THEN IPPC=1 ELSE IPPC=IPFC ENDIF * * Valeur du coefficient * IF (JTS.NE.0) THEN XTS=JTS.WELCHE(1,1,1,1,IPTS,ILTS) ENDIF * * Première coordonnée * IF (JPC.NE.0) THEN XR=JPC.WELCHE(1,1,1,1,IPPC,ILPC) ENDIF * * Copie du jacobien * CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF, $ 'COPIE ', $ JAC.XMAT,IREL,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 C SEGPRT,JAC * * Calcul de la métrique G * * Calcul de Jt CALL MAMA(JAC.XMAT,IREL,IREF, $ 'TRANSPOS',JT.XMAT,IREF,IREL, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de G=JtJ CALL MAMAMA(JT.XMAT,IREF,IREL,JAC.XMAT,IREL,IREF, $ '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 * CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,IMPR,IRET) IF (IRET.NE.0) THEN WRITE(IOIMP,*) $ 'JtJ non inversible' GOTO 9999 ENDIF SDETG=SQRT(ABS(DETG(1))) * Calcul de JPT=J (JtJ)-1 CALL MAMAMA(JAC.XMAT,IREL,IREF,IG.XMAT,IREF,IREF, $ 'FOIS ',JPT.XMAT,IREL,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * Calcul de JPT=J (JtJ)-1 CALL MAMAMA(JPT.XMAT,IREL,IREF,JT.XMAT,IREF,IREL, $ 'FOIS ',JPTJT.XMAT,IREL,IREL, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * * Calcul de la fonctionnelle * IF (NOMLOI.EQ.'TSUF ') THEN CONTRI=XTS*SDETG IF (LAXSP.EQ.1) THEN CONTRI=CONTRI*2.D0*XPI*XR ELSEIF (LAXSP.EQ.2) THEN CONTRI=CONTRI*4.D0*XPI*XR*XR ENDIF * * Calcul du résidu * ELSEIF (NOMLOI(1:4).EQ.'TSUR') THEN CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN IF (LAXSP.EQ.1) THEN CONTRI=XTS*2.D0*XPI*SDETG ELSEIF (LAXSP.EQ.2) THEN CONTRI=XTS*8.D0*XPI*XR*SDETG ELSE WRITE(IOIMP,*) 'LAXSP=',LAXSP GOTO 9999 ENDIF ELSE CONTRI=XTS*SDETG*JPT.XMAT(IDIM1,IDIM2) IF (LAXSP.EQ.1) THEN CONTRI=CONTRI*DEUPI*XR ELSEIF (LAXSP.EQ.2) THEN CONTRI=CONTRI*QUATPI*XR*XR ENDIF ENDIF * * Calcul du jacobien * ELSEIF (NOMLOI(1:4).EQ.'TSUJ') THEN CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0) $ .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN IF (LAXSP.EQ.2) THEN CONTRI=XTS*8.D0*XPI*SDETG ELSE WRITE(IOIMP,*) 'LAXSP=',LAXSP WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ENDIF ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN IF (LAXSP.GE.1) THEN CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM3,IDIM4) IF (LAXSP.EQ.2) THEN CONTRI=CONTRI*4.D0*XR ENDIF ELSE WRITE(IOIMP,*) 'LAXSP=',LAXSP WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ENDIF ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN IF (LAXSP.GE.1) THEN CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM1,IDIM2) IF (LAXSP.EQ.2) THEN CONTRI=CONTRI*4.D0*XR ENDIF ELSE WRITE(IOIMP,*) 'LAXSP=',LAXSP WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ENDIF ELSE CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4) IF (IDIM1.EQ.IDIM3) THEN CONTRI=CONTRI+IG.XMAT(IDIM2,IDIM4) ENDIF CONTRI=CONTRI- $ (JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4)) CONTRI=CONTRI- $ (JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4)) CONTRI=CONTRI*XTS*SDETG IF (LAXSP.EQ.1) THEN CONTRI=CONTRI*DEUPI*XR ELSEIF (LAXSP.EQ.2) THEN CONTRI=CONTRI*QUATPI*XR*XR ENDIF ENDIF * * Calcul d'une partie du jacobien * ELSEIF (NOMLOI(1:3).EQ.'TSU') THEN CALL CH2INT(NOMLOI(4:4),IDIM0,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0) $ .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ELSE IF (IDIM0.EQ.1) THEN CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4) ELSEIF (IDIM0.EQ.2) THEN IF (IDIM1.EQ.IDIM3) THEN CONTRI=IG.XMAT(IDIM2,IDIM4) ELSE CONTRI=0.D0 ENDIF ELSEIF (IDIM0.EQ.3) THEN CONTRI= $ -(JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4) $ ) ELSEIF (IDIM0.EQ.4) THEN CONTRI=-(JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4 $ )) ELSE WRITE(IOIMP,*) 'NOMLOI=',NOMLOI GOTO 9999 ENDIF CONTRI=CONTRI*XTS*SDETG IF (LAXSP.EQ.1) THEN CONTRI=CONTRI*DEUPI*XR ELSEIF (LAXSP.EQ.2) THEN CONTRI=CONTRI*QUATPI*XR*XR ENDIF ENDIF *! 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,JPTJT SEGSUP,JPT SEGSUP,IG SEGSUP,G SEGSUP,JT SEGSUP,JAC * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ccgtsu' RETURN * * End of subroutine CCGTSU * END