1 : C XLAP1D SOURCE MAGN 09/07/30 21:21:23 6458
2 : SUBROUTINE XLAP1D(ICOGRV,MPGRVF,MPROC,MPVITC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,
5 : $ LCLIMV,KRVIMP,
6 : $ LCLITO,KRTOIM,MPTOIM,
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 : XLAP1D
15 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
16 : C VF 2D.
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 e_t )
24 : C Les contributions sont plus "compliquées" à calculer que
25 : C les simples (cf. xlap1c) car on a à dériver des produits
26 : C de fonctions de la vitesse.
27 : C xlap1e 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
38 : C LANGAGE : ESOPE
39 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
40 : C mél : gounand@semt2.smts.cea.fr
41 : C***********************************************************************
42 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
43 : C standard) à un objet de type MATRIK.
44 : C APPELE PAR : XLAP1A : Calcul de la matrice jacobienne du
45 : C résidu du laplacien VF 2D.
46 : C***********************************************************************
47 : C ENTREES : ICOGRV (type MCHELM) : coefficients pour le
48 : C calcul du gradient de la vitesse aux
49 : C interfaces.
50 : C MPGRVF (type MPOVAL) : gradient de la vitesse
51 : C aux interfaces.
52 : C MPROC (type MPOVAL) : masse volumique par
53 : C élément.
54 : C MPVITC (type MPOVAL) : vitesse par élément.
55 : C MPVOLU (type MPOVAL) : volume des éléments.
56 : C MPNORM (type MPOVAL) : normale aux faces.
57 : C MPSURF (type MPOVAL) : surface des faces.
58 : C MELEFL (type MELEME) : connectivités face-(centre
59 : C gauche, centre droit).
60 : C KRFACE (type MLENTI) : tableau de repérage dans
61 : C le maillage des faces des éléments.
62 : C KRCENT (type MLENTI) : tableau de repérage dans
63 : C le maillage des centres des éléments.
64 : C LCLIMV (type logique) : .TRUE. => CL de Dirichlet
65 : C sur la vitesse.
66 : C KRVIMP (type MLENTI) : tableau de repérage dans
67 : C maillage des CL de Dirichlet sur la
68 : C vitesse.
69 : C LCLITO (type logique) : .TRUE. => CL de Dirichlet
70 : C sur le tenseur des contraintes.
71 : C KRTOIM (type MLENTI) : tableau de repérage dans
72 : C maillage des CL de Dirichlet sur le tenseur des
73 : C contraintes
74 : C MPTOIM (type MPOVAL) : valeurs des CL de
75 : C Dirichlet sur le tenseur des contraintes.
76 : C NOMINC (type MLMOTS) : noms des inconnues.
77 : C MPMUC (type MPOVAL) : viscosité dynamique (SI).
78 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
79 : C résidu du laplacien VF 2D.
80 : C SORTIES : -
81 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
82 : C***********************************************************************
83 : C VERSION : v1, 10/12/2001, version initiale
84 : C HISTORIQUE : v1, 10/12/2001, création
85 : C HISTORIQUE :
86 : C HISTORIQUE :
87 : C***********************************************************************
88 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
89 : C en cas de modification de ce sous-programme afin de faciliter
90 : C la maintenance !
91 : C***********************************************************************
92 : -INC CCOPTIO
93 : -INC SMCOORD
94 : -INC SMCHPOI
95 : POINTEUR MPGRVF.MPOVAL
96 : POINTEUR MPMUC.MPOVAL
97 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
98 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
99 : POINTEUR MPTOIM.MPOVAL
100 : -INC SMCHAML
101 : POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
102 : POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL
103 : -INC SMELEME
104 : POINTEUR MELEFL.MELEME
105 : POINTEUR MCOGRV.MELEME
106 : -INC SMLENTI
107 : POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
108 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
109 : -INC SMLMOTS
110 : POINTEUR NOMINC.MLMOTS
111 : -INC SMMATRIK
112 : POINTEUR IJACO.MATRIK
113 : *
114 : * Objet matrice élémentaire simplifié
115 : *
116 : SEGMENT GMATSI
117 : INTEGER POIPR1(NPP1,NEL1)
118 : INTEGER POIDU1(1,NEL1)
119 : INTEGER POIPR2(NPP2,NEL2)
120 : INTEGER POIDU2(2,NEL2)
121 : POINTEUR LMATSI(0).MATSIM
122 : ENDSEGMENT
123 : * Contributions compliquées 2 de la part du gradient de
124 : * vitesse (CCGRV2)
125 : POINTEUR CCGRV2.GMATSI
126 : SEGMENT MATSIM
127 : CHARACTER*8 NOMPRI,NOMDUA
128 : REAL*8 VALMA1(1,NPP1,NEL1)
129 : REAL*8 VALMA2(2,NPP2,NEL2)
130 : ENDSEGMENT
131 : POINTEUR RETRHO.MATSIM
132 : POINTEUR RETROU.MATSIM
133 : POINTEUR RETROV.MATSIM
134 : *
135 : REAL*8 MU
136 : *
137 : INTEGER IMPR,IRET
138 : *
139 : LOGICAL LCLIMV,LCLITO
140 : LOGICAL LMUR
141 : LOGICAL LCTR2A,LCTR2B
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,NLCLV,NLFTOI
147 : INTEGER NPTEL
148 : INTEGER ICOORX,NLCGAU,NLCDRO
149 : *
150 : REAL*8 ALPHA,UMALPH
151 : REAL*8 DRD,DRG
152 : REAL*8 DUXXF,DUXYF,DUYXF,DUYYF
153 : REAL*8 XD,XF,XG,XFMXD,XFMXG
154 : REAL*8 YD,YF,YG,YFMYD,YFMYG
155 : REAL*8 ALPHAX,ALPHAY,CNX,CNY
156 : REAL*8 SIGNOR,SURFFA,VOLUEL
157 : REAL*8 RHOP,UP,VP
158 : REAL*8 FACTOR,FACTDU,FACTDV
159 : REAL*8 DUDRHO,DUDROU,DVDRHO,DVDROV
160 : REAL*8 DSTDU
161 : REAL*8 TAUXX,TAUXY,TAUYY
162 : REAL*8 EPSIXX,EPSIXY,EPSIYX,EPSIYY
163 : REAL*8 GAMMXD,GAMMXG,GAMMYD,GAMMYG
164 : *
165 : * Executable statements
166 : *
167 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap1d.eso'
168 : * On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
169 : * prenant successivement les valeurs :
170 : * \rho, \rho u, \rho v, \rho e_t
171 : * On dérive les termes : (\tens{\tau(\grad{u})} \prod \vect{u})
172 : * \pscal \vect{n}
173 : * ce qui donne deux contributions.
174 : * Contribution 1 :
175 : * ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
176 : * Contribution 2 :
177 : * ( \tens{\tau} \prod (d \vect{u} / d var)) \pscal \vect{n}
178 : * Note :
179 : * pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
180 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
181 : * RETRHO, RETROU, RETROV
182 : IF (LCLIMV) THEN
183 : SEGACT KRVIMP
184 : ENDIF
185 : IF (LCLITO) THEN
186 : SEGACT KRTOIM
187 : SEGACT MPTOIM
188 : ENDIF
189 : SEGACT NOMINC
190 : SEGACT KRCENT
191 : SEGACT KRFACE
192 : SEGACT MELEFL
193 : SEGACT MPSURF
194 : SEGACT MPNORM
195 : SEGACT MPVOLU
196 : SEGACT MPMUC
197 : SEGACT MPROC
198 : SEGACT MPVITC
199 : SEGACT MPGRVF
200 : SEGACT ICOGRV
201 : NSOUCH=ICOGRV.IMACHE(/1)
202 : DO 1 ISOUCH=1,NSOUCH
203 : MCOGRV=ICOGRV.IMACHE(ISOUCH)
204 : JCOGRV=ICOGRV.ICHAML(ISOUCH)
205 : SEGACT JCOGRV
206 : KDUNDX=JCOGRV.IELVAL(1)
207 : KDUNDY=JCOGRV.IELVAL(2)
208 : SEGDES JCOGRV
209 : SEGACT KDUNDX
210 : SEGACT KDUNDY
211 : SEGACT MCOGRV
212 : NELEM=MCOGRV.NUM(/2)
213 : NPTEL=MCOGRV.NUM(/1)
214 : NPP1=NPTEL
215 : NPP2=NPTEL+1
216 : NEL1=NELEM
217 : NEL2=NELEM
218 : IEL1=1
219 : IEL2=1
220 : SEGINI RETRHO
221 : SEGINI RETROU
222 : SEGINI RETROV
223 : SEGINI CCGRV2
224 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
225 : RETRHO.NOMPRI(5:8)=' '
226 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
227 : RETRHO.NOMDUA(5:8)=' '
228 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
229 : RETROU.NOMPRI(5:8)=' '
230 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(4)
231 : RETROU.NOMDUA(5:8)=' '
232 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
233 : RETROV.NOMPRI(5:8)=' '
234 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(4)
235 : RETROV.NOMDUA(5:8)=' '
236 : DO 12 IELEM=1,NELEM
237 : * Le premier point du support de ICOGRV est un point FACE
238 : NGFACE=MCOGRV.NUM(1,IELEM)
239 : NLFACE=KRFACE.LECT(NGFACE)
240 : IF (NLFACE.EQ.0) THEN
241 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
242 : GOTO 9999
243 : ENDIF
244 : * On calcule la contribution 2 à la matrice jacobienne IJACO de la
245 : * face NGFACE
246 : * (points duaux : centres à gauche et à droite de la face)
247 : * (points primaux : une partie (bicoz conditions aux limites)
248 : * de ceux du stencil pour le calcul du gradient
249 : * à la face, ils doivent être des points centres
250 : * ET les centres à gauche et à droite qui servent pour le
251 : * calcul de la vitesse sur la face)
252 : * Si la vitesse sur la face est imposée par les
253 : * conditions aux limites, la contribution 2 de la face à IJACO est
254 : * nulle.
255 : LCTR2A=.TRUE.
256 : IF (LCLIMV) THEN
257 : NLCLV=KRVIMP.LECT(NGFACE)
258 : IF (NLCLV.NE.0) THEN
259 : LCTR2A=.FALSE.
260 : ENDIF
261 : ENDIF
262 : IF (LCTR2A) THEN
263 : NGCGAU=MELEFL.NUM(1,NLFACE)
264 : NGCDRO=MELEFL.NUM(3,NLFACE)
265 : NLCGAU=KRCENT.LECT(NGCGAU)
266 : NLCDRO=KRCENT.LECT(NGCDRO)
267 : LMUR=(NGCGAU.EQ.NGCDRO)
268 : * On calcule ensuite les valeurs des coefficients qui servent pour
269 : * interpoler la vitesse sur la face NGFACE à partir des vitesses
270 : * gauche (et droite si non mur) et des dérivées de la vitesse sur la
271 : * face NGFACE (de la même manière que dans xlap12.eso)
272 : IF (LMUR) THEN
273 : * Parametres geometriques
274 : ICOORX=((IDIM+1)*(NGFACE-1))+1
275 : XF=MCOORD.XCOOR(ICOORX)
276 : YF=MCOORD.XCOOR(ICOORX+1)
277 : ICOORX=((IDIM+1)*(NGCGAU-1))+1
278 : XG=MCOORD.XCOOR(ICOORX)
279 : YG=MCOORD.XCOOR(ICOORX+1)
280 : XFMXG=XF-XG
281 : YFMYG=YF-YG
282 : ALPHA=0.0D0
283 : UMALPH=1.0D0
284 : GAMMXG=1.D0
285 : EPSIXX=XFMXG
286 : EPSIXY=YFMYG
287 : GAMMYG=1.D0
288 : EPSIYX=XFMXG
289 : EPSIYY=YFMYG
290 : ELSEIF (.NOT.LMUR) THEN
291 : ICOORX=((IDIM+1)*(NGFACE-1))+1
292 : XF=MCOORD.XCOOR(ICOORX)
293 : YF=MCOORD.XCOOR(ICOORX+1)
294 : ICOORX=((IDIM+1)*(NGCGAU-1))+1
295 : XG=MCOORD.XCOOR(ICOORX)
296 : YG=MCOORD.XCOOR(ICOORX+1)
297 : XFMXG=XF-XG
298 : YFMYG=YF-YG
299 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
300 : ICOORX=((IDIM+1)*(NGCDRO-1))+1
301 : XD=MCOORD.XCOOR(ICOORX)
302 : YD=MCOORD.XCOOR(ICOORX+1)
303 : XFMXD=XF-XD
304 : YFMYD=YF-YD
305 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
306 : ALPHA=DRG/(DRG+DRD)
307 : UMALPH=1.0D0-ALPHA
308 : GAMMXG=UMALPH
309 : GAMMXD=ALPHA
310 : EPSIXX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
311 : EPSIXY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
312 : GAMMYG=UMALPH
313 : GAMMYD=ALPHA
314 : EPSIYX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
315 : EPSIYY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
316 : ELSE
317 : WRITE(IOIMP,*) 'Erreur de programmation n°2'
318 : GOTO 9999
319 : ENDIF
320 : * On distingue le cas où la face est un bord du maillage (mur)
321 : * du cas où la face est interne au maillage
322 : IF (.NOT.LMUR) THEN
323 : NPD=2
324 : NPP=NPTEL+1
325 : ELSE
326 : NPD=1
327 : NPP=NPTEL
328 : ENDIF
329 : MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) +
330 : & ALPHA*MPMUC.VPOCHA(NLCDRO,1)
331 : * On calcule tout d'abord les valeurs du tenseur des contraintes sur
332 : * la face NGFACE de la même manière que dans xlap12.eso
333 : IF (LCLITO) THEN
334 : NLFTOI=KRTOIM.LECT(NGFACE)
335 : ELSE
336 : NLFTOI=0
337 : ENDIF
338 : IF (NLFTOI.NE.0) THEN
339 : TAUXX=MPTOIM.VPOCHA(NLFTOI,1)
340 : TAUYY=MPTOIM.VPOCHA(NLFTOI,2)
341 : TAUXY=MPTOIM.VPOCHA(NLFTOI,3)
342 : ELSE
343 : DUXXF=MPGRVF.VPOCHA(NLFACE,1)
344 : DUXYF=MPGRVF.VPOCHA(NLFACE,2)
345 : DUYXF=MPGRVF.VPOCHA(NLFACE,3)
346 : DUYYF=MPGRVF.VPOCHA(NLFACE,4)
347 : DSTDU=(2.0D0/3.0D0)*(DUXXF+DUYYF)
348 : TAUXX=MU*((2.0D0*DUXXF)-DSTDU)
349 : TAUXY=MU*(DUXYF+DUYXF)
350 : TAUYY=MU*((2.0D0*DUYYF)-DSTDU)
351 : ENDIF
352 : * IPD=1 : point à gauche du point NGFACE
353 : * IPD=2 : point à droite du point NGFACE
354 : DO 122 IPD=1,NPD
355 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
356 : IF (.NOT.LMUR) THEN
357 : CCGRV2.POIDU2(IPD,IEL2)=NPDUAL
358 : ELSE
359 : CCGRV2.POIDU1(IPD,IEL1)=NPDUAL
360 : ENDIF
361 : NLCEND=KRCENT.LECT(NPDUAL)
362 : IF (NLCEND.EQ.0) THEN
363 : WRITE(IOIMP,*) 'Erreur grave n°1'
364 : GOTO 9999
365 : ENDIF
366 : DO 124 IPP=1,NPP
367 : IF (IPP.EQ.NPTEL) THEN
368 : NPPRIM=NGCGAU
369 : ELSEIF (IPP.EQ.NPTEL+1) THEN
370 : NPPRIM=NGCDRO
371 : ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
372 : NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
373 : ELSE
374 : WRITE(IOIMP,*) 'Erreur grave n°2'
375 : GOTO 9999
376 : ENDIF
377 : LCTR2B=.TRUE.
378 : IF (LCLIMV) THEN
379 : NLCLV=KRVIMP.LECT(NPPRIM)
380 : IF (NLCLV.NE.0) THEN
381 : LCTR2B=.FALSE.
382 : ENDIF
383 : ENDIF
384 : IF (.NOT.LCTR2B) THEN
385 : * Lorsque une contribution est nulle, on fixe artificiellement le
386 : * point primal égal au point dual.
387 : IF (.NOT.LMUR) THEN
388 : CCGRV2.POIPR2(IPP,IEL2)=NPDUAL
389 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
390 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
391 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
392 : ELSE
393 : CCGRV2.POIPR1(IPP,IEL1)=NPDUAL
394 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
395 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
396 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
397 : ENDIF
398 : ELSE
399 : * Les contributions 2 valent :
400 : * (d Res_{\rho e_t})_d / (d var)_p =
401 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f)
402 : * * [ [ [ ((n_x * \tau_xx) + (n_y * \tau_yx))
403 : * * ((\epsilon_xx * \alpha_x) + (\epsilon_xy * \alpha_y))
404 : * ]
405 : * * ((du)_p / (d var)_p)
406 : * ]
407 : * + [ [ ((n_x * \tau_xy) + (n_y * \tau_yy))
408 : * * ((\epsilon_yx * \alpha_x) + (\epsilon_yy * \alpha_y))
409 : * ]
410 : * * ((dv)_p / (d var)_p)
411 : * ]
412 : * ]
413 : * et (m étant le centre de gauche dans le cas mur
414 : * et les centre gauche, puis droite dans le cas normal)
415 : * (d Res_{\rho e_t})_d / (d var)_m =
416 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f)
417 : * * [ [ [ ((n_x * \tau_xx) + (n_y * \tau_yx))
418 : * * (\gamma_x)
419 : * ]
420 : * * ((du)_p / (d var)_p)
421 : * ]
422 : * + [ [ ((n_x * \tau_xy) + (n_y * \tau_yy))
423 : * * (\gamma_y)
424 : * ]
425 : * * ((dv)_p / (d var)_p)
426 : * ]
427 : * ]
428 : * sachant que l'expression donnant l'interpolation de la
429 : * vitesse (u,v) sur une face i est de la forme :
430 : * u_i = (\gamma_x,gauche * u_gauche)
431 : * + (\gamma_x,droite * u_droite) (=0 dans le cas mur)
432 : * + (\epsilon_xx,i * (du/dx)_i)
433 : * + (\epsilon_xy,i * (du/dy)_i)
434 : * v_i = (\gamma_y,gauche * v_gauche)
435 : * + (\gamma_y,droite * v_droite) (=0 dans le cas mur)
436 : * + (\epsilon_yx,i * (dv/dx)_i)
437 : * + (\epsilon_yy,i * (dv/dy)_i)
438 : * (voir dans xlap12.eso et ci-dessus pour les valeurs des coeffs
439 : * \gamma et \epsilon)
440 : NLCENP=KRCENT.LECT(NPPRIM)
441 : IF (NLCENP.EQ.0) THEN
442 : WRITE(IOIMP,*) 'Erreur grave n°3'
443 : GOTO 9999
444 : ENDIF
445 : * normale sortante pour IPD=1, rentrante pour IPD=2
446 : SIGNOR=(-1.D0)**(IPD+1)
447 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
448 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
449 : CNX =MPNORM.VPOCHA(NLFACE,1)
450 : CNY =MPNORM.VPOCHA(NLFACE,2)
451 : RHOP =MPROC.VPOCHA(NLCENP,1)
452 : UP =MPVITC.VPOCHA(NLCENP,1)
453 : VP =MPVITC.VPOCHA(NLCENP,2)
454 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA
455 : IF (IPP.EQ.NPTEL) THEN
456 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
457 : $ * GAMMXG
458 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
459 : $ * GAMMYG
460 : ELSEIF (IPP.EQ.NPTEL+1) THEN
461 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
462 : $ * GAMMXD
463 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
464 : $ * GAMMYD
465 : ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
466 : ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
467 : ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
468 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
469 : $ * ((EPSIXX*ALPHAX)+(EPSIXY*ALPHAY))
470 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
471 : $ * ((EPSIYX*ALPHAX)+(EPSIYY*ALPHAY))
472 : ELSE
473 : WRITE(IOIMP,*) 'Erreur grave n°3'
474 : GOTO 9999
475 : ENDIF
476 : DUDRHO=-UP /RHOP
477 : DUDROU=1.D0/RHOP
478 : DVDRHO=-VP /RHOP
479 : DVDROV=1.D0/RHOP
480 : IF (.NOT.LMUR) THEN
481 : CCGRV2.POIPR2(IPP,IEL2)=NPPRIM
482 : RETRHO.VALMA2(IPD,IPP,IEL2)=
483 : $ FACTOR*( (FACTDU*DUDRHO)
484 : $ +(FACTDV*DVDRHO))
485 : RETROU.VALMA2(IPD,IPP,IEL2)=
486 : $ FACTOR* (FACTDU*DUDROU)
487 : RETROV.VALMA2(IPD,IPP,IEL2)=
488 : $ FACTOR* (FACTDV*DVDROV)
489 : ELSE
490 : CCGRV2.POIPR1(IPP,IEL1)=NPPRIM
491 : RETRHO.VALMA1(IPD,IPP,IEL1)=
492 : $ FACTOR*( (FACTDU*DUDRHO)
493 : $ +(FACTDV*DVDRHO))
494 : RETROU.VALMA1(IPD,IPP,IEL1)=
495 : $ FACTOR* (FACTDU*DUDROU)
496 : RETROV.VALMA1(IPD,IPP,IEL1)=
497 : $ FACTOR* (FACTDV*DVDROV)
498 : ENDIF
499 : ENDIF
500 : 124 CONTINUE
501 : 122 CONTINUE
502 : IF (.NOT.LMUR) THEN
503 : IEL2=IEL2+1
504 : ELSE
505 : IEL1=IEL1+1
506 : ENDIF
507 : ENDIF
508 : 12 CONTINUE
509 : NPP1=NPTEL
510 : NPP2=NPTEL+1
511 : NEL1=IEL1-1
512 : NEL2=IEL2-1
513 : SEGADJ RETRHO
514 : SEGADJ RETROU
515 : SEGADJ RETROV
516 : SEGADJ CCGRV2
517 : CCGRV2.LMATSI(**)=RETRHO
518 : CCGRV2.LMATSI(**)=RETROU
519 : CCGRV2.LMATSI(**)=RETROV
520 : * On accumule les matrices résultantes dans IJACO
520 : CALL AJMTK(CCGRV2,IJACO,IMPR,IRET)
521 : IF (IRET.NE.0) GOTO 9999
522 : SEGSUP RETRHO
523 : SEGSUP RETROU
524 : SEGSUP RETROV
525 : SEGSUP CCGRV2
526 : *
527 : SEGDES MCOGRV
528 : SEGDES KDUNDY
529 : SEGDES KDUNDX
530 : 1 CONTINUE
531 : SEGDES ICOGRV
532 : SEGDES MPGRVF
533 : SEGDES MPVITC
534 : SEGDES MPMUC
535 : SEGDES MPROC
536 : SEGDES MPVOLU
537 : SEGDES MPNORM
538 : SEGDES MPSURF
539 : SEGDES MELEFL
540 : SEGDES KRFACE
541 : SEGDES KRCENT
542 : SEGDES NOMINC
543 : IF (LCLITO) THEN
544 : SEGDES KRTOIM
545 : SEGDES MPTOIM
546 : ENDIF
547 : IF (LCLIMV) THEN
548 : SEGDES KRVIMP
549 : ENDIF
550 : *
551 : * Normal termination
552 : *
553 : IRET=0
554 : RETURN
555 : *
556 : * Format handling
557 : *
558 : *
559 : * Error handling
560 : *
561 : 9999 CONTINUE
562 : IRET=1
563 : WRITE(IOIMP,*) 'An error was detected in subroutine xlap1d'
564 : RETURN
565 : *
566 : * End of subroutine XLAP1D
567 : *
568 : END
569 :
570 :
571 :
572 :
573 :
574 :
575 :
576 :
© Cast3M 2003 - All rights reserved.
Disclaimer