C CCGTDI SOURCE GOUNAND 26/01/09 21:15:12 12441 SUBROUTINE CCGTDI(LCOF, $ FC, $ IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : CCGTDI C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : C Diamètre interne d'un élément suivant une direction donnée 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, 13/09/06, version initiale C HISTORIQUE : v1, 13/09/06, 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 JDIAMA.MCHEVA POINTEUR JPC.MCHEVA POINTEUR JVIT.MCHEVA CHARACTER*8 NOMLOI INTEGER ICOF * -INC TMXMAT POINTEUR A.MXMAT POINTEUR AP.MXMAT POINTEUR JMA.MXMAT POINTEUR KJMA.MXMAT POINTEUR JM1.MXMAT,J.MXMAT POINTEUR K.MXMAT * SEGMENT MVIT POINTEUR MVCOMP(IDIM).MCHEVA ENDSEGMENT * LOGICAL LREGP,LRELP INTEGER IMPR,IRET * * Executable statements * IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgtdi' * XPETI=SQRT(XPETIT) NLFC=FC.WELCHE(/6) NPFC=FC.WELCHE(/5) ICOF=0 * * Récupération des coefficients de la metrique * SEGINI MVIT DO IIDIM=1,IDIM ICOF=ICOF+1 JVIT=LCOF.LISCHE(ICOF) IF (ICOF.EQ.1) THEN NLJV=JVIT.WELCHE(/6) NPJV=JVIT.WELCHE(/5) ELSE NLJV2=JVIT.WELCHE(/6) NPJV2=JVIT.WELCHE(/5) IF (NLJV2.NE.NLJV.OR.NPJV2.NE.NPJV) THEN WRITE(IOIMP,*) 'Erreur grave dims JVIT' GOTO 9999 ENDIF ENDIF MVIT.MVCOMP(IIDIM)=JVIT ENDDO * ICOF=ICOF+1 JMAJAC=LCOF.LISCHE(ICOF) C NLJA=JMAJAC.WELCHE(/6) C 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 JMIJAC=LCOF.LISCHE(ICOF) IF (JMIJAC.LE.0) THEN write(ioimp,*) 'Jacobien non inversible' GOTO 9999 ENDIF NLJI=JMIJAC.WELCHE(/6) NPJI=JMIJAC.WELCHE(/5) IREL2=JMIJAC.WELCHE(/4) IREF2=JMIJAC.WELCHE(/3) * IF (IREL2.NE.IREL.OR.IREF2.NE.IREF) THEN WRITE(IOIMP,*) 'Erreur dims JMIJAC' GOTO 9999 ENDIF * ICOF=ICOF+1 ICOF=ICOF+1 JMAREG=LCOF.LISCHE(ICOF) NLJR=JMAREG.WELCHE(/6) NPJR=JMAREG.WELCHE(/5) I1 =JMAREG.WELCHE(/4) I2 =JMAREG.WELCHE(/3) IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF)) $ THEN WRITE(IOIMP,*) 'Erreur dims JMAREG' GOTO 9999 ENDIF ICOF=ICOF+1 JDIAMA=LCOF.LISCHE(ICOF) NLJD=JDIAMA.WELCHE(/6) NPJD=JDIAMA.WELCHE(/5) I1 =JDIAMA.WELCHE(/4) I2 =JDIAMA.WELCHE(/3) IF ((NLJD.NE.1).OR.(NPJD.NE.1).OR.(I1.NE.1).OR.(I2.NE.1)) $ THEN WRITE(IOIMP,*) 'Erreur dims JDIAMA' GOTO 9999 ENDIF XDIAMA=JDIAMA.WELCHE(1,1,1,1,1,1) * * Objets temporaires * LDIM1=IREL LDIM2=1 SEGINI,A SEGINI,AP LDIM1=IREF LDIM2=1 SEGINI,JMA LDIM1=IREF LDIM2=1 SEGINI,KJMA LDIM1=IREL LDIM2=IREF SEGINI,J LDIM1=IREF LDIM2=IREL SEGINI,JM1 LDIM1=IREF LDIM2=IREF SEGINI,K * Copie de la matrice jacobienne ref -> regulier CALL MAMA(JMAREG.WELCHE,IREF,IREF, $ 'COPIE ',K.XMAT,IREF,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * SEGPRT,H DO ILFC=1,NLFC IF (NLJV.EQ.1) THEN ILJV=1 ELSE ILJV=ILFC ENDIF C IF (NLJA.EQ.1) THEN C ILJA=1 C ELSE C ILJA=ILFC C ENDIF IF (NLJI.EQ.1) THEN ILJI=1 ELSE ILJI=ILFC ENDIF * DO IPFC=1,NPFC IF (NPJV.EQ.1) THEN IPJV=1 ELSE IPJV=IPFC ENDIF C IF (NPJA.EQ.1) THEN C IPJA=1 C ELSE C IPJA=IPFC C ENDIF IF (NPJI.EQ.1) THEN IPJI=1 ELSE IPJI=IPFC ENDIF * * Copie de la vitesse * DO IIDIM=1,IDIM JVIT=MVIT.MVCOMP(IIDIM) A.XMAT(IIDIM,1)=JVIT.WELCHE(1,1,1,1,IPJV,ILJV) ENDDO * SEGPRT,A * * Copie de l'inverse (ou pseudo-inverse) du jacobien ref->reel * CALL MAMA(JMIJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL, $ 'COPIE ', $ JM1.XMAT,IREF,IREL, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 C SEGPRT,JM1 * * Copie du jacobien ref->reel * CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREL,IREF, $ 'COPIE ', $ J.XMAT,IREL,IREF, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 C SEGPRT,JM1 * * Vecteur A dans le repère de référence * CALL MAMAMA(JM1.XMAT,IREF,IREL,A.XMAT,IREL,1, $ 'FOIS ',JMA.XMAT,IREF,1, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * * Vecteur A' dans le repère réel (projection sur le domaine considéré) * CALL MAMAMA(J.XMAT,IREL,IREF,JMA.XMAT,IREF,1, $ 'FOIS ',AP.XMAT,IREL,1, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * * Vecteur A dans le repère de l'élément régulier * CALL MAMAMA(K.XMAT,IREF,IREF,JMA.XMAT,IREF,1, $ 'FOIS ',KJMA.XMAT,IREF,1, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 * * Normes de A * XAREG2=XZERO DO IIREF=1,IREF XAREG2=XAREG2+(KJMA.XMAT(IIREF,1)**2) ENDDO XA2=XZERO DO IIDIM=1,IDIM XA2=XA2+(A.XMAT(IIDIM,1)**2) ENDDO XAP2=XZERO DO IIDIM=1,IDIM XAP2=XAP2+(AP.XMAT(IIDIM,1)**2) ENDDO XAPA=XZERO DO IIDIM=1,IDIM XAPA=XAPA+(AP.XMAT(IIDIM,1)*A.XMAT(IIDIM,1)) ENDDO C WRITE(IOIMP,*) 'XA2=',XA2 C WRITE(IOIMP,*) 'XAP2=',XAP2 C WRITE(IOIMP,*) 'XAREG2=',XAREG2 C WRITE(IOIMP,*) 'XDIAMA=',XDIAMA * IF (XA2.LT.XPETIT) THEN WRITE(IOIMP,*) 'The given direction is 0' GOTO 9999 ENDIF IF (XAREG2.LT.XPETIT) THEN CONTRI=0.D0 ELSE CONTRI=SQRT((XAP2*XAPA)/(XAREG2*XA2))*XDIAMA * CONTRI=SQRT(XA2/XAREG2)*XDIAMA ENDIF FC.WELCHE(1,1,1,1,IPFC,ILFC)= $ FC.WELCHE(1,1,1,1,IPFC,ILFC)+ $ CONTRI ENDDO ENDDO SEGSUP,K SEGSUP,J SEGSUP,JM1 SEGSUP,KJMA SEGSUP,JMA SEGSUP,AP SEGSUP,A SEGSUP,MVIT * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ccgtdi' RETURN * * End of subroutine CCGTDI * END