C YFROT SOURCE GOUNAND 25/11/12 21:15:48 12399 SUBROUTINE YFROT IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C CE SP DISCRETISE LE TERME DE PERTE DE CHARGE DANS LES EQUATIONS DE C NAVIER STOKES C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 C L OPERATEUR EST DIAGONAL C C SYNTAXE : C C FROT K B INCO UN : C C COMMENTAIRES : C -------------- C C UN CHAMPS DE VITESSE TRANSPORTANT C VO CHAMPS DE VITESSE DE REFFERENCE C K (VECT ) Coefficient C (VECT CENTRE) C B (VECT ) exposant C (VECT CENTRE) C C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMMATRIK -INC SMCHAML -INC SMLENTI POINTEUR IZIPAD.MLENTI -INC SMELEME POINTEUR MELEM1.MELEME -INC SMCHPOI POINTEUR MCHPIN.MCHPOI POINTEUR IZG1.MCHPOI POINTEUR IZGG.MPOVAL,IZGG1.MPOVAL,IZTU1.MPOVAL POINTEUR IZGI.MCHPOI,IZGGI.MPOVAL POINTEUR IZGE.MCHPOI,IZGGE.MPOVAL POINTEUR MZK.MPOVAL,MZBETA.MPOVAL,MZV0.MPOVAL -INC SMLMOTS POINTEUR LINCO.MLMOTS CHARACTER*8 TYPE,MARG,TYPC CHARACTER*(LOCOMP) NOMP(3),NOMI,NOM,NOMZ DIMENSION IXV(3) LOGICAL LOGI PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB) c SAVE IPAS DATA LTAB/'KIZX '/,IPAS/0/ C***************************************************************************** CFROT C write(6,*)' DEBUT FROT' segact mcoord CALL LITABS(LTAB,KTAB,NTB,1,IRET) IF(IRET.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' On attend un ensemble de table soustypes' RETURN ENDIF MTABX=KTAB(1) CALL LEKTAB(MTABX,'EQEX',MTAB1) IF(MTAB1.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' On ne trouve pas l''indice EQEX ? ' RETURN ENDIF CALL LEKTAB(MTAB1,'PASDETPS',KINC) CALL ACME(KINC,'NUPASDT',nupasdt) C write(6,*)' FROT nupasdt=',nupasdt CALL LEKTAB(MTAB1,'INCO',KINC) IF(KINC.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' Il n''y a pas de table INCO ? ?.' RETURN ENDIF C***************************************************************************** C OPTIONS KIMPL=0 KFORM=0 IAXI=0 IF(IFOMOD.EQ.0)IAXI=2 KOPTI=0 TYPE=' ' CALL ACMO(MTABX,'KOPT',TYPE,IENT) IF(TYPE.EQ.'TABLE')KOPTI=IENT IF(KOPTI.NE.0)THEN TYPE=' ' CALL ACMO(KOPTI,'KIMPL',TYPE,IENT) IF(TYPE.EQ.'ENTIER')CALL ACME(KOPTI,'KIMPL',KIMPL) TYPE=' ' CALL ACMO(KOPTI,'KFORM',TYPE,IENT) IF(TYPE.EQ.'ENTIER')CALL ACME(KOPTI,'KFORM',KFORM) IF(KFORM.NE.0.AND.KIMPL.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' Seule la formulation EFM1 est autorisée' RETURN ENDIF ENDIF C***************************************************************************** CALL ACMM(MTABX,'NOMZONE',NOMZ) CALL LEKTAB(MTABX,'DOMZ',MTABZ) IF(MTABZ.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' On ne trouve pas l''indice DOMZ ? ' GO TO 90 ENDIF CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) IF(MELEME.EQ.0)GO TO 90 CALL LEKTAB(MTABZ,'XXPSOML',MCHELM) IF(MCHELM.EQ.0)GO TO 90 SEGACT MCHELM C*** TYPE='LISTMOTS' CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) SEGACT LINCO C************************************************************************* C! CALL LEKTAB(MTAB1,'DOMAINE',MTABD) C! IF(MTABD.EQ.0)THEN C! WRITE(6,*)' Operateur FROT ' C! WRITE(6,*)' On ne trouve pas l''indice DOMAINE ?' C! GO TO 90 C! ENDIF C! CALL LEKTAB(MTABD,'SOMMET',MELEM1) C! IF(MELEM1.EQ.0)THEN C! WRITE(6,*)' Operateur FROT ' C! WRITE(6,*)' On ne trouve pas l''indice SOMMET ?' C! GO TO 90 C! ENDIF CALL LEKTAB(MTAB1,'INCO',INCO) IF(INCO.EQ.0)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)'Il n''y a pas de table INCO ' RETURN ENDIF CALL LEKTAB(MTAB1,'KIZG',KIZG) IF(KIZG.EQ.0)THEN CALL CRTABL(KIZG) CALL ECMM(KIZG,'SOUSTYPE','KIZG') CALL ECMO(MTAB1,'KIZG','TABLE ',KIZG) ENDIF IF(KIMPL.EQ.1)THEN CALL LEKTAB(MTAB1,'KIZG1',KIZG1) IF(KIZG1.EQ.0)THEN CALL CRTABL(KIZG1) CALL ECMM(KIZG1,'SOUSTYPE','KIZG1') CALL ECMO(MTAB1,'KIZG1','TABLE ',KIZG1) ENDIF ENDIF C VERIFICATIONS SUR LES INCONNUES NBINC=LINCO.MOTS(/2) IF(NBINC.NE.1)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 1' RETURN ENDIF NOMI=LINCO.MOTS(1) TYPE=' ' CALL ACMO(INCO,NOMI,TYPE,MCHPOI) IF(TYPE.NE.'CHPOINT ')THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table' RETURN ELSE CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1) MCHPIN=MCHPOI NINKO = IZTU1.VPOCHA(/2) NPTI = IZTU1.VPOCHA(/1) IF (NINKO.NE.IDIM) THEN C Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes MOTERR( 1: 8) = 'INC '//NOMI MOTERR( 9:16) = 'CHPOINT ' CALL ERREUR(784) RETURN ENDIF ENDIF C************************************************************************* C Le domaine de definition est donne par le SPG de la premiere inconnue C Les inconnues suivantes devront posseder ce meme pointeur C On verifie que les points de la zone sont tous inclus dans ce SPG CALL KRIPAD(MELEM1,IZIPAD) IF(IPAS.EQ.0)THEN CALL VERPAD(IZIPAD,MELEME,IRET) IF(IRET.NE.0)THEN WRITE(6,*)' Opérateur FROT ' WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine' GO TO 90 ENDIF ENDIF C***************************************************************************** C************************************************************************* C write(6,*)' FROT KFORM=',KFORM IF(KFORM.EQ.1)THEN IHV=1 CALL YFRT('FROT ',MTABX,IHV,MCHPIN,NOMI,MATRIK,MCHPO1) CALL ECROBJ('MATRIK',MATRIK) CALL ECROBJ('CHPOINT',MCHPO1) c CALL PRLIST RETURN ENDIF C************************************************************************* C Lecture du coefficient C Type du coefficient : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur C write(6,*)' Operateur FROT lecture des coefficients' CALL ACME(MTABX,'IARG',IARG) IF(IARG.NE.2.AND.IARG.NE.3)THEN WRITE(6,*)' Operateur FROT ' WRITE(6,*)' Nombre d''arguments incorrect : ',IARG WRITE(6,*)' On attend 2 ou 3 ' RETURN ENDIF TYPE=' ' CALL ACMO(MTABZ,'CENTRE',TYPE,MELEMC) CALL ACMO(MTABZ,'SOMMET',TYPE,MELEMS) IXV(1)=-MELEMC IXV(2)=0 IXV(3)=1 CALL LEKCOF('Opérateur FROT :', & MTABX,KINC,1,IXV,MK,MZK,NPK,NC1,IKK,IRET) IF(IRET.EQ.0)RETURN IF(IKK.EQ.2)IKK=1 IXV(1)=-MELEMC IXV(2)=0 IXV(3)=1 CALL LEKCOF('Opérateur FROT :', & MTABX,KINC,2,IXV,MBETA,MZBETA,NBETA,NC2,IKB,IRET) IF(IRET.EQ.0)RETURN IF(IKB.EQ.2)IKB=1 IF(IARG.EQ.3)THEN IXV(1)=-MELEMS IXV(2)=0 IXV(3)=1 CALL LEKCOF('Opérateur FROT :', & MTABX,KINC,3,IXV,IZTG3,MZV0,NV0,NC3,IKV,IRET) IF(IRET.EQ.0)RETURN IF(IKV.EQ.2)IKV=1 ELSE Nu=3 WRITE(MARG,FMT='(A4,I1)')'ARGS',Nu CALL LEKTAB(MTABX,MARG,MCHPOI) CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM0) VPOCHA(1,1)=0.D0 VPOCHA(1,2)=0.D0 IF(IDIM.EQ.3)VPOCHA(1,3)=0.D0 MZV0=MPOVAL NV0=1 IKV=1 ENDIF C Fin lecture Arguments ************************************************ TYPE=' ' CALL ACMO(KIZG,NOMI,TYPE,IZG) IF(TYPE.NE.'CHPOINT ')THEN NC=IZTU1.VPOCHA(/2) TYPE='SOMMET' CALL CRCHPT(TYPE,MELEM1,NC,2,IZG) CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG) ENDIF CALL LICHT(IZG,IZGG,TYPC,IGEOM) IF(IGEOM.NE.MELEM1)THEN WRITE(6,*)' Opérateur FROT' WRITE(6,*)' Incompatibilité de SPG entre 1ères inconnues' RETURN ENDIF IF(KIMPL.EQ.1)THEN TYPE=' ' CALL ACMO(KIZG1,NOMI,TYPE,IZG1) IF(TYPE.NE.'CHPOINT ')THEN NC=IZTU1.VPOCHA(/2) TYPE='SOMMET' CALL CRCHPT(TYPE,MELEM1,NC,2,IZG1) CALL ECMO(KIZG1,NOMI,'CHPOINT ',IZG1) ENDIF CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) IF(IGEOM.NE.MELEM1)THEN WRITE(6,*)' Opérateur FROT' WRITE(6,*)' Incompatibilité de SPG entre 1ères inconnues' RETURN ENDIF ELSE IZGG1=IZGG ENDIF NPT=IZGG1.VPOCHA(/1) SEGACT MELEME NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 NUTOEL=0 DT=1.E30 WRITE(NOMP(1),FMT='(I1)')1 NOMP(1)=NOMP(1)(1:1)//NOMI(1:LOCOMP-1) WRITE(NOMP(2),FMT='(I1)')2 NOMP(2)=NOMP(2)(1:1)//NOMI(1:LOCOMP-1) WRITE(NOMP(3),FMT='(I1)')3 NOMP(3)=NOMP(3)(1:1)//NOMI(1:LOCOMP-1) CALL CRCHPK(TYPE,MELEM1,IDIM,IZGI,NOMP) CALL LICHT(IZGI,IZGGI,TYPC,IGEOM) CALL CRCHPK(TYPE,MELEM1,IDIM,IZGE,NOMP) CALL LICHT(IZGE,IZGGE,TYPC,IGEOM) c nbno=IZGG1.VPOCHA(/1) c nbnc=20 c write(6,*)' ==============================================' c write(6,*)' ==============================================' c write(6,*)' IZZGIIII UX avant YFRTI' c write(6,1002)(IZGGi.VPOCHA(ii,1),ii=1,nbnc) c write(6,*)' IZZGIIII UY avant YFRTI' c write(6,1002)(IZGGi.VPOCHA(ii,2),ii=1,nbnc) c c write(6,*)' IZZGEEEE UX avant YFRTI' c write(6,1002)(IZGGe.VPOCHA(ii,1),ii=1,nbnc) c write(6,*)' IZZGEEEE UY avant YFRTI' c write(6,1002)(IZGGe.VPOCHA(ii,2),ii=1,nbnc) DO 1 L=1,NBSOUS IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(L) SEGACT IPT1 MCHAML=ICHAML(L) SEGACT MCHAML MELVAL=IELVAL(1) SEGACT MELVAL IF(IMACHE(L).NE.IPT1)THEN write(*,*)'IPT1,IMACHE ',IPT1,IMACHE(L) goto 90 ENDIF NP =IPT1.NUM(/1) NBEL=IPT1.NUM(/2) IES=IDIM CALL YFRTI(VELCHE,IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,IAXI, & IZIPAD.LECT,KIMPL, & MZK.VPOCHA,NPK,IKK, & MZBETA.VPOCHA,NBETA,IKB, & MZV0.VPOCHA,NV0,IKV, & IZTU1.VPOCHA,IZGGE.VPOCHA,IZGGI.VPOCHA) SEGDES IPT1 NUTOEL=NUTOEL+NBEL 1 CONTINUE TYPE = ' ' Cµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµ Cµµ On construit un MATRIK vide si KIMPL=0 Cµµ On construit un MATRIK diagonal si KIMPL=1 C- Construction du chapeau de l'objet MATRIK IF(KIMPL.EQ.0)THEN c write(6,*)' KIMPL=',kimpl NRIGE = 7 NKID = 9 NKMT = 7 NMATRI = 0 SEGINI MATRIK SEGDES MATRIK ELSE c write(6,*)' KIMPL=',kimpl NRIGE = 7 NKID = 9 NKMT = 7 NMATRI = 1 SEGINI MATRIK IRIGEL(1,1) = MELEM1 IRIGEL(2,1) = MELEM1 C KFORM = 0 EFM1 C KFORM = 2 VF C -> matrice Diagonale IRIGEL(7,1) = 5 NBME = NINKO SEGACT MELEM1 NBSOUS = MELEM1.LISOUS(/1) IF (NBSOUS.EQ.0) NBSOUS=1 SEGINI IMATRI IRIGEL(4,1) = IMATRI KSPGP = MELEM1 KSPGD = MELEM1 DO 10 I=1,NBME WRITE(NOM,FMT='(I1)')I NOM=NOM(1:1)//NOMI(1:LOCOMP-1) LISPRI(I) = NOM LISDUA(I) = NOM NP=1 MP=1 NBEL=MELEM1.NUM(/2) SEGINI IZAFM do 13 k=1,nbel AM(K,1,1)=IZGGI.VPOCHA(K,I) 13 continue LIZAFM(1,I)=IZAFM 10 CONTINUE ENDIF c write(6,*)' ..............................................' c write(6,*)' IZZGIII UX au retour de YFRTI' c write(6,1002)(IZGGi.VPOCHA(ii,1),ii=1,nbnc) c write(6,*)' IZZGIII UY au retour de YFRTI' c write(6,1002)(IZGGi.VPOCHA(ii,2),ii=1,nbnc) c write(6,*)' IZZGEEE UX au retour de YFRTI' c write(6,1002)(IZGGe.VPOCHA(ii,1),ii=1,nbnc) c write(6,*)' IZZGEEE UY au retour de YFRTI' c write(6,1002)(IZGGe.VPOCHA(ii,2),ii=1,nbnc) c write(6,*)' ==============================================' c write(6,*)' ==============================================' c do 141 n=1,2 c do 141 i=1,nbno c IZGGE.VPOCHA(i,n)=IZGGE.VPOCHA(i,n)*(-1) c141 continue Cµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµµ SEGDES IZTU1 SEGDES LINCO SEGSUP IZIPAD CALL ECROBJ('MATRIK',MATRIK) CALL ECROBJ('CHPOINT',IZGE) c IPAS=1 RETURN 90 CONTINUE WRITE(6,*)' Interuption anormale de FROT' RETURN 1002 FORMAT(10(1X,1PE11.4)) END