OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spcompl.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!|| spcompl ../engine/source/elements/sph/spcompl.F
25!||--- called by ------------------------------------------------------
26!|| sphprep ../engine/source/elements/sph/sphprep.F
27!||--- calls -----------------------------------------------------
28!|| weight0 ../engine/source/elements/sph/weight.f
29!|| weight1 ../engine/source/elements/sph/weight.F
30!||--- uses -----------------------------------------------------
31!|| sphbox ../engine/share/modules/sphbox.F
32!||====================================================================
33 SUBROUTINE spcompl(
34 1 X ,V ,MS ,SPBUF ,ITAB ,
35 2 KXSP ,IXSP ,NOD2SP ,ISPSYM ,XSPSYM ,
36 3 VSPSYM ,IPARG ,WACOMP ,ISPCOND ,
37 4 XFRAME ,WSMCOMP ,GEO ,IPART ,IPARTSP ,
38 5 WASPACT ,ITASK ,SPH_IORD1,NUMGEO ,NCYCLE ,
39 6 MCHECK)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE sphbox
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C C o m m o n B l o c k s
50C-----------------------------------------------
51#include "sphcom.inc"
52#include "param_c.inc"
53#include "scr17_c.inc"
54#include "task_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),
59 . ISPSYM(NSPCOND,*),IPARG(NPARG,*),ISPCOND(NISPCOND,*),
60 . IPART(LIPART1,*),IPARTSP(*),WASPACT(*),
61 . ITASK
62 INTEGER, INTENT(IN) :: NUMGEO ,NCYCLE,MCHECK
63 INTEGER, INTENT(INOUT) :: SPH_IORD1
64C REAL
65 my_real
66 . x(3,*) ,v(3,*) ,ms(*) ,spbuf(nspbuf,*) ,
67 . xspsym(3,*) ,vspsym(3,*) ,wacomp(16,*),
68 . xframe(nxframe,*) ,wsmcomp(6,*),
69 . geo(npropg,*)
70C-----------------------------------------------
71C L o c a l V a r i a b l e s
72C-----------------------------------------------
73 INTEGER N,INOD,JNOD,J,NVOIS,M,NN,
74 . NVOISS,SM,JS,NC,IC,IS,
75 . I,IPRT,IGEO,IORDER,NMUN,NZERO,NUN,
76 . MWAMUN(NSPHACT),MWAUN(NSPHACT),MWAZERO(NSPHACT),
77 . IGTYP
78 my_real
79 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,dij,
80 . vj,vjx,vjy,vjz,
81 . wght,wgrad(3),
82 . wcompi,wcompxi,wcompyi,wcompzi,
83 . wgradxi,wgradyi,wgradzi,
84 . wgradxxi,wgradxyi,wgradxzi,
85 . wgradyxi,wgradyyi,wgradyzi,
86 . wgradzxi,wgradzyi,wgradzzi,det,
87 . li11,li12,li13,li21,li22,li23,li31,li32,li33,
88 . tcofa11,tcofa12,tcofa13,tcofa21,tcofa22,tcofa23,
89 . tcofa31,tcofa32,tcofa33,l(3,3),il(3,3),
90 . lxi11,lxi12,lxi13,lxi22,lxi23,lxi33,
91 . lyi11,lyi12,lyi13,lyi22,lyi23,lyi33,
92 . lzi11,lzi12,lzi13,lzi22,lzi23,lzi33,
93 . mxi11,mxi12,mxi13,mxi21,mxi22,mxi23,mxi31,mxi32,mxi33,
94 . myi11,myi12,myi13,myi21,myi22,myi23,myi31,myi32,myi33,
95 . mzi11,mzi12,mzi13,mzi21,mzi22,mzi23,mzi31,mzi32,mzi33,
96 . alphai,alphaxi,alphayi,alphazi,alphai2,
97 . betaxi,betayi,betazi,
98 . betaxxi,betayxi,betazxi,
99 . betaxyi,betayyi,betazyi,
100 . betaxzi,betayzi,betazzi,
101 . betax,wgrdx,wgrdy,wgrdz,
102 . ox,oy,oz,ux,uy,uz,vx,vy,vz,wx,wy,wz,
103 . s11,s12,s13,s22,s23,s33,
104 . alphaxj,alphayj,alphazj,
105 . betaxj,betayj,betazj
106C-----------------------------------------------
107 my_real
108 . get_u_geo
109 EXTERNAL get_u_geo
110C-----------------------------------------------
111 nmun =0
112 nzero=0
113 nun =0
114 vx = zero
115 vy = zero
116 vz = zero
117 DO i=itask+1,nsphact,nthread
118 n=waspact(i)
119 iprt =ipartsp(n)
120 igeo =ipart(2,iprt)
121 iorder=nint(get_u_geo(5,igeo))
122 IF(iorder==-1)THEN
123 nmun=nmun+1
124 mwamun(nmun)=n
125 ELSEIF(iorder== 0)THEN
126 nzero=nzero+1
127 mwazero(nzero)=n
128 ELSEIF(iorder== 1)THEN
129 nun=nun+1
130 mwaun(nun)=n
131 ENDIF
132 ENDDO
133C
134C Flag for size of SPMD communication - needed only one time
135 IF ((ncycle == 0).OR.(mcheck > 0)) THEN
136 sph_iord1 = 0
137 DO igeo=1,numgeo
138 igtyp = nint(geo(12,igeo))
139 IF (igtyp==34) THEN
140 iorder=nint(get_u_geo(5,igeo))
141 IF (iorder==1) sph_iord1 = 1
142 ENDIF
143 ENDDO
144 ENDIF
145C
146C-----------------------------------------------
147C No kernel correction.
148C-----------------------------------------------
149 DO i=1,nmun
150 n=mwamun(i)
151 wacomp(1,n)=one
152 wacomp(2,n)=zero
153 wacomp(3,n)=zero
154 wacomp(4,n)=zero
155 wacomp(5,n)=zero
156 wacomp(6,n)=zero
157 wacomp(7,n)=zero
158 wacomp( 8,n)=zero
159 wacomp( 9,n)=zero
160 wacomp(10,n)=zero
161 wacomp(11,n)=zero
162 wacomp(12,n)=zero
163 wacomp(13,n)=zero
164 wacomp(14,n)=zero
165 wacomp(15,n)=zero
166 wacomp(16,n)=zero
167 ENDDO
168C-----------------------------------------------
169C-----------------------------------------------
170C Kernel correction up to order 0.
171C-----------------------------------------------
172C-----------------------------------------------
173 DO i=1,nzero
174 n=mwazero(i)
175 inod =kxsp(3,n)
176 xi=x(1,inod)
177 yi=x(2,inod)
178 zi=x(3,inod)
179 di =spbuf(1,n)
180 rhoi=spbuf(2,n)
181C------
182 CALL weight0(xi,yi,zi,xi,yi,zi,di,wght)
183 wcompi =spbuf(12,n)*wght/max(em20,rhoi)
184 wgradxi=zero
185 wgradyi=zero
186 wgradzi=zero
187C------
188 nvois=kxsp(4,n)
189 DO j=1,nvois
190 jnod=ixsp(j,n)
191 IF(jnod>0)THEN
192 m=nod2sp(jnod)
193 xj=x(1,jnod)
194 yj=x(2,jnod)
195 zj=x(3,jnod)
196 dj =spbuf(1,m)
197 rhoj=spbuf(2,m)
198 vj=spbuf(12,m)/max(em20,rhoj)
199 ELSE
200 nn = -jnod
201 xj=xsphr(3,nn)
202 yj=xsphr(4,nn)
203 zj=xsphr(5,nn)
204 dj =xsphr(2,nn)
205 rhoj=xsphr(7,nn)
206 vj=xsphr(8,nn)/max(em20,rhoj)
207 END IF
208 dij=0.5*(di+dj)
209 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
210 wcompi=wcompi+vj*wght
211 wgradxi=wgradxi+vj*wgrad(1)
212 wgradyi=wgradyi+vj*wgrad(2)
213 wgradzi=wgradzi+vj*wgrad(3)
214 ENDDO
215C------
216C partie symetrique.
217 nvoiss=kxsp(6,n)
218 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
219 js=ixsp(j,n)
220 IF(js>0)THEN
221 sm=js/(nspcond+1)
222 nc=mod(js,nspcond+1)
223 js=ispsym(nc,sm)
224 xj =xspsym(1,js)
225 yj =xspsym(2,js)
226 zj =xspsym(3,js)
227 dj =spbuf(1,sm)
228 rhoj=spbuf(2,sm)
229 dij =half*(di+dj)
230 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
231 jnod=kxsp(3,sm)
232 vj=spbuf(12,sm)/max(em20,rhoj)
233 ELSE
234 sm=-js/(nspcond+1)
235 nc=mod(-js,nspcond+1)
236 js=ispsymr(nc,sm)
237 xj =xspsym(1,js)
238 yj =xspsym(2,js)
239 zj =xspsym(3,js)
240 dj =xsphr(2,sm)
241 rhoj=xsphr(7,sm)
242 dij =half*(di+dj)
243 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
244 vj=xsphr(8,sm)/max(em20,rhoj)
245 END IF
246 wcompi=wcompi+vj*wght
247 wgradxi=wgradxi+vj*wgrad(1)
248 wgradyi=wgradyi+vj*wgrad(2)
249 wgradzi=wgradzi+vj*wgrad(3)
250 ENDDO
251C------
252C AlphaI
253 alphai=one/max(em20,wcompi)
254 wacomp(1,n)=alphai
255C------
256 alphai2=alphai*alphai
257 alphaxi=-wgradxi*alphai2
258 alphayi=-wgradyi*alphai2
259 alphazi=-wgradzi*alphai2
260 wacomp(5,n)=alphaxi
261 wacomp(6,n)=alphayi
262 wacomp(7,n)=alphazi
263C------
264 wacomp(2,n)=zero
265 wacomp(3,n)=zero
266 wacomp(4,n)=zero
267 wacomp( 8,n)=zero
268 wacomp( 9,n)=zero
269 wacomp(10,n)=zero
270 wacomp(11,n)=zero
271 wacomp(12,n)=zero
272 wacomp(13,n)=zero
273 wacomp(14,n)=zero
274 wacomp(15,n)=zero
275 wacomp(16,n)=zero
276C------
277 ENDDO
278C-----------------------------------------------
279C-----------------------------------------------
280C Kernel correction up to order 1.
281C-----------------------------------------------
282C-----------------------------------------------
283 DO i=1,nun
284 n=mwaun(i)
285 inod =kxsp(3,n)
286 xi=x(1,inod)
287 yi=x(2,inod)
288 zi=x(3,inod)
289 di =spbuf(1,n)
290 rhoi=spbuf(2,n)
291C------
292 CALL weight0(xi,yi,zi,xi,yi,zi,di,wght)
293 wcompi =spbuf(12,n)*wght/max(em20,rhoi)
294 wgradxxi=zero
295 wgradxyi=zero
296 wgradxzi=zero
297 wgradyxi=zero
298 wgradyyi=zero
299 wgradyzi=zero
300 wgradzxi=zero
301 wgradzyi=zero
302 wgradzzi=zero
303C------
304 nvois=kxsp(4,n)
305 DO j=1,nvois
306 jnod=ixsp(j,n)
307 IF(jnod>0)THEN
308 m=nod2sp(jnod)
309 xj=x(1,jnod)
310 yj=x(2,jnod)
311 zj=x(3,jnod)
312 dj =spbuf(1,m)
313 rhoj=spbuf(2,m)
314 vj=spbuf(12,m)/max(em20,rhoj)
315 ELSE
316 nn = -jnod
317 xj=xsphr(3,nn)
318 yj=xsphr(4,nn)
319 zj=xsphr(5,nn)
320 dj =xsphr(2,nn)
321 rhoj=xsphr(7,nn)
322 vj=xsphr(8,nn)/max(em20,rhoj)
323 END IF
324 dij=0.5*(di+dj)
325 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
326 wcompi=wcompi+vj*wght
327 wcompxi=wcompxi+vj*wght*(xi-xj)
328 wcompyi=wcompyi+vj*wght*(yi-yj)
329 wcompzi=wcompzi+vj*wght*(zi-zj)
330 wgradxxi=wgradxxi+vj*wgrad(1)*(xi-xj)
331 wgradxyi=wgradxyi+vj*wgrad(1)*(yi-yj)
332 wgradxzi=wgradxzi+vj*wgrad(1)*(zi-zj)
333 wgradyxi=wgradyxi+vj*wgrad(2)*(xi-xj)
334 wgradyyi=wgradyyi+vj*wgrad(2)*(yi-yj)
335 wgradyzi=wgradyzi+vj*wgrad(2)*(zi-zj)
336 wgradzxi=wgradzxi+vj*wgrad(3)*(xi-xj)
337 wgradzyi=wgradzyi+vj*wgrad(3)*(yi-yj)
338 wgradzzi=wgradzzi+vj*wgrad(3)*(zi-zj)
339 ENDDO
340C------
341C partie symetrique.
342 nvoiss=kxsp(6,n)
343 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
344 js=ixsp(j,n)
345 IF(js>0)THEN
346 sm=js/(nspcond+1)
347 nc=mod(js,nspcond+1)
348 js=ispsym(nc,sm)
349 xj =xspsym(1,js)
350 yj =xspsym(2,js)
351 zj =xspsym(3,js)
352 dj =spbuf(1,sm)
353 rhoj=spbuf(2,sm)
354 dij =half*(di+dj)
355 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
356 jnod=kxsp(3,sm)
357 vj=spbuf(12,sm)/max(em20,rhoj)
358 ELSE
359 sm=-js/(nspcond+1)
360 nc=mod(-js,nspcond+1)
361 js=ispsymr(nc,sm)
362 xj =xspsym(1,js)
363 yj =xspsym(2,js)
364 zj =xspsym(3,js)
365 dj =xsphr(2,sm)
366 rhoj=xsphr(7,sm)
367 dij =half*(di+dj)
368 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
369 vj=xsphr(8,sm)/max(em20,rhoj)
370 END IF
371 wcompi=wcompi+vj*wght
372 wcompxi=wcompxi+vj*wght*(xi-xj)
373 wcompyi=wcompyi+vj*wght*(yi-yj)
374 wcompzi=wcompzi+vj*wght*(zi-zj)
375C
376 wgradxxi=wgradxxi+vj*wgrad(1)*(xi-xj)
377 wgradxyi=wgradxyi+vj*wgrad(1)*(yi-yj)
378 wgradxzi=wgradxzi+vj*wgrad(1)*(zi-zj)
379 wgradyxi=wgradyxi+vj*wgrad(2)*(xi-xj)
380 wgradyyi=wgradyyi+vj*wgrad(2)*(yi-yj)
381 wgradyzi=wgradyzi+vj*wgrad(2)*(zi-zj)
382 wgradzxi=wgradzxi+vj*wgrad(3)*(xi-xj)
383 wgradzyi=wgradzyi+vj*wgrad(3)*(yi-yj)
384 wgradzzi=wgradzzi+vj*wgrad(3)*(zi-zj)
385 ENDDO
386C------
387C AlphaI
388C ALPHAI=ONE/MAX(EM20,WCOMPI)
389 wacomp(2,n)=one/sign(max(em20,abs(wcompxi)),wcompxi)
390 wacomp(3,n)=one/sign(max(em20,abs(wcompyi)),wcompyi)
391 wacomp(4,n)=one/sign(max(em20,abs(wcompzi)),wcompzi)
392C------
393 li11=wgradxxi
394 li12=wgradxyi
395 li13=wgradxzi
396 li21=wgradyxi
397 li22=wgradyyi
398 li23=wgradyzi
399 li31=wgradzxi
400 li32=wgradzyi
401 li33=wgradzzi
402C Transpose of cofactors.
403 tcofa11= li22*li33-li32*li23
404 tcofa21= -li21*li33+li31*li23
405 tcofa31= li21*li32-li31*li22
406 tcofa12= -li12*li33+li32*li13
407 tcofa22= li11*li33-li31*li13
408 tcofa32= -li11*li32+li31*li12
409 tcofa13= li12*li23-li22*li13
410 tcofa23= -li11*li23+li21*li13
411 tcofa33= li11*li22-li21*li12
412C Determinant
413 det=li11*tcofa11+li12*tcofa21+li13*tcofa31
414 det=1./det
415C Inverse of correctionmatric
416 wacomp( 8,n)= -tcofa11*det
417 wacomp( 9,n)= -tcofa21*det
418 wacomp(10,n)= -tcofa31*det
419 wacomp(11,n)= -tcofa12*det
420 wacomp(12,n)= -tcofa22*det
421 wacomp(13,n)= -tcofa32*det
422 wacomp(14,n)= -tcofa13*det
423 wacomp(15,n)= -tcofa23*det
424 wacomp(16,n)= -tcofa33*det
425C------
426 wacomp( 1,n)=zero
427 wacomp( 5,n)=zero
428 wacomp( 6,n)=zero
429 wacomp( 7,n)=zero
430C------
431 ENDDO
432C-----------------------------------------------
433C-----------------------------------------------
434C Old Kernel correction up to order 1
435C Not compatible with SPMD
436C No more used
437C-----------------------------------------------
438C-----------------------------------------------
439! DO I=1,NUN_OLD
440! N=MWAUN_OLD(I)
441! INOD =KXSP(3,N)
442! XI=X(1,INOD)
443! YI=X(2,INOD)
444! ZI=X(3,INOD)
445! DI =SPBUF(1,N)
446! RHOI=SPBUF(2,N)
447C------
448! CALL WEIGHT0(XI,YI,ZI,XI,YI,ZI,DI,WGHT)
449! WCOMPI =SPBUF(12,N)*WGHT/MAX(EM20,RHOI)
450! WCOMPXI=ZERO
451! WCOMPYI=ZERO
452! WCOMPZI=ZERO
453! WGRADXI=ZERO
454! WGRADYI=ZERO
455! WGRADZI=ZERO
456! WGRADXXI=ZERO
457! WGRADXYI=ZERO
458! WGRADXZI=ZERO
459! WGRADYYI=ZERO
460! WGRADYZI=ZERO
461! WGRADZZI=ZERO
462C WGRADYXI=0.
463C WGRADZXI=0.
464C WGRADZYI=0.
465C------
466! NVOIS=KXSP(4,N)
467! DO J=1,NVOIS
468! JNOD=IXSP(J,N)
469! IF(JNOD>0)THEN
470! m=nod2sp(jnod)
471! XJ=X(1,JNOD)
472! yj=x(2,jnod)
473! ZJ=X(3,JNOD)
474! DJ =SPBUF(1,M)
475! RHOJ=SPBUF(2,M)
476! VJ=SPBUF(12,M)/MAX(EM20,RHOJ)
477! ELSE
478! NN = -JNOD
479! XJ=XSPHR(3,NN)
480! YJ=XSPHR(4,NN)
481! ZJ=XSPHR(5,NN)
482! DJ =XSPHR(2,NN)
483! RHOJ=XSPHR(7,NN)
484! VJ=XSPHR(8,NN)/MAX(EM20,RHOJ)
485! END IF
486! DIJ=HALF*(DI+DJ)
487! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
488! WCOMPI=WCOMPI+VJ*WGHT
489! WCOMPXI=WCOMPXI+VJ*WGHT*(XI-XJ)
490! WCOMPYI=WCOMPYI+VJ*WGHT*(YI-YJ)
491! WCOMPZI=WCOMPZI+VJ*WGHT*(ZI-ZJ)
492! WGRADXI=WGRADXI+VJ*WGRAD(1)
493! WGRADYI=WGRADYI+VJ*WGRAD(2)
494! WGRADZI=WGRADZI+VJ*WGRAD(3)
495! WGRADXXI=WGRADXXI+VJ*WGRAD(1)*(XI-XJ)
496! WGRADXYI=WGRADXYI+VJ*WGRAD(1)*(YI-YJ)
497! WGRADXZI=WGRADXZI+VJ*WGRAD(1)*(ZI-ZJ)
498! WGRADYYI=WGRADYYI+VJ*WGRAD(2)*(YI-YJ)
499! WGRADYZI=WGRADYZI+VJ*WGRAD(2)*(ZI-ZJ)
500! WGRADZZI=WGRADZZI+VJ*WGRAD(3)*(ZI-ZJ)
501C WGRADYXI=WGRADYXI+VJ*WGRAD(2)*(XI-XJ)
502C WGRADZXI=WGRADZXI+VJ*WGRAD(3)*(XI-XJ)
503C WGRADZYI=WGRADZYI+VJ*WGRAD(3)*(YI-YJ)
504! ENDDO
505C------
506C partie symetrique.
507! NVOISS=KXSP(6,N)
508! DO J=KXSP(5,N)+1,KXSP(5,N)+NVOISS
509! JS=IXSP(J,N)
510! IF(JS>0)THEN
511! SM=JS/(NSPCOND+1)
512! NC=MOD(JS,NSPCOND+1)
513! JS=ISPSYM(NC,SM)
514! XJ =XSPSYM(1,JS)
515! YJ =XSPSYM(2,JS)
516! ZJ =XSPSYM(3,JS)
517! DJ =SPBUF(1,SM)
518! RHOJ=SPBUF(2,SM)
519! DIJ =HALF*(DI+DJ)
520! CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
521! JNOD=KXSP(3,SM)
522! VJ=SPBUF(12,SM)/MAX(EM20,RHOJ)
523! ELSE
524! SM=-JS/(NSPCOND+1)
525! NC=MOD(-JS,NSPCOND+1)
526! JS=ISPSYMR(NC,SM)
527! XJ =XSPSYM(1,JS)
528! YJ =XSPSYM(2,JS)
529! ZJ =XSPSYM(3,JS)
530! DJ =XSPHR(2,SM)
531! RHOJ=XSPHR(7,SM)
532! DIJ =HALF*(DI+DJ)
533! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
534! VJ=XSPHR(8,SM)/MAX(EM20,RHOJ)
535! END IF
536! WCOMPI=WCOMPI+VJ*WGHT
537! WCOMPXI=WCOMPXI+VJ*WGHT*(XI-XJ)
538! WCOMPYI=WCOMPYI+VJ*WGHT*(YI-YJ)
539! WCOMPZI=WCOMPZI+VJ*WGHT*(ZI-ZJ)
540! WGRADXI=WGRADXI+VJ*WGRAD(1)
541! WGRADYI=WGRADYI+VJ*WGRAD(2)
542! WGRADZI=WGRADZI+VJ*WGRAD(3)
543! WGRADXXI=WGRADXXI+VJ*WGRAD(1)*(XI-XJ)
544! WGRADXYI=WGRADXYI+VJ*WGRAD(1)*(YI-YJ)
545! WGRADXZI=WGRADXZI+VJ*WGRAD(1)*(ZI-ZJ)
546! WGRADYYI=WGRADYYI+VJ*WGRAD(2)*(YI-YJ)
547! WGRADYZI=WGRADYZI+VJ*WGRAD(2)*(ZI-ZJ)
548! WGRADZZI=WGRADZZI+VJ*WGRAD(3)*(ZI-ZJ)
549C WGRADYXI=WGRADYXI+VJ*WGRAD(2)*(XI-XJ)
550C WGRADZXI=WGRADZXI+VJ*WGRAD(3)*(XI-XJ)
551C WGRADZYI=WGRADZYI+VJ*WGRAD(3)*(YI-YJ)
552! ENDDO
553C------
554! LI11=ZERO
555! LI12=ZERO
556! LI13=ZERO
557! LI22=ZERO
558! LI23=ZERO
559! LI33=ZERO
560C LI21=0.
561C LI31=0.
562C LI32=0.
563C------
564! LXI11=ZERO
565! LXI12=ZERO
566! LXI13=ZERO
567! LXI22=ZERO
568! LXI23=ZERO
569! LXI33=ZERO
570C LXI21=0.
571C LXI31=0.
572C LXI32=0.
573C------
574! LYI11=ZERO
575! LYI12=ZERO
576! LYI13=ZERO
577! LYI22=ZERO
578! LYI23=ZERO
579! LYI33=ZERO
580C LYI21=0.
581C LYI31=0.
582C LYI32=0.
583C------
584! LZI11=ZERO
585! LZI12=ZERO
586! LZI13=ZERO
587! LZI22=ZERO
588! LZI23=ZERO
589! LZI33=ZERO
590C LZI21=0.
591C LZI31=0.
592C LZI32=0.
593C------
594! NVOIS=KXSP(4,N)
595! DO J=1,NVOIS
596! JNOD=IXSP(J,N)
597! IF(JNOD>0)THEN
598! M=NOD2SP(JNOD)
599! XJ=X(1,JNOD)
600! YJ=X(2,JNOD)
601! ZJ=X(3,JNOD)
602! DJ =SPBUF(1,M)
603! RHOJ=SPBUF(2,M)
604! DIJ=HALF*(DI+DJ)
605! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
606C------
607! VJ=SPBUF(12,M)/MAX(EM20,RHOJ)*WGHT
608! VJX=VJ*(XI-XJ)
609! VJY=VJ*(YI-YJ)
610! VJZ=VJ*(ZI-ZJ)
611C------
612! LI11=LI11+VJX*(XI-XJ)
613! LI12=LI12+VJX*(YI-YJ)
614! li13=li13+vjx*(zi-zj)
615! LI22=LI22+VJY*(YI-YJ)
616! LI23=LI23+VJY*(ZI-ZJ)
617! LI33=LI33+VJZ*(ZI-ZJ)
618C LI21=LI21+VJY*(XI-XJ)
619C LI31=LI31+VJZ*(XI-XJ)
620C LI32=LI32+VJZ*(YI-YJ)
621C------
622! LXI11=LXI11+VJX
623! LXI12=LXI12+VJY
624! LXI13=LXI13+VJZ
625! LXI11=LXI11+VJX
626C LXI21=LXI21+VJY
627C LXI31=LXI31+VJZ
628C------
629C LYI21=LYI21+VJX
630! LYI22=LYI22+VJY
631! LYI23=LYI23+VJZ
632! LYI12=LYI12+VJX
633! LYI22=LYI22+VJY
634C LYI32=LYI32+VJZ
635C------
636C LZI31=LZI31+VJX
637C LZI32=LZI32+VJY
638! LZI33=LZI33+VJZ
639! LZI13=LZI13+VJX
640! LZI23=LZI23+VJY
641! LZI33=LZI33+VJZ
642C------
643! VJ=SPBUF(12,M)/MAX(EM20,RHOJ)*WGRAD(1)
644! vjx=vj*(xi-xj)
645! VJY=VJ*(YI-YJ)
646! VJZ=VJ*(ZI-ZJ)
647! LXI11=LXI11+VJX*(XI-XJ)
648! LXI12=LXI12+VJX*(YI-YJ)
649! LXI13=LXI13+VJX*(ZI-ZJ)
650C LXI21=LXI21+VJY*(XI-XJ)
651! LXI22=LXI22+VJY*(YI-YJ)
652! LXI23=LXI23+VJY*(ZI-ZJ)
653C LXI31=LXI31+VJZ*(XI-XJ)
654C LXI32=LXI32+VJZ*(YI-YJ)
655! LXI33=LXI33+VJZ*(ZI-ZJ)
656C------
657! VJ=SPBUF(12,M)/MAX(EM20,RHOJ)*WGRAD(2)
658! VJX=VJ*(XI-XJ)
659! VJY=VJ*(YI-YJ)
660! VJZ=VJ*(ZI-ZJ)
661! LYI11=LYI11+VJX*(XI-XJ)
662! LYI12=LYI12+VJX*(YI-YJ)
663! LYI13=LYI13+VJX*(ZI-ZJ)
664C LYI21=LYI21+VJY*(XI-XJ)
665! lyi22=lyi22+vjy*(yi-yj)
666! LYI23=LYI23+VJY*(ZI-ZJ)
667C LYI31=LYI31+VJZ*(XI-XJ)
668C LYI32=LYI32+VJZ*(YI-YJ)
669! LYI33=LYI33+VJZ*(ZI-ZJ)
670C------
671! VJ=SPBUF(12,M)/MAX(EM20,RHOJ)*WGRAD(3)
672! VJX=VJ*(XI-XJ)
673! VJY=VJ*(YI-YJ)
674! VJZ=VJ*(ZI-ZJ)
675! LZI11=LZI11+VJX*(XI-XJ)
676! LZI12=LZI12+VJX*(YI-YJ)
677! LZI13=LZI13+VJX*(ZI-ZJ)
678C LZI21=LZI21+VJY*(XI-XJ)
679! LZI22=LZI22+VJY*(YI-YJ)
680! LZI23=LZI23+VJY*(ZI-ZJ)
681C LZI31=LZI31+VJZ*(XI-XJ)
682C LZI32=LZI32+VJZ*(YI-YJ)
683! LZI33=LZI33+VJZ*(ZI-ZJ)
684! ELSE
685! NN = -JNOD
686! XJ=XSPHR(3,NN)
687! YJ=XSPHR(4,NN)
688! ZJ=XSPHR(5,NN)
689! DJ =XSPHR(2,NN)
690! RHOJ=XSPHR(7,NN)
691! DIJ=HALF*(DI+DJ)
692! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
693! VJ=XSPHR(8,NN)/MAX(EM20,RHOJ)*WGHT
694C------
695! VJX=VJ*(XI-XJ)
696! VJY=VJ*(YI-YJ)
697! VJZ=VJ*(ZI-ZJ)
698C------
699! LI11=LI11+VJX*(XI-XJ)
700! LI12=LI12+VJX*(YI-YJ)
701! LI13=LI13+VJX*(ZI-ZJ)
702! LI22=LI22+VJY*(YI-YJ)
703! LI23=LI23+VJY*(ZI-ZJ)
704! LI33=LI33+VJZ*(ZI-ZJ)
705C LI21=LI21+VJY*(XI-XJ)
706C LI31=LI31+VJZ*(XI-XJ)
707C LI32=LI32+VJZ*(YI-YJ)
708C------
709! LXI11=LXI11+VJX
710! LXI12=LXI12+VJY
711! LXI13=LXI13+VJZ
712! LXI11=LXI11+VJX
713C LXI21=LXI21+VJY
714C LXI31=LXI31+VJZ
715C------
716C LYI21=LYI21+VJX
717! LYI22=LYI22+VJY
718! LYI23=LYI23+VJZ
719! LYI12=LYI12+VJX
720! LYI22=LYI22+VJY
721C LYI32=LYI32+VJZ
722C------
723C LZI31=LZI31+VJX
724C LZI32=LZI32+VJY
725! LZI33=LZI33+VJZ
726! LZI13=LZI13+VJX
727! LZI23=LZI23+VJY
728! LZI33=LZI33+VJZ
729C------
730
731! vj=xsphr(8,nn)/max(em20,rhoj)*wgrad(1)
732! VJX=VJ*(XI-XJ)
733! VJY=VJ*(YI-YJ)
734! VJZ=VJ*(ZI-ZJ)
735! LXI11=LXI11+VJX*(XI-XJ)
736! LXI12=LXI12+VJX*(YI-YJ)
737! LXI13=LXI13+VJX*(ZI-ZJ)
738C LXI21=LXI21+VJY*(XI-XJ)
739! LXI22=LXI22+VJY*(YI-YJ)
740! LXI23=LXI23+VJY*(ZI-ZJ)
741C LXI31=LXI31+VJZ*(XI-XJ)
742C LXI32=LXI32+VJZ*(YI-YJ)
743! LXI33=LXI33+VJZ*(ZI-ZJ)
744C------
745! VJ=XSPHR(8,NN)/MAX(EM20,RHOJ)*WGRAD(2)
746! VJX=VJ*(XI-XJ)
747! VJY=VJ*(YI-YJ)
748! VJZ=VJ*(ZI-ZJ)
749! LYI11=LYI11+VJX*(XI-XJ)
750! LYI12=LYI12+VJX*(YI-YJ)
751! LYI13=LYI13+VJX*(ZI-ZJ)
752C LYI21=LYI21+VJY*(XI-XJ)
753! lyi22=lyi22+vjy*(yi-yj)
754! lyi23=lyi23+vjy*(zi-zj)
755C LYI31=LYI31+VJZ*(XI-XJ)
756C LYI32=LYI32+VJZ*(YI-YJ)
757! LYI33=LYI33+VJZ*(ZI-ZJ)
758C------
759! VJ=XSPHR(8,NN)/MAX(EM20,RHOJ)*WGRAD(3)
760! VJX=VJ*(XI-XJ)
761! VJY=VJ*(YI-YJ)
762! VJZ=VJ*(ZI-ZJ)
763! LZI11=LZI11+VJX*(XI-XJ)
764! LZI12=LZI12+VJX*(YI-YJ)
765! LZI13=LZI13+VJX*(ZI-ZJ)
766C LZI21=LZI21+VJY*(XI-XJ)
767! LZI22=LZI22+VJY*(YI-YJ)
768! LZI23=LZI23+VJY*(ZI-ZJ)
769C LZI31=LZI31+VJZ*(XI-XJ)
770C LZI32=LZI32+VJZ*(YI-YJ)
771! LZI33=LZI33+VJZ*(ZI-ZJ)
772! END IF
773! ENDDO
774C------
775C partie symetrique.
776! NVOISS=KXSP(6,N)
777! DO J=KXSP(5,N)+1,KXSP(5,N)+NVOISS
778! JS=IXSP(J,N)
779! IF(JS>0)THEN
780! SM=JS/(NSPCOND+1)
781! nc=mod(js,nspcond+1)
782! JS=ISPSYM(NC,SM)
783! XJ =XSPSYM(1,JS)
784! YJ =XSPSYM(2,JS)
785! ZJ =XSPSYM(3,JS)
786! DJ =SPBUF(1,SM)
787! RHOJ=SPBUF(2,SM)
788! DIJ =HALF*(DI+DJ)
789! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
790! JNOD=KXSP(3,SM)
791C------
792! VJ=SPBUF(12,SM)/MAX(EM20,RHOJ)*WGHT
793! VJX=VJ*(XI-XJ)
794! VJY=VJ*(YI-YJ)
795! VJZ=VJ*(ZI-ZJ)
796C------
797! LI11=LI11+VJX*(XI-XJ)
798! LI12=LI12+VJX*(YI-YJ)
799! LI13=LI13+VJX*(ZI-ZJ)
800! LI22=LI22+VJY*(YI-YJ)
801! LI23=LI23+VJY*(ZI-ZJ)
802! LI33=LI33+VJZ*(ZI-ZJ)
803C LI21=LI21+VJY*(XI-XJ)
804C LI31=LI31+VJZ*(XI-XJ)
805C LI32=LI32+VJZ*(YI-YJ)
806C------
807! LXI11=LXI11+VJX
808! LXI12=LXI12+VJY
809! LXI13=LXI13+VJZ
810! LXI11=LXI11+VJX
811C LXI21=LXI21+VJY
812C LXI31=LXI31+VJZ
813C------
814C LYI21=LYI21+VJX
815! LYI22=LYI22+VJY
816! LYI23=LYI23+VJZ
817! LYI12=LYI12+VJX
818! LYI22=LYI22+VJY
819C LYI32=LYI32+VJZ
820C------
821C LZI31=LZI31+VJX
822C LZI32=LZI32+VJY
823! LZI33=LZI33+VJZ
824! LZI13=LZI13+VJX
825! LZI23=LZI23+VJY
826! LZI33=LZI33+VJZ
827C------
828! VJ=SPBUF(12,SM)/MAX(EM20,RHOJ)*WGRAD(1)
829! VJX=VJ*(XI-XJ)
830! VJY=VJ*(YI-YJ)
831! VJZ=VJ*(ZI-ZJ)
832! lxi11=lxi11+vjx*(xi-xj)
833! LXI12=LXI12+VJX*(YI-YJ)
834! LXI13=LXI13+VJX*(ZI-ZJ)
835C LXI21=LXI21+VJY*(XI-XJ)
836! LXI22=LXI22+VJY*(YI-YJ)
837! LXI23=LXI23+VJY*(ZI-ZJ)
838C LXI31=LXI31+VJZ*(XI-XJ)
839C LXI32=LXI32+VJZ*(YI-YJ)
840! LXI33=LXI33+VJZ*(ZI-ZJ)
841C------
842! VJ=SPBUF(12,SM)/MAX(EM20,RHOJ)*WGRAD(2)
843! VJX=VJ*(XI-XJ)
844! VJY=VJ*(YI-YJ)
845! VJZ=VJ*(ZI-ZJ)
846! LYI11=LYI11+VJX*(XI-XJ)
847! LYI12=LYI12+VJX*(YI-YJ)
848! LYI13=LYI13+VJX*(ZI-ZJ)
849C LYI21=LYI21+VJY*(XI-XJ)
850! LYI22=LYI22+VJY*(YI-YJ)
851! LYI23=LYI23+VJY*(ZI-ZJ)
852C LYI31=LYI31+VJZ*(XI-XJ)
853C LYI32=LYI32+VJZ*(YI-YJ)
854! LYI33=LYI33+VJZ*(ZI-ZJ)
855C------
856! VJ=SPBUF(12,SM)/MAX(EM20,RHOJ)*WGRAD(3)
857! VJX=VJ*(XI-XJ)
858! VJY=VJ*(YI-YJ)
859! VJZ=VJ*(ZI-ZJ)
860! LZI11=LZI11+VJX*(XI-XJ)
861! LZI12=LZI12+VJX*(YI-YJ)
862! LZI13=LZI13+VJX*(ZI-ZJ)
863C LZI21=LZI21+VJY*(XI-XJ)
864! LZI22=LZI22+VJY*(YI-YJ)
865! LZI23=LZI23+VJY*(ZI-ZJ)
866C LZI31=LZI31+VJZ*(XI-XJ)
867C LZI32=LZI32+VJZ*(YI-YJ)
868! LZI33=LZI33+VJZ*(ZI-ZJ)
869! ELSE
870! SM=-JS/(NSPCOND+1)
871! NC=MOD(-JS,NSPCOND+1)
872! JS=ISPSYMR(NC,SM)
873! XJ =XSPSYM(1,JS)
874! YJ =XSPSYM(2,JS)
875! ZJ =XSPSYM(3,JS)
876! DJ =XSPHR(2,SM)
877! RHOJ=XSPHR(7,SM)
878! DIJ =HALF*(DI+DJ)
879! CALL WEIGHT1(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT,WGRAD)
880C------
881! VJ=XSPHR(8,SM)/MAX(EM20,RHOJ)*WGHT
882! VJX=VJ*(XI-XJ)
883! VJY=VJ*(YI-YJ)
884! VJZ=VJ*(ZI-ZJ)
885C------
886! LI11=LI11+VJX*(XI-XJ)
887! LI12=LI12+VJX*(YI-YJ)
888! LI13=LI13+VJX*(ZI-ZJ)
889! LI22=LI22+VJY*(YI-YJ)
890! LI23=LI23+VJY*(ZI-ZJ)
891! LI33=LI33+VJZ*(ZI-ZJ)
892C LI21=LI21+VJY*(XI-XJ)
893C LI31=LI31+VJZ*(XI-XJ)
894C LI32=LI32+VJZ*(YI-YJ)
895C------
896! LXI11=LXI11+VJX
897! LXI12=LXI12+VJY
898! LXI13=LXI13+VJZ
899! LXI11=LXI11+VJX
900C LXI21=LXI21+VJY
901C LXI31=LXI31+VJZ
902C------
903C LYI21=LYI21+VJX
904! LYI22=LYI22+VJY
905! LYI23=LYI23+VJZ
906! LYI12=LYI12+VJX
907! LYI22=LYI22+VJY
908C LYI32=LYI32+VJZ
909C------
910C LZI31=LZI31+VJX
911C LZI32=LZI32+VJY
912! LZI33=LZI33+VJZ
913! LZI13=LZI13+VJX
914! LZI23=LZI23+VJY
915! lzi33=lzi33+vjz
916C------
917! VJ=XSPHR(8,SM)/MAX(EM20,RHOJ)*WGRAD(1)
918! VJX=VJ*(XI-XJ)
919! VJY=VJ*(YI-YJ)
920! VJZ=VJ*(ZI-ZJ)
921! lxi11=lxi11+vjx*(xi-xj)
922! LXI12=LXI12+VJX*(YI-YJ)
923! LXI13=LXI13+VJX*(ZI-ZJ)
924C LXI21=LXI21+VJY*(XI-XJ)
925! LXI22=LXI22+VJY*(YI-YJ)
926! LXI23=LXI23+VJY*(ZI-ZJ)
927C LXI31=LXI31+VJZ*(XI-XJ)
928C LXI32=LXI32+VJZ*(YI-YJ)
929! LXI33=LXI33+VJZ*(ZI-ZJ)
930C------
931! VJ=XSPHR(8,SM)/MAX(EM20,RHOJ)*WGRAD(2)
932! VJX=VJ*(XI-XJ)
933! VJY=VJ*(YI-YJ)
934! VJZ=VJ*(ZI-ZJ)
935! LYI11=LYI11+VJX*(XI-XJ)
936! LYI12=LYI12+VJX*(YI-YJ)
937! LYI13=LYI13+VJX*(ZI-ZJ)
938C LYI21=LYI21+VJY*(XI-XJ)
939! LYI22=LYI22+VJY*(YI-YJ)
940! LYI23=LYI23+VJY*(ZI-ZJ)
941C LYI31=LYI31+VJZ*(XI-XJ)
942C LYI32=LYI32+VJZ*(YI-YJ)
943! LYI33=LYI33+VJZ*(ZI-ZJ)
944C------
945! VJ=XSPHR(8,SM)/MAX(EM20,RHOJ)*WGRAD(3)
946! VJX=VJ*(XI-XJ)
947! VJY=VJ*(YI-YJ)
948! VJZ=VJ*(ZI-ZJ)
949! LZI11=LZI11+VJX*(XI-XJ)
950! LZI12=LZI12+VJX*(YI-YJ)
951! LZI13=LZI13+VJX*(ZI-ZJ)
952C LZI21=LZI21+VJY*(XI-XJ)
953! LZI22=LZI22+VJY*(YI-YJ)
954! LZI23=LZI23+VJY*(ZI-ZJ)
955C LZI31=LZI31+VJZ*(XI-XJ)
956C LZI32=LZI32+VJZ*(YI-YJ)
957! LZI33=LZI33+VJZ*(ZI-ZJ)
958! END IF
959! ENDDO
960C------
961C Inverse LI.
962C------
963C Cofacteurs.
964C COFAC11= LI22*LI33-LI32*LI23
965C COFAC12= -LI21*LI33+LI31*LI23
966C COFAC13= LI21*LI32-LI31*LI22
967C COFAC21= -LI12*LI33+LI32*LI13
968C COFAC22= LI11*LI33-LI31*LI13
969C COFAC23= -LI11*LI32+LI31*LI12
970C COFAC31= LI12*LI23-LI22*LI13
971C COFAC32= -LI11*LI23+LI21*LI13
972C COFAC33= LI11*LI22-LI21*LI12
973C Transposee des cofacteurs.
974C TCOFA11= LI22*LI33-LI32*LI23
975C TCOFA21= -LI21*LI33+LI31*LI23
976C TCOFA31= LI21*LI32-LI31*LI22
977C TCOFA12= -LI12*LI33+LI32*LI13
978C TCOFA22= LI11*LI33-LI31*LI13
979C TCOFA32= -LI11*LI32+LI31*LI12
980C TCOFA13= LI12*LI23-LI22*LI13
981C TCOFA23= -LI11*LI23+LI21*LI13
982C TCOFA33= LI11*LI22-LI21*LI12
983C Determinant
984! TCOFA11= LI22*LI33-LI23*LI23
985! TCOFA12= -LI12*LI33+LI23*LI13
986! TCOFA13= LI12*LI23-LI22*LI13
987! DET=LI11*TCOFA11+LI12*TCOFA12+LI13*TCOFA13
988C------
989! IF(ABS(DET)<=EM20)THEN
990C CALL MY_LOCK
991C WRITE(ISTDO,*)' ** SPH NULL DETERMINANT.'
992C WRITE(IOUT,*) ' ** SPH NULL DETERMINANT, ID CELL=',KXSP(NISP,N)
993C CALL MY_FREE
994! L(1,1)=LI11
995! L(1,2)=LI12
996! L(1,3)=LI13
997! L(2,2)=LI22
998! L(2,3)=LI23
999! L(3,3)=LI33
1000! L(2,1)=LI12
1001! L(3,1)=LI13
1002! L(3,2)=LI23
1003C L(2,1)=LI21
1004C L(3,1)=LI31
1005C L(3,2)=LI32
1006! CALL PINVERS(L,IL)
1007! LI11=IL(1,1)
1008! LI12=IL(1,2)
1009! LI13=IL(1,3)
1010! LI22=IL(2,2)
1011! LI23=IL(2,3)
1012! LI33=IL(3,3)
1013C LI21=IL(2,1)
1014C LI31=IL(3,1)
1015C LI32=IL(3,2)
1016! ELSE
1017! DET=1./DET
1018! TCOFA22= LI11*LI33-LI13*LI13
1019! TCOFA23= -LI11*LI23+LI12*LI13
1020! TCOFA33= LI11*LI22-LI12*LI12
1021C Inverse LI
1022! LI11= TCOFA11*DET
1023! LI12= TCOFA12*DET
1024! LI13= TCOFA13*DET
1025! LI22= TCOFA22*DET
1026! LI23= TCOFA23*DET
1027! LI33= TCOFA33*DET
1028C LI21= TCOFA21*DET
1029C LI31= TCOFA31*DET
1030C LI32= TCOFA32*DET
1031! ENDIF
1032C------
1033C BetaI(3)
1034! BETAXI=-(LI11*WCOMPXI+LI12*WCOMPYI+LI13*WCOMPZI)
1035! BETAYI=-(LI12*WCOMPXI+LI22*WCOMPYI+LI23*WCOMPZI)
1036! BETAZI=-(LI13*WCOMPXI+LI23*WCOMPYI+LI33*WCOMPZI)
1037C BETAYI=-(LI21*WCOMPXI+LI22*WCOMPYI+LI23*WCOMPZI)
1038C BETAZI=-(LI31*WCOMPXI+LI32*WCOMPYI+LI33*WCOMPZI)
1039! WACOMP(2,N)=BETAXI
1040! WACOMP(3,N)=BETAYI
1041! WACOMP(4,N)=BETAZI
1042C------
1043C AlphaI
1044! ALPHAI=WCOMPI+BETAXI*WCOMPXI+BETAYI*WCOMPYI+BETAZI*WCOMPZI
1045! ALPHAI=ONE/MAX(EM20,ALPHAI)
1046! WACOMP(1,N)=ALPHAI
1047C------
1048C DBetaI(3)/dx (MXI=Inverse(LI)*LXI)
1049C MXI11=LI11*LXI11+LI12*LXI21+LI13*LXI31
1050C MXI12=LI11*LXI12+LI12*LXI22+LI13*LXI32
1051C MXI13=LI11*LXI13+LI12*LXI23+LI13*LXI33
1052C MXI21=LI21*LXI11+LI22*LXI21+LI23*LXI31
1053C MXI22=LI21*LXI12+LI22*LXI22+LI23*LXI32
1054C MXI23=LI21*LXI13+LI22*LXI23+LI23*LXI33
1055C MXI31=LI31*LXI11+LI32*LXI21+LI33*LXI31
1056C MXI32=LI31*LXI12+LI32*LXI22+LI33*LXI32
1057C MXI33=LI31*LXI13+LI32*LXI23+LI33*LXI33
1058! MXI11=LI11*LXI11+LI12*LXI12+LI13*LXI13
1059! MXI12=LI11*LXI12+LI12*LXI22+LI13*LXI23
1060! MXI13=LI11*LXI13+LI12*LXI23+LI13*LXI33
1061! MXI21=LI12*LXI11+LI22*LXI12+LI23*LXI13
1062! MXI22=LI12*LXI12+LI22*LXI22+LI23*LXI23
1063! MXI23=LI12*LXI13+LI22*LXI23+LI23*LXI33
1064! MXI31=LI13*LXI11+LI23*LXI12+LI33*LXI13
1065! MXI32=LI13*LXI12+LI23*LXI22+LI33*LXI23
1066! MXI33=LI13*LXI13+LI23*LXI23+LI33*LXI33
1067! BETAXXI=-(MXI11*BETAXI+MXI12*BETAYI+MXI13*BETAZI)
1068! . -(li11*(wcompi+wgradxxi)+li12*wgradxyi+li13*wgradxzi)
1069! BETAYXI=-(MXI21*BETAXI+MXI22*BETAYI+MXI23*BETAZI)
1070! . -(LI12*(WCOMPI+WGRADXXI)+LI22*WGRADXYI+LI23*WGRADXZI)
1071C . -(LI21*(WCOMPI+WGRADXXI)+LI22*WGRADXYI+LI23*WGRADXZI)
1072! BETAZXI=-(MXI31*BETAXI+MXI32*BETAYI+MXI33*BETAZI)
1073! . -(LI13*(WCOMPI+WGRADXXI)+LI23*WGRADXYI+LI33*WGRADXZI)
1074C . -(LI31*(WCOMPI+WGRADXXI)+LI32*WGRADXYI+LI33*WGRADXZI)
1075C------
1076C DBetaI(3)/dy (MYI=Inverse(LI)*LYI)
1077C MYI11=LI11*LYI11+LI12*LYI21+LI13*LYI31
1078C MYI12=LI11*LYI12+LI12*LYI22+LI13*LYI32
1079C MYI13=LI11*LYI13+LI12*LYI23+LI13*LYI33
1080C MYI21=LI21*LYI11+LI22*LYI21+LI23*LYI31
1081C MYI22=LI21*LYI12+LI22*LYI22+LI23*LYI32
1082C MYI23=LI21*LYI13+LI22*LYI23+LI23*LYI33
1083C MYI31=LI31*LYI11+LI32*LYI21+LI33*LYI31
1084C MYI32=LI31*LYI12+LI32*LYI22+LI33*LYI32
1085C MYI33=LI31*LYI13+LI32*LYI23+LI33*LYI33
1086! MYI11=LI11*LYI11+LI12*LYI12+LI13*LYI13
1087! MYI12=LI11*LYI12+LI12*LYI22+LI13*LYI23
1088! MYI13=LI11*LYI13+LI12*LYI23+LI13*LYI33
1089! MYI21=LI12*LYI11+LI22*LYI12+LI23*LYI13
1090! MYI22=LI12*LYI12+LI22*LYI22+LI23*LYI23
1091! MYI23=LI12*LYI13+LI22*LYI23+LI23*LYI33
1092! MYI31=LI13*LYI11+LI23*LYI12+LI33*LYI13
1093! MYI32=LI13*LYI12+LI23*LYI22+LI33*LYI23
1094! MYI33=LI13*LYI13+LI23*LYI23+LI33*LYI33
1095! BETAXYI=-(MYI11*BETAXI+MYI12*BETAYI+MYI13*BETAZI)
1096! . -(LI11*WGRADXYI+LI12*(WCOMPI+WGRADYYI)+LI13*WGRADYZI)
1097C . -(LI11*WGRADYXI+LI12*(WCOMPI+WGRADYYI)+LI13*WGRADYZI)
1098! BETAYYI=-(MYI21*BETAXI+MYI22*BETAYI+MYI23*BETAZI)
1099! . -(LI12*WGRADXYI+LI22*(WCOMPI+WGRADYYI)+LI23*WGRADYZI)
1100C . -(LI21*WGRADYXI+LI22*(WCOMPI+WGRADYYI)+LI23*WGRADYZI)
1101! BETAZYI=-(MYI31*BETAXI+MYI32*BETAYI+MYI33*BETAZI)
1102! . -(LI13*WGRADXYI+LI23*(WCOMPI+WGRADYYI)+LI33*WGRADYZI)
1103C . -(LI31*WGRADYXI+LI32*(WCOMPI+WGRADYYI)+LI33*WGRADYZI)
1104C------
1105C DBetaI(3)/dz (MZI=Inverse(LI)*LZI)
1106C MZI11=LI11*LZI11+LI12*LZI21+LI13*LZI31
1107C MZI12=LI11*LZI12+LI12*LZI22+LI13*LZI32
1108C MZI13=LI11*LZI13+LI12*LZI23+LI13*LZI33
1109C MZI21=LI21*LZI11+LI22*LZI21+LI23*LZI31
1110C MZI22=LI21*LZI12+LI22*LZI22+LI23*LZI32
1111C MZI23=LI21*LZI13+LI22*LZI23+LI23*LZI33
1112C MZI31=LI31*LZI11+LI32*LZI21+LI33*LZI31
1113C MZI32=LI31*LZI12+LI32*LZI22+LI33*LZI32
1114C MZI33=LI31*LZI13+LI32*LZI23+LI33*LZI33
1115! MZI11=LI11*LZI11+LI12*LZI12+LI13*LZI13
1116! MZI12=LI11*LZI12+LI12*LZI22+LI13*LZI23
1117! MZI13=LI11*LZI13+LI12*LZI23+LI13*LZI33
1118! MZI21=LI12*LZI11+LI22*LZI12+LI23*LZI13
1119! MZI22=LI12*LZI12+LI22*LZI22+LI23*LZI23
1120! MZI23=LI12*LZI13+LI22*LZI23+LI23*LZI33
1121! MZI31=LI13*LZI11+LI23*LZI12+LI33*LZI13
1122! MZI32=LI13*LZI12+LI23*LZI22+LI33*LZI23
1123! MZI33=LI13*LZI13+LI23*LZI23+LI33*LZI33
1124! BETAXZI=-(MZI11*BETAXI+MZI12*BETAYI+MZI13*BETAZI)
1125! . -(LI11*WGRADXZI+LI12*WGRADYZI+LI13*(WCOMPI+WGRADZZI))
1126C . -(LI11*WGRADZXI+LI12*WGRADZYI+LI13*(WCOMPI+WGRADZZI))
1127! betayzi=-(mzi21*betaxi+mzi22*betayi+mzi23*betazi)
1128! . -(LI12*WGRADXZI+LI22*WGRADYZI+LI23*(WCOMPI+WGRADZZI))
1129C . -(LI21*WGRADZXI+LI22*WGRADZYI+LI23*(WCOMPI+WGRADZZI))
1130! BETAZZI=-(MZI31*BETAXI+MZI32*BETAYI+MZI33*BETAZI)
1131! . -(LI13*WGRADXZI+LI23*WGRADYZI+LI33*(WCOMPI+WGRADZZI))
1132C . -(LI31*WGRADZXI+LI32*WGRADZYI+LI33*(WCOMPI+WGRADZZI))
1133C------
1134! WACOMP( 8,N)=BETAXXI
1135! WACOMP( 9,N)=BETAYXI
1136! WACOMP(10,N)=BETAZXI
1137! WACOMP(11,N)=BETAXYI
1138! WACOMP(12,N)=BETAYYI
1139! WACOMP(13,N)=BETAZYI
1140! WACOMP(14,N)=BETAXZI
1141! WACOMP(15,N)=BETAYZI
1142! WACOMP(16,N)=BETAZZI
1143C------
1144! ALPHAI2=ALPHAI*ALPHAI
1145! ALPHAXI= WGRADXI
1146! . +BETAXI*WGRADXXI+BETAYI*WGRADXYI+BETAZI*WGRADXZI
1147! . +BETAXXI*WCOMPXI+BETAYXI*WCOMPYI+BETAZXI*WCOMPZI
1148! . +BETAXI*WCOMPI
1149! ALPHAXI=-ALPHAXI*ALPHAI2
1150! ALPHAYI= WGRADYI
1151C . +BETAXI*WGRADYXI+BETAYI*WGRADYYI+BETAZI*WGRADYZI
1152! . +BETAXI*WGRADXYI+BETAYI*WGRADYYI+BETAZI*WGRADYZI
1153! . +BETAXYI*WCOMPXI+BETAYYI*WCOMPYI+BETAZYI*WCOMPZI
1154! . +BETAYI*WCOMPI
1155! ALPHAYI=-ALPHAYI*ALPHAI2
1156! ALPHAZI= WGRADZI
1157C . +BETAXI*WGRADZXI+BETAYI*WGRADZYI+BETAZI*WGRADZZI
1158! . +BETAXI*WGRADXZI+BETAYI*WGRADYZI+BETAZI*WGRADZZI
1159! . +BETAXZI*WCOMPXI+BETAYZI*WCOMPYI+BETAZZI*WCOMPZI
1160! . +BETAZI*WCOMPI
1161! ALPHAZI=-ALPHAZI*ALPHAI2
1162C------
1163! WACOMP(5,N)=ALPHAXI
1164! WACOMP(6,N)=ALPHAYI
1165! WACOMP(7,N)=ALPHAZI
1166C------
1167! ENDDO
1168C-----------------------------------------------
1169C Prepare les corrections sur les particules symetriques.
1170C-----------------------------------------------
1171C-------------------------------------------
1172 RETURN
1173 END
1174
1175
1176!||====================================================================
1177!|| spscomp ../engine/source/elements/sph/spcompl.F
1178!||--- called by ------------------------------------------------------
1179!|| forintp ../engine/source/elements/forintp.F
1180!||--- uses -----------------------------------------------------
1181!|| sphbox ../engine/share/modules/sphbox.F
1182!||====================================================================
1183 SUBROUTINE spscomp(
1184 1 ISPSYM ,WACOMP ,ISPCOND ,XFRAME ,WSMCOMP ,
1185 2 GEO ,IPART ,IPARTSP ,WASPACT ,ITASK )
1186C-----------------------------------------------
1187C M o d u l e s
1188C-----------------------------------------------
1189 USE sphbox
1190C-----------------------------------------------
1191C I m p l i c i t T y p e s
1192C-----------------------------------------------
1193#include "implicit_f.inc"
1194C-----------------------------------------------
1195C C o m m o n B l o c k s
1196C-----------------------------------------------
1197#include "sphcom.inc"
1198#include "param_c.inc"
1199#include "scr17_c.inc"
1200#include "task_c.inc"
1201C-----------------------------------------------
1202C D u m m y A r g u m e n t s
1203C-----------------------------------------------
1204 INTEGER ISPSYM(NSPCOND,*), ISPCOND(NISPCOND,*),
1205 . IPART(LIPART1,*), IPARTSP(*), WASPACT(*), ITASK
1206C REAL
1207 my_real
1208 . WACOMP(16,*), XFRAME(NXFRAME,*), WSMCOMP(6,*),
1209 . GEO(NPROPG,*)
1210C-----------------------------------------------
1211C L o c a l V a r i a b l e s
1212C-----------------------------------------------
1213 INTEGER MWAMUN(NSPHACT),MWAZERO(NSPHACT),MWAUN(NSPHACT)
1214 INTEGER SM,JS,NC,IC,IS,
1215 . i,n,iprt,igeo,iorder,nmun,nzero,nun
1216 my_real
1217 . ox,oy,oz,ux,uy,uz,vx,vy,vz,wx,wy,wz,
1218 . s11,s12,s13,s22,s23,s33,
1219 . alphaxi,alphayi,alphazi,
1220 . betaxi,betayi,betazi,
1221 . alphaxj,alphayj,alphazj,
1222 . betaxj,betayj,betazj
1223C-----------------------------------------------
1224 my_real
1225 . get_u_geo
1226 EXTERNAL get_u_geo
1227C-----------------------------------------------
1228 nmun =0
1229 nzero=0
1230 nun =0
1231 DO i=itask+1,nsphact,nthread
1232 n=waspact(i)
1233 iprt =ipartsp(n)
1234 igeo =ipart(2,iprt)
1235 iorder=nint(get_u_geo(5,igeo))
1236 IF(iorder==-1)THEN
1237 nmun=nmun+1
1238 mwamun(nmun)=n
1239 ELSEIF(iorder== 0)THEN
1240 nzero=nzero+1
1241 mwazero(nzero)=n
1242 ELSEIF(iorder== 1)THEN
1243 nun=nun+1
1244 mwaun(nun)=n
1245 ENDIF
1246 ENDDO
1247C-----------------------------------------------
1248C No kernel correction.
1249C-----------------------------------------------
1250 DO nc=1,nspcond
1251 ic=ispcond(2,nc)
1252 is=ispcond(3,nc)
1253 DO i=1,nmun
1254 sm=mwamun(i)
1255 js=ispsym(nc,sm)
1256 IF(js>0)THEN
1257 wsmcomp(1,js)=zero
1258 wsmcomp(2,js)=zero
1259 wsmcomp(3,js)=zero
1260 wsmcomp(4,js)=zero
1261 wsmcomp(5,js)=zero
1262 wsmcomp(6,js)=zero
1263C WACOMP( 8,JS)=0.
1264C WACOMP( 9,JS)=0.
1265C WACOMP(10,JS)=0.
1266C WACOMP(11,JS)=0.
1267C WACOMP(12,JS)=0.
1268C WACOMP(13,JS)=0.
1269C WACOMP(14,JS)=0.
1270C WACOMP(15,JS)=0.
1271C WACOMP(16,JS)=0.
1272 ENDIF
1273 ENDDO
1274 ENDDO
1275C-----------------------------------------------
1276C Kernel correction up to order 0. or 1.
1277C-----------------------------------------------
1278 DO nc=1,nspcond
1279 ic=ispcond(2,nc)
1280 is=ispcond(3,nc)
1281 IF (ic==1) THEN
1282 ox=xframe(10,is)
1283 oy=xframe(11,is)
1284 oz=xframe(12,is)
1285 ux=xframe(1,is)
1286 uy=xframe(2,is)
1287 uz=xframe(3,is)
1288 vx=xframe(4,is)
1289 vy=xframe(5,is)
1290 vz=xframe(6,is)
1291 wx=xframe(7,is)
1292 wy=xframe(8,is)
1293 wz=xframe(9,is)
1294 ELSEIF (ic==2) THEN
1295 ox=xframe(10,is)
1296 oy=xframe(11,is)
1297 oz=xframe(12,is)
1298 ux=xframe(4,is)
1299 uy=xframe(5,is)
1300 uz=xframe(6,is)
1301 vx=xframe(7,is)
1302 vy=xframe(8,is)
1303 vz=xframe(9,is)
1304 wx=xframe(1,is)
1305 wy=xframe(2,is)
1306 wz=xframe(3,is)
1307 ELSEIF (ic==3) THEN
1308 ox=xframe(10,is)
1309 oy=xframe(11,is)
1310 oz=xframe(12,is)
1311 ux=xframe(7,is)
1312 uy=xframe(8,is)
1313 uz=xframe(9,is)
1314 vx=xframe(1,is)
1315 vy=xframe(2,is)
1316 vz=xframe(3,is)
1317 wx=xframe(4,is)
1318 wy=xframe(5,is)
1319 wz=xframe(6,is)
1320 ENDIF
1321 s11=-ux*ux+vx*vx+wx*wx
1322 s12=-ux*uy+vx*vy+wx*wy
1323 s13=-ux*uz+vx*vz+wx*wz
1324 s22=-uy*uy+vy*vy+wy*wy
1325 s23=-uy*uz+vy*vz+wy*wz
1326 s33=-uz*uz+vz*vz+wz*wz
1327C S21=-UY*UX+VY*VX+WY*WX
1328C S31=-UZ*UX+VZ*VX+WZ*WX
1329C S32=-UZ*UY+VZ*VY+WZ*WY
1330C-----------------------------------------------
1331C Kernel correction up to order 0.
1332C-----------------------------------------------
1333 DO i=1,nzero
1334 sm=mwazero(i)
1335 js=ispsym(nc,sm)
1336 IF(js>0)THEN
1337C WACOMP(1,JS)=WACOMP(1,SM)
1338 wsmcomp(1,js)=zero
1339 wsmcomp(2,js)=zero
1340 wsmcomp(3,js)=zero
1341 alphaxi=wacomp( 5,sm)
1342 alphayi=wacomp( 6,sm)
1343 alphazi=wacomp( 7,sm)
1344 alphaxj=s11*alphaxi+s12*alphayi+s13*alphazi
1345 alphayj=s12*alphaxi+s22*alphayi+s23*alphazi
1346 alphazj=s13*alphaxi+s23*alphayi+s33*alphazi
1347 wsmcomp(4,js)=alphaxj
1348 wsmcomp(5,js)=alphayj
1349 wsmcomp(6,js)=alphazj
1350C WACOMP( 8,JS)=0.
1351C WACOMP( 9,JS)=0.
1352C WACOMP(10,JS)=0.
1353C WACOMP(11,JS)=0.
1354C WACOMP(12,JS)=0.
1355C WACOMP(13,JS)=0.
1356C WACOMP(14,JS)=0.
1357C WACOMP(15,JS)=0.
1358C WACOMP(16,JS)=0.
1359 ENDIF
1360 ENDDO
1361C-----------------------------------------------
1362C Kernel correction up to order 1
1363C NON COMPATIBLE SPMD SUR PLUS D ONE DOMAINE
1364C-----------------------------------------------
1365 DO i=1,nun
1366 sm=mwaun(i)
1367 js=ispsym(nc,sm)
1368 IF(js>0)THEN
1369C WACOMP(1,JS)=WACOMP(1,SM)
1370 betaxi=wacomp(2,sm)
1371 betayi=wacomp(3,sm)
1372 betazi=wacomp(4,sm)
1373 betaxj=s11*betaxi+s12*betayi+s13*betazi
1374 betayj=s12*betaxi+s22*betayi+s23*betazi
1375 betazj=s13*betaxi+s23*betayi+s33*betazi
1376 wsmcomp(1,js)=betaxj
1377 wsmcomp(2,js)=betayj
1378 wsmcomp(3,js)=betazj
1379 alphaxi=wacomp( 5,sm)
1380 alphayi=wacomp( 6,sm)
1381 alphazi=wacomp( 7,sm)
1382 alphaxj=s11*alphaxi+s12*alphayi+s13*alphazi
1383 alphayj=s12*alphaxi+s22*alphayi+s23*alphazi
1384 alphazj=s13*alphaxi+s23*alphayi+s33*alphazi
1385 wsmcomp( 4,js)=alphaxj
1386 wsmcomp( 5,js)=alphayj
1387 wsmcomp( 6,js)=alphazj
1388C WACOMP( 8,JS)=WACOMP( 8,SM)
1389C WACOMP( 9,JS)=WACOMP( 9,SM)
1390C WACOMP(10,JS)=WACOMP(10,SM)
1391C WACOMP(11,JS)=WACOMP(11,SM)
1392C WACOMP(12,JS)=WACOMP(12,SM)
1393C WACOMP(13,JS)=WACOMP(13,SM)
1394C WACOMP(14,JS)=WACOMP(14,SM)
1395C WACOMP(15,JS)=WACOMP(15,SM)
1396C WACOMP(16,JS)=WACOMP(16,SM)
1397 ENDIF
1398 ENDDO
1399C
1400C Particules symetriques de particules remotes
1401C
1402 DO sm = 1+itask,nsphr,nthread
1403 js=ispsymr(nc,sm)
1404 IF(js>0)THEN
1405 wsmcomp(1,js)=zero
1406 wsmcomp(2,js)=zero
1407 wsmcomp(3,js)=zero
1408 alphaxi=wacompr( 5,sm)
1409 alphayi=wacompr( 6,sm)
1410 alphazi=wacompr( 7,sm)
1411 alphaxj=s11*alphaxi+s12*alphayi+s13*alphazi
1412 alphayj=s12*alphaxi+s22*alphayi+s23*alphazi
1413 alphazj=s13*alphaxi+s23*alphayi+s33*alphazi
1414 wsmcomp(4,js)=alphaxj
1415 wsmcomp(5,js)=alphayj
1416 wsmcomp(6,js)=alphazj
1417 END IF
1418 END DO
1419 END DO
1420
1421C-------------------------------------------
1422 RETURN
1423 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
spscomp
subroutine spscomp(ispsym, wacomp, ispcond, xframe, wsmcomp, geo, ipart, ipartsp, waspact, itask)
Definition spcompl.F:1186
spcompl
subroutine spcompl(x, v, ms, spbuf, itab, kxsp, ixsp, nod2sp, ispsym, xspsym, vspsym, iparg, wacomp, ispcond, xframe, wsmcomp, geo, ipart, ipartsp, waspact, itask, sph_iord1, numgeo, ncycle, mcheck)
Definition spcompl.F:40
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79
weight0
subroutine weight0(xi, yi, zi, xj, yj, zj, h, w)
Definition weight.F:34