Download zchimp.eso
Back to the list
1 : C ZCHIMP SOURCE CHAT 05/01/13 04:21:24 5004
2 : SUBROUTINE ZCHIMP(MTABX,MTAB1)
3 : C=======================================================================
4 : C L'opérateur ECHIMP discrétise un échange surfacique ou volumique
5 : C avec un milieu exterieur. Le terme d'échange est de la forme
6 : * h(T-T0) où h est un coefficient d'échange, T une des inconnues du
7 : C problème traité et T0 "un champ exterieur connu".
8 : C-----------------------------------------------------------------------
9 : C La convention de signe associée à ce terme est la suivante : lorsque
10 : C h est positif et T0 En explicite :
52 : C
53 : C On suppose que les inconnues primales et duales sont identiques
54 : C Pour cela, on condense la matrice masse, quelle que soit la
55 : C formulation utilisée (donc EF -> EFM1).
56 : C 1) La matrice diagonale est stockée dans un CHPO et rangée dans la
57 : C table KIZG1 à l'indice de type MOT MOT1 (nom de l'inconnue).
58 : C 2) Le second membre est stocké dans un CHPO et rangé dans la
59 : C table KIZG à l'indice de type MOT MOT1 (nom de l'inconnue).
60 : C
61 : C
62 : C Les inconnues primales et duales peuvent etre différentes et de
63 : C masse (donc EFM1 -> EF).
64 : C 1) La matrice "masse" est stockée dans un MATRIK et rangée dans la
65 : C table KIZX à l'indice de type MOT MATELM.
66 : C 2) Le second membre est stocké dans un CHPO et assemblé dans la
67 : C table EQEX à l'indice de type MOT SMBR. Le nom de l'inconnue
68 : C duale MOT2 étant le nom de la composante du CHPO créé.
69 : C
70 : C-------------------------
71 : C Formulations acceptées :
72 : C-------------------------
73 : C 0) Défaut : EFM1 explicite.
74 : C
75 : C----------------------
76 : C Support des données :
77 : C----------------------
78 : C H : SCAL ou CHPO SCAL CENTRE
79 : C T0 : SCAL ou CHPO SCAL CENTRE ou CHPO SCAL SOMMET
80 : C
81 : C------------------------------------------
82 : C Indices de table utilisés (racine KIZX) :
83 : C------------------------------------------
84 : C 'DOMZ' : Table domaine associé à l'opérateur (domaine local)
85 : C E/ 'MAILLAGE' -> Maillage du domaine local
86 : C E/ 'CENTRE' -> Points centre du domaine local
87 : C E/ 'XXPSOML' -> Intégrale des fonctions de base par élément (MCHAML)
88 : C E/ 'IARG' : Nombre d'arguments de l'opérateur
89 : C E/ 'ARG1' : Premier argument (coefficient d'échange)
90 : C E/ 'ARG2' : Deuxième argument (champ exterieur)
91 : C E/ 'LISTINCO' : Listmot contenant le nom de l'inconnue
92 : C 'KOPT' : Table des options
93 : C E/ 'KFORM' -> Formulation spatiale
94 : C E/ 'KIMPL' -> Formulation temporelle
95 : C 'EQEX' : Table décrivant la modélisation (modèle fluide)
96 : C E/ 'INCO' -> Table des inconnues et des données
97 : C /S 'KIZG' -> Table des seconds membres (cas explicite)
98 : C /S 'KIZG1' -> Table des matrices diagonales (cas explicite)
99 : C-----------------------------------------------------------------------
100 : C
101 : IMPLICIT INTEGER(I-N)
102 : IMPLICIT REAL*8 (A-H,O-Z)
103 : C
104 : -INC CCOPTIO
105 : -INC SMCHPOI
106 : -INC SMLENTI
107 : -INC SMELEME
108 : POINTEUR MELEMC.MELEME
109 : -INC SMLMOTS
110 : C
111 : CHARACTER*8 NOMD,NOMP,TYPE,TYPC
112 : PARAMETER (NTB=1)
113 : CHARACTER*8 LTAB(NTB)
114 : DIMENSION KTAB(NTB),IXV(4)
115 : DATA LTAB/'KIZX '/
116 : C
117 : C- Récupération de la table DOMAINE associée au domaine local
118 : C
118 : CALL LEKTAB(MTABX,'DOMZ',MTABZ)
119 : IF (MTABZ.EQ.0) THEN
120 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
121 : MOTERR( 1: 8) = ' DOMZ '
122 : MOTERR( 9:16) = ' DOMZ '
123 : MOTERR(17:24) = ' KIZX '
123 : CALL ERREUR(786)
124 : RETURN
125 : ENDIF
126 : C
127 : C- Récupération de la table INCO (pointeur KINC)
128 : C
128 : CALL LEKTAB(MTAB1,'INCO',KINC)
129 : IF (KINC.EQ.0) THEN
130 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
131 : MOTERR( 1: 8) = ' INCO '
132 : MOTERR( 9:16) = ' INCO '
133 : MOTERR(17:24) = ' EQEX '
133 : CALL ERREUR(786)
134 : RETURN
135 : ENDIF
136 : C
137 : C- Identification du type d'échange :
138 : C- Surfacique (ISURF1=1) ou volumique (IVOL1=1).
139 : C
140 : ISURF1 = 0
141 : IVOL1 = 0
141 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
142 : IF (IERR.NE.0) RETURN
143 : SEGACT MELEME
144 : NBSOUS = LISOUS(/1)
145 : IF (NBSOUS.EQ.0) NBSOUS=1
146 : DO 10 I=1,NBSOUS
147 : IPT1=MELEME
148 : IF (NBSOUS.NE.1) IPT1=LISOUS(I)
149 : SEGACT IPT1
150 : C write(6,*) '***************************',ipt1.itypel
151 : IF (IDIM.EQ.2) THEN
152 : IF (IPT1.ITYPEL.EQ.2) THEN
153 : ISURF1 = 1
154 : ELSEIF (IPT1.ITYPEL.EQ.4.OR.IPT1.ITYPEL.EQ.8) THEN
155 : IVOL1 = 1
156 : ELSE
157 : C Type d'élément incorrect
157 : CALL ERREUR(16)
158 : RETURN
159 : ENDIF
160 : ELSEIF (IDIM.EQ.3) THEN
161 : IF (IPT1.ITYPEL.EQ.4.OR.IPT1.ITYPEL.EQ.8) THEN
162 : ISURF1 = 1
163 : ELSEIF (IPT1.ITYPEL.EQ.14.OR.IPT1.ITYPEL.EQ.16.OR.
164 : C pb PYR5 IPT1.ITYPEL.EQ.23.OR.IPT1.ITYPEL.EQ.25) THEN
165 : & IPT1.ITYPEL.EQ.23) THEN
166 : IVOL1 = 1
167 : ELSE
168 : C Type d'élément incorrect
168 : CALL ERREUR(16)
169 : RETURN
170 : ENDIF
171 : ENDIF
172 : SEGDES IPT1
173 : 10 CONTINUE
174 : IF (IVOL1.EQ.1.AND.ISURF1.EQ.1) THEN
175 : C Maillage incorrect : contient des éléments 1D et 2D
175 : CALL ERREUR(798)
176 : RETURN
177 : ENDIF
178 : IF (NBSOUS.NE.1) SEGDES MELEME
179 : C
180 : C- Récupération du nom des inconnues (primale et duale)
181 : C
182 : TYPE = 'LISTMOTS'
182 : CALL ACMO(MTABX,'LISTINCO',TYPE,MLMOTS)
183 : IF (IERR.NE.0) RETURN
184 : SEGACT MLMOTS
185 : NBINC = MOTS(/2)
186 : IF (NBINC.LE.0.OR.NBINC.GE.3) THEN
187 : C Indice %m1:8 : contient plus de %i1 %m9:16
188 : MOTERR( 1:8) = 'LISTINCO'
189 : INTERR(1) = 2
190 : MOTERR(9:16) = ' MOTS '
190 : CALL ERREUR(799)
191 : RETURN
192 : ENDIF
193 : NOMP = MOTS(1)
194 : IF (NBINC.EQ.1) THEN
195 : NOMD = MOTS(1)
196 : ELSE
197 : NOMD = MOTS(2)
198 : ENDIF
199 : SEGDES MLMOTS
200 : C
201 : C- Informations géométriques locales associées à l'opérateur
202 : C
202 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
203 : IF (IERR.NE.0) RETURN
203 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
204 : IF (IERR.NE.0) RETURN
205 : C
206 : C- Identification du spg de l'inconnue duale
207 : C- Trois cas sont possible (identifier par IKAS)
208 : C- 1) spg contenant melems -> formulation EF
209 : C- 2) spg contenant melemc -> formulation VF si IVOL1=1; ERREUR sinon
210 : C- 3) ni 2) ni 3) -> formulation VF si IVOL1=0; ERREUR sinon
211 : C
212 : TYPE = ' '
213 : IKAS = 0
213 : CALL ACMO(KINC,NOMD,TYPE,MCHPOI)
214 : IF (TYPE.NE.'CHPOINT ') THEN
215 : C Indice %m1:8 : ne contient pas un objet de type %m9:16
216 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
217 : MOTERR(9:16) = 'CHPOINT '
217 : CALL ERREUR(800)
218 : RETURN
219 : ELSE
219 : CALL LICHT(MCHPOI,MPOVAL,TYPC,MELEMD)
220 : NC = VPOCHA(/2)
221 : SEGDES MPOVAL
222 : IF (NC.NE.1) THEN
223 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
224 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
225 : MOTERR(9:16) = 'CHPOINT '
225 : CALL ERREUR(784)
226 : RETURN
227 : ENDIF
228 : ENDIF
228 : CALL KRIPAD(MELEMD,MLENTI)
228 : CALL VERPAD(MLENTI,MELEMS,IRET1)
228 : CALL VERPAD(MLENTI,MELEMC,IRET2)
229 : IF (IRET1.EQ.0.AND.IRET2.EQ.1) THEN
230 : IKAS = 1
231 : ELSEIF (IRET1.EQ.1.AND.IRET2.EQ.0) THEN
232 : IKAS = 2
233 : IF (IVOL1.EQ.0) THEN
234 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
235 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
236 : MOTERR(9:16) = 'CHPOINT '
236 : CALL ERREUR(788)
237 : RETURN
238 : ENDIF
239 : ELSEIF (IRET1.EQ.1.AND.IRET2.EQ.1) THEN
240 : IKAS = 3
241 : IF (IVOL1.EQ.1) THEN
242 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
243 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
244 : MOTERR(9:16) = 'CHPOINT '
244 : CALL ERREUR(788)
245 : RETURN
246 : ENDIF
247 : SEGINI,IPT4=MELEMC
248 : ELSE
249 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
250 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
251 : MOTERR(9:16) = 'CHPOINT '
251 : CALL ERREUR(788)
252 : RETURN
253 : ENDIF
254 : C
255 : C- Identification du spg de l'inconnue primale
256 : C- Deux cas sont possibles (identifier par IKP)
257 : C- 1) spg contenant melemc -> IKP=0
258 : C- 2) spg contenant melems -> IKP=4
259 : C
260 : IF (NBINC.EQ.1) THEN
261 : IF (IKAS.EQ.1) THEN
262 : IKP = 4
263 : ELSE
264 : IKP = 0
265 : ENDIF
266 : ELSE
267 : TYPE = ' '
267 : CALL ACMO(KINC,NOMP,TYPE,MCHPOI)
268 : IF (TYPE.NE.'CHPOINT ') THEN
269 : C Indice %m1:8 : ne contient pas un objet de type %m9:16
270 : MOTERR(1: 8) = 'INCO'//NOMP(1:4)
271 : MOTERR(9:16) = 'CHPOINT '
271 : CALL ERREUR(800)
272 : RETURN
273 : ELSE
273 : CALL LICHT(MCHPOI,MPOVAL,TYPC,MELEMP)
274 : NC = VPOCHA(/2)
275 : SEGDES MPOVAL
276 : IF (NC.NE.1) THEN
277 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
278 : MOTERR(1: 8) = 'INCO'//NOMP(1:4)
279 : MOTERR(9:16) = 'CHPOINT '
279 : CALL ERREUR(784)
280 : RETURN
281 : ENDIF
282 : ENDIF
282 : CALL KRIPAD(MELEMP,MLENT3)
282 : CALL VERPAD(MLENT3,MELEMS,IRET1)
282 : CALL VERPAD(MLENT3,MELEMC,IRET2)
283 : SEGSUP MLENT3
284 : IF (IRET1.EQ.0.AND.IRET2.EQ.1) THEN
285 : IKP = 4
286 : ELSEIF (IRET1.EQ.1.AND.IRET2.EQ.0) THEN
287 : IKP = 0
288 : ELSE
289 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
290 : MOTERR(1: 8) = 'INCO'//NOMP(1:4)
291 : MOTERR(9:16) = 'CHPOINT '
291 : CALL ERREUR(788)
292 : RETURN
293 : ENDIF
294 : ENDIF
295 : C
296 : C- On vérifie dans le cas surfacique que la surface est une frontière
297 : C- du domaine global dont la table DOMAINE est à l'indice 'PERE' de
298 : C- la table domaine locale.
299 : C
300 : C- Si IKAS=3, récupération des points CENTREs du maillage volumique
301 : C- associés aux points CENTREs du maillage surfacique.
302 : C
303 : IF (IVOL1.EQ.0) THEN
304 : TYPE = 'TABLE '
304 : CALL ACMO(MTABZ,'PERE',TYPE,KPERE)
305 : IF (IERR.NE.0) RETURN
305 : CALL LEKTAB(KPERE,'FACE ',IPT2)
305 : CALL KRIPAD(IPT2,MLENT2)
305 : CALL VERPAD(MLENT2,MELEMC,IRET)
306 : IF (IRET.EQ.1) THEN
307 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
308 : MOTERR(1: 8) = 'INCO'//NOMD(1:4)
309 : MOTERR(9:16) = 'CHPOINT '
309 : CALL ERREUR(788)
310 : RETURN
311 : ENDIF
311 : CALL LEKTAB(KPERE,'FACEL ',IPT3)
312 : IF (IERR.NE.0) RETURN
313 : SEGACT MELEMC
314 : SEGACT IPT3
315 : SEGACT MLENT2
316 : NBC = MELEMC.NUM(/2)
317 : DO 20 I=1,NBC
318 : J = MLENT2.LECT(MELEMC.NUM(1,I))
319 : IF (IPT3.NUM(2,J) .NE. MELEMC.NUM(1,I)) THEN
320 : C Incohérence entre tables DOMAINE
320 : CALL ERREUR(801)
321 : RETURN
322 : ELSE
323 : IF ( IPT3.NUM(1,J) .NE. IPT3.NUM(3,J)) THEN
324 : C La face %i1 n'est pas une frontière
325 : INTERR(1) = IPT3.NUM(2,J)
325 : CALL ERREUR(802)
326 : RETURN
327 : ENDIF
328 : IF (IKAS.EQ.3) THEN
329 : IPT4.NUM(1,I) = IPT3.NUM(1,J)
330 : ENDIF
331 : ENDIF
332 : 20 CONTINUE
333 : SEGDES MELEMC
334 : SEGDES IPT3
335 : IF (IVOL1.EQ.0 .AND. IKAS.EQ.3) SEGDES IPT4
336 : SEGDES MLENT2
337 : ENDIF
338 : C
339 : C- Récupération de la table des options KOPT (pointeur KOPTI)
340 : C- et initialisations des options (par défaut SI et EXPL)
341 : C- KFORM = 0 -> SI 1 -> EF 2 -> VF
342 : C- KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN
343 : C
344 : KFORM = 0
345 : KIMPL = 0
346 : TYPE = ' '
346 : CALL ACMO(MTABX,'KOPT',TYPE,KOPTI)
347 : IF (TYPE.EQ.'TABLE') THEN
348 : TYPE = ' '
348 : CALL ACMO(KOPTI,'KFORM',TYPE,IENT)
348 : IF (TYPE.EQ.'ENTIER') CALL ACME(KOPTI,'KFORM',KFORM)
349 : TYPE = ' '
349 : CALL ACMO(KOPTI,'KIMPL',TYPE,IENT)
349 : IF (TYPE.EQ.'ENTIER') CALL ACME(KOPTI,'KIMPL',KIMPL)
350 : ENDIF
351 : C
352 : C- Cohérance des options avec le support de l'inconnue duale et IKAS
353 : C Option %m1:8 incompatible avec les données
354 : IF (IKAS.EQ.1 .AND. KFORM.EQ.2) THEN
355 : MOTERR(1:8) = ' VF '
355 : CALL ERREUR(803)
356 : RETURN
357 : ENDIF
358 : IF (IKAS.NE.1 .AND. KFORM.NE.2) THEN
359 : MOTERR(1:8) = ' EF '
359 : CALL ERREUR(803)
360 : RETURN
361 : ENDIF
362 : IF (NBINC.NE.1 .AND. KIMPL.EQ.0) THEN
363 : MOTERR(1:8) = ' EXPL '
363 : CALL ERREUR(803)
364 : RETURN
365 : ENDIF
366 : C
367 : C- Récupération des arguments (2 arguments attendus) :
368 : C- 1) Coefficient d'échange
369 : C- 2) Champ exterieur connu
370 : C
370 : CALL ACME(MTABX,'IARG',IARG)
371 : IF (IERR.NE.0) RETURN
372 : IF (IARG.NE.2) THEN
373 : C Indice %m1:8 : nombre d'arguments incorrect
374 : MOTERR(1:8) = 'IARG '
374 : CALL ERREUR(804)
375 : RETURN
376 : ENDIF
377 : IXV(1) = MELEMC
378 : IXV(2) = 1
379 : IXV(3) = 0
380 : IRET = 0
380 : CALL LEKCOF('Opérateur ECHI :',
381 : & MTABX,KINC,1,IXV,MCHPO1,MPOVA1,NPT1,NC1,IKH,IRET)
382 : IF (IRET.EQ.0) RETURN
383 : IXV(1) = MELEMC
384 : IXV(2) = 1
385 : IXV(3) = 0
386 : IXV(4) = MELEMS
387 : IRET = 4
387 : CALL LEKCOF('Opérateur ECHI :',
388 : & MTABX,KINC,2,IXV,MCHPO2,MPOVA2,NPT2,NC2,IKT,IRET)
389 : IF (IRET.EQ.0) RETURN
390 : IF (IKT.EQ.0) THEN
391 : MELEME = MELEMC
392 : ELSE
393 : MELEME = MELEMS
394 : ENDIF
394 : CALL KRIPAD(MELEME,MLENT1)
395 : C
396 : C- Calcul de la discrétisation du terme d'échange
397 : C
398 : IF (KIMPL.EQ.0) THEN
398 : CALL ZECHI1(IKAS,IVOL1,MTAB1,MTABZ,MPOVA1,MPOVA2,IKH,IKT,
399 : & MELEMD,IPT4,MLENTI,MLENT1,NOMD)
400 : ELSE
400 : CALL ZECHI2(IKAS,IVOL1,MTAB1,MTABZ,MTABX,MPOVA1,MPOVA2,
401 : & IKH,IKT,IKP,IPT4,MLENT1,NOMP,NOMD)
402 : ENDIF
403 : RETURN
404 : END
405 :
406 :
407 :