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