OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
soltosph.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!|| soltosphf ../engine/source/elements/sph/soltosph.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../engine/source/output/message/message.F
29!|| arret ../engine/source/system/arret.F
30!|| my_barrier ../engine/source/system/machine.F
31!|| startimeg ../engine/source/system/timer.f
32!|| stoptimeg ../engine/source/system/timer.F
33!||--- uses -----------------------------------------------------
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!|| element_mod ../common_source/modules/elements/element_mod.F90
36!|| message_mod ../engine/share/message_module/message_mod.f
37!|| soltosph_mod ../engine/share/modules/soltosph_mod.F
38!||====================================================================
39 SUBROUTINE soltosphf(
40 1 A ,SPBUF ,IXS ,KXSP ,IPARTSP ,
41 2 NOD2SP ,IRST ,NGROUNC ,IGROUNC ,IPARG ,
42 3 STIFN ,SOL2SPH,SPH2SOL ,ELBUF_TAB,ITASK ,
43 4 NODFT ,NODLT ,ISKY ,FSKYI ,IGEO ,
44 5 SOL2SPH_TYP)
45C-----------------------------------------------
46 USE elbufdef_mod
47 USE soltosph_mod
48 USE message_mod
49 use element_mod , only : nixs
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54#include "comlock.inc"
55C-----------------------------------------------
56C C o m m o n B l o c k s
57C-----------------------------------------------
58#include "com01_c.inc"
59#include "com04_c.inc"
60#include "param_c.inc"
61#include "parit_c.inc"
62#include "sphcom.inc"
63#include "task_c.inc"
64#include "vect01_c.inc"
65C-----------------------------------------------
66C D u m m y A r g u m e n t s
67C-----------------------------------------------
68 INTEGER IXS(NIXS,*), KXSP(NISP,*), IPARTSP(*), NOD2SP(*),
69 . IRST(3,*), NGROUNC, IGROUNC(*), IPARG(NPARG,*),
70 . SOL2SPH(2,*), SPH2SOL(*), ITASK, NODFT, NODLT, ISKY(*),
71 . IGEO(NPROPGI,*),SOL2SPH_TYP(*)
72 my_real
73 . a(3,*), spbuf(nspbuf,*), stifn(*), fskyi(lskyi,nfskyi)
74 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I, J, IR, IS, IT, KP, NP, N, INOD, NG, NEL,
79 . n1, n2, n3, n4, n5, n6, n7, n8,
80 . ig , nelem, offset, nsphdir, iprtsph, iaw,
81 . my_iaw, ierror, nisk, niftsk, niltsk, niskyl, nsol
82C
83 my_real
84 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
85 . ksi, eta, zeta
86C-----
87 TYPE(g_bufel_) ,POINTER :: GBUF
88 TYPE(l_bufel_) ,POINTER :: LBUF
89 TYPE(BUF_MAT_) ,POINTER :: MBUF
90C-----------------------------------------------
91 my_real
92 . a_gauss(9,9),a_gauss_tetra(9,9)
93 DATA a_gauss /
94 1 0. ,0. ,0. ,
95 1 0. ,0. ,0. ,
96 1 0. ,0. ,0. ,
97 2 -.5 ,0.5 ,0. ,
98 2 0. ,0. ,0. ,
99 2 0. ,0. ,0. ,
100 3 -.666666666666666,0. ,0.666666666666666,
101 3 0. ,0. ,0. ,
102 3 0. ,0. ,0. ,
103 4 -.75 ,-.25 ,0.25 ,
104 4 0.75 ,0. ,0. ,
105 4 0. ,0. ,0. ,
106 5 -.8 ,-.4 ,0. ,
107 5 0.4 ,0.8 ,0. ,
108 5 0. ,0. ,0. ,
109 6 -.833333333333333,-.5 ,-.166666666666666,
110 6 0.166666666666666,0.5 ,0.833333333333333,
111 6 0. ,0. ,0. ,
112 7 -.857142857142857,-.571428571428571,-.285714285714285,
113 7 0. ,0.285714285714285,0.571428571428571,
114 7 0.857142857142857,0. ,0. ,
115 8 -.875 ,-.625 ,-.375 ,
116 8 -.125 ,0.125 ,0.375,
117 8 0.625 ,0.875 ,0. ,
118 9 -.888888888888888,-.666666666666666,-.444444444444444,
119 9 -.222222222222222,0. ,0.222222222222222,
120 9 0.444444444444444,0.666666666666666,0.888888888888888/
121C-----------------------------------------------
122 DATA a_gauss_tetra /
123 1 0.250000000000000,0.000000000000000,0.000000000000000,
124 1 0.000000000000000,0.000000000000000,0.000000000000000,
125 1 0.000000000000000,0.000000000000000,0.000000000000000,
126 2 0.166666666666667,0.500000000000000,0.000000000000000,
127 2 0.000000000000000,0.000000000000000,0.000000000000000,
128 2 0.000000000000000,0.000000000000000,0.000000000000000,
129 3 0.125000000000000,0.375000000000000,0.625000000000000,
130 3 0.000000000000000,0.000000000000000,0.000000000000000,
131 3 0.000000000000000,0.000000000000000,0.000000000000000,
132 4 0.100000000000000,0.300000000000000,0.500000000000000,
133 4 0.700000000000000,0.000000000000000,0.000000000000000,
134 4 0.000000000000000,0.000000000000000,0.000000000000000,
135 5 0.083333333333333,0.250000000000000,0.416666666666667,
136 5 0.583333333333333,0.750000000000000,0.000000000000000,
137 5 0.000000000000000,0.000000000000000,0.000000000000000,
138 6 0.071428571428571,0.214285714285714,0.357142857142857,
139 6 0.500000000000000,0.642857142857143,0.785714285714286,
140 6 0.000000000000000,0.000000000000000,0.000000000000000,
141 7 0.062500000000000,0.187500000000000,0.312500000000000,
142 7 0.437500000000000,0.562500000000000,0.687500000000000,
143 7 0.812500000000000,0.000000000000000,0.000000000000000,
144 8 0.055555555555556,0.166666666666667,0.277777777777778,
145 8 0.388888888888889,0.500000000000000,0.611111111111111,
146 8 0.722222222222222,0.833333333333333,0.000000000000000,
147 9 0.050000000000000,0.150000000000000,0.250000000000000,
148 9 0.350000000000000,0.450000000000000,0.550000000000000,
149 9 0.650000000000000,0.750000000000000,0.850000000000000/
150C-----------------------------------------------
151C MAPPING OF CONTACT FORCES APPLYING TO SLEEPING PARTICLES
152C-----------------------------------------------
153 IF(iparit==0)THEN
154 IF(itask==0)THEN
155 ALLOCATE(awork(4,nthread*numnod),stat=ierror)
156 IF(ierror/=0) THEN
157 CALL ancmsg(msgid=19,anmode=aninfo,
158 . c1='(SOLIDS to SPH)')
159 CALL arret(2)
160 ENDIF
161 END IF
162C
163 CALL my_barrier
164C
165 DO iaw=itask*numnod+1,(itask+1)*numnod
166 awork(1,iaw)=zero
167 awork(2,iaw)=zero
168 awork(3,iaw)=zero
169 awork(4,iaw)=zero
170 END DO
171C
172 CALL my_barrier
173C
174 my_iaw=itask*numnod
175!$OMP DO SCHEDULE(DYNAMIC,1)
176 DO ig = 1, ngrounc
177 ng = igrounc(ig)
178 IF(iparg(8,ng)==1)GOTO 350
179 IF (iddw>0) CALL startimeg(ng)
180 DO nelem = 1,iparg(2,ng),nvsiz
181 offset = nelem - 1
182 nel =iparg(2,ng)
183 nft =iparg(3,ng) + offset
184 iad =iparg(4,ng)
185 ity =iparg(5,ng)
186 iprtsph=iparg(69,ng)
187 lft=1
188 llt=min(nvsiz,nel-nelem+1)
189 IF(ity==1.AND.iprtsph/=0) THEN
190C-----------
191 gbuf => elbuf_tab(ng)%GBUF
192 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
193 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
194C-----
195C----TETRA---
196C-----
197 IF (iparg(28,ng)==4) THEN
198C-----
199 DO i=lft,llt
200 n=nft+i
201 IF(gbuf%OFF(i)/=zero) THEN
202C
203 n1=ixs(2,n)
204 n2=ixs(3,n)
205 n3=ixs(4,n)
206 n4=ixs(5,n)
207C
208 nsphdir=igeo(37,ixs(10,n))
209C
210C-----------------------------------------------
211C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
212 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
213 np =sol2sph(1,n)+kp
214 ir=irst(1,np-first_sphsol+1)
215 is=irst(2,np-first_sphsol+1)
216 it=irst(3,np-first_sphsol+1)
217 ksi = a_gauss_tetra(ir,nsphdir)
218 eta = a_gauss_tetra(is,nsphdir)
219 zeta = a_gauss_tetra(it,nsphdir)
220C
221 phi1=ksi
222 phi2=eta
223 phi3=zeta
224 phi4=1-ksi-eta-zeta
225C
226 inod=kxsp(3,np)
227 awork(1,my_iaw+n1)=awork(1,my_iaw+n1)+phi1*a(1,inod)
228 awork(2,my_iaw+n1)=awork(2,my_iaw+n1)+phi1*a(2,inod)
229 awork(3,my_iaw+n1)=awork(3,my_iaw+n1)+phi1*a(3,inod)
230 awork(4,my_iaw+n1)=awork(4,my_iaw+n1)+phi1*stifn(inod)
231 awork(1,my_iaw+n2)=awork(1,my_iaw+n2)+phi2*a(1,inod)
232 awork(2,my_iaw+n2)=awork(2,my_iaw+n2)+phi2*a(2,inod)
233 awork(3,my_iaw+n2)=awork(3,my_iaw+n2)+phi2*a(3,inod)
234 awork(4,my_iaw+n2)=awork(4,my_iaw+n2)+phi2*stifn(inod)
235 awork(1,my_iaw+n3)=awork(1,my_iaw+n3)+phi3*a(1,inod)
236 awork(2,my_iaw+n3)=awork(2,my_iaw+n3)+phi3*a(2,inod)
237 awork(3,my_iaw+n3)=awork(3,my_iaw+n3)+phi3*a(3,inod)
238 awork(4,my_iaw+n3)=awork(4,my_iaw+n3)+phi3*stifn(inod)
239 awork(1,my_iaw+n4)=awork(1,my_iaw+n4)+phi4*a(1,inod)
240 awork(2,my_iaw+n4)=awork(2,my_iaw+n4)+phi4*a(2,inod)
241 awork(3,my_iaw+n4)=awork(3,my_iaw+n4)+phi4*a(3,inod)
242 awork(4,my_iaw+n4)=awork(4,my_iaw+n4)+phi4*stifn(inod)
243C
244 a(1,inod)=zero
245 a(2,inod)=zero
246 a(3,inod)=zero
247 stifn(inod)=em20
248 ENDDO
249 END IF
250 ENDDO
251C-----
252 ELSE
253C-----
254C----HEXA---
255C-----
256 DO i=lft,llt
257 n=nft+i
258 IF(gbuf%OFF(i)/=zero) THEN
259C
260 n1=ixs(2,n)
261 n2=ixs(3,n)
262 n3=ixs(4,n)
263 n4=ixs(5,n)
264 n5=ixs(6,n)
265 n6=ixs(7,n)
266 n7=ixs(8,n)
267 n8=ixs(9,n)
268C
269 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
270C
271C-----------------------------------------------
272C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
273 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
274 np =sol2sph(1,n)+kp
275 ir=irst(1,np-first_sphsol+1)
276 is=irst(2,np-first_sphsol+1)
277 it=irst(3,np-first_sphsol+1)
278 ksi = a_gauss(ir,nsphdir)
279 eta = a_gauss(is,nsphdir)
280 zeta = a_gauss(it,nsphdir)
281C
282 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
283 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
284 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
285 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
286 phi5=one_over_8*(one-ksi)*(one+eta)*(one-zeta)
287 phi6=one_over_8*(one-ksi)*(one+eta)*(one+zeta)
288 phi7=one_over_8*(one+ksi)*(one+eta)*(one+zeta)
289 phi8=one_over_8*(one+ksi)*(one+eta)*(one-zeta)
290C
291 inod=kxsp(3,np)
292 awork(1,my_iaw+n1)=awork(1,my_iaw+n1)+phi1*a(1,inod)
293 awork(2,my_iaw+n1)=awork(2,my_iaw+n1)+phi1*a(2,inod)
294 awork(3,my_iaw+n1)=awork(3,my_iaw+n1)+phi1*a(3,inod)
295 awork(4,my_iaw+n1)=awork(4,my_iaw+n1)+phi1*stifn(inod)
296 awork(1,my_iaw+n2)=awork(1,my_iaw+n2)+phi2*a(1,inod)
297 awork(2,my_iaw+n2)=awork(2,my_iaw+n2)+phi2*a(2,inod)
298 awork(3,my_iaw+n2)=awork(3,my_iaw+n2)+phi2*a(3,inod)
299 awork(4,my_iaw+n2)=awork(4,my_iaw+n2)+phi2*stifn(inod)
300 awork(1,my_iaw+n3)=awork(1,my_iaw+n3)+phi3*a(1,inod)
301 awork(2,my_iaw+n3)=awork(2,my_iaw+n3)+phi3*a(2,inod)
302 awork(3,my_iaw+n3)=awork(3,my_iaw+n3)+phi3*a(3,inod)
303 awork(4,my_iaw+n3)=awork(4,my_iaw+n3)+phi3*stifn(inod)
304 awork(1,my_iaw+n4)=awork(1,my_iaw+n4)+phi4*a(1,inod)
305 awork(2,my_iaw+n4)=awork(2,my_iaw+n4)+phi4*a(2,inod)
306 awork(3,my_iaw+n4)=awork(3,my_iaw+n4)+phi4*a(3,inod)
307 awork(4,my_iaw+n4)=awork(4,my_iaw+n4)+phi4*stifn(inod)
308 awork(1,my_iaw+n5)=awork(1,my_iaw+n5)+phi5*a(1,inod)
309 awork(2,my_iaw+n5)=awork(2,my_iaw+n5)+phi5*a(2,inod)
310 awork(3,my_iaw+n5)=awork(3,my_iaw+n5)+phi5*a(3,inod)
311 awork(4,my_iaw+n5)=awork(4,my_iaw+n5)+phi5*stifn(inod)
312 awork(1,my_iaw+n6)=awork(1,my_iaw+n6)+phi6*a(1,inod)
313 awork(2,my_iaw+n6)=awork(2,my_iaw+n6)+phi6*a(2,inod)
314 awork(3,my_iaw+n6)=awork(3,my_iaw+n6)+phi6*a(3,inod)
315 awork(4,my_iaw+n6)=awork(4,my_iaw+n6)+phi6*stifn(inod)
316 awork(1,my_iaw+n7)=awork(1,my_iaw+n7)+phi7*a(1,inod)
317 awork(2,my_iaw+n7)=awork(2,my_iaw+n7)+phi7*a(2,inod)
318 awork(3,my_iaw+n7)=awork(3,my_iaw+n7)+phi7*a(3,inod)
319 awork(4,my_iaw+n7)=awork(4,my_iaw+n7)+phi7*stifn(inod)
320 awork(1,my_iaw+n8)=awork(1,my_iaw+n8)+phi8*a(1,inod)
321 awork(2,my_iaw+n8)=awork(2,my_iaw+n8)+phi8*a(2,inod)
322 awork(3,my_iaw+n8)=awork(3,my_iaw+n8)+phi8*a(3,inod)
323 awork(4,my_iaw+n8)=awork(4,my_iaw+n8)+phi8*stifn(inod)
324C
325 a(1,inod)=zero
326 a(2,inod)=zero
327 a(3,inod)=zero
328 stifn(inod)=em20
329 ENDDO
330 END IF
331 ENDDO
332C--------
333 ENDIF
334 ENDIF
335 IF (iddw>0) CALL stoptimeg(ng)
336 END DO
337C--------
338 350 CONTINUE
339 END DO
340!$OMP END DO
341C
342 CALL my_barrier
343C
344 DO it=0,nthread-1
345 DO n=nodft,nodlt
346 iaw=n+it*numnod
347 a(1,n)=a(1,n)+awork(1,iaw)
348 a(2,n)=a(2,n)+awork(2,iaw)
349 a(3,n)=a(3,n)+awork(3,iaw)
350 stifn(n)=stifn(n)+awork(4,iaw)
351 END DO
352 END DO
353C
354 CALL my_barrier
355C
356 IF(itask==0) DEALLOCATE(awork)
357 ELSE ! IPARIT==0
358C-----------------------------------------------
359 niftsk = 1+itask*nisky/ nthread
360 niltsk = (itask+1)*nisky/nthread
361C
362 CALL my_barrier
363C
364 DO nisk=niftsk,niltsk
365 inod=isky(nisk)
366 np =nod2sp(inod)
367 IF(np/=0)THEN
368 n=sph2sol(np)
369 IF(n/=0)THEN
370#include "lockon.inc"
371 niskyl = nisky
372 nisky = nisky + 8
373#include "lockoff.inc"
374 IF (niskyl+8 > lskyi) THEN
375 CALL ancmsg(msgid=243,anmode=aninfo_blind)
376 CALL arret(2)
377 ENDIF
378C
379 IF (sol2sph_typ(sph2sol(np))==4) THEN
380C---------------
381C------ Tetra --
382C---------------
383 nsol=sph2sol(n)
384C
385 n1=ixs(2,nsol)
386 n2=ixs(4,nsol)
387 n3=ixs(7,nsol)
388 n4=ixs(6,nsol)
389C
390 ir=irst(1,n-first_sphsol+1)
391 is=irst(2,n-first_sphsol+1)
392 it=irst(3,n-first_sphsol+1)
393 nsphdir=igeo(37,ixs(10,nsol))
394C
395 ksi = a_gauss_tetra(ir,nsphdir)
396 eta = a_gauss_tetra(is,nsphdir)
397 zeta = a_gauss_tetra(it,nsphdir)
398C
399 phi1=ksi
400 phi2=eta
401 phi3=zeta
402 phi4=1-ksi-eta-zeta
403C
404 niskyl=niskyl+1
405 isky(niskyl)=n1
406 DO j=1,nfskyi
407 fskyi(niskyl,j)=phi1*fskyi(nisk,j)
408 END DO
409 niskyl=niskyl+1
410 isky(niskyl)=n2
411 DO j=1,nfskyi
412 fskyi(niskyl,j)=phi2*fskyi(nisk,j)
413 END DO
414 niskyl=niskyl+1
415 isky(niskyl)=n3
416 DO j=1,nfskyi
417 fskyi(niskyl,j)=phi3*fskyi(nisk,j)
418 END DO
419 niskyl=niskyl+1
420 isky(niskyl)=n4
421 DO j=1,nfskyi
422 fskyi(niskyl,j)=phi4*fskyi(nisk,j)
423 END DO
424C
425 ELSE
426C---------------
427C------ Hexa --
428C---------------
429 n1=ixs(2,n)
430 n2=ixs(3,n)
431 n3=ixs(4,n)
432 n4=ixs(5,n)
433 n5=ixs(6,n)
434 n6=ixs(7,n)
435 n7=ixs(8,n)
436 n8=ixs(9,n)
437C
438 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
439C
440C-----------------------------------------------
441 ir=irst(1,np-first_sphsol+1)
442 is=irst(2,np-first_sphsol+1)
443 it=irst(3,np-first_sphsol+1)
444 ksi = a_gauss(ir,nsphdir)
445 eta = a_gauss(is,nsphdir)
446 zeta = a_gauss(it,nsphdir)
447C
448 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
449 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
450 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
451 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
452 phi5=one_over_8*(one-ksi)*(one+eta)*(one-zeta)
453 phi6=one_over_8*(one-ksi)*(one+eta)*(one+zeta)
454 phi7=one_over_8*(one+ksi)*(one+eta)*(one+zeta)
455 phi8=one_over_8*(one+ksi)*(one+eta)*(one-zeta)
456C
457 niskyl=niskyl+1
458 isky(niskyl)=n1
459 DO j=1,nfskyi
460 fskyi(niskyl,j)=phi1*fskyi(nisk,j)
461 END DO
462 niskyl=niskyl+1
463 isky(niskyl)=n2
464 DO j=1,nfskyi
465 fskyi(niskyl,j)=phi2*fskyi(nisk,j)
466 END DO
467 niskyl=niskyl+1
468 isky(niskyl)=n3
469 DO j=1,nfskyi
470 fskyi(niskyl,j)=phi3*fskyi(nisk,j)
471 END DO
472 niskyl=niskyl+1
473 isky(niskyl)=n4
474 DO j=1,nfskyi
475 fskyi(niskyl,j)=phi4*fskyi(nisk,j)
476 END DO
477 niskyl=niskyl+1
478 isky(niskyl)=n5
479 DO j=1,nfskyi
480 fskyi(niskyl,j)=phi5*fskyi(nisk,j)
481 END DO
482 niskyl=niskyl+1
483 isky(niskyl)=n6
484 DO j=1,nfskyi
485 fskyi(niskyl,j)=phi6*fskyi(nisk,j)
486 END DO
487 niskyl=niskyl+1
488 isky(niskyl)=n7
489 DO j=1,nfskyi
490 fskyi(niskyl,j)=phi7*fskyi(nisk,j)
491 END DO
492 niskyl=niskyl+1
493 isky(niskyl)=n8
494 DO j=1,nfskyi
495 fskyi(niskyl,j)=phi8*fskyi(nisk,j)
496 END DO
497 DO j=1,nfskyi
498 fskyi(nisk,j)=zero
499 END DO
500C
501 ENDIF
502 END IF
503 END IF
504 END DO
505 END IF
506 RETURN
507 END SUBROUTINE soltosphf
508!||====================================================================
509!|| soltosphp ../engine/source/elements/sph/soltosph.F
510!||--- called by ------------------------------------------------------
511!|| forintp ../engine/source/elements/forintp.F
512!||--- calls -----------------------------------------------------
513!|| initbuf ../engine/share/resol/initbuf.F
514!|| sig_heph1 ../engine/source/elements/sph/soltosph_hour.F
515!|| srep2glo ../engine/source/elements/sph/srep2glo.F
516!|| startimeg ../engine/source/system/timer.F
517!|| stoptimeg ../engine/source/system/timer.F
518!||--- uses -----------------------------------------------------
519!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
520!|| element_mod ../common_source/modules/elements/element_mod.F90
521!|| initbuf_mod ../engine/share/resol/initbuf.F
522!||====================================================================
523 SUBROUTINE soltosphp(
524 . X ,SPBUF ,IXS ,KXSP ,IPARTSP ,
525 . IRST ,ELBUF_TAB,IPARG ,NGROUNC ,IGROUNC ,
526 . SOL2SPH,WA ,PM)
527C-----------------------------------------------
528 USE initbuf_mod
529 USE elbufdef_mod
530 use element_mod , only : nixs
531C-----------------------------------------------
532C I m p l i c i t T y p e s
533C-----------------------------------------------
534#include "implicit_f.inc"
535C-----------------------------------------------
536C G l o b a l P a r a m e t e r s
537C-----------------------------------------------
538#include "mvsiz_p.inc"
539C-----------------------------------------------
540C C o m m o n B l o c k s
541C-----------------------------------------------
542#include "com01_c.inc"
543#include "com08_c.inc"
544#include "param_c.inc"
545#include "sphcom.inc"
546#include "task_c.inc"
547#include "vect01_c.inc"
548C-----------------------------------------------
549C D u m m y A r g u m e n t s
550C-----------------------------------------------
551 INTEGER IXS(NIXS,*), KXSP(NISP,*),
552 . IPARTSP(*), IRST(3,*), IPARG(NPARG,*), NGROUNC,
553 . IGROUNC(*), SOL2SPH(2,*)
554 my_real
555 . X(3,*), SPBUF(NSPBUF,*), WA(KWASPH,*),PM(NPROPM,*)
556 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
557C-----------------------------------------------
558C L o c a l V a r i a b l e s
559C-----------------------------------------------
560 INTEGER I, N, KP, NG, MG, J, NP, KFT, IG,
561 . NEL, OFFSET, MLW, IPLA,NELSP,K,IR,IS,IT,NSPHDIR,
562 . nptr,npts,nptt,ii(6),jj(6)
563 my_real
564 . rhon, rhoo, divv,
565 . r11(mvsiz),r12(mvsiz),r13(mvsiz),
566 . r21(mvsiz),r22(mvsiz),r23(mvsiz),
567 . r31(mvsiz),r32(mvsiz),r33(mvsiz),
568 . t11(mvsiz),t12(mvsiz),t13(mvsiz),
569 . t21(mvsiz),t22(mvsiz),t23(mvsiz),
570 . t31(mvsiz),t32(mvsiz),t33(mvsiz),
571 . rx(mvsiz),sx(mvsiz),tx(mvsiz),
572 . ry(mvsiz),sy(mvsiz),ty(mvsiz),
573 . rz(mvsiz),sz(mvsiz),tz(mvsiz),
574 . g11,g22,g33,g12,g21,g23,g32,g13,g31,
575 . s11,s22,s33,s12,s21,s23,s32,s13,s31,
576 . l11,l22,l33,l12,l23,l13,
577 . siglo(mvsiz,6), straglo(mvsiz,6), angl(mvsiz,6),
578 . dglo24(mvsiz,6),sig_heph(mvsiz,6,7),
579 . jr0(mvsiz),js0(mvsiz),jt0(mvsiz),sig_heph_glo(mvsiz,6,7),
580 . zeta,eta,ksi,sig_ha8(mvsiz,3,3,3,6)
581C
582C-----
583 TYPE(g_bufel_) ,POINTER :: GBUF, GBUFSP
584 TYPE(L_BUFEL_) ,POINTER :: LBUF, LBUFSP, LBUF2
585 TYPE(BUF_MAT_) ,POINTER :: MBUF, MBUFSP
586C-----------------------------------------------
587 my_real A_GAUSS(9,9)
588 DATA A_GAUSS /
589 1 0. ,0. ,0. ,
590 1 0. ,0. ,0. ,
591 1 0. ,0. ,0. ,
592 2 -.577350269189626,0.577350269189626,0. ,
593 2 0. ,0. ,0. ,
594 2 0. ,0. ,0. ,
595 3 -.774596669241483,0. ,0.774596669241483,
596 3 0. ,0. ,0. ,
597 3 0. ,0. ,0. ,
598 4 -.861136311594053,-.339981043584856,0.339981043584856,
599 4 0.861136311594053,0. ,0. ,
600 4 0. ,0. ,0. ,
601 5 -.906179845938664,-.538469310105683,0. ,
602 5 0.538469310105683,0.906179845938664,0. ,
603 5 0. ,0. ,0. ,
604 6 -.932469514203152,-.661209386466265,-.238619186083197,
605 6 0.238619186083197,0.661209386466265,0.932469514203152,
606 6 0. ,0. ,0. ,
607 7 -.949107912342759,-.741531185599394,-.405845151377397,
608 7 0. ,0.405845151377397,0.741531185599394,
609 7 0.949107912342759,0. ,0. ,
610 8 -.960289856497536,-.796666477413627,-.525532409916329,
611 8 -.183434642495650,0.183434642495650,0.525532409916329,
612 8 0.796666477413627,0.960289856497536,0. ,
613 9 -.968160239507626,-.836031107326636,-.613371432700590,
614 9 -.324253423403809,0. ,0.324253423403809,
615 9 0.613371432700590,0.836031107326636,0.968160239507626/
616C-----------------------------------------------
617!$OMP DO SCHEDULE(DYNAMIC,1)
618 DO ig = 1, ngrounc
619 ng = igrounc(ig)
620 IF(iparg(8,ng)==1)GOTO 300
621 IF (iddw>0) CALL startimeg(ng)
622 offset = 0
623 ity = iparg(5,ng)
624 ipartsph= iparg(69,ng)
625 IF(ity==1.AND.ipartsph/=0) THEN
626C
627C---
628 CALL initbuf(iparg ,ng ,
629 2 mlw ,nel ,nft ,iad ,ity ,
630 3 npt ,jale ,ismstr ,jeul ,jtur ,
631 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
632 5 nvaux ,jpor ,jcvt ,jclose ,ipla ,
633 6 irep ,iint ,igtyp ,israt ,isrot ,
634 7 icsen ,isorth ,isorthg ,ifailure,jsms )
635 lft = 1
636 llt = min(nvsiz,nel)
637!
638 DO i=1,6
639 ii(i) = nel*(i-1)
640 ENDDO
641!
642C-----------
643 gbuf => elbuf_tab(ng)%GBUF
644 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
645 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
646C-----------
647 CALL srep2glo(
648 1 x, ixs(1,nft+1),gbuf%GAMA, rx,
649 2 ry, rz, sx, sy,
650 3 sz, tx, ty, tz,
651 4 r11, r12, r13, r21,
652 5 r22, r23, r31, r32,
653 6 r33, t11, t12, t13,
654 7 t21, t22, t23, t31,
655 8 t32, t33, jr0, js0,
656 9 jt0, nel, lft, llt,
657 a jhbe, jcvt, isorth)
658C
659C----------- HEPH -------
660 IF (jhbe==24) THEN
661C------------------------
662 sig_heph(1:mvsiz,1:6,1:7) = zero
663 CALL sig_heph1(
664 1 jr0, js0, jt0, gbuf%SIG,
665 2 gbuf%HOURG,sig_heph, pm, ixs,
666 3 ii, nel, lft, llt)
667C
668 IF(isorth==0)THEN
669 DO j=1,7
670 DO i=lft,llt
671C hourglass stress in corotational system
672 l11 =sig_heph(i,1,j)
673 l22 =sig_heph(i,2,j)
674 l33 =sig_heph(i,3,j)
675 l12 =sig_heph(i,4,j)
676 l23 =sig_heph(i,5,j)
677 l13 =sig_heph(i,6,j)
678 s11 =l11*r11(i)+l12*r12(i)+l13*r13(i)
679 s12 =l11*r21(i)+l12*r22(i)+l13*r23(i)
680 s13 =l11*r31(i)+l12*r32(i)+l13*r33(i)
681 s21 =l12*r11(i)+l22*r12(i)+l23*r13(i)
682 s22 =l12*r21(i)+l22*r22(i)+l23*r23(i)
683 s23 =l12*r31(i)+l22*r32(i)+l23*r33(i)
684 s31 =l13*r11(i)+l23*r12(i)+l33*r13(i)
685 s32 =l13*r21(i)+l23*r22(i)+l33*r23(i)
686 s33 =l13*r31(i)+l23*r32(i)+l33*r33(i)
687 sig_heph_glo(i,1,j)=r11(i)*s11+r12(i)*s21+r13(i)*s31
688 sig_heph_glo(i,2,j)=r21(i)*s12+r22(i)*s22+r23(i)*s32
689 sig_heph_glo(i,3,j)=r31(i)*s13+r32(i)*s23+r33(i)*s33
690 sig_heph_glo(i,4,j)=r11(i)*s12+r12(i)*s22+r13(i)*s32
691 sig_heph_glo(i,5,j)=r21(i)*s13+r22(i)*s23+r23(i)*s33
692 sig_heph_glo(i,6,j)=r11(i)*s13+r12(i)*s23+r13(i)*s33
693 END DO
694 END DO
695 ELSE
696 DO j=1,7
697 DO i=lft,llt
698C hourglass stress in orthotropic system
699 l11 =sig_heph(i,1,j)
700 l22 =sig_heph(i,2,j)
701 l33 =sig_heph(i,3,j)
702 l12 =sig_heph(i,4,j)
703 l23 =sig_heph(i,5,j)
704 l13 =sig_heph(i,6,j)
705 s11 =l11*t11(i)+l12*t12(i)+l13*t13(i)
706 s12 =l11*t21(i)+l12*t22(i)+l13*t23(i)
707 s13 =l11*t31(i)+l12*t32(i)+l13*t33(i)
708 s21 =l12*t11(i)+l22*t12(i)+l23*t13(i)
709 s22 =l12*t21(i)+l22*t22(i)+l23*t23(i)
710 s23 =l12*t31(i)+l22*t32(i)+l23*t33(i)
711 s31 =l13*t11(i)+l23*t12(i)+l33*t13(i)
712 s32 =l13*t21(i)+l23*t22(i)+l33*t23(i)
713 s33 =l13*t31(i)+l23*t32(i)+l33*t33(i)
714 sig_heph_glo(i,1,j)=t11(i)*s11+t12(i)*s21+t13(i)*s31
715 sig_heph_glo(i,2,j)=t21(i)*s12+t22(i)*s22+t23(i)*s32
716 sig_heph_glo(i,3,j)=t31(i)*s13+t32(i)*s23+t33(i)*s33
717 sig_heph_glo(i,4,j)=t11(i)*s12+t12(i)*s22+t13(i)*s32
718 sig_heph_glo(i,5,j)=t21(i)*s13+t22(i)*s23+t23(i)*s33
719 sig_heph_glo(i,6,j)=t11(i)*s13+t12(i)*s23+t13(i)*s33
720 END DO
721 END DO
722 ENDIF
723C----------- HA8 -------
724 ELSEIF (jhbe==14) THEN
725C------------------------
726 nptr = elbuf_tab(ng)%NPTR
727 npts = elbuf_tab(ng)%NPTS
728 nptt = elbuf_tab(ng)%NPTT
729 IF(isorth==0)THEN
730 DO ir=1,nptr
731 DO is=1,npts
732 DO it=1,nptt
733C--------- convention HA8 : ETA-> r / ZETA -> s / KSI -> t
734C--------- convention SOL2SPH : KSI-> r / ETA -> s / ZETA -> t
735C--------- on permute KSI,ETA,ZETA
736 lbuf2 => elbuf_tab(ng)%BUFLY(1)%LBUF(it,ir,is)
737C---------
738 DO i=lft,llt
739C hourglass stress in corotational system
740 l11 =lbuf2%SIG(ii(1)+i)
741 l22 =lbuf2%SIG(ii(2)+i)
742 l33 =lbuf2%SIG(ii(3)+i)
743 l12 =lbuf2%SIG(ii(4)+i)
744 l23 =lbuf2%SIG(ii(5)+i)
745 l13 =lbuf2%SIG(ii(6)+i)
746 s11 =l11*r11(i)+l12*r12(i)+l13*r13(i)
747 s12 =l11*r21(i)+l12*r22(i)+l13*r23(i)
748 s13 =l11*r31(i)+l12*r32(i)+l13*r33(i)
749 s21 =l12*r11(i)+l22*r12(i)+l23*r13(i)
750 s22 =l12*r21(i)+l22*r22(i)+l23*r23(i)
751 s23 =l12*r31(i)+l22*r32(i)+l23*r33(i)
752 s31 =l13*r11(i)+l23*r12(i)+l33*r13(i)
753 s32 =l13*r21(i)+l23*r22(i)+l33*r23(i)
754 s33 =l13*r31(i)+l23*r32(i)+l33*r33(i)
755 sig_ha8(i,ir,is,it,1)=r11(i)*s11+r12(i)*s21+r13(i)*s31
756 sig_ha8(i,ir,is,it,2)=r21(i)*s12+r22(i)*s22+r23(i)*s32
757 sig_ha8(i,ir,is,it,3)=r31(i)*s13+r32(i)*s23+r33(i)*s33
758 sig_ha8(i,ir,is,it,4)=r11(i)*s12+r12(i)*s22+r13(i)*s32
759 sig_ha8(i,ir,is,it,5)=r21(i)*s13+r22(i)*s23+r23(i)*s33
760 sig_ha8(i,ir,is,it,6)=r11(i)*s13+r12(i)*s23+r13(i)*s33
761 END DO
762 END DO
763 END DO
764 END DO
765 ELSE
766 DO ir=1,nptr
767 DO is=1,npts
768 DO it=1,nptt
769 lbuf2 => elbuf_tab(ng)%BUFLY(1)%LBUF(ir,is,it)
770 DO i=lft,llt
771C hourglass stress in orthotropic system
772 l11 =lbuf2%SIG(ii(1)+i)
773 l22 =lbuf2%SIG(ii(2)+i)
774 l33 =lbuf2%SIG(ii(3)+i)
775 l12 =lbuf2%SIG(ii(4)+i)
776 l23 =lbuf2%SIG(ii(5)+i)
777 l13 =lbuf2%SIG(ii(6)+i)
778 s11 =l11*t11(i)+l12*t12(i)+l13*t13(i)
779 s12 =l11*t21(i)+l12*t22(i)+l13*t23(i)
780 s13 =l11*t31(i)+l12*t32(i)+l13*t33(i)
781 s21 =l12*t11(i)+l22*t12(i)+l23*t13(i)
782 s22 =l12*t21(i)+l22*t22(i)+l23*t23(i)
783 s23 =l12*t31(i)+l22*t32(i)+l23*t33(i)
784 s31 =l13*t11(i)+l23*t12(i)+l33*t13(i)
785 s32 =l13*t21(i)+l23*t22(i)+l33*t23(i)
786 s33 =l13*t31(i)+l23*t32(i)+l33*t33(i)
787 sig_ha8(i,ir,is,it,1)=t11(i)*s11+t12(i)*s21+t13(i)*s31
788 sig_ha8(i,ir,is,it,2)=t21(i)*s12+t22(i)*s22+t23(i)*s32
789 sig_ha8(i,ir,is,it,3)=t31(i)*s13+t32(i)*s23+t33(i)*s33
790 sig_ha8(i,ir,is,it,4)=t11(i)*s12+t12(i)*s22+t13(i)*s32
791 sig_ha8(i,ir,is,it,5)=t21(i)*s13+t22(i)*s23+t23(i)*s33
792 sig_ha8(i,ir,is,it,6)=t11(i)*s13+t12(i)*s23+t13(i)*s33
793 END DO
794 END DO
795 END DO
796 END DO
797 ENDIF
798C----------- Standard formulation -------
799 ELSEIF (jcvt == 0)THEN
800C------------------------------------------
801 DO i=lft,llt
802C mean stress in global system
803 siglo(i,1) =gbuf%SIG(ii(1)+i)
804 siglo(i,2) =gbuf%SIG(ii(2)+i)
805 siglo(i,3) =gbuf%SIG(ii(3)+i)
806 siglo(i,4) =gbuf%SIG(ii(4)+i)
807 siglo(i,5) =gbuf%SIG(ii(5)+i)
808 siglo(i,6) =gbuf%SIG(ii(6)+i)
809 END DO
810C----------- Corotational formulation -------
811 ELSE
812C--------------------------------------------
813C
814 IF (isorth== 0) THEN
815 DO i=lft,llt
816C mean stress in corotational system
817 l11 =gbuf%SIG(ii(1)+i)
818 l22 =gbuf%SIG(ii(2)+i)
819 l33 =gbuf%SIG(ii(3)+i)
820 l12 =gbuf%SIG(ii(4)+i)
821 l23 =gbuf%SIG(ii(5)+i)
822 l13 =gbuf%SIG(ii(6)+i)
823 s11 =l11*r11(i)+l12*r12(i)+l13*r13(i)
824 s12 =l11*r21(i)+l12*r22(i)+l13*r23(i)
825 s13 =l11*r31(i)+l12*r32(i)+l13*r33(i)
826 s21 =l12*r11(i)+l22*r12(i)+l23*r13(i)
827 s22 =l12*r21(i)+l22*r22(i)+l23*r23(i)
828 s23 =l12*r31(i)+l22*r32(i)+l23*r33(i)
829 s31 =l13*r11(i)+l23*r12(i)+l33*r13(i)
830 s32 =l13*r21(i)+l23*r22(i)+l33*r23(i)
831 s33 =l13*r31(i)+l23*r32(i)+l33*r33(i)
832 siglo(i,1)=r11(i)*s11+r12(i)*s21+r13(i)*s31
833 siglo(i,2)=r21(i)*s12+r22(i)*s22+r23(i)*s32
834 siglo(i,3)=r31(i)*s13+r32(i)*s23+r33(i)*s33
835 siglo(i,4)=r11(i)*s12+r12(i)*s22+r13(i)*s32
836 siglo(i,5)=r21(i)*s13+r22(i)*s23+r23(i)*s33
837 siglo(i,6)=r11(i)*s13+r12(i)*s23+r13(i)*s33
838 END DO
839 ELSE
840 DO i=lft,llt
841C mean stress in orthotropic system
842 l11 =gbuf%SIG(ii(1)+i)
843 l22 =gbuf%SIG(ii(2)+i)
844 l33 =gbuf%SIG(ii(3)+i)
845 l12 =gbuf%SIG(ii(4)+i)
846 l23 =gbuf%SIG(ii(5)+i)
847 l13 =gbuf%SIG(ii(6)+i)
848 s11 =l11*t11(i)+l12*t12(i)+l13*t13(i)
849 s12 =l11*t21(i)+l12*t22(i)+l13*t23(i)
850 s13 =l11*t31(i)+l12*t32(i)+l13*t33(i)
851 s21 =l12*t11(i)+l22*t12(i)+l23*t13(i)
852 s22 =l12*t21(i)+l22*t22(i)+l23*t23(i)
853 s23 =l12*t31(i)+l22*t32(i)+l23*t33(i)
854 s31 =l13*t11(i)+l23*t12(i)+l33*t13(i)
855 s32 =l13*t21(i)+l23*t22(i)+l33*t23(i)
856 s33 =l13*t31(i)+l23*t32(i)+l33*t33(i)
857 siglo(i,1)=t11(i)*s11+t12(i)*s21+t13(i)*s31
858 siglo(i,2)=t21(i)*s12+t22(i)*s22+t23(i)*s32
859 siglo(i,3)=t31(i)*s13+t32(i)*s23+t33(i)*s33
860 siglo(i,4)=t11(i)*s12+t12(i)*s22+t13(i)*s32
861 siglo(i,5)=t21(i)*s13+t22(i)*s23+t23(i)*s33
862 siglo(i,6)=t11(i)*s13+t12(i)*s23+t13(i)*s33
863 END DO
864 END IF
865C------------------------
866 ENDIF
867C-----------------------
868 IF(elbuf_tab(ng)%BUFLY(1)%L_STRA > 0)THEN
869 IF(jcvt == 0)THEN
870 DO i=lft,llt
871 straglo(i,1)=lbuf%STRA(ii(1)+i)
872 straglo(i,2)=lbuf%STRA(ii(2)+i)
873 straglo(i,3)=lbuf%STRA(ii(3)+i)
874 straglo(i,4)=lbuf%STRA(ii(4)+i)
875 straglo(i,5)=lbuf%STRA(ii(5)+i)
876 straglo(i,6)=lbuf%STRA(ii(6)+i)
877 END DO
878 ELSEIF(isorth==0)THEN
879 DO i=lft,llt
880C
881C strain in corotational system
882 l11 =lbuf%STRA(ii(1)+i)
883 l22 =lbuf%STRA(ii(2)+i)
884 l33 =lbuf%STRA(ii(3)+i)
885 l12 =half*lbuf%STRA(ii(4)+i)
886 l23 =half*lbuf%STRA(ii(5)+i)
887 l13 =half*lbuf%STRA(ii(6)+i)
888 s11 =l11*r11(i)+l12*r12(i)+l13*r13(i)
889 s12 =l11*r21(i)+l12*r22(i)+l13*r23(i)
890 s13 =l11*r31(i)+l12*r32(i)+l13*r33(i)
891 s21 =l12*r11(i)+l22*r12(i)+l23*r13(i)
892 s22 =l12*r21(i)+l22*r22(i)+l23*r23(i)
893 s23 =l12*r31(i)+l22*r32(i)+l23*r33(i)
894 s31 =l13*r11(i)+l23*r12(i)+l33*r13(i)
895 s32 =l13*r21(i)+l23*r22(i)+l33*r23(i)
896 s33 =l13*r31(i)+l23*r32(i)+l33*r33(i)
897 straglo(i,1)=r11(i)*s11+r12(i)*s21+r13(i)*s31
898 straglo(i,2)=r21(i)*s12+r22(i)*s22+r23(i)*s32
899 straglo(i,3)=r31(i)*s13+r32(i)*s23+r33(i)*s33
900 straglo(i,4)=two*(r11(i)*s12+r12(i)*s22+r13(i)*s32)
901 straglo(i,5)=two*(r21(i)*s13+r22(i)*s23+r23(i)*s33)
902 straglo(i,6)=two*(r11(i)*s13+r12(i)*s23+r13(i)*s33)
903 END DO
904 ELSE
905 DO i=lft,llt
906C
907C strain in orthotropic system
908 l11 =lbuf%STRA(ii(1)+i)
909 l22 =lbuf%STRA(ii(2)+i)
910 l33 =lbuf%STRA(ii(3)+i)
911 l12 =half*lbuf%STRA(ii(4)+i)
912 l23 =half*lbuf%STRA(ii(5)+i)
913 l13 =half*lbuf%STRA(ii(6)+i)
914 s11 =l11*t11(i)+l12*t12(i)+l13*t13(i)
915 s12 =l11*t21(i)+l12*t22(i)+l13*t23(i)
916 s13 =l11*t31(i)+l12*t32(i)+l13*t33(i)
917 s21 =l12*t11(i)+l22*t12(i)+l23*t13(i)
918 s22 =l12*t21(i)+l22*t22(i)+l23*t23(i)
919 s23 =l12*t31(i)+l22*t32(i)+l23*t33(i)
920 s31 =l13*t11(i)+l23*t12(i)+l33*t13(i)
921 s32 =l13*t21(i)+l23*t22(i)+l33*t23(i)
922 s33 =l13*t31(i)+l23*t32(i)+l33*t33(i)
923 straglo(i,1)=t11(i)*s11+t12(i)*s21+t13(i)*s31
924 straglo(i,2)=t21(i)*s12+t22(i)*s22+t23(i)*s32
925 straglo(i,3)=t31(i)*s13+t32(i)*s23+t33(i)*s33
926 straglo(i,4)=two*(t11(i)*s12+t12(i)*s22+t13(i)*s32)
927 straglo(i,5)=two*(t21(i)*s13+t22(i)*s23+t23(i)*s33)
928 straglo(i,6)=two*(t11(i)*s13+t12(i)*s23+t13(i)*s33)
929 END DO
930 END IF
931 END IF
932C-----------
933
934 IF(elbuf_tab(ng)%BUFLY(1)%L_ANG > 0)THEN
935 IF(jcvt == 0 .AND. isorth == 0)THEN
936 DO i=lft,llt
937 g11=lbuf%ANG(ii(1)+i)
938 g21=lbuf%ANG(ii(2)+i)
939 g31=lbuf%ANG(ii(3)+i)
940 g12=lbuf%ANG(ii(4)+i)
941 g22=lbuf%ANG(ii(5)+i)
942 g32=lbuf%ANG(ii(6)+i)
943 g13=g21*g32-g31*g22
944 g23=g31*g12-g11*g32
945 g33=g11*g22-g21*g12
946C MATRICE DE PASSAGE GLOBAL -> ANG
947 s11=rx(i)*g11+sx(i)*g21+tx(i)*g31
948 s12=rx(i)*g12+sx(i)*g22+tx(i)*g32
949 s13=rx(i)*g13+sx(i)*g23+tx(i)*g33
950 s21=ry(i)*g11+sy(i)*g21+ty(i)*g31
951 s22=ry(i)*g12+sy(i)*g22+ty(i)*g32
952 s23=ry(i)*g13+sy(i)*g23+ty(i)*g33
953 s31=rz(i)*g11+sz(i)*g21+tz(i)*g31
954 s32=rz(i)*g12+sz(i)*g22+tz(i)*g32
955 s33=rz(i)*g13+sz(i)*g23+tz(i)*g33
956 angl(i,1)=s11
957 angl(i,2)=s21
958 angl(i,3)=s31
959 angl(i,4)=s12
960 angl(i,5)=s22
961 angl(i,6)=s32
962 END DO
963 ELSEIF(jcvt /=0 .AND. isorth == 0)THEN
964 DO i=lft,llt
965 g11=lbuf%ANG(ii(1)+i)
966 g21=lbuf%ANG(ii(2)+i)
967 g31=lbuf%ANG(ii(3)+i)
968 g12=lbuf%ANG(ii(4)+i)
969 g22=lbuf%ANG(ii(5)+i)
970 g32=lbuf%ANG(ii(6)+i)
971 g13=g21*g32-g31*g22
972 g23=g31*g12-g11*g32
973 g33=g11*g22-g21*g12
974C MATRICE DE PASSAGE GLOBAL -> ANG
975 s11=r11(i)*g11+r12(i)*g21+r13(i)*g31
976 s12=r11(i)*g12+r12(i)*g22+r13(i)*g32
977 s13=r11(i)*g13+r12(i)*g23+r13(i)*g33
978 s21=r21(i)*g11+r22(i)*g21+r23(i)*g31
979 s22=r21(i)*g12+r22(i)*g22+r23(i)*g32
980 s23=r21(i)*g13+r22(i)*g23+r23(i)*g33
981 s31=r31(i)*g11+r32(i)*g21+r33(i)*g31
982 s32=r31(i)*g12+r32(i)*g22+r33(i)*g32
983 s33=r31(i)*g13+r32(i)*g23+r33(i)*g33
984 angl(i,1)=s11
985 angl(i,2)=s21
986 angl(i,3)=s31
987 angl(i,4)=s12
988 angl(i,5)=s22
989 angl(i,6)=s32
990 END DO
991 ELSE
992 DO i=lft,llt
993C
994C ISorth /=0 (ANG is given both for solid & sph wrt orthotropic system)
995C MATRICE DE PASSAGE ORTHOTROPE -> ANG
996 angl(i,1)=lbuf%ANG(ii(1)+i)
997 angl(i,2)=lbuf%ANG(ii(2)+i)
998 angl(i,3)=lbuf%ANG(ii(3)+i)
999 angl(i,4)=lbuf%ANG(ii(4)+i)
1000 angl(i,5)=lbuf%ANG(ii(5)+i)
1001 angl(i,6)=lbuf%ANG(ii(6)+i)
1002 END DO
1003 END IF
1004 END IF
1005C-----------
1006 IF(elbuf_tab(ng)%BUFLY(1)%L_DGLO > 0)THEN
1007
1008 IF(jcvt == 0 .AND. isorth == 0)THEN
1009 DO i=lft,llt
1010C TENSOR wrt ISOPARAMETRIC SYSTEM
1011 g11=lbuf%DGLO(ii(1)+i)
1012 g22=lbuf%DGLO(ii(2)+i)
1013 g33=lbuf%DGLO(ii(3)+i)
1014 g12=lbuf%DGLO(ii(4)+i)
1015 g23=lbuf%DGLO(ii(5)+i)
1016 g13=lbuf%DGLO(ii(6)+i)
1017 s11=g11*rx(i)+g12*sx(i)+g13*tx(i)
1018 s12=g11*ry(i)+g12*sy(i)+g13*ty(i)
1019 s13=g11*rz(i)+g12*sz(i)+g13*tz(i)
1020 s21=g12*rx(i)+g22*sx(i)+g23*tx(i)
1021 s22=g12*ry(i)+g22*sy(i)+g23*ty(i)
1022 s23=g12*rz(i)+g22*sz(i)+g23*tz(i)
1023 s31=g13*rx(i)+g23*sx(i)+g33*tx(i)
1024 s32=g13*ry(i)+g23*sy(i)+g33*ty(i)
1025 s33=g13*rz(i)+g23*sz(i)+g33*tz(i)
1026C TENSOR wrt GLOBAL SYSTEM
1027 dglo24(i,1)=rx(i)*s11+sx(i)*s21+tx(i)*s31
1028 dglo24(i,2)=ry(i)*s12+sy(i)*s22+ty(i)*s32
1029 dglo24(i,3)=rz(i)*s13+sz(i)*s23+tz(i)*s33
1030 dglo24(i,4)=rx(i)*s12+sx(i)*s22+tx(i)*s32
1031 dglo24(i,5)=ry(i)*s13+sy(i)*s23+ty(i)*s33
1032 dglo24(i,6)=rx(i)*s13+sx(i)*s23+tx(i)*s33
1033 END DO
1034 ELSEIF(jcvt /=0 .AND. isorth == 0)THEN
1035 DO i=lft,llt
1036C TENSOR wrt COROTATIONAL SYSTEM
1037 g11=lbuf%DGLO(ii(1)+i)
1038 g22=lbuf%DGLO(ii(2)+i)
1039 g33=lbuf%DGLO(ii(3)+i)
1040 g12=lbuf%DGLO(ii(4)+i)
1041 g23=lbuf%DGLO(ii(5)+i)
1042 g13=lbuf%DGLO(ii(6)+i)
1043 s11=g11*r11(i)+g12*r12(i)+g13*r13(i)
1044 s12=g11*r21(i)+g12*r22(i)+g13*r23(i)
1045 s13=g11*r31(i)+g12*r32(i)+g13*r33(i)
1046 s21=g12*r11(i)+g22*r12(i)+g23*r13(i)
1047 s22=g12*r21(i)+g22*r22(i)+g23*r23(i)
1048 s23=g12*r31(i)+g22*r32(i)+g23*r33(i)
1049 s31=g13*r11(i)+g23*r12(i)+g33*r13(i)
1050 s32=g13*r21(i)+g23*r22(i)+g33*r23(i)
1051 s33=g13*r31(i)+g23*r32(i)+g33*r33(i)
1052C TENSOR wrt GLOBAL SYSTEM
1053 dglo24(i,1)=r11(i)*s11+r12(i)*s21+r13(i)*s31
1054 dglo24(i,2)=r21(i)*s12+r22(i)*s22+r23(i)*s32
1055 dglo24(i,3)=r31(i)*s13+r32(i)*s23+r33(i)*s33
1056 dglo24(i,4)=r11(i)*s12+r12(i)*s22+r13(i)*s32
1057 dglo24(i,5)=r21(i)*s13+r22(i)*s23+r23(i)*s33
1058 dglo24(i,6)=r11(i)*s13+r12(i)*s23+r13(i)*s33
1059 END DO
1060 ELSE
1061C
1062C ISorth /=0 (orthotropic system is the same for solid & sph)
1063C TENSOR wrt ORTHOTROPIC SYSTEM
1064 DO i=lft,llt
1065 dglo24(i,1)=lbuf%DGLO(ii(1)+i)
1066 dglo24(i,2)=lbuf%DGLO(ii(2)+i)
1067 dglo24(i,3)=lbuf%DGLO(ii(3)+i)
1068 dglo24(i,4)=lbuf%DGLO(ii(4)+i)
1069 dglo24(i,5)=lbuf%DGLO(ii(5)+i)
1070 dglo24(i,6)=lbuf%DGLO(ii(6)+i)
1071 END DO
1072 END IF
1073 END IF
1074C-----------
1075 DO i=lft,llt
1076 IF(gbuf%OFF(i)==zero) cycle
1077 n=nft+i
1078C
1079C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
1080 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
1081 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
1082C
1083 np=sol2sph(1,n)+kp
1084 IF (kxsp(2,np)/=0) THEN
1085C
1086 mg =mod(-kxsp(2,np),ngroup+1)
1087 nelsp=iparg(2,mg)
1088 kft=iparg(3,mg)
1089 gbufsp => elbuf_tab(mg)%GBUF
1090 lbufsp => elbuf_tab(mg)%BUFLY(1)%LBUF(1,1,1)
1091 mbufsp => elbuf_tab(mg)%BUFLY(1)%MAT(1,1,1)
1092 j=np-kft
1093 rhon = gbuf%RHO(i)
1094 rhoo = wa(10,np)
1095 divv = (rhoo-rhon)/max(em30,rhoo*dt1)
1096 wa(13,np) = divv
1097 wa(14,np) = zero
1098 spbuf(2,np) = rhon
1099 gbufsp%RHO(j) = rhon
1100C
1101C group was not computed
1102C (there is no cloud active particle within the group)
1103C IF(IPARG(8,MG)==1)CYCLE
1104C
1105 gbufsp%EINT(j) =gbuf%EINT(i)
1106C
1107!
1108 DO k=1,6
1109 jj(k) = nelsp*(k-1)
1110 ENDDO
1111!
1112C-----------
1113 IF (jhbe==14) THEN
1114C HA8 stress
1115 ir=irst(1,np-first_sphsol+1)
1116 is=irst(2,np-first_sphsol+1)
1117 it=irst(3,np-first_sphsol+1)
1118 DO k=1,6
1119 gbufsp%SIG(jj(k)+j)=sig_ha8(i,ir,is,it,k)
1120 ENDDO
1121 ELSEIF (jhbe==24) THEN
1122C Hourglass stress contribution for HEPH
1123 ir=irst(1,np-first_sphsol+1)
1124 is=irst(2,np-first_sphsol+1)
1125 it=irst(3,np-first_sphsol+1)
1126C--- Permutation (KSI-> ETA;ETA->ZETA;ZETA->KSI) for consistencvy with HEPH convention
1127 eta = a_gauss(ir,nsphdir)
1128 zeta = a_gauss(is,nsphdir)
1129 ksi = a_gauss(it,nsphdir)
1130C-----------
1131 DO k=1,6
1132 gbufsp%SIG(jj(k)+j) = sig_heph_glo(i,k,1)
1133 . +zeta*sig_heph_glo(i,k,2)
1134 . +eta*sig_heph_glo(i,k,3)
1135 . +ksi*sig_heph_glo(i,k,4)
1136 . +zeta*eta*sig_heph_glo(i,k,5)
1137 . +zeta*ksi*sig_heph_glo(i,k,6)
1138 . +eta*ksi*sig_heph_glo(i,k,7)
1139 END DO
1140 ELSE
1141 gbufsp%SIG(jj(1)+j) = siglo(i,1)
1142 gbufsp%SIG(jj(2)+j) = siglo(i,2)
1143 gbufsp%SIG(jj(3)+j) = siglo(i,3)
1144 gbufsp%SIG(jj(4)+j) = siglo(i,4)
1145 gbufsp%SIG(jj(5)+j) = siglo(i,5)
1146 gbufsp%SIG(jj(6)+j) = siglo(i,6)
1147 ENDIF
1148C-----------
1149 wa(1,np)=gbufsp%SIG(jj(1)+j)
1150 wa(2,np)=gbufsp%SIG(jj(2)+j)
1151 wa(3,np)=gbufsp%SIG(jj(3)+j)
1152 wa(4,np)=gbufsp%SIG(jj(4)+j)
1153 wa(5,np)=gbufsp%SIG(jj(5)+j)
1154 wa(6,np)=gbufsp%SIG(jj(6)+j)
1155C
1156C (particles ARE computed with rho_old, Eold...)
1157C-----------
1158 IF(gbuf%G_PLA > 0) gbufsp%PLA(j) = gbuf%PLA(i)
1159 IF(gbuf%G_EPSD> 0) gbufsp%EPSD(j)= gbuf%EPSD(i)
1160 IF(gbuf%G_EPSQ> 0) gbufsp%EPSQ(j)= gbuf%EPSQ(i)
1161C-----------
1162 IF(gbuf%G_GAMA > 0)THEN
1163C
1164C TIJ stores global to orthotropic system
1165 gbufsp%GAMA(jj(1)+j)=t11(i)
1166 gbufsp%GAMA(jj(2)+j)=t21(i)
1167 gbufsp%GAMA(jj(3)+j)=t31(i)
1168 gbufsp%GAMA(jj(4)+j)=t12(i)
1169 gbufsp%GAMA(jj(5)+j)=t22(i)
1170 gbufsp%GAMA(jj(6)+j)=t32(i)
1171 END IF
1172C-----------
1173 IF(elbuf_tab(ng)%BUFLY(1)%L_STRA > 0.AND.
1174 . elbuf_tab(mg)%BUFLY(1)%L_STRA > 0)THEN
1175 lbufsp%STRA(jj(1)+j)=straglo(i,1)
1176 lbufsp%STRA(jj(2)+j)=straglo(i,2)
1177 lbufsp%STRA(jj(3)+j)=straglo(i,3)
1178 lbufsp%STRA(jj(4)+j)=straglo(i,4)
1179 lbufsp%STRA(jj(5)+j)=straglo(i,5)
1180 lbufsp%STRA(jj(6)+j)=straglo(i,6)
1181 END IF
1182C-----------
1183 IF(elbuf_tab(ng)%BUFLY(1)%L_ANG > 0)THEN
1184 lbufsp%ANG(jj(1)+j)=angl(i,1)
1185 lbufsp%ANG(jj(2)+j)=angl(i,2)
1186 lbufsp%ANG(jj(3)+j)=angl(i,3)
1187 lbufsp%ANG(jj(4)+j)=angl(i,4)
1188 lbufsp%ANG(jj(5)+j)=angl(i,5)
1189 lbufsp%ANG(jj(6)+j)=angl(i,6)
1190 END IF
1191C-----------
1192 IF(elbuf_tab(ng)%BUFLY(1)%L_SF > 0)THEN
1193 lbufsp%SF(jj(1)+j)=lbuf%SF(ii(1)+i)
1194 lbufsp%SF(jj(2)+j)=lbuf%SF(ii(2)+i)
1195 lbufsp%SF(jj(3)+j)=lbuf%SF(ii(3)+i)
1196 END IF
1197C-----------
1198 IF(elbuf_tab(ng)%BUFLY(1)%L_DAM > 0)THEN
1199 DO k=1,elbuf_tab(ng)%BUFLY(1)%L_DAM
1200 lbufsp%DAM(jj(k)+j)=lbuf%DAM(ii(k)+i)
1201 ENDDO
1202 END IF
1203C-----------
1204 IF(elbuf_tab(ng)%BUFLY(1)%L_DSUM > 0)
1205 . lbufsp%DSUM(j)=lbuf%DSUM(i)
1206C-----------
1207 IF(elbuf_tab(ng)%BUFLY(1)%L_DGLO > 0)THEN
1208 lbufsp%DGLO(jj(1)+j)=dglo24(i,1)
1209 lbufsp%DGLO(jj(2)+j)=dglo24(i,2)
1210 lbufsp%DGLO(jj(3)+j)=dglo24(i,3)
1211 lbufsp%DGLO(jj(4)+j)=dglo24(i,4)
1212 lbufsp%DGLO(jj(5)+j)=dglo24(i,5)
1213 lbufsp%DGLO(jj(6)+j)=dglo24(i,6)
1214 END IF
1215C-----------
1216 IF(elbuf_tab(ng)%BUFLY(1)%L_ROB > 0)
1217 . lbufsp%ROB(j)=lbuf%ROB(i)
1218C-----------
1219 IF(elbuf_tab(ng)%BUFLY(1)%L_SIGC > 0)THEN
1220C
1221C Law24 / Stress in crack frame (same for SPH and solids)
1222 lbufsp%SIGC(jj(1)+j)=lbuf%SIGC(ii(1)+i)
1223 lbufsp%SIGC(jj(2)+j)=lbuf%SIGC(ii(2)+i)
1224 lbufsp%SIGC(jj(3)+j)=lbuf%SIGC(ii(3)+i)
1225 lbufsp%SIGC(jj(4)+j)=lbuf%SIGC(ii(4)+i)
1226 lbufsp%SIGC(jj(5)+j)=lbuf%SIGC(ii(5)+i)
1227 lbufsp%SIGC(jj(6)+j)=lbuf%SIGC(ii(6)+i)
1228 END IF
1229C-----------
1230 IF(elbuf_tab(ng)%BUFLY(1)%L_CRAK > 0)THEN
1231 lbufsp%CRAK(jj(1)+j)=lbuf%CRAK(ii(1)+i)
1232 lbufsp%CRAK(jj(2)+j)=lbuf%CRAK(ii(2)+i)
1233 lbufsp%CRAK(jj(3)+j)=lbuf%CRAK(ii(3)+i)
1234 END IF
1235C-----------
1236 IF(elbuf_tab(ng)%BUFLY(1)%L_EPSA > 0)THEN
1237 lbufsp%EPSA(jj(1)+j)=lbuf%EPSA(ii(1)+i)
1238 lbufsp%EPSA(jj(2)+j)=lbuf%EPSA(ii(2)+i)
1239 lbufsp%EPSA(jj(3)+j)=lbuf%EPSA(ii(3)+i)
1240 END IF
1241C-----------
1242 IF(elbuf_tab(ng)%BUFLY(1)%L_SIGA > 0)THEN
1243 lbufsp%SIGA(jj(1)+j)=lbuf%SIGA(ii(1)+i)
1244 lbufsp%SIGA(jj(2)+j)=lbuf%SIGA(ii(2)+i)
1245 lbufsp%SIGA(jj(3)+j)=lbuf%SIGA(ii(3)+i)
1246 END IF
1247C-----------
1248C Stress in orthotropic skew for mulaw
1249 IF(elbuf_tab(ng)%BUFLY(1)%L_SIGL > 0)THEN
1250 lbufsp%SIGL(jj(1)+j)=lbuf%SIGL(ii(1)+i)
1251 lbufsp%SIGL(jj(2)+j)=lbuf%SIGL(ii(2)+i)
1252 lbufsp%SIGL(jj(3)+j)=lbuf%SIGL(ii(3)+i)
1253 lbufsp%SIGL(jj(4)+j)=lbuf%SIGL(ii(4)+i)
1254 lbufsp%SIGL(jj(5)+j)=lbuf%SIGL(ii(5)+i)
1255 lbufsp%SIGL(jj(6)+j)=lbuf%SIGL(ii(6)+i)
1256 END IF
1257C-----------
1258C Volume fraction
1259 IF(elbuf_tab(ng)%BUFLY(1)%L_BFRAC > 0)THEN
1260 lbufsp%BFRAC(jj(1)+j)=lbuf%BFRAC(ii(1)+i)
1261 END IF
1262C-----------
1263C Afterburning param
1264 IF(elbuf_tab(ng)%BUFLY(1)%L_ABURN > 0)THEN
1265 lbufsp%ABURN(jj(1)+j)=lbuf%ABURN(ii(1)+i)
1266 END IF
1267C-----------
1268C User variables (same for SPH and solids)
1269 IF(elbuf_tab(ng)%BUFLY(1)%NVAR_MAT > 0)THEN
1270 DO k=1,elbuf_tab(ng)%BUFLY(1)%NVAR_MAT
1271 mbufsp%VAR(nelsp*(k-1)+j) = mbuf%VAR(nel*(k-1)+i)
1272 END DO
1273 ENDIF
1274C
1275 ENDIF
1276C-----------
1277 ENDDO
1278 ENDDO
1279 END IF
1280 IF (iddw>0) CALL stoptimeg(ng)
1281C--------
1282 300 CONTINUE
1283 END DO
1284!$OMP END DO
1285C-----------------------------------------------
1286 RETURN
1287 END SUBROUTINE soltosphp
soltosphp
subroutine soltosphp(x, spbuf, ixs, kxsp, ipartsp, irst, elbuf_tab, iparg, ngrounc, igrounc, sol2sph, wa, pm)
Definition soltosph.F:527
soltosphf
subroutine soltosphf(a, spbuf, ixs, kxsp, ipartsp, nod2sp, irst, ngrounc, igrounc, iparg, stifn, sol2sph, sph2sol, elbuf_tab, itask, nodft, nodlt, isky, fskyi, igeo, sol2sph_typ)
Definition soltosph.F:45
startimeg
subroutine startimeg(ng)
Definition timer.F:1371
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1419
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
initbuf_mod
Definition initbuf.F:160
initbuf_mod::initbuf
subroutine initbuf(iparg, ng, mtn, llt, nft, iad, ity, npt, jale, ismstr, jeul, jtur, jthe, jlag, jmult, jhbe, jivf, mid, jpor, jcvt, jclose, jpla, irep, iint, igtyp, israt, isrot, icsen, isorth, isorthg, ifailure, jsms)
Definition initbuf.F:261
message_mod
Definition message_mod.F:1257
soltosph_mod
Definition soltosph_mod.F:32
sig_heph1
subroutine sig_heph1(jr0, js0, jt0, gsig, fhour, sig_heph, pm, ixs, ii, nel, lft, llt)
Definition soltosph_hour.F:34
srep2glo
subroutine srep2glo(x, ixs, gama, rx, ry, rz, sx, sy, sz, tx, ty, tz, r11, r12, r13, r21, r22, r23, r31, r32, r33, t11, t12, t13, t21, t22, t23, t31, t32, t33, jr0, js0, jt0, nel, lft, llt, jhbe, jcvt, isorth)
Definition srep2glo.F:45
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895
arret
subroutine arret(nn)
Definition arret.F:86
my_barrier
subroutine my_barrier
Definition machine.F:31