1 : C YLAP2E SOURCE MAGN 09/07/30 21:22:37 6458
2 : SUBROUTINE YLAP2E(ICOGRV,MPGRVF,MPROC,MPVITC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,
5 : $ LCLIMV,KRVIMP,MPVIMP,
6 : $ LCLITO,KRTOIM,
7 : $ NOMINC,
8 : $ MU,
9 : $ IJACO,
10 : $ IMPR,IRET)
11 : IMPLICIT INTEGER(I-N)
12 : IMPLICIT REAL*8 (A-H,O-Z)
13 : C***********************************************************************
14 : C NOM : YLAP2E
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. ylap1c) car on a à dériver des produits
26 : C de fonctions de la vitesse.
27 : C ylap2d 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 : YLAP2A : 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 MU (type réel) : 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, 28/08/2001, version initiale
83 : C HISTORIQUE : v1, 28/08/2001, 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 MPROC.MPOVAL ,MPVITC.MPOVAL
96 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
97 : POINTEUR MPVIMP.MPOVAL
98 : -INC SMCHAML
99 : POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
100 : POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL,KDUNDZ.MELVAL
101 : -INC SMELEME
102 : POINTEUR MELEFL.MELEME
103 : POINTEUR MCOGRV.MELEME
104 : -INC SMLENTI
105 : POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
106 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
107 : -INC SMLMOTS
108 : POINTEUR NOMINC.MLMOTS
109 : -INC SMMATRIK
110 : POINTEUR IJACO.MATRIK
111 : *
112 : * Objet matrice élémentaire simplifié
113 : *
114 : SEGMENT GMATSI
115 : INTEGER POIPR1(NPP1,NEL1)
116 : INTEGER POIDU1(1,NEL1)
117 : INTEGER POIPR2(NPP2,NEL2)
118 : INTEGER POIDU2(2,NEL2)
119 : POINTEUR LMATSI(0).MATSIM
120 : ENDSEGMENT
121 : * Contributions compliquées 1 de la part du gradient de
122 : * vitesse (CCGRV1)
123 : POINTEUR CCGRV1.GMATSI
124 : SEGMENT MATSIM
125 : CHARACTER*8 NOMPRI,NOMDUA
126 : REAL*8 VALMA1(1,NPP1,NEL1)
127 : REAL*8 VALMA2(2,NPP2,NEL2)
128 : ENDSEGMENT
129 : POINTEUR RETRHO.MATSIM
130 : POINTEUR RETROU.MATSIM
131 : POINTEUR RETROV.MATSIM
132 : POINTEUR RETROW.MATSIM
133 : *
134 : REAL*8 MU
135 : *
136 : INTEGER IMPR,IRET
137 : *
138 : LOGICAL LCLIMV,LCLITO
139 : LOGICAL LMUR
140 : LOGICAL LCTR1A,LCTR1B
141 : *
142 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
143 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
144 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
145 : INTEGER NLCENP,NLCEND,NLFACE,NLCLTO,NLCLV
146 : INTEGER NLCED,NLCEG,NLFVI
147 : INTEGER NPTEL
148 : INTEGER ICOORX
149 : *
150 : REAL*8 ALPHA,UMALPH
151 : REAL*8 DRD,DRG
152 : REAL*8 DUXXF,DUXYF,DUXZF
153 : REAL*8 DUYXF,DUYYF,DUYZF
154 : REAL*8 DUZXF,DUZYF,DUZZF
155 : REAL*8 UXD,UXF,UXG
156 : REAL*8 UYD,UYF,UYG
157 : REAL*8 UZD,UZF,UZG
158 : REAL*8 XD,XF,XG,XFMXD,XFMXG
159 : REAL*8 YD,YF,YG,YFMYD,YFMYG
160 : REAL*8 ZD,ZF,ZG,ZFMZD,ZFMZG
161 : REAL*8 ALPHAX,ALPHAY,ALPHAZ
162 : REAL*8 CNX,CNY,CNZ
163 : REAL*8 SIGNOR,SURFFA,VOLUEL
164 : REAL*8 RHOP,UP,VP,WP
165 : REAL*8 FACTOR,FACTDU,FACTDV,FACTDW
166 : REAL*8 DUDRHO,DUDROU
167 : REAL*8 DVDRHO,DVDROV
168 : REAL*8 DWDRHO,DWDROW
169 : *
170 : * Executable statements
171 : *
172 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap2e.eso'
173 : * On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
174 : * prenant successivement les valeurs :
175 : * \rho, \rho u, \rho v, \rho w, \rho e_t
176 : * On dérive les termes : (\tens{\tau(\grad{u})} \prod \vect{u})
177 : * \pscal \vect{n}
178 : * ce qui donne deux contributions.
179 : * Contribution 1 :
180 : * ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
181 : * Contribution 2 :
182 : * ( \tens{\tau} \prod (d \vect{u} / d var)) \pscal \vect{n}
183 : * Note :
184 : * pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
185 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
186 : * RETRHO, RETROU, RETROV, RETROW
187 : IF (LCLIMV) THEN
188 : SEGACT KRVIMP
189 : SEGACT MPVIMP
190 : ENDIF
191 : IF (LCLITO) THEN
192 : SEGACT KRTOIM
193 : ENDIF
194 : SEGACT NOMINC
195 : SEGACT KRCENT
196 : SEGACT KRFACE
197 : SEGACT MELEFL
198 : SEGACT MPSURF
199 : SEGACT MPNORM
200 : SEGACT MPVOLU
201 : SEGACT MPROC
202 : SEGACT MPVITC
203 : SEGACT MPGRVF
204 : SEGACT ICOGRV
205 : NSOUCH=ICOGRV.IMACHE(/1)
206 : DO 1 ISOUCH=1,NSOUCH
207 : MCOGRV=ICOGRV.IMACHE(ISOUCH)
208 : JCOGRV=ICOGRV.ICHAML(ISOUCH)
209 : SEGACT JCOGRV
210 : KDUNDX=JCOGRV.IELVAL(1)
211 : KDUNDY=JCOGRV.IELVAL(2)
212 : KDUNDZ=JCOGRV.IELVAL(3)
213 : SEGDES JCOGRV
214 : SEGACT KDUNDX
215 : SEGACT KDUNDY
216 : SEGACT KDUNDZ
217 : SEGACT MCOGRV
218 : NELEM=MCOGRV.NUM(/2)
219 : NPTEL=MCOGRV.NUM(/1)
220 : NPP1=NPTEL-1
221 : NPP2=NPTEL-1
222 : NEL1=NELEM
223 : NEL2=NELEM
224 : IEL1=1
225 : IEL2=1
226 : SEGINI RETRHO
227 : SEGINI RETROU
228 : SEGINI RETROV
229 : SEGINI RETROW
230 : SEGINI CCGRV1
231 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
232 : RETRHO.NOMPRI(5:8)=' '
233 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
234 : RETRHO.NOMDUA(5:8)=' '
235 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
236 : RETROU.NOMPRI(5:8)=' '
237 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
238 : RETROU.NOMDUA(5:8)=' '
239 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
240 : RETROV.NOMPRI(5:8)=' '
241 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
242 : RETROV.NOMDUA(5:8)=' '
243 : RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
244 : RETROW.NOMPRI(5:8)=' '
245 : RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
246 : RETROW.NOMDUA(5:8)=' '
247 : DO 12 IELEM=1,NELEM
248 : * Le premier point du support de ICOGRV est un point FACE
249 : NGFACE=MCOGRV.NUM(1,IELEM)
250 : NLFACE=KRFACE.LECT(NGFACE)
251 : IF (NLFACE.EQ.0) THEN
252 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
253 : GOTO 9999
254 : ENDIF
255 : * On calcule la contribution 1 à la matrice jacobienne IJACO de la
256 : * face NGFACE
257 : * (points duaux : centres à gauche et à droite de la face)
258 : * (points primaux : une partie (bicoz conditions aux limites)
259 : * de ceux du stencil pour le calcul du gradient
260 : * à la face, ils doivent être des points centres)
261 : * Si le tenseur des contraintes sur la face est imposé par les
262 : * conditions aux limites, la contribution 1 de la face à IJACO est
263 : * nulle.
264 : LCTR1A=.TRUE.
265 : IF (LCLITO) THEN
266 : NLCLTO=KRTOIM.LECT(NGFACE)
267 : IF (NLCLTO.NE.0) THEN
268 : LCTR1A=.FALSE.
269 : ENDIF
270 : ENDIF
271 : IF (LCTR1A) THEN
272 : NGCGAU=MELEFL.NUM(1,NLFACE)
273 : NGCDRO=MELEFL.NUM(3,NLFACE)
274 : LMUR=(NGCGAU.EQ.NGCDRO)
275 : * On calcule tout d'abord une interpolation de la vitesse sur la
276 : * face NGFACE de la même manière que dans ylap13.eso.
277 : IF (LCLIMV) THEN
278 : NLFVI=KRVIMP.LECT(NGFACE)
279 : ELSE
280 : NLFVI=0
281 : ENDIF
282 : * La vitesse est imposée sur la face, rien à calculer
283 : IF (NLFVI.NE.0) THEN
284 : UXF=MPVIMP.VPOCHA(NLFVI,1)
285 : UYF=MPVIMP.VPOCHA(NLFVI,2)
286 : UZF=MPVIMP.VPOCHA(NLFVI,3)
287 : ELSE
288 : NLCEG=KRCENT.LECT(NGCGAU)
289 : NLCED=KRCENT.LECT(NGCDRO)
290 : * Cas non au bord
291 : IF (.NOT.LMUR) THEN
292 : * Paramètres géométriques
293 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
294 : XF = MCOORD.XCOOR(ICOORX)
295 : YF = MCOORD.XCOOR(ICOORX+1)
296 : ZF = MCOORD.XCOOR(ICOORX+2)
297 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
298 : XG = MCOORD.XCOOR(ICOORX)
299 : YG = MCOORD.XCOOR(ICOORX+1)
300 : ZG = MCOORD.XCOOR(ICOORX+2)
301 : XFMXG = XF - XG
302 : YFMYG = YF - YG
303 : ZFMZG = ZF - ZG
304 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
305 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
306 : XD = MCOORD.XCOOR(ICOORX)
307 : YD = MCOORD.XCOOR(ICOORX+1)
308 : ZD = MCOORD.XCOOR(ICOORX+2)
309 : XFMXD = XF - XD
310 : YFMYD = YF - YD
311 : ZFMZD = ZF - ZD
312 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
313 : ALPHA=DRG/(DRG+DRD)
314 : UMALPH= 1.0D0 - ALPHA
315 : DUXXF = MPGRVF.VPOCHA(NLFACE,1)
316 : DUXYF = MPGRVF.VPOCHA(NLFACE,2)
317 : DUXZF = MPGRVF.VPOCHA(NLFACE,3)
318 : DUYXF = MPGRVF.VPOCHA(NLFACE,4)
319 : DUYYF = MPGRVF.VPOCHA(NLFACE,5)
320 : DUYZF = MPGRVF.VPOCHA(NLFACE,6)
321 : DUZXF = MPGRVF.VPOCHA(NLFACE,7)
322 : DUZYF = MPGRVF.VPOCHA(NLFACE,8)
323 : DUZZF = MPGRVF.VPOCHA(NLFACE,9)
324 : * Calcul de la vitesse
325 : UXG = MPVITC.VPOCHA(NLCEG,1)
326 : UYG = MPVITC.VPOCHA(NLCEG,2)
327 : UZG = MPVITC.VPOCHA(NLCEG,3)
328 : UXD = MPVITC.VPOCHA(NLCED,1)
329 : UYD = MPVITC.VPOCHA(NLCED,2)
330 : UZD = MPVITC.VPOCHA(NLCED,3)
331 : UXF = UMALPH * UXG + ALPHA * UXD
332 : UYF = UMALPH * UYG + ALPHA * UYD
333 : UZF = UMALPH * UZG + ALPHA * UZD
334 : * Correction de la vitesse lineaire exacte
335 : UXF = UXF +
336 : & (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
337 : & (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
338 : & (DUXZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
339 : UYF = UYF +
340 : & (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
341 : & (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
342 : & (DUYZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
343 : UZF = UZF +
344 : & (DUZXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
345 : & (DUZYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
346 : & (DUZZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
347 :
348 : * Cas au bord
349 : ELSE
350 : * Parametres geometriques
351 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
352 : XF = MCOORD.XCOOR(ICOORX)
353 : YF = MCOORD.XCOOR(ICOORX+1)
354 : ZF = MCOORD.XCOOR(ICOORX+2)
355 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
356 : XG = MCOORD.XCOOR(ICOORX)
357 : YG = MCOORD.XCOOR(ICOORX+1)
358 : ZG = MCOORD.XCOOR(ICOORX+2)
359 : XFMXG = XF - XG
360 : YFMYG = YF - YG
361 : ZFMZG = ZF - ZG
362 : DUXXF = MPGRVF.VPOCHA(NLFACE,1)
363 : DUXYF = MPGRVF.VPOCHA(NLFACE,2)
364 : DUXZF = MPGRVF.VPOCHA(NLFACE,3)
365 : DUYXF = MPGRVF.VPOCHA(NLFACE,4)
366 : DUYYF = MPGRVF.VPOCHA(NLFACE,5)
367 : DUYZF = MPGRVF.VPOCHA(NLFACE,6)
368 : DUZXF = MPGRVF.VPOCHA(NLFACE,7)
369 : DUZYF = MPGRVF.VPOCHA(NLFACE,8)
370 : DUZZF = MPGRVF.VPOCHA(NLFACE,9)
371 : * Calcul de la vitesse
372 : UXF = MPVITC.VPOCHA(NLCEG,1)
373 : UYF = MPVITC.VPOCHA(NLCEG,2)
374 : UZF = MPVITC.VPOCHA(NLCEG,3)
375 : * Correction de la vitesse lineaire exacte
376 : UXF = UXF +
377 : & (DUXXF * XFMXG ) +
378 : & (DUXYF * YFMYG ) +
379 : & (DUXZF * ZFMZG )
380 : UYF = UYF +
381 : & (DUYXF * XFMXG ) +
382 : & (DUYYF * YFMYG ) +
383 : & (DUYZF * ZFMZG )
384 : UZF = UZF +
385 : & (DUZXF * XFMXG ) +
386 : & (DUZYF * YFMYG ) +
387 : & (DUZZF * ZFMZG )
388 : ENDIF
389 : ENDIF
390 : * On distingue le cas où la face est un bord du maillage (mur)
391 : * du cas où la face est interne au maillage
392 : IF (.NOT.LMUR) THEN
393 : NPD=2
394 : ELSE
395 : NPD=1
396 : ENDIF
397 : NPP=NPTEL-1
398 : * IPD=1 : point à gauche du point NGFACE
399 : * IPD=2 : point à droite du point NGFACE
400 : DO 122 IPD=1,NPD
401 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
402 : IF (.NOT.LMUR) THEN
403 : CCGRV1.POIDU2(IPD,IEL2)=NPDUAL
404 : ELSE
405 : CCGRV1.POIDU1(IPD,IEL1)=NPDUAL
406 : ENDIF
407 : NLCEND=KRCENT.LECT(NPDUAL)
408 : IF (NLCEND.EQ.0) THEN
409 : WRITE(IOIMP,*) 'Erreur grave n°1'
410 : GOTO 9999
411 : ENDIF
412 : DO 124 IPP=1,NPP
413 : NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
414 : NLCENP=KRCENT.LECT(NPPRIM)
415 : IF(NLCENP .EQ. 0)THEN
416 : * Lorsque une contribution est nulle, on fixe artificiellement le
417 : * point primal égal au point dual.
418 : IF (.NOT.LMUR) THEN
419 : CCGRV1.POIPR2(IPP,IEL2)=NPDUAL
420 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
421 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
422 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
423 : RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
424 : ELSE
425 : CCGRV1.POIPR1(IPP,IEL1)=NPDUAL
426 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
427 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
428 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
429 : RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
430 : ENDIF
431 : ELSE
432 : * Les contributions 1 valent :
433 : * (d Res_{\rho e_t})_d / (d var)_p =
434 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \mu
435 : * * [ [ [ (n_x * (( 4/3 * u_f * \alpha_x) + (v_f * \alpha_y)
436 : * + (w_f * \alpha_z)))
437 : * + (n_y * ((u_f * \alpha_y) + (-2/3 * v_f * \alpha_x)))
438 : * + (n_z * ((u_f * \alpha_z) + (-2/3 * w_f * \alpha_x)))
439 : * ]
440 : * * ((du)_p / (d var)_p)
441 : * ]
442 : * + [ [ (n_x * ((-2/3 * u_f * \alpha_y) + (v_f * \alpha_x)))
443 : * + (n_y * ((u_f * \alpha_x) + ( 4/3 * v_f * \alpha_y)
444 : * + (w_f * \alpha_z)))
445 : * + (n_z * ((v_f * \alpha_z) + (-2/3 * w_f * \alpha_y)))
446 : * ]
447 : * * ((dv)_p / (d var)_p)
448 : * ]
449 : * + [ [ (n_x * ((-2/3 * u_f * \alpha_z) + (w_f * \alpha_x)))
450 : * + (n_y * ((-2/3 * v_f * \alpha_z) + (w_f * \alpha_y)))
451 : * + (n_z * ((u_f * \alpha_x) + (v_f * \alpha_y)
452 : * + ( 4/3 * w_f * \alpha_z)))
453 : * ]
454 : * * ((dw)_p / (d var)_p)
455 : * ]
456 : * ]
457 : *
458 : * normale sortante pour IPD=1, rentrante pour IPD=2
459 : SIGNOR=(-1.D0)**(IPD+1)
460 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
461 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
462 : CNX =MPNORM.VPOCHA(NLFACE,1)
463 : CNY =MPNORM.VPOCHA(NLFACE,2)
464 : CNZ =MPNORM.VPOCHA(NLFACE,3)
465 : ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
466 : ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
467 : ALPHAZ=KDUNDZ.VELCHE(IPP+1,IELEM)
468 : RHOP =MPROC.VPOCHA(NLCENP,1)
469 : UP =MPVITC.VPOCHA(NLCENP,1)
470 : VP =MPVITC.VPOCHA(NLCENP,2)
471 : WP =MPVITC.VPOCHA(NLCENP,3)
472 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU
473 : FACTDU=(CNX*((( 4.D0/3.D0)*UXF*ALPHAX)
474 : $ +(UYF*ALPHAY)
475 : $ +(UZF*ALPHAZ)))
476 : $ + (CNY*((UXF*ALPHAY)
477 : $ +((-2.D0/3.D0)*UYF*ALPHAX)))
478 : $ + (CNZ*((UXF*ALPHAZ)
479 : $ +((-2.D0/3.D0)*UZF*ALPHAX)))
480 : FACTDV=(CNX*(((-2.D0/3.D0)*UXF*ALPHAY)
481 : $ +(UYF*ALPHAX)))
482 : $ + (CNY*((UXF*ALPHAX)
483 : $ +(( 4.D0/3.D0)*UYF*ALPHAY)
484 : $ +(UZF*ALPHAZ)))
485 : $ + (CNZ*((UYF*ALPHAZ)
486 : $ +((-2.D0/3.D0)*UZF*ALPHAY)))
487 : FACTDW=(CNX*(((-2.D0/3.D0)*UXF*ALPHAZ)
488 : $ +(UZF*ALPHAX)))
489 : $ + (CNY*(((-2.D0/3.D0)*UYF*ALPHAZ)
490 : $ +(UZF*ALPHAY)))
491 : $ + (CNZ*((UXF*ALPHAX)
492 : $ +(UYF*ALPHAY)
493 : $ +(( 4.D0/3.D0)*UZF*ALPHAZ)))
494 : DUDRHO=-UP /RHOP
495 : DUDROU=1.D0/RHOP
496 : DVDRHO=-VP /RHOP
497 : DVDROV=1.D0/RHOP
498 : DWDRHO=-WP /RHOP
499 : DWDROW=1.D0/RHOP
500 : IF (.NOT.LMUR) THEN
501 : CCGRV1.POIPR2(IPP,IEL2)=NPPRIM
502 : RETRHO.VALMA2(IPD,IPP,IEL2)=
503 : $ FACTOR*( (FACTDU*DUDRHO)
504 : $ +(FACTDV*DVDRHO)
505 : $ +(FACTDW*DWDRHO))
506 : RETROU.VALMA2(IPD,IPP,IEL2)=
507 : $ FACTOR* (FACTDU*DUDROU)
508 : RETROV.VALMA2(IPD,IPP,IEL2)=
509 : $ FACTOR* (FACTDV*DVDROV)
510 : RETROW.VALMA2(IPD,IPP,IEL2)=
511 : $ FACTOR* (FACTDW*DWDROW)
512 : ELSE
513 : CCGRV1.POIPR1(IPP,IEL1)=NPPRIM
514 : RETRHO.VALMA1(IPD,IPP,IEL1)=
515 : $ FACTOR*( (FACTDU*DUDRHO)
516 : $ +(FACTDV*DVDRHO)
517 : $ +(FACTDW*DWDRHO))
518 : RETROU.VALMA1(IPD,IPP,IEL1)=
519 : $ FACTOR* (FACTDU*DUDROU)
520 : RETROV.VALMA1(IPD,IPP,IEL1)=
521 : $ FACTOR* (FACTDV*DVDROV)
522 : RETROW.VALMA1(IPD,IPP,IEL1)=
523 : $ FACTOR* (FACTDW*DWDROW)
524 : ENDIF
525 : ENDIF
526 : 124 CONTINUE
527 : 122 CONTINUE
528 : IF (.NOT.LMUR) THEN
529 : IEL2=IEL2+1
530 : ELSE
531 : IEL1=IEL1+1
532 : ENDIF
533 : ENDIF
534 : 12 CONTINUE
535 : NPP1=NPTEL-1
536 : NPP2=NPTEL-1
537 : NEL1=IEL1-1
538 : NEL2=IEL2-1
539 : SEGADJ RETRHO
540 : SEGADJ RETROU
541 : SEGADJ RETROV
542 : SEGADJ RETROW
543 : SEGADJ CCGRV1
544 : CCGRV1.LMATSI(**)=RETRHO
545 : CCGRV1.LMATSI(**)=RETROU
546 : CCGRV1.LMATSI(**)=RETROV
547 : CCGRV1.LMATSI(**)=RETROW
548 : * On accumule les matrices résultantes dans IJACO
548 : CALL AJMTK(CCGRV1,IJACO,IMPR,IRET)
549 : IF (IRET.NE.0) GOTO 9999
550 : SEGSUP RETRHO
551 : SEGSUP RETROU
552 : SEGSUP RETROV
553 : SEGSUP RETROW
554 : SEGSUP CCGRV1
555 : *
556 : SEGDES MCOGRV
557 : SEGDES KDUNDZ
558 : SEGDES KDUNDY
559 : SEGDES KDUNDX
560 : 1 CONTINUE
561 : SEGDES ICOGRV
562 : SEGDES MPGRVF
563 : SEGDES MPVITC
564 : SEGDES MPROC
565 : SEGDES MPVOLU
566 : SEGDES MPNORM
567 : SEGDES MPSURF
568 : SEGDES MELEFL
569 : SEGDES KRFACE
570 : SEGDES KRCENT
571 : SEGDES NOMINC
572 : IF (LCLITO) THEN
573 : SEGDES KRTOIM
574 : ENDIF
575 : IF (LCLIMV) THEN
576 : SEGDES KRVIMP
577 : SEGDES MPVIMP
578 : ENDIF
579 : *
580 : * Normal termination
581 : *
582 : IRET=0
583 : RETURN
584 : *
585 : * Format handling
586 : *
587 : *
588 : * Error handling
589 : *
590 : 9999 CONTINUE
591 : IRET=1
592 : WRITE(IOIMP,*) 'An error was detected in subroutine ylap2e'
593 : RETURN
594 : *
595 : * End of subroutine YLAP2E
596 : *
597 : END
598 :
599 :
600 :
601 :
602 :
603 :
604 :
605 :
606 :
© Cast3M 2003 - All rights reserved.
Disclaimer