OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sptemp.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| spgradt ../engine/source/elements/sph/sptemp.F
25!||--- called by ------------------------------------------------------
26!|| forintp ../engine/source/elements/forintp.F
27!||--- calls -----------------------------------------------------
28!|| weight1 ../engine/source/elements/sph/weight.F
29!||--- uses -----------------------------------------------------
30!|| sphbox ../engine/share/modules/sphbox.F
31!||====================================================================
32 SUBROUTINE spgradt(
33 1 X, MS, SPBUF, KXSP,
34 2 IXSP, NOD2SP, ISPSYM, XSPSYM,
35 3 WA, WACOMP, WTEMP, WTR,
36 4 WGRADT, LFT, LLT, NFT)
37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE sphbox
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "sphcom.inc"
49C-----------------------------------------------
50C D u m m y A r g u m e n t s
51C-----------------------------------------------
52 INTEGER, INTENT(INOUT) :: LFT
53 INTEGER, INTENT(INOUT) :: LLT
54 INTEGER, INTENT(INOUT) :: NFT
55 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
56 . ispsym(nspcond,*)
57 my_real
58 . x(3,*) ,ms(*) ,
59 . spbuf(nspbuf,*) ,xspsym(3,*) ,
60 . wa(kwasph,*) ,wacomp(16,*),
61 . wtemp(*), wtr(*), wgradt(3,*)
62C-----------------------------------------------
63C L o c a l V a r i a b l e s
64C-----------------------------------------------
65 INTEGER I,N,INOD,JNOD,J,NVOIS,M,
66 . NVOISS,SM,JS,NC,NS,NN
67 my_real
68 . XI,YI,ZI,DI,RHOI,XJ,YJ,ZJ,DJ,RHOJ,DIJ,
69 . vj,wght,wgrad(3),
70 . alphai,alphaxi,alphayi,alphazi,
71 . betaxi,betayi,betazi,
72 . betaxxi,betayxi,betazxi,
73 . betaxyi,betayyi,betazyi,
74 . betaxzi,betayzi,betazzi,
75 . betax,wgrdx,wgrdy,wgrdz,
76 . ti,tj
77C-------------------------------------------
78 DO i=lft,llt
79 n =nft+i
80 wgradt(1,n)=zero
81 wgradt(2,n)=zero
82 wgradt(3,n)=zero
83 ENDDO
84C-----------------------------------------------
85C Calcul du gradient de temperature.
86C-----------------------------------------------
87 DO 10 i=lft,llt
88 n =nft+i
89 IF(kxsp(2,n)<=0)GOTO 10
90 inod =kxsp(3,n)
91 xi=x(1,inod)
92 yi=x(2,inod)
93 zi=x(3,inod)
94 di=spbuf(1,n)
95 ti=wtemp(n)
96 rhoi =wa(10,n)
97C-----
98 alphai=wacomp(1,n)
99C BETAXI=WACOMP(2,N)
100C BETAYI=WACOMP(3,N)
101C BETAZI=WACOMP(4,N)
102 alphaxi=wacomp( 5,n)
103 alphayi=wacomp( 6,n)
104 alphazi=wacomp( 7,n)
105 betaxxi=wacomp( 8,n)
106 betayxi=wacomp( 9,n)
107 betazxi=wacomp(10,n)
108 betaxyi=wacomp(11,n)
109 betayyi=wacomp(12,n)
110 betazyi=wacomp(13,n)
111 betaxzi=wacomp(14,n)
112 betayzi=wacomp(15,n)
113 betazzi=wacomp(16,n)
114C------
115 nvois=kxsp(4,n)
116 DO j=1,nvois
117 jnod=ixsp(j,n)
118 IF(jnod>0)THEN
119 m=nod2sp(jnod)
120 xj=x(1,jnod)
121 yj=x(2,jnod)
122 zj=x(3,jnod)
123 dj=spbuf(1,m)
124 dij=0.5*(di+dj)
125 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
126 rhoj=wa(10,m)
127 vj=spbuf(12,m)/max(em20,rhoj)
128 tj=wtemp(m)
129 ELSE
130 nn = -jnod
131 xj=xsphr(3,nn)
132 yj=xsphr(4,nn)
133 zj=xsphr(5,nn)
134 dj=xsphr(2,nn)
135 dij=0.5*(di+dj)
136 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
137 rhoj=xsphr(7,nn)
138 vj=xsphr(8,nn)/max(em20,rhoj)
139 tj=wtr(nn)
140 END IF
141 wgrdx=wgrad(1)*alphai+wght*alphaxi
142 . +wgrad(1)*betaxxi+wgrad(2)*betaxyi+wgrad(3)*betaxzi
143 wgrdy=wgrad(2)*alphai+wght*alphayi
144 . +wgrad(1)*betayxi+wgrad(2)*betayyi+wgrad(3)*betayzi
145 wgrdz=wgrad(3)*alphai+wght*alphazi
146 . +wgrad(1)*betazxi+wgrad(2)*betazyi+wgrad(3)*betazzi
147! Old order1 correction
148! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
149! WGRDX=WGRAD(1)*ALPHAI*BETAX
150! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
151! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
152! WGRDY=WGRAD(2)*ALPHAI*BETAX
153! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
154! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
155! WGRDZ=WGRAD(3)*ALPHAI*BETAX
156! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
157! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
158 wgrad(1)=wgrdx
159 wgrad(2)=wgrdy
160 wgrad(3)=wgrdz
161C
162 wgradt(1,n)=wgradt(1,n)+vj*(tj-ti)*wgrad(1)
163 wgradt(2,n)=wgradt(2,n)+vj*(tj-ti)*wgrad(2)
164 wgradt(3,n)=wgradt(3,n)+vj*(tj-ti)*wgrad(3)
165C--------
166 END DO
167C------
168C partie symetrique.
169 nvoiss=kxsp(6,n)
170 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
171 js=ixsp(j,n)
172 IF(js>0)THEN
173 sm=js/(nspcond+1)
174 nc=mod(js,nspcond+1)
175 js=ispsym(nc,sm)
176 xj =xspsym(1,js)
177 yj =xspsym(2,js)
178 zj =xspsym(3,js)
179 dj =spbuf(1,sm)
180 dij =half*(di+dj)
181 rhoj=wa(10,sm)
182 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
183 jnod=kxsp(3,sm)
184 vj=spbuf(12,sm)/max(em20,rhoj)
185 tj=wtemp(sm)
186 ELSE
187 sm=-js/(nspcond+1)
188 nc=mod(-js,nspcond+1)
189 js=ispsymr(nc,sm)
190 xj =xspsym(1,js)
191 yj =xspsym(2,js)
192 zj =xspsym(3,js)
193 dj =xsphr(2,sm)
194 dij =half*(di+dj)
195 rhoj=xsphr(7,sm)
196 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
197 jnod=kxsp(3,sm)
198 vj=xsphr(8,sm)/max(em20,rhoj)
199 tj=wtr(sm)
200 END IF
201C
202 wgrdx=wgrad(1)*alphai+wght*alphaxi
203 . +wgrad(1)*betaxxi+wgrad(2)*betaxyi+wgrad(3)*betaxzi
204 wgrdy=wgrad(2)*alphai+wght*alphayi
205 . +wgrad(1)*betayxi+wgrad(2)*betayyi+wgrad(3)*betayzi
206 wgrdz=wgrad(3)*alphai+wght*alphazi
207 . +wgrad(1)*betazxi+wgrad(2)*betazyi+wgrad(3)*betazzi
208! Old order1 correction
209! BETAX=ONE + BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
210! WGRDX=WGRAD(1)*ALPHAI*BETAX
211! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
212! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
213! WGRDY=WGRAD(2)*ALPHAI*BETAX
214! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
215! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
216! WGRDZ=WGRAD(3)*ALPHAI*BETAX
217! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
218! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
219 wgrad(1)=wgrdx
220 wgrad(2)=wgrdy
221 wgrad(3)=wgrdz
222C
223 wgradt(1,n)=wgradt(1,n)+vj*(tj-ti)*wgrad(1)
224 wgradt(2,n)=wgradt(2,n)+vj*(tj-ti)*wgrad(2)
225 wgradt(3,n)=wgradt(3,n)+vj*(tj-ti)*wgrad(3)
226 END DO
227C------
228 10 CONTINUE
229C-----------------------------------------------
230 RETURN
231 END
232!||====================================================================
233!|| splaplt ../engine/source/elements/sph/sptemp.F
234!||--- called by ------------------------------------------------------
235!|| forintp ../engine/source/elements/forintp.F
236!||--- calls -----------------------------------------------------
237!|| weight1 ../engine/source/elements/sph/weight.F
238!||--- uses -----------------------------------------------------
239!|| sphbox ../engine/share/modules/sphbox.F
240!||====================================================================
241 SUBROUTINE splaplt(
242 1 X, MS, SPBUF, KXSP,
243 2 IXSP, NOD2SP, ISPSYM, XSPSYM,
244 3 WA, WACOMP, WGRADT, WGR,
245 4 WGRADTSM,WLAPLT, WSMCOMP, LAMBDA,
246 5 LAMBDR, LFT, LLT, NFT)
247C-----------------------------------------------
248C M o d u l e s
249C-----------------------------------------------
250 USE sphbox
251C-----------------------------------------------
252C I m p l i c i t T y p e s
253C-----------------------------------------------
254#include "implicit_f.inc"
255C-----------------------------------------------
256C C o m m o n B l o c k s
257C-----------------------------------------------
258#include "sphcom.inc"
259C-----------------------------------------------
260C D u m m y A r g u m e n t s
261C-----------------------------------------------
262 INTEGER, INTENT(INOUT) :: LFT
263 INTEGER, INTENT(INOUT) :: LLT
264 INTEGER, INTENT(INOUT) :: NFT
265 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
266 . ISPSYM(NSPCOND,*)
267 my_real
268 . X(3,*) ,MS(*) ,
269 . SPBUF(NSPBUF,*) ,XSPSYM(3,*) ,
270 . WA(KWASPH,*) ,WACOMP(16,*),
271 . wgradt(3,*),wgr(3,*),wgradtsm(3,*),
272 . wlaplt(*),wsmcomp(6,*),lambda(*),lambdr(*)
273C-----------------------------------------------
274C L o c a l V a r i a b l e s
275C-----------------------------------------------
276 INTEGER I,N,INOD,JNOD,J,NVOIS,M,
277 . NVOISS,SM,JS,NC,NS,NN
278 my_real
279 . XI,YI,ZI,DI,RHOI,XJ,YJ,ZJ,DJ,RHOJ,DIJ,
280 . VJ,WGHT,WGRAD(3),
281 . ALPHAI,ALPHAXI,ALPHAYI,ALPHAZI,
282 . BETAXI,BETAYI,BETAZI,
283 . BETAXXI,BETAYXI,BETAZXI,
284 . BETAXYI,BETAYYI,BETAZYI,
285 . betaxzi,betayzi,betazzi,
286 . alphaj,alphaxj,alphayj,alphazj,
287 . betaxj,betayj,betazj,
288 . betaxxj,betayxj,betazxj,
289 . betaxyj,betayyj,betazyj,
290 . betaxzj,betayzj,betazzj,
291 . betax,wgrdx,wgrdy,wgrdz,wgrd(3),
292 . gradtxi,gradtyi,gradtzi,
293 . gradtxj,gradtyj,gradtzj
294C-------------------------------------------
295 DO i=lft,llt
296 n =nft+i
297 wlaplt(n)=zero
298 ENDDO
299C-----------------------------------------------
300C Calcul du Laplacien(T).
301C-----------------------------------------------
302 DO 10 i=lft,llt
303 n =nft+i
304 IF(kxsp(2,n)<=0)GOTO 10
305 inod =kxsp(3,n)
306 xi=x(1,inod)
307 yi=x(2,inod)
308 zi=x(3,inod)
309 di=spbuf(1,n)
310 gradtxi=wgradt(1,n)
311 gradtyi=wgradt(2,n)
312 gradtzi=wgradt(3,n)
313 rhoi =wa(10,n)
314C-----
315 alphai=wacomp(1,n)
316C BETAXI=WACOMP(2,N)
317C BETAYI=WACOMP(3,N)
318C BETAZI=WACOMP(4,N)
319 alphaxi=wacomp( 5,n)
320 alphayi=wacomp( 6,n)
321 alphazi=wacomp( 7,n)
322 betaxxi=wacomp( 8,n)
323 betayxi=wacomp( 9,n)
324 betazxi=wacomp(10,n)
325 betaxyi=wacomp(11,n)
326 betayyi=wacomp(12,n)
327 betazyi=wacomp(13,n)
328 betaxzi=wacomp(14,n)
329 betayzi=wacomp(15,n)
330 betazzi=wacomp(16,n)
331C------
332 nvois=kxsp(4,n)
333 DO j=1,nvois
334 jnod=ixsp(j,n)
335 IF(jnod>0)THEN
336 m=nod2sp(jnod)
337 xj=x(1,jnod)
338 yj=x(2,jnod)
339 zj=x(3,jnod)
340 dj=spbuf(1,m)
341 dij=0.5*(di+dj)
342 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
343 rhoj=wa(10,m)
344 vj=spbuf(12,m)/max(em20,rhoj)
345 gradtxj=wgradt(1,m)
346 gradtyj=wgradt(2,m)
347 gradtzj=wgradt(3,m)
348 wgrdx=wgrad(1)
349 wgrdy=wgrad(2)
350 wgrdz=wgrad(3)
351 wgrad(1)=wgrdx*alphai+wght*alphaxi
352 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
353 wgrad(2)=wgrdy*alphai+wght*alphayi
354 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
355 wgrad(3)=wgrdz*alphai+wght*alphazi
356 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
357! Old order1 correction
358! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
359! WGRAD(1)=WGRDX*ALPHAI*BETAX
360! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
361! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
362! WGRAD(2)=WGRDY*ALPHAI*BETAX
363! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
364! . (betaxyi*(xi-xj)+betayyi*(yi-yj)+betazyi*(zi-zj)+betayi))
365! WGRAD(3)=WGRDZ*ALPHAI*BETAX
366! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
367! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
368C----------
369C noyau conjugue Grad[Wa(b)]
370 alphaj=wacomp(1,m)
371C BETAXJ=WACOMP(2,M)
372C BETAYJ=WACOMP(3,M)
373C BETAZJ=WACOMP(4,M)
374 alphaxj=wacomp( 5,m)
375 alphayj=wacomp( 6,m)
376 alphazj=wacomp( 7,m)
377 betaxxj=wacomp( 8,m)
378 betayxj=wacomp( 9,m)
379 betazxj=wacomp(10,m)
380 betaxyj=wacomp(11,m)
381 betayyj=wacomp(12,m)
382 betazyj=wacomp(13,m)
383 betaxzj=wacomp(14,m)
384 betayzj=wacomp(15,m)
385 betazzj=wacomp(16,m)
386C
387 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
388 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
389 wgrd(2)=-wgrdy*alphaj+wght*alphayj
390 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
391 wgrd(3)=-wgrdz*alphaj+wght*alphazj
392 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
393! Old order1 correction
394! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
395! WGRD(1)=-WGRDX*ALPHAJ*BETAX
396! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
397! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
398! WGRD(2)=-WGRDY*ALPHAJ*BETAX
399! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
400! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
401! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
402! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
403! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
404C
405 wlaplt(n)=wlaplt(n)+vj*(
406 . -lambda(m)*(gradtxj*wgrd(1)+gradtyj*wgrd(2)+gradtzj*wgrd(3))
407 . +lambda(n)*(gradtxi*wgrad(1)+gradtyi*wgrad(2)+gradtzi*wgrad(3)))
408C--------
409 ELSE
410 nn = -jnod
411 xj=xsphr(3,nn)
412 yj=xsphr(4,nn)
413 zj=xsphr(5,nn)
414 dj=xsphr(2,nn)
415 dij=0.5*(di+dj)
416 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
417 rhoj=xsphr(7,nn)
418 vj=xsphr(8,nn)/max(em20,rhoj)
419 gradtxj=wgr(1,nn)
420 gradtyj=wgr(2,nn)
421 gradtzj=wgr(3,nn)
422 wgrdx=wgrad(1)
423 wgrdy=wgrad(2)
424 wgrdz=wgrad(3)
425 wgrad(1)=wgrdx*alphai+wght*alphaxi
426 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
427 wgrad(2)=wgrdy*alphai+wght*alphayi
428 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
429 wgrad(3)=wgrdz*alphai+wght*alphazi
430 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
431! Old order1 correction
432! BETAX=1.+BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
433! WGRAD(1)=WGRDX*ALPHAI*BETAX
434! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
435! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
436! WGRAD(2)=WGRDY*ALPHAI*BETAX
437! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
438! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
439! WGRAD(3)=WGRDZ*ALPHAI*BETAX
440! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
441! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
442C----------
443C noyau conjugue Grad[Wa(b)]
444 alphaj=wacompr(1,nn)
445C BETAXJ=WACOMPR(2,NN)
446C BETAYJ=WACOMPR(3,NN)
447C BETAZJ=WACOMPR(4,NN)
448 alphaxj=wacompr( 5,nn)
449 alphayj=wacompr( 6,nn)
450 alphazj=wacompr( 7,nn)
451 betaxxj=wacompr( 8,nn)
452 betayxj=wacompr( 9,nn)
453 betazxj=wacompr(10,nn)
454 betaxyj=wacompr(11,nn)
455 betayyj=wacompr(12,nn)
456 betazyj=wacompr(13,nn)
457 betaxzj=wacompr(14,nn)
458 betayzj=wacompr(15,nn)
459 betazzj=wacompr(16,nn)
460C
461 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
462 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
463 wgrd(2)=-wgrdy*alphaj+wght*alphayj
464 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
465 wgrd(3)=-wgrdz*alphaj+wght*alphazj
466 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
467! Old order1 correction
468! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
469! WGRD(1)=-WGRDX*ALPHAJ*BETAX
470! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
471! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
472! WGRD(2)=-WGRDY*ALPHAJ*BETAX
473! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
474! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
475! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
476! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
477! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
478C
479 wlaplt(n)=wlaplt(n)+vj*(
480 . -lambdr(nn)*(gradtxj*wgrd(1)+gradtyj*wgrd(2)+gradtzj*wgrd(3))
481 . +lambda(n)*(gradtxi*wgrad(1)+gradtyi*wgrad(2)+gradtzi*wgrad(3)))
482 END IF
483C--------
484 END DO
485C------
486C partie symetrique.
487 nvoiss=kxsp(6,n)
488 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
489 js=ixsp(j,n)
490 IF(js>0)THEN
491 sm=js/(nspcond+1)
492 nc=mod(js,nspcond+1)
493 js=ispsym(nc,sm)
494 xj =xspsym(1,js)
495 yj =xspsym(2,js)
496 zj =xspsym(3,js)
497 dj =spbuf(1,sm)
498 dij =half*(di+dj)
499 rhoj=wa(10,sm)
500 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
501 jnod=kxsp(3,sm)
502 vj=spbuf(12,sm)/max(em20,rhoj)
503 gradtxj=wgradtsm(1,js)
504 gradtyj=wgradtsm(2,js)
505 gradtzj=wgradtsm(3,js)
506 wgrdx=wgrad(1)
507 wgrdy=wgrad(2)
508 wgrdz=wgrad(3)
509 wgrad(1)=wgrdx*alphai+wght*alphaxi
510 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
511 wgrad(2)=wgrdy*alphai+wght*alphayi
512 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
513 wgrad(3)=wgrdz*alphai+wght*alphazi
514 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
515! Old order1 correction
516! BETAX=ONE + BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
517! WGRAD(1)=WGRDX*ALPHAI*BETAX
518! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
519! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
520! WGRAD(2)=WGRDY*ALPHAI*BETAX
521! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
522! . (betaxyi*(xi-xj)+betayyi*(yi-yj)+betazyi*(zi-zj)+betayi))
523! WGRAD(3)=WGRDZ*ALPHAI*BETAX
524! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
525! . (betaxzi*(xi-xj)+betayzi*(yi-yj)+betazzi*(zi-zj)+betazi))
526C----------
527C noyau conjugue.
528 alphaj=wacomp(1,sm)
529C BETAXJ=WSMCOMP(1,JS)
530C BETAYJ=WSMCOMP(2,JS)
531C BETAZJ=WSMCOMP(3,JS)
532 alphaxj=wsmcomp( 4,js)
533 alphayj=wsmcomp( 5,js)
534 alphazj=wsmcomp( 6,js)
535 betaxxj=wacomp( 8,sm)
536 betayxj=wacomp( 9,sm)
537 betazxj=wacomp(10,sm)
538 betaxyj=wacomp(11,sm)
539 betayyj=wacomp(12,sm)
540 betazyj=wacomp(13,sm)
541 betaxzj=wacomp(14,sm)
542 betayzj=wacomp(15,sm)
543 betazzj=wacomp(16,sm)
544C
545 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
546 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
547 wgrd(2)=-wgrdy*alphaj+wght*alphayj
548 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
549 wgrd(3)=-wgrdz*alphaj+wght*alphazj
550 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
551! Old order1 correction
552! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
553! wgrd(1)=-wgrdx*alphaj*betax
554! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
555! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
556! WGRD(2)=-WGRDY*ALPHAJ*BETAX
557! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
558! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
559! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
560! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
561! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
562C
563 wlaplt(n)=wlaplt(n)+vj*(
564 . -lambda(sm)*(gradtxj*wgrd(1)+gradtyj*wgrd(2)+gradtzj*wgrd(3))
565 . +lambda(n)*(gradtxi*wgrad(1)+gradtyi*wgrad(2)+gradtzi*wgrad(3)))
566 ELSE
567 sm=-js/(nspcond+1)
568 nc=mod(-js,nspcond+1)
569 js=ispsymr(nc,sm)
570 xj =xspsym(1,js)
571 yj =xspsym(2,js)
572 zj =xspsym(3,js)
573 dj =xsphr(2,sm)
574 dij =half*(di+dj)
575 rhoj=xsphr(7,sm)
576 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
577 jnod=kxsp(3,sm)
578 vj=xsphr(8,sm)/max(em20,rhoj)
579 gradtxj=wgradtsm(1,js)
580 gradtyj=wgradtsm(2,js)
581 gradtzj=wgradtsm(3,js)
582C
583 wgrdx=wgrad(1)
584 wgrdy=wgrad(2)
585 wgrdz=wgrad(3)
586 wgrad(1)=wgrdx*alphai+wght*alphaxi
587 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
588 wgrad(2)=wgrdy*alphai+wght*alphayi
589 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
590 wgrad(3)=wgrdz*alphai+wght*alphazi
591 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
592! Old order1 correction
593! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
594! WGRAD(1)=WGRDX*ALPHAI*BETAX
595! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
596! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
597! WGRAD(2)=WGRDY*ALPHAI*BETAX
598! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
599! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
600! WGRAD(3)=WGRDZ*ALPHAI*BETAX
601! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
602! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
603C----------
604C noyau conjugue.
605 alphaj=wacompr(1,sm)
606C BETAXJ=WSMCOMP(1,JS)
607C BETAYJ=WSMCOMP(2,JS)
608C BETAZJ=WSMCOMP(3,JS)
609 alphaxj=wsmcomp( 4,js)
610 alphayj=wsmcomp( 5,js)
611 alphazj=wsmcomp( 6,js)
612 betaxxj=wacompr( 8,sm)
613 betayxj=wacompr( 9,sm)
614 betazxj=wacompr(10,sm)
615 betaxyj=wacompr(11,sm)
616 betayyj=wacompr(12,sm)
617 betazyj=wacompr(13,sm)
618 betaxzj=wacompr(14,sm)
619 betayzj=wacompr(15,sm)
620 betazzj=wacompr(16,sm)
621C
622 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
623 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
624 wgrd(2)=-wgrdy*alphaj+wght*alphayj
625 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
626 wgrd(3)=-wgrdz*alphaj+wght*alphazj
627 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
628! Old order1 correction
629! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
630! WGRD(1)=-WGRDX*ALPHAJ*BETAX
631! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
632! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
633! WGRD(2)=-WGRDY*ALPHAJ*BETAX
634! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
635! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
636! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
637! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
638! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
639C
640 wlaplt(n)=wlaplt(n)+vj*(
641 . -lambdr(sm)*(gradtxj*wgrd(1)+gradtyj*wgrd(2)+gradtzj*wgrd(3))
642 . +lambda(n)*(gradtxi*wgrad(1)+gradtyi*wgrad(2)+gradtzi*wgrad(3)))
643 END IF
644 END DO
645C------
646 10 CONTINUE
647C-----------------------------------------------
648 RETURN
649 END
650!||====================================================================
651!|| spgtsym ../engine/source/elements/sph/sptemp.F
652!||--- called by ------------------------------------------------------
653!|| forintp ../engine/source/elements/forintp.F
654!||--- uses -----------------------------------------------------
655!|| sphbox ../engine/share/modules/sphbox.F
656!||====================================================================
657 SUBROUTINE spgtsym(
658 1 ISPCOND, XFRAME, ISPSYM, XSPSYM,
659 2 WGRADT, WGRADTSM,WASPACT, WGR,
660 3 LFT, LLT, NFT)
661C-----------------------------------------------
662C M o d u l e s
663C-----------------------------------------------
664 USE sphbox
665C-----------------------------------------------
666C I m p l i c i t T y p e s
667C-----------------------------------------------
668#include "implicit_f.inc"
669C-----------------------------------------------
670C C o m m o n B l o c k s
671C-----------------------------------------------
672#include "sphcom.inc"
673#include "param_c.inc"
674C-----------------------------------------------
675C D u m m y A r g u m e n t s
676C-----------------------------------------------
677 INTEGER, INTENT(INOUT) :: LFT
678 INTEGER, INTENT(INOUT) :: LLT
679 INTEGER, INTENT(INOUT) :: NFT
680 INTEGER ISPCOND(NISPCOND,*), ISPSYM(NSPCOND,*), WASPACT(*)
681 my_real
682 . XFRAME(NXFRAME,*) ,XSPSYM(3,*) , WGRADT(3,*),
683 . WGRADTSM(3,*), WGR(3,*)
684C-----------------------------------------------
685C L o c a l V a r i a b l e s
686C-----------------------------------------------
687 INTEGER IC,NC,IS,SM,JS,ISLIDE,SS
688 my_real
689 . sx,sy,sz,
690 . nx,ny,nz,tx,ty,tz,nn,
691 . ox,oy,oz,ux,uy,uz,vx,vy,vz,wx,wy,wz
692C-----------------------------------------------
693C Prepare le gradient de temperature sur les particules symetriques.
694C-----------------------------------------------
695 ox = zero
696 oy = zero
697 oz = zero
698 ux = zero
699 uy = zero
700 uz = zero
701 vx = zero
702 vy = zero
703 vz = zero
704 wx = zero
705 wy = zero
706 wz = zero
707 DO nc=1,nspcond
708 ic=ispcond(2,nc)
709 is=ispcond(3,nc)
710 islide=ispcond(5,nc)
711 IF (ic==1) THEN
712 ox=xframe(10,is)
713 oy=xframe(11,is)
714 oz=xframe(12,is)
715 ux=xframe(1,is)
716 uy=xframe(2,is)
717 uz=xframe(3,is)
718 ELSEIF (ic==2) THEN
719 ox=xframe(10,is)
720 oy=xframe(11,is)
721 oz=xframe(12,is)
722 ux=xframe(4,is)
723 uy=xframe(5,is)
724 uz=xframe(6,is)
725 ELSEIF (ic==3) THEN
726 ox=xframe(10,is)
727 oy=xframe(11,is)
728 oz=xframe(12,is)
729 ux=xframe(7,is)
730 uy=xframe(8,is)
731 uz=xframe(9,is)
732 ENDIF
733 DO ss=1,nsphact
734 sm=waspact(ss)
735 js=ispsym(nc,sm)
736 IF(js>0)THEN
737 sx=wgradt(1,sm)
738 sy=wgradt(2,sm)
739 sz=wgradt(3,sm)
740C IF(ISLIDE==0)THEN
741C----------
742 nn=sx*ux+sy*uy+sz*uz
743 nx=nn*ux
744 ny=nn*uy
745 nz=nn*uz
746 tx=sx-nx
747 ty=sy-ny
748 tz=sz-nz
749 wgradtsm(1,js)=tx-nx
750 wgradtsm(2,js)=ty-ny
751 wgradtsm(3,js)=tz-nz
752C ELSE
753C ENDIF
754 ENDIF
755 ENDDO
756C
757C Particules symetriques de particules remotes
758C
759 DO ss=1,nsphr
760 js=ispsymr(nc,ss)
761 IF(js>0)THEN
762 sx=wgr(1,ss)
763 sy=wgr(2,ss)
764 sz=wgr(3,ss)
765C IF(ISLIDE==0)THEN
766C----------
767 nn=sx*ux+sy*uy+sz*uz
768 nx=nn*ux
769 ny=nn*uy
770 nz=nn*uz
771 tx=sx-nx
772 ty=sy-ny
773 tz=sz-nz
774 wgradtsm(1,js)=tx-nx
775 wgradtsm(2,js)=ty-ny
776 wgradtsm(3,js)=tz-nz
777C ELSE
778C ENDIF
779 ENDIF
780 ENDDO
781C----------------------------------
782 ENDDO
783 RETURN
784 END
785!||====================================================================
786!|| sptempel ../engine/source/elements/sph/sptemp.F
787!||--- called by ------------------------------------------------------
788!|| spstres ../engine/source/elements/sph/spstres.F
789!||====================================================================
790 SUBROUTINE sptempel(
791 1 KXSP, TEMP, TEMPEL, LFT,
792 2 LLT, NFT)
793C-----------------------------------------------
794C I m p l i c i t T y p e s
795C-----------------------------------------------
796#include "implicit_f.inc"
797C-----------------------------------------------
798C C o m m o n B l o c k s
799C-----------------------------------------------
800#include "sphcom.inc"
801C-----------------------------------------------
802C D u m m y A r g u m e n t s
803C-----------------------------------------------
804 INTEGER, INTENT(INOUT) :: LFT
805 INTEGER, INTENT(INOUT) :: LLT
806 INTEGER, INTENT(INOUT) :: NFT
807 INTEGER KXSP(NISP,*)
808 my_real TEMP(*),TEMPEL(*)
809C-----------------------------------------------
810C L o c a l V a r i a b l e s
811C-----------------------------------------------
812 INTEGER I,N,INOD
813C-----------------------------------------------
814C Temperature in element is equivalent to Nodal temperature for SPH
815C-----------------------------------------------
816 DO I=lft,llt
817 n = nft+i
818 IF(kxsp(2,n)>0)THEN
819 inod = kxsp(3,n)
820 tempel(i)=temp(inod)
821 ENDIF
822 ENDDO
823C
824 RETURN
825 END
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
spgtsym
subroutine spgtsym(ispcond, xframe, ispsym, xspsym, wgradt, wgradtsm, waspact, wgr, lft, llt, nft)
Definition sptemp.F:661
spgradt
subroutine spgradt(x, ms, spbuf, kxsp, ixsp, nod2sp, ispsym, xspsym, wa, wacomp, wtemp, wtr, wgradt, lft, llt, nft)
Definition sptemp.F:37
splaplt
subroutine splaplt(x, ms, spbuf, kxsp, ixsp, nod2sp, ispsym, xspsym, wa, wacomp, wgradt, wgr, wgradtsm, wlaplt, wsmcomp, lambda, lambdr, lft, llt, nft)
Definition sptemp.F:247
sptempel
subroutine sptempel(kxsp, temp, tempel, lft, llt, nft)
Definition sptemp.F:793
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79