C ZZFIMP SOURCE MAGN 10/05/18 21:18:09 6675 SUBROUTINE ZZFIMP(MTABX,MTAB1) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C SYNTAXE : C C FIMP coef C C EN 2D C elements SEG2 -> Flux C elements TRI3 -> Source volumique C elements QUA4 -> Source volumique C EN 3D C elements SEG2 -> Pas de sens !! C elements TRI3 -> Flux C elements QUA4 -> Flux C elements CUB8 -> Source volumique C elements PRI6 -> Source volumique C elements TET4 -> Source volumique C C C MTAB1 : Table type EQEX -> RV C MTABZ : Table type DOMAINE -> Zone definition opérateur C MTABD : Table type DOMAINE -> Zone Totale apres assemblage C MTABX : Table type KIZX -> Description opérateur C C les cartes correspondantes sont commentées C¤. C Gare à l'explicite !!! C C C C C*********************************************************************** -INC CCOPTIO -INC CCGEOME -INC SMCHAML -INC SMCOORD -INC SMLENTI -INC SMELEME POINTEUR MELEM1.MELEME,SPGID.MELEME,SPGZ.MELEME POINTEUR MELEMP.MELEME -INC SMCHPOI POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL POINTEUR IZTU1.MPOVAL POINTEUR MZFLU.MPOVAL POINTEUR IZVOL.MPOVAL -INC SIZFFB POINTEUR IZF1.IZFFM -INC SMLMOTS POINTEUR LINCO.MLMOTS CHARACTER*8 NOMZ,NOMI,NOMA,TYPE,NOM0,TYPC,MTERR LOGICAL LOGI PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB),IXV(4) SAVE IPAS DATA LTAB/'KIZX '/,IPAS/0/ C***************************************************************************** CFIMP C? write(6,*)' Debut FIMP' C C- Récupération de la table INCO (pointeur KINC) C CALL LEKTAB(MTAB1,'INCO',KINC) IF(KINC.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' INCO ' MOTERR( 9:16) = ' INCO ' MOTERR(17:24) = ' EQEX ' CALL ERREUR(786) RETURN ENDIF C***************************************************************************** C OPTIONS C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC C IDCEN = 0->rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG C KPOIN = 0->SOMMET 1-> FACE 2-> CENTRE 3-> CENTREP0 4-> CENTREP1 IAXI=0 IF(IFOMOD.EQ.0)IAXI=2 C C- Récupération de la table des options KOPT (pointeur KOPTI) C CALL LEKTAB(MTABX,'KOPT',KOPTI) IF (KOPTI.EQ.0) THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' KOPT ' MOTERR( 9:16) = ' KOPT ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL ACME(KOPTI,'MTRMASS ',MMPG) IPG=0 IF(MMPG.EQ.3)IPG=1 CALL ACME(KOPTI,'IDCEN',IDCEN) KDIM=1 IF(IDCEN.EQ.2)KDIM=IDIM CALL ACME(KOPTI,'IKOMP',IKOMP) CALL ACME(KOPTI,'KIMPL',KIMPL) CALL ACME(KOPTI,'KFORM',KFORM) CALL ACME(KOPTI,'KPOIND',KPOIND) C write(6,*)' INCOD=',KPOIND IF(KFORM.GE.2)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = 'EF/EFM1 ' CALL ERREUR(803) RETURN ENDIF CALL ACME(KOPTI,'KMACO',KMACO) CALL ACMF(KOPTI,'AIMPL',AIMPL) IF (IERR.NE.0) RETURN C write(6,*)' Apres les options ' C***************************************************************************** C C- Récupération de la table DOMAINE associée au domaine local C CALL ACMM(MTABX,'NOMZONE',NOMZ) CALL LEKTAB(MTABX,'DOMZ',MTABZ) IF(MTABZ.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' DOMZ ' MOTERR( 9:16) = ' DOMZ ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF C C- Récupération des indices CENTRE, FACE, SOMMET et MAILLAGE C- de la table DOMAINE C C? CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) CALL LEKTAB(MTABZ,'MAILLAGE',MELEMZ) CALL LEKTAB(MTABZ,'SOMMET',MELEMS) CALL LEKTAB(MTABZ,'CENTRE',MELEMC) CALL LEKTAB(MTABZ,'XXPSOML',MCHELM) CALL LEKTAB(MTABZ,'XXVOLUM',MCHVOL) CALL LEKTAB(MTABZ,'MACRO',MACRO) C? IF(CALL LEKTAB(MTABZ,'MACRO1',MACRO1) CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) IF (IERR.NE.0) RETURN CALL LICHT(MCHVOL,IZVOL,TYPC,IGEOM) SEGACT MCHELM C C- Vérification des compatiblités Formulation/SPG C- Identification du spg de l'inconnue C- SPGZ=spg inco duale de la zone; MELEME=connectivité associée ; C C write(6,*)' KPOIND,KFORM=',KPOIND,KFORM IF(KFORM.EQ.0)THEN IF(KPOIND.EQ.99)THEN KPOIND=0 SPGZ =MELEMS C MELEME=MELEMS MELEME=MELEMZ ELSEIF(KPOIND.EQ.2)THEN SPGZ =MELEMC MELEME=MELEMC MELEMP=MELEMC ELSEIF(KPOIND.EQ.0)THEN SPGZ =MELEMS C MELEME=MELEMS MELEME=MELEMZ ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = ' EF ' CALL ERREUR(803) RETURN ENDIF ELSEIF(KFORM.EQ.1)THEN C write(6,*)' KFORM=',kform,' KPOIND=',kpoind IF(KPOIND.EQ.99)THEN KPOIND=0 SPGZ =MELEMS MELEME=MELEMZ ELSEIF(KPOIND.EQ.0)THEN SPGZ =MELEMS MELEME=MELEMZ ELSEIF(KPOIND.EQ.2)THEN SPGZ =MELEMC MELEME=MELEMZ MELEMP=MELEMC ELSEIF(KPOIND.EQ.3)THEN MTERR='EF CTRP0' IF(MACRO.EQ.0.AND.MQUAD.EQ.0)GO TO 90 CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) SPGZ =MELEMC MELEME=MELEMZ IF(MACRO.NE.0)CALL LEKTAB(MTABZ,'MACRO1',MELEME) MELEMP=MELEMC ELSEIF(KPOIND.EQ.4)THEN MTERR='EF CTRP1' IF(MACRO.EQ.0.AND.MQUAD.EQ.0)GO TO 90 CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) SPGZ =MELEMC MELEME=MELEMZ IF(MACRO.NE.0)CALL LEKTAB(MTABZ,'MACRO1',MELEME) CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP) ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0.AND.KPOIND.NE.3 & .AND.KPOIND.NE.4)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = ' EF ' CALL ERREUR(803) RETURN ENDIF ELSEIF(KFORM.EQ.2)THEN IF(KPOIND.EQ.99)THEN KPOIND=2 SPGZ =MELEMC MELEME=MELEMC MELEMP=MELEMC ELSEIF(KPOIND.EQ.2)THEN SPGZ =MELEMC MELEME=MELEMC MELEMP=MELEMC ELSEIF(KPOIND.NE.2)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = ' VF ' CALL ERREUR(803) RETURN ENDIF ELSEIF(KFORM.EQ.3)THEN MTERR='EFMC' GO TO 90 IF(KPOIND.EQ.99.OR.KPOIND.EQ.2)THEN KPOIND=2 SPGZ =MELEMC MELEME=MELEMC ELSEIF(KPOIND.EQ.3.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) SPGZ =MELEMC MELEME=MELEMC ELSEIF(KPOIND.EQ.4.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) CALL LEKTAB(MTABZ,'ELTP1NC',MELEMQ) SPGZ =MELEMC MELEME=MELEMQ IF(MACRO.NE.0)MELEMO=MACRO1 IF(MQUAD.NE.0)MELEMO=MELEMZ C ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.3.AND.KPOIND.NE.4)THEN ELSE C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = ' EFMC ' CALL ERREUR(803) RETURN ENDIF ENDIF C************************************************************************* C VERIFICATIONS SUR LES INCONNUES C C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI C TYPE='LISTMOTS' CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) IF (IERR.NE.0) RETURN SEGACT LINCO NBINC=LINCO.MOTS(/2) IF(NBINC.NE.1)THEN C Indice %m1:8 : contient plus de %i1 %m9:16 MOTERR( 1:8) = 'LISTINCO' INTERR(1) = 1 MOTERR(9:16) = ' MOTS ' CALL ERREUR(799) RETURN ENDIF NOMI=LINCO.MOTS(1) IF(KIMPL.EQ.0)THEN 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 ENDIF C --> 1 ere Inconnue NOMI=LINCO.MOTS(1) TYPE=' ' CALL ACMO(KINC,NOMI,TYPE,MCHPOI) IF(TYPE.NE.'CHPOINT ')THEN C Indice %m1:8 : ne contient pas un objet de type %m9:16 MOTERR( 1: 8) = 'INC '//NOMI MOTERR( 9:16) = 'CHPOINT ' CALL ERREUR(800) RETURN ELSE CALL LICHT(MCHPOI,IZTU1,TYPC,SPGID) 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(SPGID,MLENTI) IF(IPAS.EQ.0)THEN CALL VERPAD(MLENTI,SPGZ,IRET) IF(IRET.NE.0)THEN C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique MOTERR(1: 8) = 'INC '//NOMI MOTERR(9:16) = 'CHPOINT ' CALL ERREUR(788) IPAS=0 RETURN ENDIF ENDIF C************************************************************************* C Lecture du coefficient C Type du coefficient : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur C write(6,*)' Lecture des coefficients ' CALL ACME(MTABX,'IARG',IARG) IF(IARG.NE.1)THEN C Indice %m1:8 : nombre d'arguments incorrect MOTERR(1:8) = 'IARG ' CALL ERREUR(804) RETURN ENDIF IRET=3 IXV(1)=MELEMC IXV(2)=1 IXV(3)=0 C? IXV(4)=MELEMS C write(6,*)' MELEMQ=',melemq CALL LEKCOF('Opérateur FIMP :', & MTABX,KINC,1,IXV,MFLU,MZFLU,NPT1,NC1,IK1,IRET) IF(IRET.EQ.0)RETURN C Fin lecture Arguments ************************************************** IF(KIMPL.EQ.0)THEN TYPE=' ' CALL ACMO(KIZG,NOMI,TYPE,IZG1) IF(TYPE.NE.'CHPOINT ')THEN NC=IZTU1.VPOCHA(/2) TYPE='SOMMET' CALL CRCHPT(TYPE,SPGID,NC,IZG1) C SEGACT IZG1 C MSOUPO=IZG1.IPCHP(1) C SEGACT MSOUPO C NOCOMP(1)=NOMI(1:4) C SEGDES MSOUPO CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG1) ENDIF ELSE NC=IZTU1.VPOCHA(/2) TYPE='SOMMET' C? pour plutard CALL CRCHPT(TYPE,SPGZ,NC,IZG1) CALL CRCHPT(TYPE,SPGID,NC,IZG1) SEGACT IZG1 MSOUPO=IZG1.IPCHP(1) SEGACT MSOUPO*MOD NOCOMP(1)=NOMI(1:4) ENDIF CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) IF(IGEOM.NE.SPGID)THEN C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique MOTERR(1: 8) = 'INC '//NOMI MOTERR(9:16) = 'CHPOINT ' CALL ERREUR(788) RETURN ENDIF SEGACT MELEME NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 NUTOEL=0 NPTD=IZTU1.VPOCHA(/1) IES=IDIM 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 NP =IPT1.NUM(/1) NBEL=IPT1.NUM(/2) IF(KPOIND.EQ.0)THEN CALL ZXFIMP(NBEL,NUTOEL,NP,LECT,IPT1.NUM, & VELCHE,IZGG1.VPOCHA,MZFLU.VPOCHA,IK1) ELSE IPT2=MELEMP IF(NBSOUS.NE.1)IPT2=LISOUS(L) SEGACT IPT2 IF(MQUAD.NE.0)THEN IF(KPOIND.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' IF(KPOIND.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' IF(KPOIND.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1' ELSEIF(MACRO.NE.0)THEN IF(KPOIND.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' IF(KPOIND.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0' IF(KPOIND.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1' ELSE IF(KPOIND.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' ENDIF CALL KALPBG(NOM0,'FONFORM ',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=GR(/1) NPG=GR(/3) IZF1=KTP(1) SEGACT IZF1*MOD MP1=IZF1.FN(/1) NP = IPT1.NUM(/1) NK=NUTOEL DO 21 K=1,NBEL NK=NK+1 JC=(1-IK1)*(NK-1)+1 DO 109 I=1,NP J=IPT1.NUM(I,K) DO 109 N=1,IDIM XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 109 CONTINUE CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE) DO 39 M=1,MP1 M1=LECT(IPT2.NUM(M,K)) U=0.D0 DO 33 LL=1,NPG U=U+IZF1.FN(M,LL)*MZFLU.VPOCHA(JC,1)*PGSQ(LL) 33 CONTINUE IZGG1.VPOCHA(M1,1)=IZGG1.VPOCHA(M1,1)-U 39 CONTINUE 21 CONTINUE ENDIF SEGDES IPT1,MCHAML,MELVAL NUTOEL=NUTOEL+NBEL 1 CONTINUE C SEGDES MZFLU,IZVOL SEGDES MZFLU IF(KIMPL.NE.0)THEN TYPE=' ' CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2) IF(TYPE.NE.'CHPOINT')THEN CALL ECMO(MTAB1,'SMBR','CHPOINT',IZG1) ELSE CALL ECROBJ('CHPOINT',MCHPO2) CALL ECROBJ('CHPOINT',IZG1) CALL PRFUSE CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI) ENDIF ENDIF SEGDES MELEME SEGDES IZTU1 SEGDES IZG1,IZGG1 SEGDES LINCO SEGSUP MLENTI IPAS=1 RETURN 90 CONTINUE C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = MTERR CALL ERREUR(803) RETURN 1002 FORMAT(10(1X,1PE11.4)) END