1 : C XLAP2E SOURCE MAGN 09/07/30 21:21:36 6458
2 : SUBROUTINE XLAP2E(ICOGRV,MPGRVF,MPROC,MPVITC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,
5 : $ LCLIMV,KRVIMP,MPVIMP,
6 : $ LCLITO,KRTOIM,
7 : $ NOMINC,
8 : $ MPMUC,
9 : $ IJACO,
10 : $ IMPR,IRET)
11 : IMPLICIT INTEGER(I-N)
12 : IMPLICIT REAL*8 (A-H,O-Z)
13 : C***********************************************************************
14 : C NOM : XLAP2E
15 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
16 : C VF 3D.
17 : C Ici, on ne calcule que des contributions 'compliquées' à
18 : C la matrice jacobienne faisant intervenir les
19 : C coefficients pour le calcul des gradients de vitesse
20 : C (ICOGRV)
21 : C (contributions à (d Res_{\rho e_t} / d var)
22 : C var prenant successivement les valeurs :
23 : C \rho, \rho u, \rho v, \rho w, \rho e_t )
24 : C Les contributions sont plus "compliquées" à calculer que
25 : C les simples (cf. xlap2c) car on a à dériver des produits
26 : C de fonctions de la vitesse.
27 : C xlap2d calcule aussi une partie des contributions
28 : C 'compliquées'.
29 : C
30 : C NOTE PERSO : On pourrait programmer ça plus lisiblement en stockant
31 : C les contributions dans un tableau de pointeurs (2
32 : C indices, c'est possible ?) et en effectuant des produits
33 : C matriciels (coeff. x matrices de dérivées).
34 : C On n'y coupera pas si on veut regrouper 2D et 3D...
35 : C On ne va pas le faire.
36 : C
37 : C LANGAGE : ESOPE
38 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
39 : C mél : gounand@semt2.smts.cea.fr
40 : C***********************************************************************
41 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
42 : C standard) à un objet de type MATRIK.
43 : C APPELE PAR : XLAP2A : Calcul de la matrice jacobienne du
44 : C résidu du laplacien VF 3D.
45 : C***********************************************************************
46 : C ENTREES : ICOGRV (type MCHELM) : coefficients pour le
47 : C calcul du gradient de la vitesse aux
48 : C interfaces.
49 : C MPGRVF (type MPOVAL) : gradient de la vitesse
50 : C aux interfaces.
51 : C MPROC (type MPOVAL) : masse volumique par
52 : C élément.
53 : C MPVITC (type MPOVAL) : vitesse par élément.
54 : C MPVOLU (type MPOVAL) : volume des éléments.
55 : C MPNORM (type MPOVAL) : normale aux faces.
56 : C MPSURF (type MPOVAL) : surface des faces.
57 : C MELEFL (type MELEME) : connectivités face-(centre
58 : C gauche, centre droit).
59 : C KRFACE (type MLENTI) : tableau de repérage dans
60 : C le maillage des faces des éléments.
61 : C KRCENT (type MLENTI) : tableau de repérage dans
62 : C le maillage des centres des éléments.
63 : C LCLIMV (type logique) : .TRUE. => CL de Dirichlet
64 : C sur la vitesse.
65 : C KRVIMP (type MLENTI) : tableau de repérage dans
66 : C maillage des CL de Dirichlet sur la
67 : C vitesse.
68 : C MPVIMP (type MPOVAL) : valeurs des CL de
69 : C Dirichlet sur la vitesse.
70 : C LCLITO (type logique) : .TRUE. => CL de Dirichlet
71 : C sur le tenseur des contraintes.
72 : C KRTOIM (type MLENTI) : tableau de repérage dans
73 : C maillage des CL de Dirichlet sur le tenseur des
74 : C contraintes
75 : C NOMINC (type MLMOTS) : noms des inconnues.
76 : C MPMUC (type MPOVAL) : viscosité dynamique (SI).
77 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
78 : C résidu du laplacien VF 3D.
79 : C SORTIES : -
80 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
81 : C***********************************************************************
82 : C VERSION : v1, 08/03/2002, version initiale
83 : C HISTORIQUE : v1, 08/03/2002, création
84 : C HISTORIQUE :
85 : C HISTORIQUE :
86 : C***********************************************************************
87 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
88 : C en cas de modification de ce sous-programme afin de faciliter
89 : C la maintenance !
90 : C***********************************************************************
91 : -INC CCOPTIO
92 : -INC SMCOORD
93 : -INC SMCHPOI
94 : POINTEUR MPGRVF.MPOVAL
95 : POINTEUR MPMUC.MPOVAL
96 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
97 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
98 : POINTEUR MPVIMP.MPOVAL
99 : -INC SMCHAML
100 : POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
101 : POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL,KDUNDZ.MELVAL
102 : -INC SMELEME
103 : POINTEUR MELEFL.MELEME
104 : POINTEUR MCOGRV.MELEME
105 : -INC SMLENTI
106 : POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
107 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
108 : -INC SMLMOTS
109 : POINTEUR NOMINC.MLMOTS
110 : -INC SMMATRIK
111 : POINTEUR IJACO.MATRIK
112 : *
113 : * Objet matrice élémentaire simplifié
114 : *
115 : SEGMENT GMATSI
116 : INTEGER POIPR1(NPP1,NEL1)
117 : INTEGER POIDU1(1,NEL1)
118 : INTEGER POIPR2(NPP2,NEL2)
119 : INTEGER POIDU2(2,NEL2)
120 : POINTEUR LMATSI(0).MATSIM
121 : ENDSEGMENT
122 : * Contributions compliquées 1 de la part du gradient de
123 : * vitesse (CCGRV1)
124 : POINTEUR CCGRV1.GMATSI
125 : SEGMENT MATSIM
126 : CHARACTER*8 NOMPRI,NOMDUA
127 : REAL*8 VALMA1(1,NPP1,NEL1)
128 : REAL*8 VALMA2(2,NPP2,NEL2)
129 : ENDSEGMENT
130 : POINTEUR RETRHO.MATSIM
131 : POINTEUR RETROU.MATSIM
132 : POINTEUR RETROV.MATSIM
133 : POINTEUR RETROW.MATSIM
134 : *
135 : REAL*8 MU
136 : *
137 : INTEGER IMPR,IRET
138 : *
139 : LOGICAL LCLIMV,LCLITO
140 : LOGICAL LMUR
141 : LOGICAL LCTR1A,LCTR1B
142 : *
143 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
144 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
145 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
146 : INTEGER NLCENP,NLCEND,NLFACE,NLCLTO,NLCLV
147 : INTEGER NLCED,NLCEG,NLFVI
148 : INTEGER NPTEL
149 : INTEGER ICOORX,NLCGAU,NLCDRO
150 : *
151 : REAL*8 ALPHA,UMALPH
152 : REAL*8 DRD,DRG
153 : REAL*8 DUXXF,DUXYF,DUXZF
154 : REAL*8 DUYXF,DUYYF,DUYZF
155 : REAL*8 DUZXF,DUZYF,DUZZF
156 : REAL*8 UXD,UXF,UXG
157 : REAL*8 UYD,UYF,UYG
158 : REAL*8 UZD,UZF,UZG
159 : REAL*8 XD,XF,XG,XFMXD,XFMXG
160 : REAL*8 YD,YF,YG,YFMYD,YFMYG
161 : REAL*8 ZD,ZF,ZG,ZFMZD,ZFMZG
162 : REAL*8 ALPHAX,ALPHAY,ALPHAZ
163 : REAL*8 CNX,CNY,CNZ
164 : REAL*8 SIGNOR,SURFFA,VOLUEL
165 : REAL*8 RHOP,UP,VP,WP
166 : REAL*8 FACTOR,FACTDU,FACTDV,FACTDW
167 : REAL*8 DUDRHO,DUDROU
168 : REAL*8 DVDRHO,DVDROV
169 : REAL*8 DWDRHO,DWDROW
170 : *
171 : * Executable statements
172 : *
173 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap2e.eso'
174 : * On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
175 : * prenant successivement les valeurs :
176 : * \rho, \rho u, \rho v, \rho w, \rho e_t
177 : * On dérive les termes : (\tens{\tau(\grad{u})} \prod \vect{u})
178 : * \pscal \vect{n}
179 : * ce qui donne deux contributions.
180 : * Contribution 1 :
181 : * ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
182 : * Contribution 2 :
183 : * ( \tens{\tau} \prod (d \vect{u} / d var)) \pscal \vect{n}
184 : * Note :
185 : * pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
186 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
187 : * RETRHO, RETROU, RETROV, RETROW
188 : IF (LCLIMV) THEN
189 : SEGACT KRVIMP
190 : SEGACT MPVIMP
191 : ENDIF
192 : IF (LCLITO) THEN
193 : SEGACT KRTOIM
194 : ENDIF
195 : SEGACT NOMINC
196 : SEGACT KRCENT
197 : SEGACT KRFACE
198 : SEGACT MELEFL
199 : SEGACT MPSURF
200 : SEGACT MPNORM
201 : SEGACT MPVOLU
202 : SEGACT MPMUC
203 : SEGACT MPROC
204 : SEGACT MPVITC
205 : SEGACT MPGRVF
206 : SEGACT ICOGRV
207 : NSOUCH=ICOGRV.IMACHE(/1)
208 : DO 1 ISOUCH=1,NSOUCH
209 : MCOGRV=ICOGRV.IMACHE(ISOUCH)
210 : JCOGRV=ICOGRV.ICHAML(ISOUCH)
211 : SEGACT JCOGRV
212 : KDUNDX=JCOGRV.IELVAL(1)
213 : KDUNDY=JCOGRV.IELVAL(2)
214 : KDUNDZ=JCOGRV.IELVAL(3)
215 : SEGDES JCOGRV
216 : SEGACT KDUNDX
217 : SEGACT KDUNDY
218 : SEGACT KDUNDZ
219 : SEGACT MCOGRV
220 : NELEM=MCOGRV.NUM(/2)
221 : NPTEL=MCOGRV.NUM(/1)
222 : NPP1=NPTEL-1
223 : NPP2=NPTEL-1
224 : NEL1=NELEM
225 : NEL2=NELEM
226 : IEL1=1
227 : IEL2=1
228 : SEGINI RETRHO
229 : SEGINI RETROU
230 : SEGINI RETROV
231 : SEGINI RETROW
232 : SEGINI CCGRV1
233 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
234 : RETRHO.NOMPRI(5:8)=' '
235 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
236 : RETRHO.NOMDUA(5:8)=' '
237 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
238 : RETROU.NOMPRI(5:8)=' '
239 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
240 : RETROU.NOMDUA(5:8)=' '
241 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
242 : RETROV.NOMPRI(5:8)=' '
243 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
244 : RETROV.NOMDUA(5:8)=' '
245 : RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
246 : RETROW.NOMPRI(5:8)=' '
247 : RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
248 : RETROW.NOMDUA(5:8)=' '
249 : DO 12 IELEM=1,NELEM
250 : * Le premier point du support de ICOGRV est un point FACE
251 : NGFACE=MCOGRV.NUM(1,IELEM)
252 : NLFACE=KRFACE.LECT(NGFACE)
253 : IF (NLFACE.EQ.0) THEN
254 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
255 : GOTO 9999
256 : ENDIF
257 : * On calcule la contribution 1 à la matrice jacobienne IJACO de la
258 : * face NGFACE
259 : * (points duaux : centres à gauche et à droite de la face)
260 : * (points primaux : une partie (bicoz conditions aux limites)
261 : * de ceux du stencil pour le calcul du gradient
262 : * à la face, ils doivent être des points centres)
263 : * Si le tenseur des contraintes sur la face est imposé par les
264 : * conditions aux limites, la contribution 1 de la face à IJACO est
265 : * nulle.
266 : LCTR1A=.TRUE.
267 : IF (LCLITO) THEN
268 : NLCLTO=KRTOIM.LECT(NGFACE)
269 : IF (NLCLTO.NE.0) THEN
270 : LCTR1A=.FALSE.
271 : ENDIF
272 : ENDIF
273 : IF (LCTR1A) THEN
274 : NGCGAU=MELEFL.NUM(1,NLFACE)
275 : NGCDRO=MELEFL.NUM(3,NLFACE)
276 : NLCGAU=KRCENT.LECT(NGCGAU)
277 : NLCDRO=KRCENT.LECT(NGCDRO)
278 : LMUR=(NGCGAU.EQ.NGCDRO)
279 : * On distingue le cas où la face est un bord du maillage (mur)
280 : * du cas où la face est interne au maillage
281 : IF (.NOT.LMUR) THEN
282 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
283 : XF = MCOORD.XCOOR(ICOORX)
284 : YF = MCOORD.XCOOR(ICOORX+1)
285 : ZF = MCOORD.XCOOR(ICOORX+2)
286 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
287 : XG = MCOORD.XCOOR(ICOORX)
288 : YG = MCOORD.XCOOR(ICOORX+1)
289 : ZG = MCOORD.XCOOR(ICOORX+2)
290 : XFMXG = XF - XG
291 : YFMYG = YF - YG
292 : ZFMZG = ZF - ZG
293 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
294 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
295 : XD = MCOORD.XCOOR(ICOORX)
296 : YD = MCOORD.XCOOR(ICOORX+1)
297 : ZD = MCOORD.XCOOR(ICOORX+2)
298 : XFMXD = XF - XD
299 : YFMYD = YF - YD
300 : ZFMZD = ZF - ZD
301 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
302 : ALPHA=DRG/(DRG+DRD)
303 : UMALPH= 1.0D0 - ALPHA
304 : ELSE
305 : ALPHA=0.0D0
306 : UMALPH=1.0D0
307 : ENDIF
308 : MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) +
309 : & ALPHA*MPMUC.VPOCHA(NLCDRO,1)
310 : * On calcule tout d'abord une interpolation de la vitesse sur la
311 : * face NGFACE de la même manière que dans xlap12.eso.
312 : IF (LCLIMV) THEN
313 : NLFVI=KRVIMP.LECT(NGFACE)
314 : ELSE
315 : NLFVI=0
316 : ENDIF
317 : * La vitesse est imposée sur la face, rien à calculer
318 : IF (NLFVI.NE.0) THEN
319 : UXF=MPVIMP.VPOCHA(NLFVI,1)
320 : UYF=MPVIMP.VPOCHA(NLFVI,2)
321 : UZF=MPVIMP.VPOCHA(NLFVI,3)
322 : ELSE
323 : NLCEG=KRCENT.LECT(NGCGAU)
324 : NLCED=KRCENT.LECT(NGCDRO)
325 : * Cas non au bord
326 : IF (.NOT.LMUR) THEN
327 : * Paramètres géométriques
328 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
329 : XF = MCOORD.XCOOR(ICOORX)
330 : YF = MCOORD.XCOOR(ICOORX+1)
331 : ZF = MCOORD.XCOOR(ICOORX+2)
332 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
333 : XG = MCOORD.XCOOR(ICOORX)
334 : YG = MCOORD.XCOOR(ICOORX+1)
335 : ZG = MCOORD.XCOOR(ICOORX+2)
336 : XFMXG = XF - XG
337 : YFMYG = YF - YG
338 : ZFMZG = ZF - ZG
339 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
340 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
341 : XD = MCOORD.XCOOR(ICOORX)
342 : YD = MCOORD.XCOOR(ICOORX+1)
343 : ZD = MCOORD.XCOOR(ICOORX+2)
344 : XFMXD = XF - XD
345 : YFMYD = YF - YD
346 : ZFMZD = ZF - ZD
347 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
348 : ALPHA=DRG/(DRG+DRD)
349 : UMALPH= 1.0D0 - ALPHA
350 : DUXXF = MPGRVF.VPOCHA(NLFACE,1)
351 : DUXYF = MPGRVF.VPOCHA(NLFACE,2)
352 : DUXZF = MPGRVF.VPOCHA(NLFACE,3)
353 : DUYXF = MPGRVF.VPOCHA(NLFACE,4)
354 : DUYYF = MPGRVF.VPOCHA(NLFACE,5)
355 : DUYZF = MPGRVF.VPOCHA(NLFACE,6)
356 : DUZXF = MPGRVF.VPOCHA(NLFACE,7)
357 : DUZYF = MPGRVF.VPOCHA(NLFACE,8)
358 : DUZZF = MPGRVF.VPOCHA(NLFACE,9)
359 : * Calcul de la vitesse
360 : UXG = MPVITC.VPOCHA(NLCEG,1)
361 : UYG = MPVITC.VPOCHA(NLCEG,2)
362 : UZG = MPVITC.VPOCHA(NLCEG,3)
363 : UXD = MPVITC.VPOCHA(NLCED,1)
364 : UYD = MPVITC.VPOCHA(NLCED,2)
365 : UZD = MPVITC.VPOCHA(NLCED,3)
366 : UXF = UMALPH * UXG + ALPHA * UXD
367 : UYF = UMALPH * UYG + ALPHA * UYD
368 : UZF = UMALPH * UZG + ALPHA * UZD
369 : * Correction de la vitesse lineaire exacte
370 : UXF = UXF +
371 : & (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
372 : & (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
373 : & (DUXZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
374 : UYF = UYF +
375 : & (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
376 : & (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
377 : & (DUYZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
378 : UZF = UZF +
379 : & (DUZXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
380 : & (DUZYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
381 : & (DUZZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
382 :
383 : * Cas au bord
384 : ELSE
385 : * Parametres geometriques
386 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
387 : XF = MCOORD.XCOOR(ICOORX)
388 : YF = MCOORD.XCOOR(ICOORX+1)
389 : ZF = MCOORD.XCOOR(ICOORX+2)
390 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
391 : XG = MCOORD.XCOOR(ICOORX)
392 : YG = MCOORD.XCOOR(ICOORX+1)
393 : ZG = MCOORD.XCOOR(ICOORX+2)
394 : XFMXG = XF - XG
395 : YFMYG = YF - YG
396 : ZFMZG = ZF - ZG
397 : DUXXF = MPGRVF.VPOCHA(NLFACE,1)
398 : DUXYF = MPGRVF.VPOCHA(NLFACE,2)
399 : DUXZF = MPGRVF.VPOCHA(NLFACE,3)
400 : DUYXF = MPGRVF.VPOCHA(NLFACE,4)
401 : DUYYF = MPGRVF.VPOCHA(NLFACE,5)
402 : DUYZF = MPGRVF.VPOCHA(NLFACE,6)
403 : DUZXF = MPGRVF.VPOCHA(NLFACE,7)
404 : DUZYF = MPGRVF.VPOCHA(NLFACE,8)
405 : DUZZF = MPGRVF.VPOCHA(NLFACE,9)
406 : * Calcul de la vitesse
407 : UXF = MPVITC.VPOCHA(NLCEG,1)
408 : UYF = MPVITC.VPOCHA(NLCEG,2)
409 : UZF = MPVITC.VPOCHA(NLCEG,3)
410 : * Correction de la vitesse lineaire exacte
411 : UXF = UXF +
412 : & (DUXXF * XFMXG ) +
413 : & (DUXYF * YFMYG ) +
414 : & (DUXZF * ZFMZG )
415 : UYF = UYF +
416 : & (DUYXF * XFMXG ) +
417 : & (DUYYF * YFMYG ) +
418 : & (DUYZF * ZFMZG )
419 : UZF = UZF +
420 : & (DUZXF * XFMXG ) +
421 : & (DUZYF * YFMYG ) +
422 : & (DUZZF * ZFMZG )
423 : ENDIF
424 : ENDIF
425 : * On distingue le cas où la face est un bord du maillage (mur)
426 : * du cas où la face est interne au maillage
427 : IF (.NOT.LMUR) THEN
428 : NPD=2
429 : ELSE
430 : NPD=1
431 : ENDIF
432 : NPP=NPTEL-1
433 : * IPD=1 : point à gauche du point NGFACE
434 : * IPD=2 : point à droite du point NGFACE
435 : DO 122 IPD=1,NPD
436 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
437 : IF (.NOT.LMUR) THEN
438 : CCGRV1.POIDU2(IPD,IEL2)=NPDUAL
439 : ELSE
440 : CCGRV1.POIDU1(IPD,IEL1)=NPDUAL
441 : ENDIF
442 : NLCEND=KRCENT.LECT(NPDUAL)
443 : IF (NLCEND.EQ.0) THEN
444 : WRITE(IOIMP,*) 'Erreur grave n°1'
445 : GOTO 9999
446 : ENDIF
447 : DO 124 IPP=1,NPP
448 : NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
449 : LCTR1B=.TRUE.
450 : IF (LCLIMV) THEN
451 : NLCLV=KRVIMP.LECT(NPPRIM)
452 : IF (NLCLV.NE.0) THEN
453 : LCTR1B=.FALSE.
454 : ENDIF
455 : ENDIF
456 : IF (.NOT.LCTR1B) THEN
457 : * Lorsque une contribution est nulle, on fixe artificiellement le
458 : * point primal égal au point dual.
459 : IF (.NOT.LMUR) THEN
460 : CCGRV1.POIPR2(IPP,IEL2)=NPDUAL
461 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
462 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
463 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
464 : RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
465 : ELSE
466 : CCGRV1.POIPR1(IPP,IEL1)=NPDUAL
467 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
468 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
469 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
470 : RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
471 : ENDIF
472 : ELSE
473 : * Les contributions 1 valent :
474 : * (d Res_{\rho e_t})_d / (d var)_p =
475 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \mu
476 : * * [ [ [ (n_x * (( 4/3 * u_f * \alpha_x) + (v_f * \alpha_y)
477 : * + (w_f * \alpha_z)))
478 : * + (n_y * ((u_f * \alpha_y) + (-2/3 * v_f * \alpha_x)))
479 : * + (n_z * ((u_f * \alpha_z) + (-2/3 * w_f * \alpha_x)))
480 : * ]
481 : * * ((du)_p / (d var)_p)
482 : * ]
483 : * + [ [ (n_x * ((-2/3 * u_f * \alpha_y) + (v_f * \alpha_x)))
484 : * + (n_y * ((u_f * \alpha_x) + ( 4/3 * v_f * \alpha_y)
485 : * + (w_f * \alpha_z)))
486 : * + (n_z * ((v_f * \alpha_z) + (-2/3 * w_f * \alpha_y)))
487 : * ]
488 : * * ((dv)_p / (d var)_p)
489 : * ]
490 : * + [ [ (n_x * ((-2/3 * u_f * \alpha_z) + (w_f * \alpha_x)))
491 : * + (n_y * ((-2/3 * v_f * \alpha_z) + (w_f * \alpha_y)))
492 : * + (n_z * ((u_f * \alpha_x) + (v_f * \alpha_y)
493 : * + ( 4/3 * w_f * \alpha_z)))
494 : * ]
495 : * * ((dw)_p / (d var)_p)
496 : * ]
497 : * ]
498 : NLCENP=KRCENT.LECT(NPPRIM)
499 : IF (NLCENP.EQ.0) THEN
500 : WRITE(IOIMP,*) 'Erreur grave n°2'
501 : GOTO 9999
502 : ENDIF
503 : * normale sortante pour IPD=1, rentrante pour IPD=2
504 : SIGNOR=(-1.D0)**(IPD+1)
505 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
506 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
507 : CNX =MPNORM.VPOCHA(NLFACE,1)
508 : CNY =MPNORM.VPOCHA(NLFACE,2)
509 : CNZ =MPNORM.VPOCHA(NLFACE,3)
510 : ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
511 : ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
512 : ALPHAZ=KDUNDZ.VELCHE(IPP+1,IELEM)
513 : RHOP =MPROC.VPOCHA(NLCENP,1)
514 : UP =MPVITC.VPOCHA(NLCENP,1)
515 : VP =MPVITC.VPOCHA(NLCENP,2)
516 : WP =MPVITC.VPOCHA(NLCENP,3)
517 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU
518 : FACTDU=(CNX*((( 4.D0/3.D0)*UXF*ALPHAX)
519 : $ +(UYF*ALPHAY)
520 : $ +(UZF*ALPHAZ)))
521 : $ + (CNY*((UXF*ALPHAY)
522 : $ +((-2.D0/3.D0)*UYF*ALPHAX)))
523 : $ + (CNZ*((UXF*ALPHAZ)
524 : $ +((-2.D0/3.D0)*UZF*ALPHAX)))
525 : FACTDV=(CNX*(((-2.D0/3.D0)*UXF*ALPHAY)
526 : $ +(UYF*ALPHAX)))
527 : $ + (CNY*((UXF*ALPHAX)
528 : $ +(( 4.D0/3.D0)*UYF*ALPHAY)
529 : $ +(UZF*ALPHAZ)))
530 : $ + (CNZ*((UYF*ALPHAZ)
531 : $ +((-2.D0/3.D0)*UZF*ALPHAY)))
532 : FACTDW=(CNX*(((-2.D0/3.D0)*UXF*ALPHAZ)
533 : $ +(UZF*ALPHAX)))
534 : $ + (CNY*(((-2.D0/3.D0)*UYF*ALPHAZ)
535 : $ +(UZF*ALPHAY)))
536 : $ + (CNZ*((UXF*ALPHAX)
537 : $ +(UYF*ALPHAY)
538 : $ +(( 4.D0/3.D0)*UZF*ALPHAZ)))
539 : DUDRHO=-UP /RHOP
540 : DUDROU=1.D0/RHOP
541 : DVDRHO=-VP /RHOP
542 : DVDROV=1.D0/RHOP
543 : DWDRHO=-WP /RHOP
544 : DWDROW=1.D0/RHOP
545 : IF (.NOT.LMUR) THEN
546 : CCGRV1.POIPR2(IPP,IEL2)=NPPRIM
547 : RETRHO.VALMA2(IPD,IPP,IEL2)=
548 : $ FACTOR*( (FACTDU*DUDRHO)
549 : $ +(FACTDV*DVDRHO)
550 : $ +(FACTDW*DWDRHO))
551 : RETROU.VALMA2(IPD,IPP,IEL2)=
552 : $ FACTOR* (FACTDU*DUDROU)
553 : RETROV.VALMA2(IPD,IPP,IEL2)=
554 : $ FACTOR* (FACTDV*DVDROV)
555 : RETROW.VALMA2(IPD,IPP,IEL2)=
556 : $ FACTOR* (FACTDW*DWDROW)
557 : ELSE
558 : CCGRV1.POIPR1(IPP,IEL1)=NPPRIM
559 : RETRHO.VALMA1(IPD,IPP,IEL1)=
560 : $ FACTOR*( (FACTDU*DUDRHO)
561 : $ +(FACTDV*DVDRHO)
562 : $ +(FACTDW*DWDRHO))
563 : RETROU.VALMA1(IPD,IPP,IEL1)=
564 : $ FACTOR* (FACTDU*DUDROU)
565 : RETROV.VALMA1(IPD,IPP,IEL1)=
566 : $ FACTOR* (FACTDV*DVDROV)
567 : RETROW.VALMA1(IPD,IPP,IEL1)=
568 : $ FACTOR* (FACTDW*DWDROW)
569 : ENDIF
570 : ENDIF
571 : 124 CONTINUE
572 : 122 CONTINUE
573 : IF (.NOT.LMUR) THEN
574 : IEL2=IEL2+1
575 : ELSE
576 : IEL1=IEL1+1
577 : ENDIF
578 : ENDIF
579 : 12 CONTINUE
580 : NPP1=NPTEL-1
581 : NPP2=NPTEL-1
582 : NEL1=IEL1-1
583 : NEL2=IEL2-1
584 : SEGADJ RETRHO
585 : SEGADJ RETROU
586 : SEGADJ RETROV
587 : SEGADJ RETROW
588 : SEGADJ CCGRV1
589 : CCGRV1.LMATSI(**)=RETRHO
590 : CCGRV1.LMATSI(**)=RETROU
591 : CCGRV1.LMATSI(**)=RETROV
592 : CCGRV1.LMATSI(**)=RETROW
593 : * On accumule les matrices résultantes dans IJACO
593 : CALL AJMTK(CCGRV1,IJACO,IMPR,IRET)
594 : IF (IRET.NE.0) GOTO 9999
595 : SEGSUP RETRHO
596 : SEGSUP RETROU
597 : SEGSUP RETROV
598 : SEGSUP RETROW
599 : SEGSUP CCGRV1
600 : *
601 : SEGDES MCOGRV
602 : SEGDES KDUNDZ
603 : SEGDES KDUNDY
604 : SEGDES KDUNDX
605 : 1 CONTINUE
606 : SEGDES ICOGRV
607 : SEGDES MPGRVF
608 : SEGDES MPVITC
609 : SEGDES MPMUC
610 : SEGDES MPROC
611 : SEGDES MPVOLU
612 : SEGDES MPNORM
613 : SEGDES MPSURF
614 : SEGDES MELEFL
615 : SEGDES KRFACE
616 : SEGDES KRCENT
617 : SEGDES NOMINC
618 : IF (LCLITO) THEN
619 : SEGDES KRTOIM
620 : ENDIF
621 : IF (LCLIMV) THEN
622 : SEGDES KRVIMP
623 : SEGDES MPVIMP
624 : ENDIF
625 : *
626 : * Normal termination
627 : *
628 : IRET=0
629 : RETURN
630 : *
631 : * Format handling
632 : *
633 : *
634 : * Error handling
635 : *
636 : 9999 CONTINUE
637 : IRET=1
638 : WRITE(IOIMP,*) 'An error was detected in subroutine xlap2e'
639 : RETURN
640 : *
641 : * End of subroutine XLAP2E
642 : *
643 : END
644 :
645 :
646 :
647 :
648 :
649 :
650 :
651 :
© Cast3M 2003 - All rights reserved.
Disclaimer