C YFROT SOURCE PV 09/03/12 21:37:49 6325 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 CCOPTIO -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 MZK.MPOVAL,MZBETA.MPOVAL,MZV0.MPOVAL -INC SMLMOTS POINTEUR LINCO.MLMOTS CHARACTER*8 NOMZ,NOMI,TYPE,MARG,TYPC DIMENSION IXV(3) LOGICAL LOGI PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB) SAVE IPAS DATA LTAB/'KIZX '/,IPAS/0/ C***************************************************************************** CFROT C write(6,*)' DEBUT FROT' 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,'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************************************************************************* IF(KFORM.EQ.1)THEN IHV=1 CALL YFRT('FROT ',MTABX,IHV,MCHPIN,NOMI,MATRIK,MCHPO1) CALL ECROBJ('MATRIK',MATRIK) CALL ECROBJ('CHPOINT',MCHPO1) 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,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,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 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,IZGG.VPOCHA,IZGG1.VPOCHA) SEGDES IPT1 NUTOEL=NUTOEL+NBEL 1 CONTINUE SEGDES IZTU1 SEGDES LINCO SEGSUP IZIPAD IPAS=1 RETURN 90 CONTINUE WRITE(6,*)' Interuption anormale de FROT' RETURN 1002 FORMAT(10(1X,1PE11.4)) END