C YDFDT SOURCE FANDEUR 22/01/03 21:15:57 11136 SUBROUTINE YDFDT C--------------------------------------------------------------------- C Cet opérateur discrétise le terme de dérivée temporelle. C--------------------------------------------------------------------- C Les formulations numériques sont EF Galerkin, EF Petrov-Galerkin et C VF, la matrice masse pouvant etre condensée en EF. C Les schema en temps disponibles CRANK NICOLSON (DT**2) CNG (DT**4) C Les schema en temps disponibles BDF2 (DT**2) BDF4 (DT**4) C--------------------------------------------------------------------- C C------------------------- C Phrase d'appel GIBIANE : C------------------------- C C 'DFDT' TAB1 ; C C TAB1 : Table de sous-type 'KIZX' contenant les données pour DFDT C C-------------------------------- C Construction de TAB1 via EQEX : C-------------------------------- C C 'ZONE' TAB2 'OPER' 'DFDT' CHP1 CHP2 FLOT1 (CHP3) (CHP4) 'INCO' MOT1 C C TAB2 : TABLE DOMAINE associée à la zone d'action de l'opérateur C CHP1 : Coefficient multiplicateur de la dérivée temporelle C (densité dans le cas de la qdm ou rho*cp en thermique) C (FLOTTANT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) C CHP2 : Inconnue au pas de temps précédant C (MOT ou CHPO SCAL ou CHPO VECT) - (spg : CENTRE ou SOMMET) C FLOT1 : Valeur du pas de temps (MOT ou REEL ou ENTIER) C CHP3 : Champ de vitesse - Optionnel (utilisé en Petrov-Galerkin) C (POINT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) C CHP4 : Diffusion - Optionnel (utilisé en Petrov-Galerkin) C (FLOTTANT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) C MOT1 : Nom de l'inconnue sur laquelle agit l'opérateur (MOT) C C------------ C Résultats : C------------ C C Les matrices élémentaires sont stockés dans un objet MATRIK qui C est rangé à l'indice de type MOT MATELM de la table KIZX. C Le second membre est stocké dans un CHPO et assemblé dans la table C EQEX à l'indice de type MOT SMBR. C--------------------------------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C -INC PPARAM -INC CCOPTIO -INC CCREEL *- -INC CCGEOME -INC SIZFFB POINTEUR IZF1.IZFFM PARAMETER (LRV=64) SEGMENT PETROV REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG) REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV) REAL*8 AJK(LRV,IDIM,IDIM,NPG) ENDSEGMENT -INC SMCOORD -INC SMLENTI POINTEUR IPADU.MLENTI,IPADS.MLENTI,IPADI.MLENTI -INC SMELEME POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEMZ.MELEME POINTEUR MELEMO.MELEME -INC SMCHAML -INC SMCHPOI POINTEUR IZRO.MPOVAL,MZDT.MPOVAL,IZTGG2.MPOVAL,IZTU1.MPOVAL POINTEUR IZTGN2.MPOVAL POINTEUR IZTCO.MPOVAL,MZUN.MPOVAL,MZMU.MPOVAL -INC SMMATRIK -INC SMLMOTS C CHARACTER*8 NOMI,TYPE,TYPC,NOM,NOM0,NOMDT,MTYP,CHAI PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB),IXV(4) DATA LTAB/'KIZX '/ C C DFDT C write(6,*)' Opérateur DFDT ' segact mcoord C Nouvelle directive EQUA de EQEX MTYP=' ' CALL QUETYP(MTYP,0,IRET) IF(IRET.EQ.0)THEN C% On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40 MOTERR( 1: 8) = 'CHAI ' MOTERR( 9:16) = 'MMODEL ' MOTERR(17:24) = 'TABLE ' CALL ERREUR(472) RETURN ENDIF IF(MTYP.EQ.'MMODEL')THEN CALL YTCLSF('T ','DFDT ') RETURN ELSEIF(MTYP.EQ.'MOT ')THEN CALL LIRCHA(CHAI,1,IRET) CALL YTCLSF(CHAI,'DFDT ') RETURN ENDIF C Fin Nouvelle directive EQUA de EQEX DTM1 = 0.D0 DTX = 0.D0 MLENTI = 0 IAXI = 0 NPTU = 0 IF (IFOMOD.EQ.0) IAXI=2 DEUPI=1.D0 IF(IAXI.NE.0)DEUPI=2.D0*XPI C C- Lecture de la table KIZX (pointeur MTABX) associée à DFDT C c CALL LITABS(LTAB,KTAB,NTB,1,IRET) c IF (IERR.NE.0) RETURN c MTABX = KTAB(1) C C- Récupération de la table xNS type KIZX (pointeur MTABX) C C write(6,*)'Récupération de la table xNS type KIZX pointeur MTABX' CALL LIROBJ('TABLE',MTABX,1,IRET) IF(IRET.EQ.0)THEN C% On ne trouve pas d'objet de type %m1:8 MOTERR( 1: 8) = 'TABLE' CALL ERREUR(37) RETURN ENDIF C A tout hazard on regarde si une table ne peut en cacher une autre C write(6,*)' A tout hazard on regarde si une table ne peut en', C &' cacher une autre !!!!!' MTYP='TABLE' CALL QUETYP(MTYP,0,IRET) IF(IRET.EQ.1)THEN CALL LIROBJ('TABLE',MTABX,1,IRET) IF(IRET.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' MTABX ' MOTERR( 9:16) = ' MTABX ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF C write(6,*)' EH BEN C EST LE CAS une KIZX // !!!!!!!!!!' C call pplist(MTABX) C C- Récupération de la table DOMAINE associée au domaine local C TYPE=' ' CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) IF(TYPE.EQ.'MMODEL')THEN CALL LEKMOD(MMDZ,MTABZ,INEFMD) TYPE='TABLE' CALL ECMO(MTABX,'TDOMZ',TYPE,MTABZ) ELSEIF(TYPE.EQ.'TABLE')THEN CALL ECMO(MTABX,'TDOMZ',TYPE,MMDZ) ELSE 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 ELSE TYPE=' ' CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) IF(TYPE.EQ.'MMODEL')THEN CALL LEKMOD(MMDZ,MTABZ,INEFMD) TYPE='TABLE' CALL ECMO(MTABX,'TDOMZ',TYPE,MTABZ) ELSEIF(TYPE.EQ.'TABLE')THEN CALL ECMO(MTABX,'TDOMZ',TYPE,MMDZ) ELSE 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 ENDIF C....................................................................... C C- Récupération de la table EQEX (pointeur MTAB1) C CALL LEKTAB(MTABX,'EQEX',MTAB1) IF (MTAB1.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' EQEX ' MOTERR( 9:16) = ' EQEX ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL ACME(MTAB1,'NAVISTOK',NASTOK) IF(NASTOK.EQ.0)THEN CALL ZDFDT(MTABX,MTAB1) RETURN ENDIF C C- Récupération de la table DOMAINE associée au domaine local C MTYP='MMODEL' CALL QUETYP(MTYP,0,IRET) IF(IRET.EQ.1)THEN CALL LIROBJ('MMODEL',MMDZ,0,IRET) ELSE TYPE=' ' CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) IF(TYPE.NE.'MMODEL')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 ENDIF CALL LEKMOD(MMDZ,MTABZ,INEFMD) 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- Récupération de la table PASDETPS (pointeur KINT) C CALL LEKTAB(MTAB1,'PASDETPS',KINT) C?2 IF(KINT.EQ.0)KINT=KINC 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 ELSE IKOMP=0 CALL ACME(KOPTI,'ISCHT',ISCHT) C- Initialisation des constantes en fonction du schema CDT=1.D0 CN=1.D0 CNM=1.D0 CNMM=0.D0 IF(ISCHT.EQ.1)THEN CDT=2.D0 CN=3.D0 CNM=4.D0 CNMM=-1.D0 ELSEIF(ISCHT.EQ.2)THEN CDT=2.D0 CN=6.D0 CNM=0.D0 CNMM=6.D0 ENDIF KDIM=1 CALL ACME(KOPTI,'IDCEN',IDCEN) IF(IDCEN.EQ.2)KDIM=IDIM CALL ACME(KOPTI,'KFORM',KFORM) CALL ACME(KOPTI,'KPOIND',KPOIND) CALL ACME(KOPTI,'KMACO',KMACO) CALL ACMF(KOPTI,'CMD',CMD) ENDIF IDCENV = IDCEN IF (IDCEN.EQ.5.OR.IDCEN.EQ.4) IDCENV=1 C C- Récupération des indices CENTRE, FACE, SOMMET et MAILLAGE C- de la table DOMAINE C CALL LEKTAB(MTABZ,'SOMMET',MELEMS) CALL LEKTAB(MTABZ,'FACE',MELEMF) CALL LEKTAB(MTABZ,'CENTRE',MELEMC) CALL LEKTAB(MTABZ,'MAILLAGE',MELEMZ) CALL LEKTAB(MTABZ,'MACRO1',MELEMI) CALL LEKTAB(MTABZ,'MACRO',MACRO) CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) IF (IERR.NE.0) RETURN C C- Vérification des compatiblités Formulation/SPG C- Identification du spg de l'inconnue C- MELEM1=spg ; MELEME=connectivité associée ; C IF(KFORM.EQ.0)THEN IF(KPOIND.EQ.99)THEN KPOIND=0 MELEM1=MELEMS MELEME=MELEMS ELSEIF(KPOIND.EQ.2)THEN MELEM1=MELEMC MELEME=MELEMC ELSEIF(KPOIND.EQ.0)THEN MELEM1=MELEMS MELEME=MELEMS 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 IF(KPOIND.EQ.99)THEN KPOIND=0 MELEM1=MELEMS MELEME=MELEMZ ELSEIF(KPOIND.EQ.0)THEN MELEM1=MELEMS MELEME=MELEMZ ELSEIF(KPOIND.EQ.2)THEN MELEM1=MELEMC MELEME=MELEMC 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.2)THEN IF(KPOIND.EQ.99)THEN KPOIND=2 MELEM1=MELEMC MELEME=MELEMC ELSEIF(KPOIND.EQ.2)THEN MELEM1=MELEMC MELEME=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 IF(KPOIND.EQ.99.OR.KPOIND.EQ.2)THEN KPOIND=2 MELEM1=MELEMC MELEME=MELEMC ELSEIF(KPOIND.EQ.3.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) MELEM1=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) MELEM1=MELEMC MELEME=MELEMQ IF(MACRO.NE.0)MELEMO=MELEMI 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- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI C TYPE = 'LISTMOTS' CALL ACMO(MTABX,'LISTINCO',TYPE,MLMOTS) IF (IERR.NE.0) RETURN SEGACT MLMOTS NBINC = MOTS(/2) IF (NBINC.LE.0.OR.NBINC.GE.3) THEN C Indice %m1:8 : contient plus de %i1 %m9:16 MOTERR( 1:8) = 'LISTINCO' INTERR(1) = 2 MOTERR(9:16) = ' MOTS ' CALL ERREUR(799) RETURN ENDIF NOMI = MOTS(1) IF (NBINC.EQ.2) THEN IF (MOTS(1).NE.MOTS(2)) THEN C Indice %m1:8 : contient plus de %i1 %m9:16 MOTERR( 1:8) = 'LISTINCO' INTERR(1) = 1 MOTERR(9:16) = ' MOT ' CALL ERREUR(799) RETURN ENDIF ENDIF SEGDES MLMOTS C C- Récupération de l'inconnue C 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 LICHTL(MCHPOI,IZTU1,TYPC,MELEMD) NINKO = IZTU1.VPOCHA(/2) IF (NINKO.NE.1.AND.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************************************************************************* 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(MELEMD,IPADI) IPADS=IPADI IPADU=IPADI IF(MELEMD.NE.MELEMS)CALL KRIPAD(MELEMS,IPADS) MLENTI=IPADI IF(MELEM1.NE.MELEMD)CALL KRIPAD(MELEM1,MLENTI) CALL VERPAD(IPADI,MELEM1,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) RETURN ENDIF C C- Test du nombre d'arguments C ICD=0 IF(ISCHT.NE.0)ICD=1 CALL ACME(MTABX,'IARG',IARG) IF (IDCENV.EQ.1) THEN IF (IARG.NE.(3+ICD).AND.IARG.NE.(5+ICD)) THEN C Indice %m1:8 : nombre d'arguments incorrect MOTERR(1:8) = 'IARG ' CALL ERREUR(804) RETURN ENDIF ELSE IF (IARG.NE.(5+ICD)) THEN C Indice %m1:8 : nombre d'arguments incorrect MOTERR(1:8) = 'IARG ' CALL ERREUR(804) RETURN ENDIF ENDIF C C- Lecture du coefficient multiplicateur SCAL CENTRE ou SCAL SOMMET C IXV(1) = MELEMC IXV(2) = 1 IXV(3) = 0 IXV(4) = MELEMS IRET = 4 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,1,IXV,MRO,IZRO,NPT1,NC1,IK1,IRET) IF (IRET.EQ.0) RETURN C C- Lecture de l'inconnue au pas de temps précédant C IF (NINKO.EQ.1) THEN IXV(1) = MELEMD ELSE IXV(1) =-MELEMD ENDIF IXV(2) = 0 IXV(3) = 0 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,2,IXV,IZTG2,IZTGG2,NPT2,NC2,IK2,IRET) IF (IRET.EQ.0) RETURN C C- Lecture de l'inconnue au pas de temps n-2 (ISCHT NE 0) C IF (NINKO.EQ.1) THEN IXV(1) = MELEMD ELSE IXV(1) =-MELEMD ENDIF IXV(2) = 0 IXV(3) = 0 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,(2+ICD),IXV,IZTG22,IZTGN2,NPT21,NC21, & IK21,IRET) IF (IRET.EQ.0) RETURN C C- Lecture du pas de temps C TYPE=' ' CALL ACMO(MTABX,'ARG3',TYPE,KK) IF(TYPE.EQ.'FLOTTANT')THEN CALL ACMF(MTABX,'ARG3',DTR) IF(KINT.NE.0)CALL ECMF(KINT,'DELTAT',DTR) ELSEIF(TYPE.EQ.'MOT')THEN CALL ACMM(MTABX,'ARG3',NOMDT) IF(NOMDT.NE.'DELTAT')THEN IXV(1) = 0 IXV(2) = 1 IXV(3) = 0 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,(3+ICD),IXV,IZTG3,MZDT,NPT3,NC3,IK3,IRET) IF (IRET.EQ.0) RETURN DTR = MZDT.VPOCHA(1,1) IF(KINT.NE.0)CALL ECMF(KINT,'DELTAT',DTR) ELSE IF(KINT.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = 'PASDETPS' MOTERR( 9:16) = 'PASDETPS' MOTERR(17:24) = ' EQEX ' CALL ERREUR(786) RETURN ENDIF CALL ACMF(KINT,'DELTAT',DTR) CALL ACMF(MTAB1,'ALFA',ALFA) C write(6,*)' DFDT nomdt=',nomdt,' DTR=',dtr,alfa,(DTR*ALFA) DTR=DTR*ALFA ENDIF ENDIF DTR = DTR * CDT C C- Lecture des données pour le décentrement C- Champ de vitesse et viscosité C IKU=1 MZUN=IZTU1 MZMU=IZRO IKMU=1 IF (IARG.EQ.5) THEN IXV(1) =-MELEMS IXV(2) = 0 IXV(3) = 1 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,(4+ICD),IXV,MUN,MZUN,NPTU,NCU,IKU,IRET) IF (IRET.EQ.0) RETURN IPADU=IPADS IF (IKU.EQ.2) IKU=1 C IXV(1) = MELEMC IXV(2) = 1 IXV(3) = 0 IXV(4) = MELEMS IRET = 4 CALL LEKCOF(' Operateur DFDT :', & MTABX,KINC,(5+ICD),IXV,MMU,MZMU,NPTMU,NCMU,IKMU,IRET) IF (IRET.EQ.0) RETURN ENDIF C C- Transformation éventuelle des arguments du SOMMET vers le CENTRE C IF (IK1.EQ.4) THEN CALL COSOCE(MELEMC,MELEMS,MELEMZ,IZRO,MPOVAL) IZRO = MPOVAL IK1 = 0 IK1DETR = 1 ELSE IK1DETR = 0 ENDIF IF (IKMU.EQ.4) THEN CALL COSOCE(MELEMC,MELEMS,MELEMZ,MZMU,MPOVAL) MZMU = MPOVAL IKMU = 0 IKMUDETR = 1 ELSE IKMUDETR = 0 ENDIF C C- Récupération de caractéristiques géométriques dans le cas EF C- ou du volume de chaque élément dans le cas VF C IF (KPOIND.EQ.0) THEN ELSEIF (KPOIND.EQ.2) THEN CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI) IF (IERR.NE.0) RETURN CALL LICHTL(MCHPOI,MPOVAL,TYPC,IGEOM) ELSEIF (KPOIND.EQ.3.AND.MACRO.EQ.0) THEN CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI) IF (IERR.NE.0) RETURN CALL LICHTL(MCHPOI,MPOVAL,TYPC,IGEOM) ELSEIF (KPOIND.EQ.3.AND.MACRO.NE.0) THEN CALL LEKTAB(MTABZ,'VOLUMAC',MCHPOI) IF (IERR.NE.0) RETURN CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM) ENDIF C C- Récupération/Construction du chapeau de l'objet MATRIK C TYPE = ' ' CALL ACMO(MTABX,'MATELM',TYPE,MATRIK) IF (TYPE.EQ.'MATRIK'.AND.KMACO.NE.0) THEN SEGACT MATRIK NMATRI = IRIGEL(/2) MELEME = IRIGEL(1,1) SEGACT MELEME IMATRI = IRIGEL(4,1) SEGACT IMATRI NBME = LIZAFM(/2) NINKO = NBME NBSOUS = LIZAFM(/1) MELEM1 = KSPGP ELSE NRIGE = 7 NKID = 9 NKMT = 7 NMATRI = 1 SEGINI MATRIK IRIGEL(1,1) = MELEME IRIGEL(2,1) = MELEME IF (KFORM.EQ.0.OR.KFORM.EQ.2) THEN C EFM1 ou VF -> Diagonale IRIGEL(7,1) = 5 ELSEIF (KFORM.EQ.1) THEN C EF matrice pleine sym si centree non sinon IF (IDCENV.EQ.1) THEN IRIGEL(7,1) = 0 ELSE IRIGEL(7,1) = 2 ENDIF ELSEIF (KFORM.EQ.3) THEN IF(KPOIND.EQ.2.OR.KPOIND.EQ.3)THEN IRIGEL(7,1) = 5 ELSEIF(KPOIND.EQ.4)THEN IRIGEL(7,1) = 0 ENDIF ENDIF NBME = NINKO SEGACT MELEME NBSOUS = LISOUS(/1) IF (NBSOUS.EQ.0) NBSOUS=1 SEGINI IMATRI IRIGEL(4,1) = IMATRI KSPGP = MELEM1 KSPGD = MELEM1 IF (NBME.EQ.1) THEN LISPRI(1) = NOMI(1:4)//' ' LISDUA(1) = NOMI(1:4)//' ' ELSE DO 10 I=1,NBME WRITE(NOM,FMT='(I1,A7)')I,NOMI(1:7) LISPRI(I) = NOM(1:4)//' ' LISDUA(I) = NOM(1:4)//' ' 10 CONTINUE ENDIF ENDIF C C- Construction des segments IZAFM contenant les matrices élémentaires C IF (KMACO.NE.2) THEN DO 30 L=1,NBSOUS IPT1 = MELEME IF (NBSOUS.NE.1) IPT1=LISOUS(L) SEGACT IPT1 NP = IPT1.NUM(/1) MP = NP NBEL = IPT1.NUM(/2) SEGDES IPT1 DO 20 I=1,NBME SEGINI IZAFM LIZAFM(L,I) = IZAFM 20 CONTINUE 30 CONTINUE ENDIF C C--------------------- C Calcul des matrices C--------------------- C IF (KMACO.NE.2) THEN C C-------------------------------------------------- C-- Cas ou l'inconnue est au SOMMET (EF ou EFM1) -- C-------------------------------------------------- C IF (KPOIND.EQ.0) THEN K0 = 0 IF (KFORM.EQ.0) THEN MELEME = MELEMZ SEGACT MELEME NBSOUS = LISOUS(/1) IF (NBSOUS.EQ.0) NBSOUS=1 ENDIF DO 140 L=1,NBSOUS IPT1 = MELEME IF (NBSOUS.NE.1) IPT1=LISOUS(L) SEGACT IPT1 NP = IPT1.NUM(/1) NBEL = IPT1.NUM(/2) NOM0 = NOMS(IPT1.ITYPEL)//' ' CALL KALPBG(NOM0,'FONFORM ',IZFFM) IF (IZFFM.EQ.0) THEN C Echec lors de la lecture des fonctions de forme d'un élément CALL ERREUR(662) RETURN ENDIF SEGACT IZFFM IZHR = KZHR(1) SEGACT IZHR*MOD NPG = GR(/3) NES = GR(/1) C C- Allocation des tableaux de travail pour PETROV-GALERKIN C SEGINI PETROV C C- Traitement des éléments par paquets de LRV éléments C NNN = MOD(NBEL,LRV) IF (NNN.EQ.0) THEN NPACK = NBEL / LRV ELSE NPACK = 1 + (NBEL-NNN)/LRV ENDIF KPACKD = 1 KPACKF = NPACK DO 130 KPACK=KPACKD,KPACKF KDEB = 1 + (KPACK-1)*LRV KFIN = MIN(NBEL,KDEB+LRV-1) C C- Réduction des coefficients sur les éléments du paquet C DO 40 K=KDEB,KFIN KP = K - KDEB + 1 NK = K + K0 NKR = (1-IK1)*(NK-1) + 1 COEFK(KP) = IZRO.VPOCHA(NKR,1) / DTR 40 CONTINUE IF(IDCENV.NE.1)THEN DO 41 K=KDEB,KFIN KP = K - KDEB + 1 NK = K + K0 NKR = (1-IK1)*(NK-1) + 1 NKM = (1-IKMU)*(NK-1) + 1 ANUK(KP) = MZMU.VPOCHA(NKM,1)/(IZRO.VPOCHA(NKR,1)+XPETIT) 41 CONTINUE ENDIF C C- Récupération de la fonction de projection suivant IDCEN C C CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN DT=0.D0 IF(IDCENV.EQ.2)THEN DO 119 N=1,NBME CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, & NES,NP,NPG,IAXI,XCOOR, & IPT1.NUM,KDEB,KFIN,LRV,IDCENV,CMD,IKOMP, & IZTU1.VPOCHA(1,N),IPADI.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU, & ANUK,WT(1,1,1,N),WS(1,1,1,N),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 119 CONTINUE ELSE CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, & NES,NP,NPG,IAXI,XCOOR, & IPT1.NUM,KDEB,KFIN,LRV,IDCENV,CMD,IKOMP, & IZTU1.VPOCHA,IPADI.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU, & ANUK,WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) ENDIF C DO 120 N=1,NBME N1=1 IF(IDCEN.EQ.2)N1=N C C- Matrice masse consistante C IF (KFORM.EQ.1) THEN IZAFM = LIZAFM(L,N) DO 80 K=KDEB,KFIN KP = K - KDEB + 1 DO 70 I=1,NP DO 60 J=1,NP UU = 0.D0 DO 50 LG=1,NPG UU=UU+FN(I,LG)*WT(KP,J,LG,N1)*PGSK(KP,LG) & *DEUPI*RPGK(KP,LG) 50 CONTINUE AM(K,I,J)=UU*COEFK(KP)*CN 60 CONTINUE 70 CONTINUE 80 CONTINUE C C- Matrice masse condensée C ELSE IZAFM = LIZAFM(1,N) DO 110 K=KDEB,KFIN KP = K - KDEB + 1 DO 100 J=1,NP UU = 0.D0 DO 90 LG=1,NPG UU = UU + WT(KP,J,LG,N1)*PGSK(KP,LG) & *DEUPI*RPGK(KP,LG) 90 CONTINUE KJ = MLENTI.LECT(IPT1.NUM(J,K)) AM(KJ,1,1)=AM(KJ,1,1) + (UU*COEFK(KP)*CN) 100 CONTINUE 110 CONTINUE ENDIF 120 CONTINUE 130 CONTINUE K0 = K0 + NBEL SEGSUP PETROV 140 CONTINUE C C-------------------------------------------------- C-- Cas ou l'inconnue est au CENTRE (VF ou EFM1) -- C-------------------------------------------------- C ELSEIF (KPOIND.EQ.2) THEN DO 160 N=1,NBME IZAFM = LIZAFM(1,N) DO 150 I=1,NBEL NKR = (1-IK1)*(I-1) + 1 AM(I,1,1) = (VPOCHA(I,1) * IZRO.VPOCHA(NKR,1)*CN)/DTR 150 CONTINUE 160 CONTINUE C C------------------------------------------------------------ C-- Cas ou l'inconnue est au CENTREP0 (EFMC MACRO ou QUAD) -- C------------------------------------------------------------ C ELSEIF (KPOIND.EQ.3) THEN DO 163 N=1,NBME IZAFM = LIZAFM(1,N) DO 153 I=1,NBEL NKR = (1-IK1)*(I-1) + 1 AM(I,1,1) = (VPOCHA(I,1) * IZRO.VPOCHA(NKR,1)*CN)/DTR 153 CONTINUE 163 CONTINUE C C------------------------------------------------------------ C-- Cas ou l'inconnue est au CENTREP1 (EFMC MACRO ou QUAD) -- C------------------------------------------------------------ C ELSEIF (KPOIND.EQ.4) THEN DO N=1,NBME IZAFM = LIZAFM(1,N) SEGACT MELEMO NBSOUS=MELEMO.LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 DO LS=1,NBSOUS IPT1=MELEMO IF(NBSOUS.NE.1)IPT1=MELEMO.LISOUS(LS) SEGACT IPT1 NP1=IPT1.NUM(/1) NBEL1=IPT1.NUM(/2) IF(MQUAD.NE.0)THEN WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'PRP1' ELSEIF(MACRO.NE.0)THEN WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'MCP1' 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 MP=IZF1.FN(/1) NKR=1 DO K=1,NBEL DO I1=1,NP1 J1=IPT1.NUM(I1,K) DO ND=1,IDIM XYZ(ND,I1)=XCOOR((J1-1)*(IDIM+1)+ND) ENDDO ENDDO CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES, & IDIM,NP1,NPG,IAXI,AIR1) DO 154 I=1,MP DO 155 J=1,MP U=0.D0 DO 157 L=1,NPG U=U+IZF1.FN(I,L)*IZF1.FN(J,L)*PGSQ(L) 157 CONTINUE AM(K,I,J)=(U*IZRO.VPOCHA(NKR,1) * CN) / DTR 155 CONTINUE NKR = NKR + (1-IK1) 154 CONTINUE ENDDO ENDDO ENDDO ENDIF ENDIF C C-------------------------- C Création du second-membre C-------------------------- C Saturation de la matrice calculée ci-dessus par l'inconnue C au pas de temps précédent C C C- Construction du chapeau de l'objet CHAMPOIN qui contiendra le RHS C NAT = 2 NSOUPO = 1 SEGINI MCHPO1 MCHPO1.IFOPOI = IFOUR MCHPO1.MOCHDE = TITREE MCHPO1.MTYPOI = 'SMBR' MCHPO1.JATTRI(1) = 2 NC = NINKO SEGINI MSOUP1 MCHPO1.IPCHP(1) = MSOUP1 DO 170 N=1,NINKO MSOUP1.NOCOMP(N) = LISDUA(N) 170 CONTINUE MSOUP1.IGEOC = MELEM1 SEGACT MELEM1 N = MELEM1.NUM(/2) SEGINI MPOVA1 MSOUP1.IPOVAL = MPOVA1 SEGACT MELEM1 NBSOUS = LIZAFM(/1) C MELEME = IRIGEL(1,1) IF(ISCHT.NE.2)THEN DO 220 L=1,NBSOUS IPT1 = MELEME IF (NBSOUS.NE.1) IPT1=LISOUS(L) SEGACT IPT1 NP = IPT1.NUM(/1) NBEL = IPT1.NUM(/2) DO 210 N=1,NINKO IZAFM = LIZAFM(L,N) SEGACT IZAFM DO 200 K=1,NBEL DO 190 J=1,NP UU = 0.D0 IU = LECT(IPT1.NUM(J,K)) DO 180 I=1,NP IK = IPADI.LECT(IPT1.NUM(I,K)) UU = UU + AM(K,I,J)*IZTGG2.VPOCHA(IK,N)/CN*CNM 180 CONTINUE MPOVA1.VPOCHA(IU,N) = MPOVA1.VPOCHA(IU,N) + UU 190 CONTINUE 200 CONTINUE 210 CONTINUE 220 CONTINUE ENDIF C Cas BDF2 M Tn-2 IF(ISCHT.EQ.1)THEN DO 222 L=1,NBSOUS IPT1 = MELEME IF (NBSOUS.NE.1) IPT1=LISOUS(L) SEGACT IPT1 NP = IPT1.NUM(/1) NBEL = IPT1.NUM(/2) DO 212 N=1,NINKO IZAFM = LIZAFM(L,N) SEGACT IZAFM DO 202 K=1,NBEL DO 192 J=1,NP UU = 0.D0 IU = LECT(IPT1.NUM(J,K)) DO 182 I=1,NP IK = IPADI.LECT(IPT1.NUM(I,K)) UU = UU + AM(K,I,J)*IZTGN2.VPOCHA(IK,N)/CN*CNMM 182 CONTINUE MPOVA1.VPOCHA(IU,N) = MPOVA1.VPOCHA(IU,N) + UU 192 CONTINUE 202 CONTINUE 212 CONTINUE 222 CONTINUE ENDIF C C- Sauvegarde de l'objet MATRIK à l'indice MATELM de la table KIZX et du C- second-membre à l'indice SMBR de la table EQEX (assemblage éventuel) C CALL ECROBJ('MATRIK',MATRIK) CALL ECROBJ('CHPOINT ',MCHPO1) C C- Désactivation et ménage C SEGDES MCHPO1,MSOUP1,MPOVA1 SEGDES IZTGG2 SEGDES IMATRI,MATRIK*NOMOD IF (IK1DETR.EQ.1) SEGSUP IZRO IF (IKMUDETR.EQ.1) SEGSUP MZMU IF(IPADS.NE.IPADI)SEGSUP IPADS IF(MLENTI.NE.IPADI)SEGSUP IPADI SEGSUP IPADI RETURN 1002 FORMAT(10(1X,1PE11.4)) END