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