Download ytclsf.eso
Back to the list
1 : C YTCLSF SOURCE FANDEUR 13/02/07 21:15:17 7694
2 : SUBROUTINE YTCLSF(CHHH,NOMPR)
3 : IMPLICIT INTEGER(I-N)
4 : IMPLICIT REAL*8 (A-H,O-Z)
5 : C***********************************************************************
6 : C
7 : C CET OPERATEUR DISCRETISE TOUT
8 : C
9 : C SYNTAXE :
10 : C ---------
11 : C
12 : C s m = OPER $mod LSMOT1 LSMOT2
13 : C DT SCHT TN >
14 : C TN
15 : C DT SCHT TN >
16 : C >>
17 : C
18 : C
19 : C OPTIONs LINE MACRO QUAF LINB ISOQ
20 : C CENTREP0 CENTREP1 CENTRE MSOMMET RIGIDITE MDIA
21 : C CENTREE SUPG SUPGDC CMD a CONS STABP
22 : C xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx
23 : C
24 : C SCHT EUL_EXPL EUL_IMPL CN CNG BDF2 SEMI b
25 : C
26 : C
27 : C
28 : C COEFFICIENTS :
29 : C --------------
30 : C
31 : C
32 : C ALF (SCAL DOMA) DIFFUSIVITE THERMIQUE
33 : C (SCAL ELEM)
34 : C (SCAL NOEU)
35 : C
36 : C INCONNUES :
37 : C -----------
38 : C
39 : C TN CHAMP DE TEMPERATURE
40 : C UN CHAMP DE VITESSE
41 : C PN CHAMP DE PRESSION
42 : C
43 : C INCOD Type inconnue duale : VITEsse TEMPerature ou PRESsion
44 : C MPRE Type pression : CENTRE CENTREP1 CENTREP0 MSOMMET
45 : C
46 : C IHV=0 SCALAIRE
47 : C IHV=1 VITESSE
48 : C
49 : C IKAS = 1 KMAB calcul B (Div U)
50 : C IKAS = 2 KMBT calcul Bt (Grad p)
51 : C IKAS = 3 KBBT calcul B assemble B et Bt
52 : C
53 : C TTRAN Transitoire/permanent VRAI/FAUX
54 : C TDFDT on discretise dfdt => TTRAN vrai sinon erreur
55 : C (en permanent pour NS TDFDT faux)
56 : C (en transitoire pour LAPN TDFDT faux)
57 : C TPROJ Projection => TTRAN vrai
58 : C
59 : C
60 : C************************************************************************
61 :
62 : -INC CCOPTIO
63 : -INC CCREEL
64 : -INC CCGEOME
65 : -INC SMCOORD
66 : -INC SMLENTI
67 : POINTEUR MLENT4.MLENTI
68 : -INC SMLMOTS
69 : -INC SMMODEL
70 : -INC SMELEME
71 : POINTEUR IGEOM.MELEME,MELEMS.MELEME,MELEMC.MELEME,MELEMP.MELEME
72 : POINTEUR MELEM2.MELEME
73 : POINTEUR MELEMA.MELEME,MELSTB.MELEME
74 : -INC SMCHPOI
75 : POINTEUR MTETA1.MPOVAL,MTETA2.MPOVAL,MTETA3.MPOVAL,MTETA4.MPOVAL
76 : -INC SMTABLE
77 : POINTEUR MTABZ.MTABLE
78 : -INC SIZFFB
79 : POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR
80 : SEGMENT SAJT
81 : REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM)
82 : REAL*8 TN1(NP,IDIM),TN2(NP,IDIM)
83 : REAL*8 WT(NP,NPG),WS(NP,NPG)
84 : ENDSEGMENT
85 : -INC SMRIGID
86 : C-INC SMMATRIK
87 : C*******************************************************************
88 : C
89 : SEGMENT MATRIK
90 : REAL*8 COEMTK(NMATRI)
91 : INTEGER JRIGEL(NRIGE,NMATRI)
92 : INTEGER KSYM,KMINC,KMINCP,KMINCD,KIZM
93 : INTEGER KISPGT,KISPGP,KISPGD
94 : INTEGER KNTTT,KNTTP,KNTTD
95 : INTEGER KIDMAT(NKID)
96 : INTEGER KKMMT(NKMT)
97 : ENDSEGMENT
98 :
99 : SEGMENT JMATRI
100 : CHARACTER*8 LISPRI(NBMF),LJSDUA(NBMF)
101 : INTEGER LIZAFM(NBSOUS,NBMF)
102 : INTEGER KSPGP,KSPGD
103 : ENDSEGMENT
104 : POINTEUR JMATRS.JMATRI,JMATR1.JMATRI,JMATR2.JMATRI,JMATR3.JMATRI
105 :
106 : C Stokage matrices elementaires non assemblees (valeurs)
107 : SEGMENT IZAFM
108 : REAL*8 AM(NBEL,NP,MP)
109 : ENDSEGMENT
110 : POINTEUR IPM1.IZAFM,IPM2.IZAFM,IPM3.IZAFM,IPM4.IZAFM
111 :
112 : C*******************************************************************
113 : POINTEUR IPM.IZAFM
114 : -INC SMCHAML
115 :
116 : LOGICAL TLAPN,TDFDT,TCONV,TSOUR,TECHI,TMDIA,TMASS,TKBBT,TVNIP
117 : LOGICAL ICAL2,ICAL3,TVITP,TTRAN,TPROJ,TLINCO
118 : LOGICAL XDIAG,XPG,XTV,XTG,XBDF
119 : LOGICAL XRIG,XCONS
120 : CHARACTER*8 CHHH,TYPE,NOM,NOM0,CHAI,SCHT,NOMPER,NOMPR,MPRE
121 : CHARACTER*4 INCOD
122 : CHARACTER*4 NOMI(3),NOMA(3),NOMP(3),NOMD(3)
123 : PARAMETER (NBH=7,NBO=24)
124 : CHARACTER*8 LSCHE(NBH),LSOPT(NBO)
125 : DATA LSOPT /'LINE ','MACRO ','QUAF ','LINB ',
126 : & 'ISOQ ','MMDIAGO ','RGDT ','CONS ','TRANS ',
127 : & 'PERM ','PROJ ','STABP ',' ',' ',
128 : & 'CENTRE ','CENTREP0','CENTREP1','MSOMMET ',' ',
129 : & ' ','CMD ','CENTREE ','SUPGDC ','SUPG '/
130 : DATA LSCHE /'EUL_EXPL','EUL_IMPL','TVISQ ','SEMI ',
131 : & 'CN ','CNG ','BDF2 '/
132 : DIMENSION ICOEF(4),SU(3)
133 : C*****************************************************************************
134 : CTCLSF
135 : NOMPER=NOMPR
136 : write(6,*)' DEBUT YTCLSF appele par ',NOMPER
137 : NBMF=1
138 : IDCEN=1
139 :
140 : C KPOIN = 0->SOMMET 1-> FACE 2-> CENTRE 3-> CENTREP0 4-> CENTREP1 5-> MSOMMET
141 : C KPOIND support dual du schéma d'intégration
142 : KPOIND=0
143 :
144 : TVNIP=.FALSE.
145 : TPROJ=.FALSE.
146 : TTRAN=.FALSE.
147 : TLAPN=.FALSE.
148 : TDFDT=.FALSE.
149 : TCONV=.FALSE.
150 : TSOUR=.FALSE.
151 : TECHI=.FALSE.
152 : TMDIA=.FALSE.
153 : TMASS=.FALSE.
154 : TKBBT=.FALSE.
155 : TVITP=.FALSE.
156 :
157 : L=0
158 : 21 CONTINUE
159 : L=L+1
160 : IF(CHHH(L:L).EQ.'T')TDFDT=.TRUE.
161 : IF(CHHH(L:L).EQ.'L')TLAPN=.TRUE.
162 : IF(CHHH(L:L).EQ.'C')TCONV=.TRUE.
163 : IF(CHHH(L:L).EQ.'S')TSOUR=.TRUE.
164 : IF(CHHH(L:L).EQ.'V')TVNIP=.TRUE.
165 : IF(CHHH(L:L).EQ.'K')THEN
166 : TKBBT=.TRUE.
167 : L=L+1
168 : IKAS=0
169 : IF(CHHH(L:L).EQ.'1')IKAS=1
170 : IF(CHHH(L:L).EQ.'2')IKAS=2
171 : IF(CHHH(L:L).EQ.'3')IKAS=3
172 : ENDIF
173 : C E et M sont exclusifs
174 : IF(CHHH(L:L).EQ.'E')TECHI=.TRUE.
175 : IF(CHHH(L:L).EQ.'M')TMDIA=.TRUE.
176 : IF(L.LT.8)GO TO 21
177 :
178 : IF(TECHI)THEN
179 : IF(TDFDT.OR.TLAPN.OR.TCONV.OR.TSOUR.OR.TMDIA)THEN
180 : C Données incompatibles
181 : WRITE(IOIMP,*)'Operateur : ',NOMPER
181 : CALL ERREUR (21)
182 : RETURN
183 : ENDIF
184 : TMASS=.TRUE.
185 : TSOUR=.TRUE.
186 : ENDIF
187 :
188 : IF(TMDIA)THEN
189 : IF(TDFDT.OR.TLAPN.OR.TCONV.OR.TSOUR.OR.TECHI)THEN
190 : C Données incompatibles
191 : WRITE(IOIMP,*)'Operateur : ',NOMPER
191 : CALL ERREUR (21)
192 : RETURN
193 : ENDIF
194 : TMASS=.TRUE.
195 : ENDIF
196 :
197 : IF(TKBBT.OR.TVNIP)TVITP=.TRUE.
198 : IF(TVNIP)NOMPER='VNIMP'
199 :
200 : C
201 : C Lecture du model et verification OK NAVIER_STOKES
202 : C
203 : TYPE='MMODEL'
203 : CALL LIROBJ (TYPE,MMODEL,1,IRET)
204 :
205 : IF(IRET.EQ.0)THEN
206 : C On ne trouve pas d'objet de type %m1:8
207 : MOTERR( 1: 8) = 'MMODEL '
208 : WRITE(IOIMP,*)'Operateur : ',NOMPER
208 : CALL ERREUR (37)
209 : RETURN
210 : ENDIF
211 :
211 : CALL LEKMOD (MMODEL,MTABZ,INEFMD)
212 : IF(MTABZ.EQ.0)THEN
213 : C% Le modele %m1:8 %m9:16 est obligatoire.
214 : MOTERR( 1: 8) = 'NAVIER_S'
215 : MOTERR( 9:16) = 'TOKES '
216 : WRITE(IOIMP,*)'Operateur : ',NOMPER
216 : CALL ERREUR (935)
217 : RETURN
218 : ENDIF
219 : INEFM0=INEFMD
220 :
221 : C /S INEFMD : Type formulation
222 : C INEFMD = 1 LINE,
223 : C = 2 MACRO,
224 : C = 3 QUADRATIQUE,
225 : C = 4 LINB.
226 : C = 5 ISOQ.
227 : IF(INEFMD.NE.1.AND.INEFMD.NE.2.AND.INEFMD.NE.3.AND.
228 : & INEFMD.NE.4.AND.INEFMD.NE.5)THEN
229 : * Le type d'élément fini ne convient pas
230 : MOTERR( 1: 8) =' '
231 : WRITE(IOIMP,*)'Operateur : ',NOMPER
231 : CALL ERREUR (926)
232 : RETURN
233 : ENDIF
234 : C*************************************************************************
235 : C Lecture des Noms d'Inco(s) et type d'Inco(s)
236 : TYPE='LISTMOTS'
236 : CALL LIROBJ (TYPE,MLMOT3,1,IRET)
237 : IF(IRET.EQ.0)THEN
238 : C On ne trouve pas d'objet de type %m1:8
239 : MOTERR( 1: 8) = 'LISTMOTS'
240 : WRITE(IOIMP,*)'Operateur : ',NOMPER
240 : CALL ERREUR (37)
241 : RETURN
242 : ENDIF
243 :
244 : TYPE='LISTMOTS'
244 : CALL LIROBJ (TYPE,MLMOT4,1,IRET)
245 : IF(IRET.EQ.0)THEN
246 : C On ne trouve pas d'objet de type %m1:8
247 : MOTERR( 1: 8) = 'LISTMOTS'
248 : WRITE(IOIMP,*)'Operateur : ',NOMPER
248 : CALL ERREUR (37)
249 : RETURN
250 : ENDIF
251 :
252 : SEGACT MLMOT3,MLMOT4
253 : NINKO=MLMOT3.MOTS(/2)
254 : c write(6,*)'MLMOT3=',MLMOT3.MOTS(1),MLMOT3.MOTS(/2)
255 : c write(6,*)'MLMOT4=',MLMOT4.MOTS(1),MLMOT4.MOTS(/2)
256 :
257 : IHV=-1
258 : IF(MLMOT4.MOTS(1)(1:4).EQ.'PRES')IHV=0
259 : IF(MLMOT4.MOTS(1)(1:4).EQ.'TEMP')IHV=0
260 : IF(MLMOT4.MOTS(1)(1:4).EQ.'VITE')IHV=1
261 : INCOD=MLMOT4.MOTS(1)
262 : IF(NINKO.EQ.2)INCOD=MLMOT4.MOTS(2)
263 : IF(IHV.EQ.-1)THEN
264 : C Le type d'inconnue %m1:8 ne convient pas.
265 : MOTERR( 1: 8) = MLMOT4.MOTS(1)
266 : WRITE(IOIMP,*)'Operateur : ',NOMPER
266 : CALL ERREUR (927)
267 : RETURN
268 : ENDIF
269 :
270 : IHP=0
271 :
272 : IF(IHV.EQ.0.AND.NINKO.NE.1)THEN
273 : C% Le nombre d'inconnue %i1 ne convient pas.
274 : INTERR(1) = NINKO
275 : WRITE(IOIMP,*)'Operateur : ',NOMPER
275 : CALL ERREUR (928)
276 : RETURN
277 : ENDIF
278 :
279 : IF (TKBBT) THEN
280 : IF (NINKO.NE.2) THEN
281 : C% Le nombre d'inconnues %i1 ne convient pas.
282 : INTERR(1) = NINKO
283 : WRITE(IOIMP,*)'Operateur : ',NOMPER
283 : CALL ERREUR (928)
284 : RETURN
285 : ENDIF
286 : IF ((MLMOT4.MOTS(1).NE.'VITESSE').OR.
287 : & (MLMOT4.MOTS(2).NE.'PRESSION')) THEN
288 : C% Les types d'inconnue ne conviennent pas %m1:8 %m9:16 %m17:24 .
289 : MOTERR( 1: 8) = MLMOT4.MOTS(1)
290 : MOTERR( 9:16) = MLMOT4.MOTS(2)
291 : WRITE(IOIMP,*)'Operateur : ',NOMPER
291 : CALL ERREUR (932)
292 : RETURN
293 : ENDIF
294 : ENDIF
295 :
296 : DO 48 N=1,NINKO
297 : NOMI(N)=MLMOT3.MOTS(N)
298 : NOMA(N)=MLMOT3.MOTS(N)
299 : 48 CONTINUE
300 :
301 : SEGDES MLMOT3,MLMOT4
302 : c write(6,*)' INCOD=',INCOD
303 :
304 : C*****************************************************************************
305 : C OPTIONS
306 :
307 : IKOMP=0
308 : IAXI=0
309 : IF(IFOMOD.EQ.0)IAXI=2
310 : DEUPI=1.D0
311 : IF(IAXI.NE.0)DEUPI=2.D0*XPI
312 :
313 : C Initialisation par defaut de quelques options
314 : CSTAB=1.D-2
315 : MPRE ='????????'
316 : XCONS=.FALSE.
317 : XRIG =.FALSE.
318 : XDIAG=.FALSE.
319 :
320 : 1950 CONTINUE
320 : CALL QUETYP (TYPE,0,IRET)
321 : IF(IRET.EQ.0)THEN
322 : C% Il manque une donnée
322 : CALL ERREUR (641)
323 : RETURN
324 : ENDIF
325 : c write(6,*)' TYPE=',type
326 : IF(TYPE.NE.'MOT ')GO TO 1949
326 : CALL LIRCHA (CHAI,1,LCHAR)
327 : c write(6,*)' OPtion lue ',CHAI
328 :
329 : IF(LCHAR.EQ.0)THEN
330 : C% Directive %m1:8 on ne trouve pas d'objet de type %m9:16 .
331 : MOTERR( 1: 8) = 'OPTION '
332 : MOTERR( 9:16) = 'MOT '
333 : WRITE(IOIMP,*)'Operateur : ',NOMPER
333 : CALL ERREUR (928)
334 : RETURN
335 : ENDIF
336 :
337 : NOM=' '
338 : NOM(1:LCHAR)=CHAI
339 :
339 : CALL OPTLI (IP,LSOPT,NOM,NBO)
340 : IF(IP.EQ.0)THEN
341 : C% On a lu : %m1:8 : , alors qu'on attend un des mots clés suivant :
342 : C% %m9:16 %m17:24 %m25:32 %m33:40 Consulter la notice
343 : MOTERR( 1: 8) = NOM
344 : MOTERR( 9:16) = LSOPT(1)
345 : MOTERR(17:24) = LSOPT(3)
346 : MOTERR(25:32) = LSOPT(5)
347 : MOTERR(33:40) = '... etc '
348 : WRITE(IOIMP,*)'Operateur : ',NOMPER
348 : CALL ERREUR (930)
349 : RETURN
350 : ENDIF
351 :
352 : GO TO (1951,1952,1953,1954,1955,1956,1957,1958,1959,1960,
353 : & 1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,
354 : & 1971,1972,1973,1974),IP
355 :
356 : C LINE
357 : 1951 CONTINUE
358 : INEFMD=1
359 : GO TO 1950
360 :
361 : C MACRO
362 : 1952 CONTINUE
363 : INEFMD=2
364 : GO TO 1950
365 :
366 : C QUAF
367 : 1953 CONTINUE
368 : INEFMD=3
369 : GO TO 1950
370 :
371 : C LINB
372 : 1954 CONTINUE
373 : INEFMD=4
374 : GO TO 1950
375 :
376 : C ISOQ
377 : 1955 CONTINUE
378 : INEFMD=5
379 : GO TO 1950
380 :
381 : C MMDIAGO
382 : 1956 CONTINUE
383 : XDIAG=.TRUE.
384 : GO TO 1950
385 :
386 : C RGDT
387 : 1957 CONTINUE
388 : XRIG=.TRUE.
389 : c write(6,*)' XRIG=',XRIG
390 : GO TO 1950
391 :
392 : C CONS
393 : 1958 CONTINUE
394 : XCONS=.TRUE.
395 : GO TO 1950
396 :
397 : C TRANS
398 : 1959 CONTINUE
399 : TTRAN=.TRUE.
400 : GO TO 1949
401 :
402 : C PERM
403 : 1960 CONTINUE
404 : TTRAN=.FALSE.
405 : TDFDT=.FALSE.
406 : IF(NOMPER.EQ.'DFDT')THEN
407 : C% L'operateur %m1:8 ne peut etre appele dans un algorithme Permanent
408 : MOTERR( 1: 8) = NOMPER
408 : CALL ERREUR (936)
409 : RETURN
410 : ENDIF
411 : GO TO 1949
412 :
413 : C PROJ
414 : 1961 CONTINUE
415 : TPROJ=.TRUE.
416 : TTRAN=.TRUE.
417 : GO TO 1949
418 :
419 : C STABP
420 : 1962 CONTINUE
420 : CALL LIRREE (CSTAB,1,IRET)
421 : IF(IRET.EQ.0)THEN
422 : C On ne trouve pas d'objet de type %m1:8
423 : MOTERR( 1: 8) = 'FLOTTANT'
424 : WRITE(IOIMP,*)'Operateur : ',NOMPER
424 : CALL ERREUR (37)
425 : RETURN
426 : ENDIF
427 : GO TO 1950
428 :
429 : C Emplacements libres
430 : 1963 CONTINUE
431 : 1964 CONTINUE
432 : GO TO 1950
433 :
434 : C CENTRE
435 : 1965 CONTINUE
436 : C CENTREP0
437 : 1966 CONTINUE
438 : C CENTREP1
439 : 1967 CONTINUE
440 : C MSOMMET
441 : 1968 CONTINUE
442 : MPRE=LSOPT(IP)
443 : GO TO 1950
444 :
445 : C Emplacements libres
446 : 1969 CONTINUE
447 : 1970 CONTINUE
448 : GO TO 1950
449 :
450 : C CMD
451 : 1971 CONTINUE
451 : CALL LIRREE (CMD,1,IRET)
452 : IF(IRET.EQ.0)THEN
453 : C On ne trouve pas d'objet de type %m1:8
454 : MOTERR( 1: 8) = 'FLOTTANT'
455 : WRITE(IOIMP,*)'Operateur : ',NOMPER
455 : CALL ERREUR (37)
456 : RETURN
457 : ENDIF
458 : GO TO 1950
459 :
460 : C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG
461 : C CENTREE
462 : 1972 CONTINUE
463 : IDCEN=1
464 : c write(6,*)'LSOPT 1972 ',LSOPT(IP),IP,' IDCEN=',IDCEN
465 : GO TO 1950
466 :
467 : C SUPGDC
468 : 1973 CONTINUE
469 : IDCEN=2
470 : c write(6,*)'LSOPT 1973 ',LSOPT(IP),IP,' IDCEN=',IDCEN
471 : GO TO 1950
472 :
473 : C SUPG
474 : 1974 CONTINUE
475 : IDCEN=3
476 : c write(6,*)'LSOPT 1974 ',LSOPT(IP),IP,' IDCEN=',IDCEN
477 : GO TO 1950
478 :
479 :
480 : C*********** Fin Lecture des options ***********************************
481 :
482 : 1949 CONTINUE
483 :
484 : IF(TKBBT.AND.MPRE.EQ.'????????')THEN
485 : C% L'option %m1:8 n'a pas ete defini, elle est necessaire dans ce contexte.
486 : MOTERR( 1: 8) = 'Pression'
486 : CALL ERREUR (934)
487 : RETURN
488 : ENDIF
489 :
490 :
491 : IF(XRIG)THEN
492 : IF(IHV.EQ.0)THEN
493 : NOMP(1)='T '
494 : NOMD(1)='Q '
495 : ELSEIF(IHV.EQ.1)THEN
496 : NOMP(1)='UX '
497 : NOMP(2)='UY '
498 : NOMP(3)='UZ '
499 : NOMD(1)='FX '
500 : NOMD(2)='FY '
501 : NOMD(3)='FZ '
502 : ENDIF
503 : ELSE
504 : IF(IHV.EQ.0)THEN
505 : NOMP(1)=NOMI(1)
506 : NOMD(1)=NOMI(1)
507 : ELSEIF(IHV.EQ.1)THEN
508 : WRITE(NOMP(1),FMT='(I1,A3)')1,NOMI(1)(1:3)
509 : WRITE(NOMD(1),FMT='(I1,A3)')1,NOMI(1)(1:3)
510 : WRITE(NOMP(2),FMT='(I1,A3)')2,NOMI(1)(1:3)
511 : WRITE(NOMD(2),FMT='(I1,A3)')2,NOMI(1)(1:3)
512 : WRITE(NOMP(3),FMT='(I1,A3)')3,NOMI(1)(1:3)
513 : WRITE(NOMD(3),FMT='(I1,A3)')3,NOMI(1)(1:3)
514 : IF(TVITP)THEN
515 : DO 1948 I=1,IDIM
516 : NOMD(I)=NOMI(2)(1:4)
517 : 1948 CONTINUE
518 : ENDIF
519 : ENDIF
520 : ENDIF
521 : c write(6,*)' NOMP=',nomp
522 : c write(6,*)' NOMD=',nomd
523 :
524 :
525 : c write(6,*)'DECENTREMENT',IDCEN,CMD
526 :
527 : AIMPL=1.D0
528 :
529 : CMT=CMD
530 : XPG=.FALSE.
531 : XTG=.FALSE.
532 : XTV=.FALSE.
533 : XBDF=.FALSE.
534 : C***********************************************************************
535 : C*************** Lecture des caracteristiques du schema en temps *******
536 : C***********************************************************************
537 :
538 : c write(6,*)' TTRAN,TPROJ=',TTRAN,TPROJ
539 : IF(TTRAN)THEN
540 :
541 : C On lit le pas de temps
542 : c write(6,*)'On lit le pas de temps'
542 : CALL LIRREE (DELTAT,1,IRET)
543 : IF(IRET.EQ.0)THEN
544 : C On ne trouve pas d'objet de type %m1:8
545 : MOTERR( 1: 8) = 'FLOTTANT'
546 : WRITE(IOIMP,*)'Operateur : ',NOMPER
546 : CALL ERREUR (37)
547 : RETURN
548 : ENDIF
549 :
550 : C On lit le schema en temps
551 : c write(6,*)'On lit le schema en temps '
551 : CALL LIRCHA (SCHT,1,LCHAR)
552 : IF(LCHAR.EQ.0)THEN
553 : C On ne trouve pas d'objet de type %m1:8
554 : MOTERR( 1: 8) = 'MOT '
555 : WRITE(IOIMP,*)'Operateur : ',NOMPER
555 : CALL ERREUR (37)
556 : RETURN
557 : ENDIF
558 : c write(6,*)'SCHEMA en TEMPS ',SCHT
558 : CALL OPTLI (IPST,LSCHE,SCHT,NBH)
559 : IF(IPST.EQ.0)THEN
560 : C% On a lu : %m1:8 : , alors qu'on attend un des mots clés suivant :
561 : C% %m9:16 %m17:24 %m25:32 %m33:40 Consulter la notice
562 : MOTERR( 1: 8) = 'SCHT '
563 : MOTERR( 9:16) = LSCHE(2)
564 : MOTERR(17:24) = LSCHE(4)
565 : MOTERR(25:32) = LSCHE(7)
566 : MOTERR(33:40) = '... etc '
567 : WRITE(IOIMP,*)'Operateur : ',NOMPER
567 : CALL ERREUR (930)
568 : RETURN
569 : ENDIF
570 :
571 : IF(SCHT.EQ.'SEMI')THEN
571 : CALL LIRREE (XVAL,1,IRET)
572 : IF(IRET.EQ.0)THEN
573 : C On ne trouve pas d'objet de type %m1:8
574 : MOTERR( 1: 8) = 'FLOTTANT'
575 : WRITE(IOIMP,*)'Operateur : ',NOMPER
575 : CALL ERREUR (37)
576 : RETURN
577 : ENDIF
578 : AIMPL=XVAL
579 : ENDIF
580 :
581 : IF(SCHT.EQ.'EUL_EXPL')AIMPL=0.D0
582 : IF(SCHT.EQ.'EUL_IMPL')AIMPL=1.D0
583 : IF(SCHT.EQ.'CN ')AIMPL=0.5D0
584 :
585 : IF(SCHT.EQ.'EUL_EXPL'.OR.
586 : & SCHT.EQ.'EUL_IMPL'.OR.
587 : & SCHT.EQ.'CN '.OR.
588 : & SCHT.EQ.'BDF2 '.OR.
589 : & SCHT.EQ.'SEMI ')XPG=.TRUE.
590 : IF(.NOT.TCONV)XPG=.FALSE.
591 :
592 : IF(XPG.AND.IDCEN.NE.1.AND.IDCEN.NE.2.AND.IDCEN.NE.3)THEN
593 : C% L'option PETROV-GALERKIN est incompatible ce schema en temps
593 : CALL ERREUR (937)
594 : RETURN
595 : ENDIF
596 :
597 : IF(.NOT.TCONV)IDCEN=1
598 :
599 : IF(SCHT.EQ.'BDF2 ')THEN
600 : AIMPL=1.D0
601 : XBDF=.TRUE.
602 : CMT=CMD
603 : ENDIF
604 :
605 : IF(SCHT.EQ.'CNG ')THEN
606 : AIMPL=0.5D0
607 : XTG=.TRUE.
608 : CMT=DELTAT
609 : ENDIF
610 :
611 : IF(SCHT.EQ.'TVISQ')THEN
612 : AIMPL=0.5D0
613 : XTV=.TRUE.
614 : CMT=DELTAT
615 : ENDIF
616 : ENDIF
617 :
618 : C*********** Fin Lecture des caracteristiques du schema en temps *******
619 :
620 : c write(6,*)'TCLSF appele par ',NOMPER,' AIMPL=',AIMPL,' IHV=',IHV
621 : c write(6,*)'XBDF =',XBDF
622 : c write(6,*)' TLAPN=',TLAPN,' TTRAN=',TTRAN,' TCONV=',TCONV
623 : c write(6,*)' TSOUR=',TSOUR,' TMASS=',TMASS,' TECHI=',TECHI
624 : c write(6,*)' TMDIA=',TMDIA,' XDIAG=',XDIAG,' TKBBT=',TKBBT
625 :
626 : C*****************************************************************************
627 : C
627 : CALL MEKTAB (MTABZ,INEFMD,'MAILLAGE',MELEME)
628 : IF(TKBBT.AND.(MPRE.NE.'CENTRE').AND.(INEFMD.EQ.2))THEN
628 : CALL MEKTAB (MTABZ,INEFMD,'MACRO1',MELEME)
629 : ENDIF
630 :
630 : CALL MEKTAB (MTABZ,INEFMD,'SOMMET',MELEMS)
631 : IF(MPRE.NE.'????????')THEN
631 : CALL MEKTAB (MTABZ,INEFMD,MPRE,MELEMP)
632 : ENDIF
632 : CALL MEKTAB (MTABZ,INEFMD,'CENTRE',MELEMC)
633 : IF (IERR.NE.0) RETURN
634 :
635 : c write(6,*)'INCOD=',INCOD,' MPRE=',MPRE
636 : IF(INCOD.EQ.'VITE'.OR.INCOD.EQ.'TEMP')THEN
637 : MELEM2=MELEME
638 : ELSE
638 : IF(MPRE.EQ.'CENTREP1')CALL MEKTAB (MTABZ,INEFMD,'ELTP1NC ',MELEM2)
638 : IF(MPRE.EQ.'MSOMMET ')CALL MEKTAB (MTABZ,INEFMD,'MMAIL ',MELEM2)
638 : IF(MPRE.EQ.'CENTRE ')CALL MEKTAB (MTABZ,INEFMD,'CENTRE ',MELEM2)
638 : IF(MPRE.EQ.'CENTREP0')CALL MEKTAB (MTABZ,INEFMD,'CENTREP0',MELEM2)
639 :
640 : C /S INEFMD : Type formulation
641 : C INEFMD = 1 LINE,
642 : C = 2 MACRO,
643 : C = 3 QUADRATIQUE,
644 : C = 4 LINB,
645 : C = 5 ISOQ.
646 : NOM0='????????'
647 : IF(INEFMD.EQ.1)THEN
648 : IF(MPRE.EQ.'CENTRE ')NOM0='OK'
649 : IF(MPRE.EQ.'MSOMMET ')NOM0='OK'
650 : ELSEIF(INEFMD.EQ.2)THEN
651 : IF(MPRE.EQ.'CENTRE ')NOM0='OK'
652 : c ces elements ne marchent pas
653 : IF(MPRE.EQ.'CENTREP0')NOM0='OK'
654 : IF(MPRE.EQ.'CENTREP1')NOM0='OK'
655 : IF(MPRE.EQ.'MSOMMET ')NOM0='OK'
656 : ELSEIF(INEFMD.EQ.3)THEN
657 : IF(MPRE.EQ.'CENTREP0')NOM0='OK'
658 : IF(MPRE.EQ.'CENTREP1')NOM0='OK'
659 : IF(MPRE.EQ.'MSOMMET ')NOM0='OK'
660 : ELSEIF(INEFMD.EQ.4)THEN
661 : IF(MPRE.EQ.'MSOMMET ')NOM0='OK'
662 : ENDIF
663 : IF(NOM0.EQ.'????????')THEN
664 : C% Le type d'element fini Vitesse/pression ne convient pas : %m1:8 %m9:16 .
665 : MOTERR( 1: 8) =LSOPT(INEFMD)
666 : MOTERR( 9:16) =MPRE
667 : WRITE(IOIMP,*)'Operateur : ',NOMPER
667 : CALL ERREUR (933)
668 : RETURN
669 : ENDIF
670 :
671 : ENDIF
672 :
673 : C*************************************************************************
674 :
675 : SEGACT MELEME,MELEM2
676 : c write(6,*)'MELEME=',MELEME,'MELEM2=',MELEM2
677 : SEGDES MTABZ
678 :
679 : C*************************************************************************
680 : C Lecture des Inco(s) aux temps precedents si Transitoire
681 : C Lecture des Inco(s) a l'iteration precedente si Permanent
682 : c write(6,*)' Lecture des Inco(s)'
683 :
683 : CALL QUETYP (TYPE,0,IRET)
684 : IF(IRET.EQ.0)THEN
685 : C% Il manque une donnée
685 : CALL ERREUR (641)
686 : RETURN
687 : ENDIF
688 :
689 : TLINCO=.FALSE.
690 : IF(IDCEN.EQ.2)TLINCO=.TRUE.
691 : IF(TYPE.EQ.'MOT ')THEN
692 : TLINCO=.TRUE.
692 : CALL LIRCHA (CHAI,1,LCHAR)
693 : IF(CHAI.NE.'INCO') THEN
694 : C% On a lu : %m1:8 : , alors qu'on attend un des mots clés suivant :
695 : C% %m9:16 %m17:24 %m25:32 %m33:40 Consulter la notice
696 : MOTERR( 1: 8) = CHAI
697 : MOTERR( 9:16) = 'INCO'
698 : MOTERR(17:24) = ' '
699 : MOTERR(25:32) = ' '
700 : MOTERR(33:40) = ' '
701 : WRITE(IOIMP,*)'Operateur : ',NOMPER
701 : CALL ERREUR (930)
702 : RETURN
703 : ENDIF
704 : ENDIF
705 : IF(.NOT.TLINCO.AND.TTRAN)THEN
706 : C% Opération illicite dans ce contexte
707 : WRITE(IOIMP,*)'Operateur : ',NOMPER,' On attend des inconnues'
707 : CALL ERREUR (153)
708 : RETURN
709 : ENDIF
710 :
711 : C 1ere inconnue
712 : IF(NINKO.GE.1.AND.TLINCO)THEN
713 : TYPE='CHPOINT'
713 : CALL LIROBJ (TYPE,MCHPOI,1,IRET)
714 :
715 : IF(IRET.EQ.0)THEN
716 : C% Rubrique %m1:8 on ne trouve pas d'objet de type %m9:16 .
717 : MOTERR( 1: 8) = 'INCO '
718 : MOTERR( 9:16) = 'CHPOINT '
719 : WRITE(IOIMP,*)'Operateur : ',NOMPER
719 : CALL ERREUR (929)
720 : RETURN
721 : ENDIF
722 :
722 : CALL LRCHT (MCHPOI,MTETA1,TYPE,IGEOM)
723 : MTETA2=MTETA1
723 : CALL KRIPAD (IGEOM,MLENT1)
723 : CALL VERPAD (MLENT1,MELEMS,IRET)
724 : IF(IRET.NE.0)THEN
725 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
726 : MOTERR(1: 8) = 'TETA1'
727 : MOTERR(9:16) = 'CHPOINT '
728 : WRITE(IOIMP,*)'Operateur : ',NOMPER
728 : CALL ERREUR (788)
729 : RETURN
730 : ENDIF
731 :
732 : NC=MTETA1.VPOCHA(/2)
733 :
734 : IF((IHV.EQ.0.AND.NC.NE.1).OR.
735 : & (IHV.EQ.1.AND.NC.NE.IDIM))THEN
736 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
737 : MOTERR( 1: 8) = 'Inconnue'
738 : MOTERR( 9:16) = 'CHPOINT'
739 : WRITE(IOIMP,*)'Operateur : ',NOMPER
739 : CALL ERREUR (784)
740 : RETURN
741 : ENDIF
742 :
743 : SEGDES IGEOM,MELEMS
744 :
745 : IF(XBDF)THEN
746 : TYPE='CHPOINT'
746 : CALL LIROBJ (TYPE,MCHPOI,1,IRET)
747 :
748 : IF(IRET.EQ.0)THEN
749 : C% On ne trouve pas d'objet de type %m1:8
750 : MOTERR( 1: 8) = 'CHPOINT '
751 : WRITE(IOIMP,*)'Operateur : ',NOMPER
751 : CALL ERREUR (37)
752 : RETURN
753 : ENDIF
754 :
754 : CALL LRCHT (MCHPOI,MTETA2,TYPE,IGEOM)
754 : CALL KRIPAD (IGEOM,MLENT2)
754 : CALL VERPAD (MLENT2,MELEMS,IRET)
755 : IF(IRET.NE.0)THEN
756 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
757 : MOTERR(1: 8) = 'TETA2'
758 : MOTERR(9:16) = 'CHPOINT '
759 : WRITE(IOIMP,*)'Operateur : ',NOMPER
759 : CALL ERREUR (788)
760 : RETURN
761 : ENDIF
762 :
763 : NC=MTETA2.VPOCHA(/2)
764 :
765 : IF((IHV.EQ.0.AND.NC.NE.1).OR.
766 : & (IHV.EQ.1.AND.NC.NE.IDIM))THEN
767 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
768 : MOTERR( 1: 8) = 'Inconnue'
769 : MOTERR( 9:16) = 'CHPOINT'
770 : WRITE(IOIMP,*)'Operateur : ',NOMPER
770 : CALL ERREUR (785)
771 : RETURN
772 : ENDIF
773 :
774 : SEGDES IGEOM,MELEMS
775 :
776 : ENDIF
777 : ENDIF
778 :
779 : C 2eme inconnue
780 : IF(NINKO.GE.2.AND.TLINCO)THEN
781 : TYPE='CHPOINT'
781 : CALL LIROBJ (TYPE,MCHPOI,1,IRET)
782 :
783 : IF(IRET.EQ.0)THEN
784 : C% Rubrique %m1:8 on ne trouve pas d'objet de type %m9:16 .
785 : MOTERR( 1: 8) = 'INCO2 '
786 : MOTERR( 9:16) = 'CHPOINT '
787 : WRITE(IOIMP,*)'Operateur : ',NOMPER
787 : CALL ERREUR (929)
788 : RETURN
789 : ENDIF
790 :
790 : CALL LRCHT (MCHPOI,MTETA3,TYPE,IGEOM)
791 : MTETA4=MTETA3
791 : CALL KRIPAD (IGEOM,MLENT3)
792 : IRET=0
792 : IF(.NOT.TVNIP)CALL VERPAD (MLENT3,MELEMP,IRET)
793 : IF(IRET.NE.0)THEN
794 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
795 : MOTERR(1: 8) = 'TETA3'
796 : MOTERR(9:16) = 'CHPOINT '
797 : WRITE(IOIMP,*)'Operateur : ',NOMPER
797 : CALL ERREUR (788)
798 : RETURN
799 : ENDIF
800 :
801 : NC=MTETA3.VPOCHA(/2)
802 :
803 : IF(IHP.EQ.0.AND.NC.NE.1)THEN
804 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
805 : MOTERR( 1: 8) = 'Inconnue'
806 : MOTERR( 9:16) = 'CHPOINT'
807 : WRITE(IOIMP,*)'Operateur : ',NOMPER
807 : CALL ERREUR (784)
808 : RETURN
809 : ENDIF
810 :
811 : SEGDES IGEOM,MELEMP
812 :
813 : IF(XBDF)THEN
814 : TYPE='CHPOINT'
814 : CALL LIROBJ (TYPE,MCHPOI,1,IRET)
815 :
816 : IF(IRET.EQ.0)THEN
817 : C% On ne trouve pas d'objet de type %m1:8
818 : MOTERR( 1: 8) = 'CHPOINT '
819 : WRITE(IOIMP,*)'Operateur : ',NOMPER
819 : CALL ERREUR (37)
820 : RETURN
821 : ENDIF
822 :
822 : CALL LRCHT (MCHPOI,MTETA4,TYPE,IGEOM)
822 : CALL KRIPAD (IGEOM,MLENT4)
823 : IRET=0
823 : IF(.NOT.TVNIP)CALL VERPAD (MLENT4,MELEMP,IRET)
824 : IF(IRET.NE.0)THEN
825 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
826 : MOTERR(1: 8) = 'TETA4'
827 : MOTERR(9:16) = 'CHPOINT '
828 : WRITE(IOIMP,*)'Operateur : ',NOMPER
828 : CALL ERREUR (788)
829 : RETURN
830 : ENDIF
831 :
832 : NC=MTETA4.VPOCHA(/2)
833 :
834 : IF(IHP.EQ.0.AND.NC.NE.1)THEN
835 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
836 : MOTERR( 1: 8) = 'Inconnue'
837 : MOTERR( 9:16) = 'CHPOINT'
838 : WRITE(IOIMP,*)'Operateur : ',NOMPER
838 : CALL ERREUR (785)
839 : RETURN
840 : ENDIF
841 :
842 : SEGDES IGEOM,MELEMP
843 :
844 : ENDIF
845 : ENDIF
846 : C On cree un second membre non vide
847 : SEGACT MELEMS
848 : N=MELEMS.NUM(/2)
849 : IF(IHV.EQ.0)NC=1
850 : IF(IHV.EQ.1)NC=IDIM
851 : IF(TVITP)NC=1
852 : SEGDES MELEMS
853 : NSOUPO=1
854 : NAT=2
855 : SEGINI MCHPO1,MSOUP1,MPOVA1
856 : MCHPO1.JATTRI(1)=2
857 : MCHPO1.IFOPOI=IFOMOD
858 : MCHPO1.MTYPOI=' '
859 : MCHPO1.MOCHDE=' '
860 : MCHPO1.IPCHP(1)=MSOUP1
861 : MSOUP1.IGEOC=MELEMS
862 : MSOUP1.IPOVAL=MPOVA1
863 : IF(IHV.EQ.0.OR.TVITP)THEN
864 : MSOUP1.NOCOMP(1)=NOMD(1)
865 : ELSE
866 : DO 91 N=1,NC
867 : MSOUP1.NOCOMP(N)=NOMD(N)
868 : 91 CONTINUE
869 : ENDIF
870 : SEGDES MSOUP1,MCHPO1
871 :
872 : C*****************************************************************************
873 : C Lecture du coefficient
874 :
874 : CALL LIRENT (IARG,1,IRET)
875 :
876 : IF(IRET.EQ.0)THEN
877 : WRITE(IOIMP,*)'Operateur : ',NOMPER
877 : CALL ERREUR (766)
878 : RETURN
879 : ENDIF
880 : c write(6,*)' IARG=',iarg,' TDFDT=',TDFDT
881 :
881 : CALL INITI (ICOEF,4,0)
882 : IF(TDFDT)ICOEF(1)=1
883 : IF(TMASS)ICOEF(1)=1
884 : IF(TCONV)THEN
885 : ICOEF(1)=1
886 : ICOEF(2)=1
887 : ENDIF
888 : IF(TLAPN)ICOEF(3)=1
889 : IF(TSOUR)ICOEF(4)=1
890 : IF(TKBBT)ICOEF(1)=1
891 : ICT=ICOEF(1)+ICOEF(2)+ICOEF(3)+ICOEF(4)
892 : IF(ICT.NE.IARG)THEN
893 : C% Données incompatibles
894 : WRITE(6,*)'Nombre d arguments incorrect'
895 : WRITE(IOIMP,*)'Operateur : ',NOMPER
895 : CALL ERREUR (22)
896 : RETURN
897 : ENDIF
898 :
899 : NUCOEF=0
900 :
901 : IF(TDFDT.OR.TMASS.OR.TCONV)THEN
902 : NUCOEF=NUCOEF+1
902 : CALL LEKMIF (NUCOEF,MTABZ,0,MCHEL1,KPOIND)
903 : IF (IERR.NE.0) RETURN
904 : ENDIF
905 :
906 : c write(6,*)' Avt TLAPN LEKMIF '
907 : IF(TLAPN)THEN
908 : NUCOEF=NUCOEF+1
908 : CALL LEKMIF (NUCOEF,MTABZ,0,MCHEL2,KPOIND)
909 : IF (IERR.NE.0) RETURN
910 : ELSEIF(TCONV)THEN
911 : c write(6,*)'Cas decentrement'
912 : SEGACT MELEME
913 : L1=72
914 : N1=MAX(1,LISOUS(/1))
915 : N2=1
916 : N3=6
917 : SEGINI MCHEL2
918 : SEGACT MCHEL1
919 : DO 71 L=1,N1
920 : MCHAM1=MCHEL1.ICHAML(L)
921 : SEGACT MCHAM1
922 : MELVA1=MCHAM1.IELVAL(1)
923 : SEGACT MELVA1
924 : N1PTEL=MELVA1.VELCHE(/1)
925 : N1EL=1
926 : N2PTEL=0
927 : N2EL=0
928 : c fauteDO 71 L=1,N1
929 : SEGINI MCHAM2,MELVA2
930 : MCHEL2.ICHAML(L)=MCHAM2
931 : MCHAM2.IELVAL(1)=MELVA2
932 : DO 72 LG=1,N1PTEL
933 : MELVA2.VELCHE(LG,1)=1.D-30
934 : 72 CONTINUE
935 : 71 CONTINUE
936 : ENDIF
937 :
938 : c write(6,*)' Avt TCONV LEKMIF XPG,XTV,XTG,N.TDFDT=',
939 : c & XPG,XTV,XTG,(.NOT.TDFDT)
940 : IF(TCONV)THEN
941 : c? IF(XPG)THEN
942 : NUCOEF=NUCOEF+1
942 : CALL LEKMIF (NUCOEF,MTABZ,1,MCHEL3,KPOIND)
943 : IF (IERR.NE.0) RETURN
944 : c? ELSEIF(XTV)THEN
945 : c? NUCOEF=NUCOEF+1
945 : c? CALL LEKMIF (NUCOEF,MTABZ,1,MCHEL3,KPOIND)
946 : c? & 4 ,DELTAT,MCHEL1,MCHEL1)
947 : c? IF (IERR.NE.0) RETURN
948 : c? ELSEIF(XTG)THEN
949 : c? NUCOEF=NUCOEF+1
949 : c? CALL LEKMIF (NUCOEF,MTABZ,1,MCHEL3,KPOIND)
950 : c? & 5 ,DELTAT,MCHEL1,MCHEL1)
951 : c? IF (IERR.NE.0) RETURN
952 : C ELSEIF(.NOT.TDFDT)THEN
953 : c? ELSE
954 : c write(6,*)' On est bien ds ce cas'
955 : c? NUCOEF=NUCOEF+1
955 : c? CALL LEKMIF (NUCOEF,MTABZ,1,MCHEL3,KPOIND)
956 : c? IF (IERR.NE.0) RETURN
957 : c? ENDIF
958 : ENDIF
959 :
960 : c write(6,*)' Avt LEKMIF TSOUR'
961 : IF(TSOUR)THEN
962 : NUCOEF=NUCOEF+1
962 : CALL LEKMIF (NUCOEF,MTABZ,IHV,MCHEL4,KPOIND)
963 : c write(6,*)' Apr LEKMIF TSOUR'
964 : IF (IERR.NE.0) RETURN
965 : ENDIF
966 :
967 : IF(TKBBT.OR.TVNIP)THEN
968 : NUCOEF=NUCOEF+1
968 : CALL LEKMIF (NUCOEF,MTABZ,0,MCHEL1,KPOIND)
969 : c write(6,*)' meleme,melem2=',meleme,melem2
970 : IF (IERR.NE.0) RETURN
971 : ENDIF
972 :
973 : C*************************************************************************
974 : C******* CALCUL **********************************************************
975 : C*************************************************************************
976 :
977 : C Initialisation RIGIDITE Vide
978 : IF(XRIG)THEN
979 : NRIGE=8
980 : NRIGEL=0
981 : SEGINI MRIGID
982 :
983 : ELSE
984 : C Initialisation MATRIK Vide
985 : NRIGE=7
986 : NKID =9
987 : NKMT =7
988 : NMATRI=0
989 : SEGINI MATRIK
990 : ENDIF
991 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
992 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
993 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
994 : IF(.NOT.TVITP)THEN
995 :
995 : CALL KRIPAD (MELEMS,MLENTI)
996 : SEGDES MELEMS
997 :
998 : NUTOEL=0
999 : SEGACT MELEME,MELEM2
1000 : NBSOUS=LISOUS(/1)
1001 : IF(NBSOUS.EQ.0)NBSOUS=1
1002 :
1003 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN)THEN
1004 : C Cas RIGIDITE
1005 : IF(XRIG)THEN
1006 : NRIGE=8
1007 : NRI=IRIGEL(/2)
1008 : NRIGEL=NBSOUS+NRI
1009 : SEGADJ MRIGID
1010 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
1011 :
1012 : ELSE
1013 : C Cas MATRIK Pleine
1014 : NRIGE=7
1015 : NKID =9
1016 : NKMT =7
1017 : NMATR0=JRIGEL(/2)
1018 :
1019 : NBME=1
1020 : IF(IHV.EQ.1)NBME=IDIM
1021 : NMATRI=NMATR0+NBME
1022 : SEGADJ MATRIK
1023 :
1024 : DO 41 M=1,NBME
1025 : JRIGEL(1,NMATR0+M)=MELEME
1026 : JRIGEL(2,NMATR0+M)=MELEM2
1027 : JRIGEL(7,NMATR0+M)=0
1028 : IF(TCONV)JRIGEL(7,NMATR0+M)=2
1029 : SEGINI JMATRI
1030 : JRIGEL(4,NMATR0+M)=JMATRI
1031 : KSPGP=MELEMS
1032 : KSPGD=MELEMS
1033 : LISPRI(1)=NOMP(M)//' '
1034 : LJSDUA(1)=NOMD(M)//' '
1035 : 41 CONTINUE
1036 : ENDIF
1037 : ENDIF
1038 :
1039 : IF(TDFDT.OR.TMASS.OR.TCONV)SEGACT MCHEL1
1040 : IF(TLAPN)SEGACT MCHEL2
1041 : IF(TCONV)SEGACT MCHEL3
1042 : IF(TSOUR)SEGACT MCHEL4
1043 : IF(TLINCO)THEN
1044 : SEGACT MTETA1
1045 : IF(XBDF)SEGACT MTETA2
1046 : ENDIF
1047 :
1048 : IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN
1049 : WRITE(6,*)' Geometries incompatibles dans ',nomper
1050 : C% Données incompatibles
1050 : CALL ERREUR (22)
1051 : RETURN
1052 : ENDIF
1053 :
1054 : DO 101 L=1,MAX(1,LISOUS(/1))
1055 : IPT1=MELEME
1056 : IPT2=MELEM2
1057 : IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
1058 : IF(MELEM2.LISOUS(/1).NE.0)IPT2=MELEM2.LISOUS(L)
1059 : SEGACT IPT1,IPT2
1060 :
1061 : NOM0 = NOMS(IPT1.ITYPEL)//' '
1062 : c write(6,*)' 1er KALPBG NOM0=',NOM0,IPT1
1062 : CALL KALPBG (NOM0,'FONFORM ',IZFFM)
1063 : IF(IZFFM.EQ.0)RETURN
1064 : SEGACT IZFFM*MOD
1065 : IZHR=KZHR(1)
1066 : SEGACT IZHR*MOD
1067 : IZF1 = KTP(1)
1068 : IZH2 = KZHR(2)
1069 : IZW = IZFFM
1070 : IZWH = IZHR
1071 :
1072 : IF(INCOD.EQ.'PRES')THEN
1073 : IZW = IZF1
1074 : IZWH = IZH2
1075 : IF(IPT2.NUM(/1).NE.IZF1.FN(/1))THEN
1076 : * Le type d'élément fini ne convient pas
1077 : MOTERR( 1: 8) =' '
1078 : WRITE(IOIMP,*)'Operateur : ',NOMPER
1078 : CALL ERREUR (926)
1079 : RETURN
1080 : ENDIF
1081 : ENDIF
1082 :
1083 : NES=GR(/1)
1084 : NPG=GR(/3)
1085 :
1086 : NP = IPT1.NUM(/1)
1087 : MP = IPT2.NUM(/1)
1088 : C? MP = IZW.FN(/1) ceci doit etre identique
1089 :
1090 : NBEL=IPT1.NUM(/2)
1091 :
1092 : SEGINI SAJT
1093 :
1094 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN)THEN
1095 : C Cas RIGIDITE
1096 : IF(XRIG)THEN
1097 : IRIGEL(1,NRI+L)=MELEME
1098 : COERIG(L)=1.D0
1099 :
1100 : IRIGEL(7,NRI+L)=0
1101 : IF(TCONV)IRIGEL(7,NRI+L)=2
1102 :
1103 : NBME=1
1104 : IF(IHV.EQ.1)NBME=IDIM
1105 : NLIGRP=NP
1106 : NLIGRD=MP
1107 : SEGINI DESCR
1108 : IRIGEL(3,NRI+L)=DESCR
1109 : IF(NBME.EQ.1)THEN
1110 : DO 102 I=1,NLIGRP
1111 : LISINC(I)=NOMP(1)//' '
1112 : NOELEP(I)=I
1113 : 102 CONTINUE
1114 : DO 103 I=1,NLIGRD
1115 : LISDUA(I)=NOMD(1)//' '
1116 : NOELED(I)=I
1117 : 103 CONTINUE
1118 : ELSE
1119 : ENDIF
1120 : SEGDES DESCR
1121 :
1122 : NELRIG=NBEL
1123 : SEGINI xMATRI
1124 : IRIGEL(4,NRI+L)=xMATRI
1125 : c write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI
1126 :
1127 : * DO 104 K=1,NELRIG
1128 : * SEGINI XMATRI
1129 : c write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI
1130 : * IMATTT(K)=XMATRI
1131 : * 104 CONTINUE
1132 :
1133 : ELSE
1134 : C Cas MATRIK
1135 : SEGINI IPM1
1136 : JMATR1=JRIGEL(4,NMATR0+1)
1137 : JMATR1.LIZAFM(L,1)=IPM1
1138 : IPM2=IPM1
1139 : IPM3=IPM1
1140 : IF(NBME.GE.2)THEN
1141 : JMATR2=JRIGEL(4,NMATR0+2)
1142 :
1143 : IF(IAXI.NE.0)THEN
1144 : SEGINI IPM2
1145 : JMATR2.LIZAFM(L,1)=IPM2
1146 : ICAL2=.TRUE.
1147 : ELSE
1148 : IPM2=IPM1
1149 : JMATR2.LIZAFM(L,1)=IPM2
1150 : ICAL2=.FALSE.
1151 : ENDIF
1152 :
1153 : ENDIF
1154 : IF(NBME.GE.3)THEN
1155 : JMATR3=JRIGEL(4,NMATR0+3)
1156 :
1157 : IPM3=IPM1
1158 : JMATR3.LIZAFM(L,1)=IPM3
1159 : ICAL3=.FALSE.
1160 :
1161 : ENDIF
1162 : ENDIF
1163 : ENDIF
1164 :
1165 : C----Ro
1166 : IK1=1
1167 : IF(TDFDT.OR.TMASS.OR.TCONV)THEN
1168 : MCHAM1=MCHEL1.ICHAML(L)
1169 : SEGACT MCHAM1
1170 :
1171 : MELVA1=MCHAM1.IELVAL(1)
1172 : SEGACT MELVA1
1173 : N1PTEL=MELVA1.VELCHE(/1)
1174 : N1EL=MELVA1.VELCHE(/2)
1175 : IF(N1EL.EQ.1)THEN
1176 : IK1=1
1177 : ELSEIF(N1EL.EQ.NBEL)THEN
1178 : IK1=0
1179 : ENDIF
1180 : ENDIF
1181 :
1182 : C----Lambda
1183 : IK2=1
1184 : IF(TLAPN)THEN
1185 : MCHAM2=MCHEL2.ICHAML(L)
1186 : SEGACT MCHAM2
1187 :
1188 : MELVA2=MCHAM2.IELVAL(1)
1189 : SEGACT MELVA2
1190 : N1PTEL=MELVA2.VELCHE(/1)
1191 : N1EL=MELVA2.VELCHE(/2)
1192 : IF(N1EL.EQ.1)THEN
1193 : IK2=1
1194 : ELSEIF(N1EL.EQ.NBEL)THEN
1195 : IK2=0
1196 : ENDIF
1197 : ENDIF
1198 :
1199 : C----U
1200 : IK3=1
1201 : IF(TCONV)THEN
1202 : MCHAM3=MCHEL3.ICHAML(L)
1203 : SEGACT MCHAM3
1204 :
1205 : MELVA3=MCHAM3.IELVAL(1)
1206 : SEGACT MELVA3
1207 : N1PTEL=MELVA3.VELCHE(/1)
1208 : N1EL=MELVA3.VELCHE(/2)
1209 : IF(N1EL.EQ.1)THEN
1210 : IK3=1
1211 : ELSEIF(N1EL.EQ.NBEL)THEN
1212 : IK3=0
1213 : ENDIF
1214 : ENDIF
1215 :
1216 : C----source
1217 : IK4=1
1218 : IF(TSOUR)THEN
1219 : MCHAM4=MCHEL4.ICHAML(L)
1220 : SEGACT MCHAM4
1221 :
1222 : MELVA4=MCHAM4.IELVAL(1)
1223 : SEGACT MELVA4
1224 : N1PTEL=MELVA4.VELCHE(/1)
1225 : N1EL=MELVA4.VELCHE(/2)
1226 : IF(N1EL.EQ.1)THEN
1227 : IK4=1
1228 : ELSEIF(N1EL.EQ.NBEL)THEN
1229 : IK4=0
1230 : ENDIF
1231 : ENDIF
1232 :
1233 : c write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
1234 : C===============================================
1235 : DO 108 KE=1,NBEL
1236 :
1237 : NK1=KE + IK1*(1 - KE)
1238 : NK2=KE + IK2*(1 - KE)
1239 : NK3=KE + IK3*(1 - KE)
1240 : NK4=KE + IK4*(1 - KE)
1241 :
1242 : DO 109 I=1,NP
1243 : J=IPT1.NUM(I,KE)
1244 : DO 109 N=1,IDIM
1245 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
1246 : 109 CONTINUE
1247 :
1248 : IF(TLINCO)THEN
1249 : DO 107 I=1,NP
1250 : I1=MLENT1.LECT(IPT1.NUM(I,KE))
1251 : DO 107 N=1,NC
1252 : TN1(I,N)=MTETA1.VPOCHA(I1,N)
1253 : 107 CONTINUE
1254 : ENDIF
1255 :
1256 : IF(XBDF)THEN
1257 : DO 106 I=1,NP
1258 : I2=MLENT2.LECT(IPT1.NUM(I,KE))
1259 : DO 106 N=1,NC
1260 : TN2(I,N)=MTETA2.VPOCHA(I2,N)
1261 : 106 CONTINUE
1262 : ENDIF
1263 :
1263 : CALL CALJBR (FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
1264 : * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
1265 :
1266 : IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN
1266 : CALL RSETD (WT,FN,(NP*NPG))
1266 : CALL RSETD (WS,FN,(NP*NPG))
1267 : ELSE
1267 : CALL CALDEC (WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT,
1268 : & IDCEN,CMT,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2),
1269 : & MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,AIRE,KE)
1270 : ENDIF
1271 :
1271 : CALL INITD (RF1,(NP*MP*IDIM),0.D0)
1271 : CALL INITD (SM1,(MP*IDIM),0.D0)
1272 :
1273 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1274 : IF(TDFDT) THEN
1275 : C...... Terme Transitoire
1276 :
1277 : c? write(6,*)' IDCEN TRAN=',IDCEN,CMD,DELTAT,' NC=',nc
1278 : c? write(6,*)' XBDF=',XBDF,' XDIAG=',XDIAG,' MP=',mp
1279 : DO 110 I=1,MP
1279 : CALL INITD (SU,NC,0.D0)
1280 : UD=0.D0
1281 :
1282 : DO 111 J=1,NP
1283 : U1=0.D0
1284 : TV=0.D0
1285 : DO 112 LG=1,NPG
1286 :
1287 : C1=MELVA1.VELCHE(LG,NK1)
1288 :
1289 : c? DO 114 N=1,IDIM
1290 : c? IF(XTV)THEN
1291 : c? C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)
1292 : c? TV=TV+MELVA3.VELCHE((IDIM+I-1)*NPG+LG,NK3)*HR(N,J,LG)
1293 : c? & *PGSQ(LG)*C3*C1*DEUPI*RPG(LG)
1294 : c? ELSEIF(XTG)THEN
1295 : c? C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)
1296 : c? TV=TV+MELVA3.VELCHE((IDIM+I-1)*NPG+LG,NK3)*HR(N,J,LG)
1297 : c? & *PGSQ(LG)*C3*C1*DEUPI*RPG(LG)
1298 : c? ENDIF
1299 : c? 114 CONTINUE
1300 :
1301 : U1=U1+WT(I,LG)*FN(J,LG)*PGSQ(LG)*C1*DEUPI*RPG(LG)
1302 :
1303 : 112 CONTINUE
1304 :
1305 : IF(XDIAG)THEN
1306 : UD=UD+U1
1307 : U1=0.D0
1308 : ENDIF
1309 :
1310 : IF(XBDF.AND.(.NOT.XDIAG))THEN
1311 : U1=U1/DELTAT
1312 : DO 116 N=1,NC
1313 : SU(N) = SU(N) + (2.D0*TN1(J,N) - 0.5D0*TN2(J,N))*U1
1314 : 116 continue
1315 : U1=U1*1.5D0
1316 : ELSEIF(.NOT.XBDF)THEN
1317 : U1=U1/DELTAT
1318 : do 1161 n=1,nc
1319 : SU(N) = SU(N) + (U1 + (AIMPL-1.D0)*(-TV))*TN1(J,N)
1320 : 1161 CONTINUE
1321 : ENDIF
1322 :
1323 : RF1(J,I,1)=RF1(J,I,1)+U1+(AIMPL*TV)
1324 :
1325 : 111 CONTINUE
1326 :
1327 : IF(XDIAG)THEN
1328 : UD=UD/DELTAT
1329 : IF(XBDF)THEN
1330 : do 119 n=1,nc
1331 : SU(N) = SU(N) + (2.D0*TN1(I,N) - 0.5D0*TN2(I,N))*UD
1332 : 119 continue
1333 : UD=UD*1.5D0
1334 : ELSE
1335 : do 118 n=1,nc
1336 : SU(N) = SU(N) + UD*TN1(I,N)
1337 : 118 continue
1338 : ENDIF
1339 : RF1(I,I,1)=RF1(I,I,1)+UD
1340 : ENDIF
1341 :
1342 : DO 117 N=1,NC
1343 : SM1(I,N)=SM1(I,N)+ SU(N)
1344 : 117 CONTINUE
1345 :
1346 : 110 CONTINUE
1347 : c write(6,*)' apres SM1(I,N)='
1348 : c write(6,1002)(SM1(I,1),I=1,np)
1349 : c write(6,1002)(SM1(I,2),I=1,np)
1350 : C...... Transitoire Fin
1351 : ENDIF
1352 : C=======================================================================
1353 :
1354 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1355 : IF(TCONV) THEN
1356 : C...... Convection
1357 : DO 410 I=1,MP
1357 : CALL INITD (SU,NC,0.D0)
1358 :
1359 : DO 411 J=1,NP
1360 : U3=0.D0
1361 : DO 412 LG=1,NPG
1362 :
1363 : C1=MELVA1.VELCHE(LG,NK1)
1364 : DO 414 N=1,IDIM
1365 : C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)
1366 : U3=U3+WT(I,LG)*HR(N,J,LG)*PGSQ(LG)*C3*C1*DEUPI*RPG(LG)
1367 : 414 CONTINUE
1368 :
1369 : 412 CONTINUE
1370 :
1371 : DO 416 N=1,NC
1372 : SU(N) = SU(N) + (AIMPL-1.D0)*U3*TN1(J,N)
1373 : 416 CONTINUE
1374 :
1375 : RF1(J,I,1)=RF1(J,I,1)+(AIMPL*U3)
1376 :
1377 : 411 CONTINUE
1378 :
1379 : DO 417 N=1,NC
1380 : SM1(I,N)=SM1(I,N)+ SU(N)
1381 : 417 CONTINUE
1382 :
1383 : 410 CONTINUE
1384 : C...... Convection Fin
1385 : ENDIF
1386 : C=======================================================================
1387 :
1388 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1389 : IF(TLAPN) THEN
1390 : C...... Laplacien
1391 : DO 310 I=1,MP
1391 : CALL INITD (SU,NC,0.D0)
1392 : DO 311 J=1,NP
1393 : U2=0.D0
1394 : UR=0.D0
1395 : DO 312 LG=1,NPG
1396 :
1397 : C2=MELVA2.VELCHE(LG,NK2)
1398 :
1399 : DO 313 N=1,IDIM
1400 : U2=U2+IZWH.HR(N,I,LG)*HR(N,J,LG)*PGSQ(LG)*C2*DEUPI*RPG(LG)
1401 : 313 CONTINUE
1402 :
1403 : IF(IHV.EQ.1.AND.IAXI.EQ.2)THEN
1404 : C cas 2D axi Laplacien Champ vitesse
1405 : UR=UR+WT(I,LG)*FN(J,LG)/RPG(LG)*PGSQ(LG)*C2*DEUPI
1406 : C cas 2D axi Laplacien Champ vitesse Fin
1407 : ENDIF
1408 :
1409 : 312 CONTINUE
1410 :
1411 : DO 316 N=1,NC
1412 : SU(N) = SU(N) + (AIMPL-1.D0)*(U2+UR)*TN1(J,N)
1413 : 316 CONTINUE
1414 :
1415 : RF1(J,I,1)=RF1(J,I,1)+(AIMPL*(U2+UR))
1416 :
1417 : 311 CONTINUE
1418 :
1419 : DO 317 N=1,NC
1420 : SM1(I,N)=SM1(I,N) + SU(N)
1421 : 317 CONTINUE
1422 :
1423 : 310 CONTINUE
1424 : C...... Laplacien Fin
1425 : ENDIF
1426 : C=======================================================================
1427 :
1428 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1429 : IF(TMASS) THEN
1430 : C...... Masse
1431 : DO 210 I=1,MP
1431 : CALL INITD (SU,NC,0.D0)
1432 : UD=0.D0
1433 : DO 211 J=1,NP
1434 : U1=0.D0
1435 : DO 212 LG=1,NPG
1436 :
1437 : C1=MELVA1.VELCHE(LG,NK1)
1438 : U1=U1+WT(I,LG)*FN(J,LG)*PGSQ(LG)*C1*DEUPI*RPG(LG)
1439 :
1440 :
1441 : 212 CONTINUE
1442 :
1443 : IF(XDIAG)THEN
1444 : UD=UD+U1
1445 : U1=0.D0
1446 : ELSE
1447 : DO 216 N=1,NC
1448 : SU(N) = SU(N) + (AIMPL-1.D0)*U1*TN1(J,N)
1449 : 216 CONTINUE
1450 : RF1(J,I,1)=RF1(J,I,1)+AIMPL*U1
1451 : ENDIF
1452 :
1453 : 211 CONTINUE
1454 :
1455 : IF(XDIAG)THEN
1456 : RF1(I,I,1)=RF1(I,I,1)+AIMPL*UD
1457 : DO 218 N=1,NC
1458 : SM1(I,N)=SM1(I,N)+ (AIMPL-1.D0)*UD*TN1(I,N)
1459 : 218 CONTINUE
1460 : ELSE
1461 : DO 217 N=1,NC
1462 : SM1(I,N)=SM1(I,N)+ SU(N)
1463 : 217 CONTINUE
1464 : ENDIF
1465 :
1466 : 210 CONTINUE
1467 : C...... Masse Fin
1468 : ENDIF
1469 : C=======================================================================
1470 :
1471 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1472 : IF(TSOUR) THEN
1473 : C...... Source
1474 : DO 610 I=1,MP
1475 : C1=1.D0
1476 : DO 617 N=1,NC
1477 : U4=0.D0
1478 : DO 615 LG=1,NPG
1479 : IF(TECHI)THEN
1480 : C1=MELVA1.VELCHE(LG,NK1)
1481 : ENDIF
1482 : C4=MELVA4.VELCHE((N-1)*NPG+LG,NK4)
1483 : U4=U4+WS(I,LG)*PGSQ(LG)*C4*C1*DEUPI*RPG(LG)
1484 : 615 CONTINUE
1485 : SM1(I,N)=SM1(I,N)+ U4
1486 : 617 CONTINUE
1487 : 610 CONTINUE
1488 : C...... Source Fin
1489 : ENDIF
1490 : C=======================================================================
1491 :
1492 :
1493 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1494 : C ...... Chargement Rigidite ou Matrik
1495 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN)THEN
1496 : C Cas RIGIDITE
1497 : IF(XRIG)THEN
1498 : * XMATRI=IMATTT(KE)
1499 : DO 901 I=1,MP
1500 : DO 901 J=1,NP
1501 : RE(I,J,ke)=RF1(J,I,1)
1502 : 901 CONTINUE
1503 : * SEGDES XMATRI
1504 : ELSE
1505 : C Cas MATRIK
1506 : DO 902 N=1,NBME
1507 : JMATR1=JRIGEL(4,NMATR0+N)
1508 : IPM4=JMATR1.LIZAFM(L,1)
1509 : DO 902 I=1,NP
1510 : DO 902 J=1,NP
1511 : IPM4.AM(KE,J,I)=RF1(J,I,1)
1512 : 902 CONTINUE
1513 : ENDIF
1514 : ENDIF
1515 : C ...... Chargement Second membre
1516 : DO 910 I=1,NP
1517 : I1=LECT(IPT1.NUM(I,KE))
1518 : DO 910 N=1,NC
1519 : MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N)
1520 : 910 CONTINUE
1521 : C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1522 :
1523 : 108 CONTINUE
1524 :
1525 : IF(TDFDT.OR.TMASS.OR.TCONV)THEN
1526 : SEGSUP MCHAM1,MELVA1
1527 : ENDIF
1528 :
1529 : IF(TSOUR)THEN
1530 : SEGSUP MCHAM4,MELVA4
1531 : ENDIF
1532 :
1533 : IF(TLAPN.OR.TCONV)THEN
1534 : SEGSUP MCHAM2,MELVA2
1535 : ENDIF
1536 :
1537 : IF(TCONV)THEN
1538 : SEGSUP MCHAM3,MELVA3
1539 : ENDIF
1540 :
1541 : NUTOEL=NUTOEL+NBEL
1542 :
1543 : SEGDES IPT1,IPT2
1544 :
1545 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN.OR.TKBBT)THEN
1546 : C Cas RIGIDITE
1547 : IF(XRIG)THEN
1548 : SEGDES xMATRI
1549 : ELSE
1550 : C Cas MATRIK
1551 : SEGDES IPM1
1552 : IF(NBME.GE.2)SEGDES IPM2
1553 : IF(NBME.GE.3)SEGDES IPM3
1554 : ENDIF
1555 : ENDIF
1556 :
1557 : SEGSUP IZFFM,IZHR,IZF1,IZH2
1558 : SEGSUP SAJT
1559 :
1560 : 101 CONTINUE
1561 :
1562 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN)THEN
1563 : IF(.NOT.XRIG)THEN
1564 : DO 141 M=1,NBME
1565 : JMATRI=JRIGEL(4,NMATR0+M)
1566 : SEGDES JMATRI
1567 : 141 CONTINUE
1568 : ENDIF
1569 : ENDIF
1570 :
1571 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
1572 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
1573 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
1574 : ELSEIF(TVITP)THEN
1575 :
1575 : CALL KRIPAD (MELEMS,MLENTI)
1576 : SEGDES MELEMS
1577 :
1578 : NUTOEL=0
1579 : SEGACT MELEME,MELEM2
1580 : NBSOUS=LISOUS(/1)
1581 : IF(NBSOUS.EQ.0)NBSOUS=1
1582 :
1583 : C Cas RIGIDITE Pleine
1584 : IF(XRIG)THEN
1585 : NRIGE=8
1586 : NRI=IRIGEL(/2)
1587 : NRIGEL=NBSOUS+NRI
1588 : SEGADJ MRIGID
1589 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
1590 :
1591 : ELSE
1592 : C Cas MATRIK Pleine
1593 : NRIGE=7
1594 : NKID =9
1595 : NKMT =7
1596 : NMATR0=JRIGEL(/2)
1597 : NMATRI=NMATR0+1
1598 : SEGADJ MATRIK
1599 :
1600 : NBMF=1
1601 : IF(IHV.EQ.1)NBMF=IDIM
1602 : SEGINI JMATRI
1603 : c write(6,*)' KBBT JMATRI=',JMATRI,MELEMS,MELEMP
1604 : JRIGEL(4,NMATR0+1)=JMATRI
1605 :
1606 : IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
1607 : JRIGEL(7,NMATR0+1)=-3
1608 : IF(IKAS.EQ.3)JRIGEL(7,NMATR0+1)=4
1609 : JRIGEL(1,NMATR0+1)=MELEME
1610 : JRIGEL(2,NMATR0+1)=MELEM2
1611 : KSPGP=MELEMS
1612 : KSPGD=MELEMP
1613 : IF(NBMF.EQ.1)THEN
1614 : LISPRI(1)=NOMP(1)//' '
1615 : LJSDUA(1)=NOMD(1)//' '
1616 : ELSE
1617 : DO 42 I=1,NBMF
1618 : LISPRI(I)=NOMP(I)//' '
1619 : LJSDUA(I)=NOMD(I)//' '
1620 : 42 CONTINUE
1621 : ENDIF
1622 : ELSEIF(IKAS.EQ.2)THEN
1623 : JRIGEL(7,NMATR0+1)=3
1624 : JRIGEL(1,NMATR0+1)=MELEM2
1625 : JRIGEL(2,NMATR0+1)=MELEME
1626 : KSPGP=MELEMP
1627 : KSPGD=MELEMS
1628 : IF(NBMF.EQ.1)THEN
1629 : LISPRI(1)=NOMD(1)//' '
1630 : LJSDUA(1)=NOMP(1)//' '
1631 : ELSE
1632 : DO 43 I=1,NBMF
1633 : LISPRI(I)=NOMD(I)//' '
1634 : LJSDUA(I)=NOMP(I)//' '
1635 : 43 CONTINUE
1636 : ENDIF
1637 : ELSE
1638 : RETURN
1639 : ENDIF
1640 : ENDIF
1641 :
1642 : IF(TPROJ)JRIGEL(7,NMATR0+1)=4
1643 :
1644 : IF(TKBBT)SEGACT MCHEL1
1645 : IF(TLINCO)THEN
1646 : SEGACT MTETA1
1647 : IF(XBDF)SEGACT MTETA2
1648 : ENDIF
1649 :
1650 : IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN
1651 : WRITE(6,*)' Geometries incompatibles dans ',nomper
1652 : C% Données incompatibles
1652 : CALL ERREUR (22)
1653 : RETURN
1654 : ENDIF
1655 :
1656 : IF(INEFMD.EQ.2.AND.MPRE.EQ.'CENTRE'.AND.IKAS.NE.2)THEN
1657 : C CAS Stabilisation via MACRO CENTRE
1658 : c write(6,*)' CAS Stabilisation via MACRO CENTRE '
1659 :
1660 : TYPE=' '
1660 : CALL ACMO (MTABZ,'MELSTB',TYPE,MELSTB)
1661 : TYPE=' '
1661 : CALL ACMO (MTABZ,'MCHPOC',TYPE,MCHPOC)
1662 :
1663 : SEGACT MELSTB
1664 : IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4
1665 : IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8
1666 : NBNN=MELSTB.NUM(/1)
1667 : NBSOUS=0
1668 : NBREF=0
1669 : SEGINI MELEMA
1670 : MELEMA.ITYPEL=MELSTB.ITYPEL
1671 :
1672 : NKPE=4
1673 : IF(IDIM.EQ.3)NKPE=8
1674 : DO 4878 k=1,NBELEM
1675 : MI=(K-1)*NKPE+1
1676 : DO 4879 I=1,NBNN
1677 : MELEMA.NUM(I,K)=MELSTB.NUM(I,MI)
1678 : 4879 CONTINUE
1679 : 4878 CONTINUE
1680 :
1681 : BETA0=-ABS(CSTAB)
1682 : IF(TPROJ)BETA0=ABS(CSTAB)
1683 : c write(6,*)'TPROJ BETA0=',TPROJ,BETA0
1684 : NK=0
1685 : NMATR0=JRIGEL(/2)
1686 : NMATRI=NMATR0+1
1687 : SEGADJ MATRIK
1688 : NBMF=1
1689 : NBSOUS=1
1690 : SEGINI JMATRS
1691 : JRIGEL(4,NMATR0+1)=JMATRS
1692 : JMATRS.KSPGP=MELEMC
1693 : JMATRS.KSPGD=MELEMC
1694 : JRIGEL(1,NMATR0+1)=MELEMA
1695 : JRIGEL(2,NMATR0+1)=MELEMA
1696 : JRIGEL(7,NMATR0+1)=0
1696 : CALL LRCHT (MCHPOC,MPOVAL,TYPE,IGEOM)
1697 :
1698 :
1699 : NBSOUS=MELSTB.LISOUS(/1)
1700 : IF(NBSOUS.NE.0)THEN
1700 : CALL ERREUR (5)
1701 : ENDIF
1702 :
1703 : NBEL=MELEMA.NUM(/2)
1704 : NBCI=MELSTB.NUM(/2)
1705 : NP =MELSTB.NUM(/1)
1706 : MP =NP
1707 :
1708 : SEGINI IZAFM
1709 : JMATRS.LIZAFM(1,1)=IZAFM
1710 : JMATRS.LISPRI(1)=NOMD(1)
1711 : JMATRS.LJSDUA(1)=NOMD(1)
1712 :
1712 : CALL KRIPAD (MELEMC,MLENT1)
1713 :
1714 : DO 33 K=1,NBEL
1715 :
1716 : NK=NK+1
1717 :
1718 : DO 32 J=1,NP
1719 : K1=MLENT1.LECT(MELEMA.NUM(J,K))
1720 : II=J
1721 : DO 34 I=1,NP
1722 : U=VPOCHA(K1,I)*BETA0
1723 : IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0
1724 : IF(II.LE.NP)THEN
1725 : AM(K,II,J)=U
1726 : ELSE
1727 : AM(K,II-NP,J)=U
1728 : ENDIF
1729 : II=II+1
1730 : 34 CONTINUE
1731 : 32 CONTINUE
1732 : 33 CONTINUE
1733 :
1734 : SEGDES MELEMA,MELSTB,MELEMC,IZAFM,MPOVAL
1735 : SEGSUP MLENT1
1736 :
1737 : c write(6,*)' Fin Stab '
1738 : ENDIF
1739 : SEGACT MELEME,MELEM2
1740 :
1741 : DO 201 L=1,MAX(1,LISOUS(/1))
1742 : IPT1=MELEME
1743 : IPT2=MELEM2
1744 : IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
1745 : IF(MELEM2.LISOUS(/1).NE.0)IPT2=MELEM2.LISOUS(L)
1746 : SEGACT IPT1,IPT2
1747 :
1748 : C /S INEFMD : Type formulation
1749 : C INEFMD = 1 LINE,
1750 : C = 2 MACRO,
1751 : C = 3 QUADRATIQUE,
1752 : C = 4 LINB.
1753 : NOM0='????????'
1754 : IF(INEFMD.EQ.1)THEN
1755 : IF(MPRE.EQ.'CENTRE ')NOM0=NOMS(IPT1.ITYPEL)//' '
1756 : IF(MPRE.EQ.'MSOMMET ')NOM0=NOMS(IPT1.ITYPEL)//'P1P1'
1757 : ELSEIF(INEFMD.EQ.2)THEN
1758 : IF(MPRE.EQ.'CENTRE ')NOM0=NOMS(IPT1.ITYPEL)//' '
1759 :
1760 : c ces elements ne marchent pas
1761 : IF(MPRE.EQ.'CENTREP0')NOM0=NOMS(IPT1.ITYPEL)//'MCP0'
1762 : IF(MPRE.EQ.'CENTREP1')NOM0=NOMS(IPT1.ITYPEL)//'MCP1'
1763 :
1764 : IF(MPRE.EQ.'MSOMMET ')NOM0=NOMS(IPT1.ITYPEL)//'MCF1'
1765 : ELSEIF(INEFMD.EQ.3)THEN
1766 : IF(MPRE.EQ.'CENTREP0')NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
1767 : IF(MPRE.EQ.'CENTREP1')NOM0=NOMS(IPT1.ITYPEL)//'PRP1'
1768 : IF(MPRE.EQ.'MSOMMET ')NOM0=NOMS(IPT1.ITYPEL)//'PFP1'
1769 : ELSEIF(INEFMD.EQ.4)THEN
1770 : IF(MPRE.EQ.'MSOMMET ')NOM0=NOMS(IPT1.ITYPEL)//' '
1771 : ENDIF
1772 : IF(NOM0.EQ.'????????')THEN
1773 : C% Le type d'element fini Vitesse/pression ne convient pas : %m1:8 %m9:16 .
1774 : MOTERR( 1: 8) =LSOPT(INEFMD)
1775 : MOTERR( 9:16) =MPRE
1776 : WRITE(IOIMP,*)'Operateur : ',NOMPER
1776 : CALL ERREUR (933)
1777 : RETURN
1778 : ENDIF
1779 :
1780 : c write(6,*)nomper,' NOM0=',nom0,
1781 : c write(6,*)' 2eme KALPBG NOM0=',NOM0,IPT1
1781 : CALL KALPBG (NOM0,'FONFORM ',IZFFM)
1782 : if(IZFFM.eq.0)return
1783 :
1784 : SEGACT IZFFM*MOD
1785 : IZHR=KZHR(1)
1786 : IZH2=KZHR(2)
1787 :
1788 : SEGACT IZHR*MOD,IZH2*MOD
1789 :
1790 : NES=GR(/1)
1791 : NPG=GR(/3)
1792 : IZF1=KTP(1)
1793 : SEGACT IZF1*MOD
1794 : NP = IPT1.NUM(/1)
1795 : MP = IPT2.NUM(/1)
1796 : IF(IKAS.EQ.2)THEN
1797 : MP = IPT1.NUM(/1)
1798 : NP = IPT2.NUM(/1)
1799 : ENDIF
1800 :
1801 : NBEL=IPT1.NUM(/2)
1802 :
1803 : SEGINI SAJT
1804 :
1805 : C Cas RIGIDITE
1806 : IF(XRIG)THEN
1807 : IRIGEL(1,NRI+L)=MELEME
1808 : COERIG(L)=1.D0
1809 :
1810 : IRIGEL(7,NRI+L)=0
1811 : IF(TCONV)IRIGEL(7,NRI+L)=2
1812 :
1813 : NBME=1
1814 : IF(IHV.EQ.1)NBME=IDIM
1815 : NLIGRP=NP
1816 : NLIGRD=MP
1817 : SEGINI DESCR
1818 : IRIGEL(3,NRI+L)=DESCR
1819 : IF(NBME.EQ.1)THEN
1820 : DO 202 I=1,NLIGRP
1821 : LISINC(I)=NOMP(1)//' '
1822 : NOELEP(I)=I
1823 : 202 CONTINUE
1824 : DO 203 I=1,NLIGRD
1825 : LISDUA(I)=NOMD(1)//' '
1826 : NOELED(I)=I
1827 : 203 CONTINUE
1828 : ELSE
1829 : ENDIF
1830 : SEGDES DESCR
1831 :
1832 : NELRIG=NBEL
1833 : SEGINI xMATRI
1834 : IRIGEL(4,NRI+L)=xMATRI
1835 : c write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI
1836 :
1837 : * DO 204 K=1,NELRIG
1838 : * SEGINI XMATRI
1839 : c write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI
1840 : * IMATTT(K)=XMATRI
1841 : * 204 CONTINUE
1842 :
1843 : ELSE
1844 : C Cas MATRIK
1845 : NBMF=LIZAFM(/2)
1846 : SEGINI IPM1
1847 : LIZAFM(L,1)=IPM1
1848 : IPM2=IPM1
1849 : IPM3=IPM1
1850 : IF(NBMF.GE.2)THEN
1851 : SEGINI IPM2
1852 : LIZAFM(L,2)=IPM2
1853 : ENDIF
1854 : IF(NBMF.GE.3)THEN
1855 : SEGINI IPM3
1856 : LIZAFM(L,3)=IPM3
1857 : ENDIF
1858 : ENDIF
1859 :
1860 : C----Ro
1861 : IK1=1
1862 : IF(TKBBT)THEN
1863 : MCHAM1=MCHEL1.ICHAML(L)
1864 : SEGACT MCHAM1
1865 :
1866 : MELVA1=MCHAM1.IELVAL(1)
1867 : SEGACT MELVA1
1868 : N1PTEL=MELVA1.VELCHE(/1)
1869 : N1EL=MELVA1.VELCHE(/2)
1870 : IF(N1EL.EQ.1)THEN
1871 : IK1=1
1872 : ELSEIF(N1EL.EQ.NBEL)THEN
1873 : IK1=0
1874 : ENDIF
1875 : ENDIF
1876 :
1877 : c write(6,*)'AVT208 NC=',NC,'IK1=',IK1,'NP=',NP,'MP=',MP,'MP1=',MP1
1878 : C===============================================
1879 : DO 208 KE=1,NBEL
1880 :
1881 : NK1=KE + IK1*(1 - KE)
1882 :
1883 : DO 209 I=1,NP
1884 : J=IPT1.NUM(I,KE)
1885 : DO 209 N=1,IDIM
1886 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
1887 : 209 CONTINUE
1888 :
1888 : CALL CALJBR (FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
1889 : * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
1890 :
1890 : CALL INITD (RF1,(NP*MP*IDIM),0.D0)
1890 : c CALL INITD (SM1,NP,0.D0)
1891 :
1892 : c DO 207 I=1,NP
1893 : c I1=MLENT1.LECT(IPT1.NUM(I,KE))
1894 : c DO 207 N=1,NC
1895 : c TN1(I,N)=MTETA1.VPOCHA(I1,N)
1896 : c207 CONTINUE
1897 :
1898 : c IF(XBDF)THEN
1899 : c DO 206 I=1,NP
1900 : c I2=MLENT2.LECT(IPT1.NUM(I,KE))
1901 : c DO 206 N=1,NC
1902 : c TN2(I,N)=MTETA2.VPOCHA(I2,N)
1903 : c206 CONTINUE
1904 : c ENDIF
1905 :
1906 :
1907 : C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1908 : IF(TKBBT) THEN
1909 : C...... Matrice C Ct
1910 :
1911 : c write(6,*)' Calcul C Ct'
1912 : DO 510 N=1,IDIM
1913 : DO 510 I=1,MP
1914 : DO 511 J=1,NP
1915 : U1=0.D0
1916 : DO 512 LG=1,NPG
1917 : C1=MELVA1.VELCHE(LG,NK1)
1918 : U1=U1+IZF1.FN(I,LG)*HR(N,J,LG)*PGSQ(LG)*DEUPI*RPG(LG)*C1
1919 : 512 CONTINUE
1920 : RF1(J,I,N)=RF1(J,I,N)+U1
1921 : 511 CONTINUE
1922 : 510 CONTINUE
1923 : C...... Matrice C Ct Fin
1924 : ENDIF
1925 : C=======================================================================
1926 :
1927 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1928 : C ...... Chargement Rigidite ou Matrik
1929 : c write(6,*)' Chargement Rigidite ou Matrik'
1930 : C Cas RIGIDITE
1931 : IF(XRIG)THEN
1932 : * XMATRI=IMATTT(KE)
1933 : DO 921 I=1,MP
1934 : DO 921 J=1,NP
1935 : RE(I,J,ke)=RF1(J,I,1)
1936 : 921 CONTINUE
1937 : * SEGDES XMATRI
1938 : ELSE
1939 : C Cas MATRIK
1940 : DO 923 N=1,NBMF
1941 : IPM4=LIZAFM(L,N)
1942 : c write(6,*)'IPM4=',IPM4
1943 : DO 922 I=1,MP
1944 : DO 922 J=1,NP
1945 : IPM4.AM(KE,J,I)=RF1(J,I,N)
1946 : 922 CONTINUE
1947 : 923 CONTINUE
1948 : c write(6,*)' Fin Chargement'
1949 : ENDIF
1950 : C ...... Chargement Second membre
1951 : c DO 920 I=1,NP
1952 : c I1=LECT(IPT1.NUM(I,KE))
1953 : c DO 920 N=1,NC
1954 : c MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N)
1955 : c920 CONTINUE
1956 : C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1957 :
1958 :
1959 :
1960 : 208 CONTINUE
1961 :
1962 : IF(TKBBT)THEN
1963 : SEGSUP MCHAM1,MELVA1
1964 : ENDIF
1965 :
1966 : NUTOEL=NUTOEL+NBEL
1967 :
1968 : SEGDES IPT1,IPT2
1969 :
1970 : IF(TDFDT.OR.TMASS.OR.TCONV.OR.TLAPN.OR.TKBBT)THEN
1971 : C Cas RIGIDITE
1972 : IF(XRIG)THEN
1973 : SEGDES xMATRI
1974 : ELSE
1975 : C Cas MATRIK
1976 : SEGDES IPM1
1977 : IF(NBMF.GE.2)SEGDES IPM2
1978 : IF(NBMF.GE.3)SEGDES IPM3
1979 : ENDIF
1980 : ENDIF
1981 :
1982 : SEGSUP IZFFM,IZHR,IZF1,IZH2
1983 : SEGSUP SAJT
1984 :
1985 : 201 CONTINUE
1986 :
1987 : C Cas RIGIDITE
1988 : IF(XRIG)THEN
1989 : ELSE
1990 : C Cas MATRIK
1991 : NMATRI=JRIGEL(/2)
1992 : DO 205 I=1,NMATRI
1993 : JMATRI=JRIGEL(4,I)
1994 : SEGDES JMATRI
1995 : 205 CONTINUE
1996 : ENDIF
1997 :
1998 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
1999 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
2000 : C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
2001 : ENDIF
2002 :
2003 :
2004 : IF(TDFDT)THEN
2005 : SEGSUP MCHEL1
2006 : ENDIF
2007 :
2008 : IF(TSOUR)THEN
2009 : SEGSUP MCHEL4
2010 : ENDIF
2011 :
2012 : IF(TLAPN.OR.TCONV)THEN
2013 : SEGSUP MCHEL2
2014 : ENDIF
2015 :
2016 : IF(TCONV)THEN
2017 : SEGSUP MCHEL3
2018 : ENDIF
2019 :
2020 :
2021 : SEGDES MCHPO1,MPOVA1
2022 : SEGDES MELEME,MELEM2
2023 :
2024 : SEGSUP MLENTI
2025 : IF(TLINCO)THEN
2026 : SEGSUP MLENT1
2027 : SEGDES MTETA1
2028 : IF(XBDF)THEN
2029 : SEGSUP MLENT2
2030 : SEGDES MTETA2
2031 : ENDIF
2032 : ENDIF
2033 :
2034 : c write(6,*)' NINKO=',NINKO
2035 : IF(NINKO.GE.2.AND.TLINCO)THEN
2036 : SEGSUP MLENT3
2037 : SEGDES MTETA3
2038 : IF(XBDF)THEN
2039 : SEGSUP MLENT4
2040 : SEGDES MTETA4
2041 : ENDIF
2042 : ENDIF
2043 :
2044 : C Cas RIGIDITE
2045 : IF(XRIG)THEN
2046 : c write(6,*)' On ecrit bien une rigidite ',MRIGID
2047 : SEGDES MRIGID
2047 : CALL ECROBJ ('RIGIDITE',MRIGID)
2048 : ELSE
2049 : C Cas MATRIK
2050 : c write(6,*)' On ecrit bien un MATRIK ',MATRIK
2051 : SEGDES MATRIK
2051 : CALL ECROBJ ('MATRIK',MATRIK)
2052 : ENDIF
2053 :
2054 : c write(6,*)'On ecrit bien un CHPOINT',MCHPO1
2054 : CALL ECROBJ ('CHPOINT',MCHPO1)
2055 : C*************************************************************************
2056 :
2057 : c write(6,*)' FIN YTCLSF'
2058 : RETURN
2059 : 1001 FORMAT(20(1X,I5))
2060 : 1002 FORMAT(10(1X,1PE11.4))
2061 : END
2062 :
2063 :
2064 :
2065 :
2066 :
2067 :
2068 :
2069 :
2070 :
2071 :
2072 :
2073 :
2074 :
2075 :
2076 :
2077 :
2078 :
2079 :
2080 :
2081 :
2082 :
2083 :
2084 :
2085 :