1 : C DIVU SOURCE CHAT 11/03/16 21:19:15 6902 2 : SUBROUTINE DIVU 3 : C----------------------------------------------------------------------- 4 : C Calcul la divergence d'un champ de vecteur (en général la vitesse) 5 : C discretise dans l'espace de Raviart Thomas de plus petit degre. 6 : C----------------------------------------------------------------------- 7 : C 8 : C On obtient la DIVERGENCE de la vitesse en sommant les débits sur 9 : C chaque élément, chaque débit étant de support la normale sortante. 10 : C C'est pourquoi nous avons besoin du MCHAML d'orientation des normales 11 : C de référence afin de modifier la valeur algébrique du débit. 12 : C 13 : C--------------------------- 14 : C Phrase d'appel (GIBIANE) : 15 : C--------------------------- 16 : C 17 : C CHP1 = DIVU MMODEL TABLE1 CHP2 CHAM1 (CHP3 (CHP4)) 18 : C 19 : C------------------------ 20 : C Opérandes et résultat : 21 : C------------------------ 22 : C 23 : C CHP1 : CHAMPOINT centre contenant la divergence de la vitesse 24 : C MMODEL : Objet MMODEL spécifiant la formulation 25 : C TABLE1 : Objet table DOMAINE regroupant les infos géometriques 26 : C CHP2 : CHAMPOINT contenant le debit aux faces 27 : C CHAM1 : MCHAML sens de la normale à la face (1=out , -1=in) 28 : C CHP3 : CHAMPOIN contenant une fonction aux centres des elements 29 : C CHP4 : CHAMPOIN contenant les valeurs de cette fonction sur les 30 : C limites du domaine(le support est contitué par des points face) 31 : C 32 : C----------------------------------------------------------------------- 33 : C 34 : C Langage : ESOPE + FORTRAN77 35 : C 36 : C Auteurs : F.DABBENE 01/94 37 : C Complements : C. LE POTIER ET F. AURIOL 20/00 38 : C----------------------------------------------------------------------- 39 : IMPLICIT INTEGER(I-N) 40 : IMPLICIT REAL*8 (A-H,O-Z) 41 : * 42 : -INC CCOPTIO 43 : -INC SMCHAML 44 : -INC SMCHPOI 45 : -INC SMELEME 46 : -INC SMMODEL 47 : -INC SMTABLE 48 : * 49 : SEGMENT IPMAHY 50 : INTEGER MAHYBR(NSOUS) 51 : ENDSEGMENT 52 : * 53 : REAL*8 XVALIN,XVALRE 54 : LOGICAL LOGRE,LOGIN 55 : CHARACTER*4 NOMTOT(1) 56 : CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE 57 : * 58 : * Initialisations 59 : * 60 : IVALIN = 0 61 : XVALIN = 0.D0 62 : LOGIN = .TRUE. 63 : IOBIN = 0 64 : TAPIND = 'MOT ' 65 : CHARIN = 'MAILLAGE' 66 : TYPOBJ = 'MAILLAGE' 67 : * 68 : * Lecture du MMODEL 69 : * 69 : CALL LIROBJ('MMODEL',IPMODE,1,IRET1) 70 : IF (IERR.NE.0) RETURN 71 : * 72 : * Lecture de la TABLE domaine 73 : * 74 : IPTABL = 0 74 : CALL LEKMOD(IPMODE,IPTABL,IRET) 75 : IF (IERR.NE.0) RETURN 75 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 76 : IF (IERR.NE.0) RETURN 77 : IPGEOM = IOBRE 77 : CALL LEKTAB(IPTABL,'ELTFA',IOBRE) 78 : IF (IERR.NE.0) RETURN 79 : IELTFA = IOBRE 79 : CALL LEKTAB(IPTABL,'CENTRE',IOBRE) 80 : IF (IERR.NE.0) RETURN 81 : IPCENT = IOBRE 81 : CALL LEKTAB(IPTABL,'FACE',IOBRE) 82 : IF (IERR.NE.0) RETURN 83 : IPFACE = IOBRE 83 : CALL LEKTAB(IPTABL,'FACEL',IOBRE) 84 : IF (IERR.NE.0) RETURN 85 : IFACEL = IOBRE 86 : * 87 : * Lecture du CHAMPOINT de debit (nom de composantes FLUX) 88 : * 88 : CALL LIROBJ('CHPOINT',ICHP1,1,IRET1) 89 : IF (IERR.NE.0) RETURN 90 : * 91 : * Lecture du CHAMELEM 92 : * 92 : CALL LIROBJ('MCHAML',IPCHEL,1,IRET1) 92 : C CALL LEKTAB(IPTABL,'XXNORMAE',IPCHEL) 93 : IF (IERR.NE.0) RETURN 94 : * 95 : * Lecture du CHAMPOINT de la fonction aux centres 96 : * 97 : IPFONC=0 97 : CALL LIROBJ('CHPOINT',ICHP2,0,IRET2) 98 : IF (IERR.NE.0) RETURN 99 : IF(IRET2.EQ.1)THEN 100 : INDIC = 1 101 : NBCOMP = 1 102 : NOMTOT(1) = 'SCAL' 102 : CALL QUEPOI(ICHP2,IPCENT,INDIC,NBCOMP,NOMTOT) 103 : IF (IERR.NE.0) RETURN 104 : ICLIM=0 104 : CALL LIROBJ('CHPOINT',ICLIM,0,IRLIM) 104 : CALL DIVUFN(ICHP2,ICLIM,IPFACE,IFACEL,ICHP1,IPCHEL,IPFONC) 105 : 106 : * VERRUE 106 : * CALL ECCHPO(ipfonc) 107 : ENDIF 108 : * 109 : * Test CHAMPOINT de FLUX 110 : * 111 : INDIC = 1 112 : NBCOMP = 1 113 : NOMTOT(1) = 'FLUX' 113 : CALL QUEPOI(ICHP1,IPFACE,INDIC,NBCOMP,NOMTOT) 114 : IF (IERR.NE.0) RETURN 115 : * 116 : * Test de la formulation contenu dans MMODEL 117 : * 118 : MMODEL = IPMODE 119 : SEGACT MMODEL 120 : NSOUS = KMODEL(/1) 121 : SEGINI IPMAHY 122 : IDARCY = 0 123 : IPT1=IPGEOM 124 : IPT2=IPGEOM 125 : DO 10 ISOUS=1,NSOUS 126 : IF(NSOUS.GT.1)THEN 127 : SEGACT IPT2 128 : IPT1=IPT2.LISOUS(ISOUS) 129 : SEGDES IPT2 130 : ENDIF 131 : IMODEL = KMODEL(ISOUS) 132 : SEGACT IMODEL 133 : LETYPE = FORMOD(1) 134 : IF (LETYPE.EQ.'DARCY') THEN 135 : IDARCY = IDARCY + 1 136 : MAHYBR(ISOUS) = IPT1 137 : ENDIF 138 : SEGDES IMODEL 139 : 10 CONTINUE 140 : SEGDES MMODEL 141 : IF (IDARCY.EQ.0) THEN 142 : SEGSUP IPMAHY 143 : MOTERR = LETYPE 143 : CALL ERREUR(193) 144 : RETURN 145 : ENDIF 146 : * 147 : * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini 148 : * 149 : MELEME = IELTFA 150 : SEGACT MELEME 151 : LZONES = LISOUS(/1) 152 : IF (LZONES.EQ.0) LZONES = 1 153 : IPT1 = IPGEOM 154 : SEGACT IPT1 155 : DO 30 ISOUS=1,NSOUS 156 : IMAMEL = MAHYBR(ISOUS) 157 : IF (IMAMEL.NE.0) THEN 158 : DO 20 ISZ=1,LZONES 159 : IF (LZONES.EQ.1) THEN 160 : IPT2 = IPT1 161 : IPT3 = MELEME 162 : ELSE 163 : IPT2 = IPT1.LISOUS(ISZ) 164 : IPT3 = LISOUS(ISZ) 165 : ENDIF 166 : IF (IPT2.EQ.IMAMEL) THEN 167 : MAHYBR(ISOUS) = IPT3 168 : GOTO 30 169 : ENDIF 170 : 20 CONTINUE 171 : IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN 172 : MOTERR(1:8) = ' MODELE ' 173 : MOTERR(1:16) = ' ELTFA ' 174 : INTERR(1) = ISOUS 174 : CALL ERREUR(698) 175 : SEGDES IPT1 176 : SEGDES MELEME 177 : SEGSUP IPMAHY 178 : RETURN 179 : ENDIF 180 : ENDIF 181 : 30 CONTINUE 182 : SEGDES IPT1 183 : SEGDES MELEME 184 : * 185 : * Test du SUPPORT du MCHAML d'orientation : On teste qu'à chaque zone 186 : * du MODELE ou Darcy est défini correspond une zone du MCHAML 187 : * 188 : MCHELM = IPCHEL 189 : SEGACT MCHELM 190 : NCH1 = IMACHE(/1) 191 : NBOK = 0 192 : DO 50 ISOUS=1,NSOUS 193 : IMAMEL = MAHYBR(ISOUS) 194 : IF (IMAMEL.NE.0) THEN 195 : DO 40 JSOUS=1,NCH1 196 : IF (IMACHE(JSOUS).EQ.IMAMEL) NBOK=NBOK+1 197 : 40 CONTINUE 198 : ENDIF 199 : 50 CONTINUE 200 : SEGDES MCHELM 201 : SEGSUP IPMAHY 202 : IF (NBOK.NE.NCH1) THEN 203 : MOTERR(1:8) = ' ORIENT ' 204 : MOTERR(9:16)= ' ELTFA ' 204 : CALL ERREUR(135) 205 : RETURN 206 : ENDIF 207 : * 208 : * Calcul de la divergence du champ de vecteur connu via ses flux 209 : * 209 : CALL DIVU1(ICHP1,IPCHEL,IPFACE,IPCENT,IPFONC,IRET) 209 : CALL ECROBJ('CHPOINT',IRET) 210 : IF(IPFONC.NE.0)THEN 210 : CALL ECROBJ('CHPOINT',IPFONC) 210 : * CALL ECCHPO(ipfonc) 210 : * CALL ECCHPO(iret) 211 : MCHPOI=IPFONC 212 : SEGACT MCHPOI 213 : MSOUPO=IPCHP(1) 214 : SEGACT MSOUPO 215 : MPOVAL=IPOVAL 216 : SEGDES MPOVAL,MSOUPO,MCHPOI 217 : ENDIF 218 : * 219 : RETURN 220 : END 221 : 222 : 223 : 224 : 225 : 226 : 227 : 228 :
© Cast3M 2003 - All rights reserved.
Disclaimer