OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sponfprs.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!|| sponfprs ../engine/source/elements/sph/sponfprs.F
25!||--- called by ------------------------------------------------------
26!|| forintp ../engine/source/elements/forintp.F
27!||--- calls -----------------------------------------------------
28!|| initbuf ../engine/share/resol/initbuf.F
29!||--- uses -----------------------------------------------------
30!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
31!|| initbuf_mod ../engine/share/resol/initbuf.F
32!|| sphbox ../engine/share/modules/sphbox.F
33!||====================================================================
34 SUBROUTINE sponfprs(X ,V ,A ,MS ,
35 2 SPBUF ,ITAB ,KXSP ,IXSP ,NOD2SP ,
36 3 ISPHIO ,VSPHIO ,NPC ,PLD ,PM ,
37 4 IPARG ,ELBUF_TAB,IPART ,IPARTSP ,WASPACT ,
38 5 VNORMAL ,WA ,SPHVELN ,WAR ,WFEXT)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE initbuf_mod
43 USE sphbox
44 USE elbufdef_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49#include "comlock.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "vect01_c.inc"
54#include "com01_c.inc"
55#include "com06_c.inc"
56#include "com08_c.inc"
57#include "param_c.inc"
58#include "sphcom.inc"
59#include "scr17_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),
64 . ISPHIO(NISPHIO,*),IPART(LIPART1,*),IPARTSP(*),WASPACT(*),
65 . NPC(*),IPARG(NPARG,*)
66 my_real
67 . x(3,*) ,v(3,*) ,a(3,*) ,ms(*) ,spbuf(nspbuf,*) ,vsphio(*) ,
68 . pld(*) ,pm(npropm,*),vnormal(3,*),wa(kwasph,*), sphveln(2,*),
69 . war(10,*)
70 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
71 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER I, ITYPE, IVAD,
76 . NEL,KAD,NG,K,
77 . II,IPT,JJ,NPF,
78 . NS,N,INOD,
79 . iprt,ippv,j,m,jnod,impose,jmpose,
80 . imat,ifpres,np,jmpose2,nn,kk(6)
81 my_real
82 . pentp,pn,pshft,px,
83 . xi,yi,zi,xj,yj,zj,dmin,dd,vx,vy,vz,
84 . pmin,vn,vo,po,rhoi,sspi,pinfini,lc,alp,
85 . wfextt,dvol
86 TYPE(g_bufel_) ,POINTER :: GBUF
87 LOGICAL :: lBOOL
88C=======================================================================
89 WFEXTT= zero
90 pinfini = zero
91 pn = zero
92C-----------------------------------------------
93C 1. GENERAL OUTLET : REIMPOSE P=f(t) ou Pvoisin.
94C-----------------------------------------------
95 DO ns=1,nsphact
96 n=waspact(ns)
97 impose=kxsp(2,n)/(ngroup+1)
98 IF(impose/=0)THEN
99 itype=isphio(1,impose)
100 IF(itype==2)THEN
101C------
102C general outlet.
103 iprt=isphio(2,impose)
104 ivad=isphio(4,impose)
105C
106 ifpres=isphio(6,impose)
107 IF(ifpres/=0)THEN
108C------
109C P=f(t)
110 npf = (npc(ifpres+1)-npc(ifpres))/2
111 ii = npc(ifpres)
112 IF (tt<=pld(ii)) THEN
113 pentp=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
114 pn =pld(ii+1)+pentp*(tt-pld(ii))
115 ELSEIF (tt>=pld(ii+2*(npf-1))) THEN
116 jj=ii+2*(npf-1)
117 pentp=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
118 pn =pld(jj+1)+max(-pld(jj+1),pentp*(tt-pld(jj)))
119 ELSE
120 DO ipt=1,npf-1
121 IF (pld(ii)<=tt.AND.tt<=pld(ii+2)) THEN
122 pentp=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
123 pn =pld(ii+1)+pentp*(tt-pld(ii))
124 GOTO 140
125 ENDIF
126 ii=ii+2
127 ENDDO
128 ENDIF
129 140 CONTINUE
130 pn=pn*vsphio(ivad+1)
131 iprt =ipartsp(n)
132 imat =ipart(1,iprt)
133 pmin= pm(37,imat)
134 pn=max(pn,pmin)
135 ELSE
136C------
137C continuite.
138 inod=kxsp(3,n)
139 xi=x(1,inod)
140 yi=x(2,inod)
141 zi=x(3,inod)
142C
143 ippv=0
144 dmin=1.e+20
145 DO j=1,kxsp(4,n)
146 jnod=ixsp(j,n)
147 IF(jnod>0)THEN
148 m=nod2sp(jnod)
149 IF(kxsp(2,m)>=0)THEN
150 jmpose=kxsp(2,m)/(ngroup+1)
151 lbool=.false.
152 IF(jmpose==0)THEN
153 lbool=.true.
154 ELSE
155 IF(isphio(1,jmpose)==1)lbool=.true.
156 ENDIF
157 IF(lbool)THEN
158 xj =x(1,jnod)
159 yj =x(2,jnod)
160 zj =x(3,jnod)
161 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
162 IF(dd<dmin)THEN
163 ippv=jnod
164 dmin=dd
165 ENDIF
166 ENDIF
167 ENDIF
168 ELSE
169 nn = -jnod
170 jmpose = nint(xsphr(12,nn))
171 IF(jmpose>0)THEN
172 jmpose2=isphio(1,jmpose)
173 ELSE
174 jmpose2=0
175 ENDIF
176 IF(jmpose2==0.OR.jmpose2==1)THEN
177 xj =xsphr(3,nn)
178 yj =xsphr(4,nn)
179 zj =xsphr(5,nn)
180 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
181 IF(dd<dmin)THEN
182 ippv=jnod
183 dmin=dd
184 ENDIF
185 ENDIF
186 ENDIF
187 ENDDO
188C
189 IF(ippv>0)THEN
190 np=nod2sp(ippv)
191 ng=mod(kxsp(2,np),ngroup+1)
192 CALL initbuf(iparg ,ng ,
193 2 mtn ,nel ,nft ,kad ,ity ,
194 3 npt ,jale ,ismstr ,jeul ,jtur ,
195 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
196 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
197 6 irep ,iint ,igtyp ,israt ,isrot ,
198 7 icsen ,isorth ,isorthg ,ifailure,jsms )
199!
200 DO i=1,6
201 kk(i) = nel*(i-1)
202 ENDDO
203!
204 gbuf => elbuf_tab(ng)%GBUF
205 k=np-nft
206!
207 pn=-( gbuf%SIG(kk(1)+k)
208 . +gbuf%SIG(kk(2)+k)
209 . +gbuf%SIG(kk(3)+k))*third
210
211 ELSEIF(ippv<0)THEN !cas IPPV negatif on se sert des infos recuperes ds la routine de com
212
213 pn=-( war(1,-ippv)
214 . +war(2,-ippv)
215 . +war(3,-ippv))*third
216 ENDIF
217 ENDIF
218C
219 ng=mod(kxsp(2,n),ngroup+1)
220 CALL initbuf(iparg ,ng ,
221 2 mtn ,nel ,nft ,kad ,ity ,
222 3 npt ,jale ,ismstr ,jeul ,jtur ,
223 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
224 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
225 6 irep ,iint ,igtyp ,israt ,isrot ,
226 7 icsen ,isorth ,isorthg ,ifailure,jsms )
227!
228 DO i=1,6
229 kk(i) = nel*(i-1)
230 ENDDO
231!
232 gbuf => elbuf_tab(ng)%GBUF
233 k=n-nft
234!
235 px=-( gbuf%SIG(kk(1)+k)
236 . +gbuf%SIG(kk(2)+k)
237 . +gbuf%SIG(kk(3)+k))*third
238 gbuf%SIG(kk(1)+k)=gbuf%SIG(kk(1)+k)+px-pn
239 gbuf%SIG(kk(2)+k)=gbuf%SIG(kk(2)+k)+px-pn
240 gbuf%SIG(kk(3)+k)=gbuf%SIG(kk(3)+k)+px-pn
241 wa(1,n)=gbuf%SIG(kk(1)+k)
242 wa(2,n)=gbuf%SIG(kk(2)+k)
243 wa(3,n)=gbuf%SIG(kk(3)+k)
244C
245 inod=kxsp(3,n)
246 dvol=spbuf(12,n)/max(em20,spbuf(2,n))
247 . -spbuf(12,n)/max(em20,wa(10,n))
248 wfextt=wfextt+(pn-px)*dvol
249 ENDIF
250 ENDIF
251 ENDDO
252C-----------------------------------------------
253C 2. SILENT BOUNDARY.
254C-----------------------------------------------
255 DO ns=1,nsphact
256 n=waspact(ns)
257 impose=kxsp(2,n)/(ngroup+1)
258 IF(impose/=0)THEN
259 itype=isphio(1,impose)
260 IF(itype==3)THEN
261 ivad=isphio(4,impose)
262C-------
263 inod=kxsp(3,n)
264 xi=x(1,inod)
265 yi=x(2,inod)
266 zi=x(3,inod)
267C-------
268 ippv=0
269 dmin=1.e+20
270 DO j=1,kxsp(4,n)
271 jnod=ixsp(j,n)
272 IF(jnod>0)THEN
273 m=nod2sp(jnod)
274 IF(kxsp(2,m)>=0)THEN
275 jmpose=kxsp(2,m)/(ngroup+1)
276 lbool=.false.
277 IF(jmpose == 0)THEN
278 lbool=.true.
279 ELSE
280 IF(isphio(1,jmpose) == 1)lbool=.true.
281 ENDIF
282 IF(lbool)THEN
283 xj =x(1,jnod)
284 yj =x(2,jnod)
285 zj =x(3,jnod)
286 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
287 IF(dd<dmin)THEN
288 ippv=jnod
289 dmin=dd
290 ENDIF
291 ENDIF
292 ENDIF
293 ELSE
294 nn = -jnod
295 jmpose = nint(xsphr(12,nn))
296 IF(jmpose>0)THEN
297 jmpose2=isphio(1,jmpose)
298 ELSE
299 jmpose2=0
300 ENDIF
301 IF(jmpose2==0.OR.jmpose2==1)THEN
302 xj =xsphr(3,nn)
303 yj =xsphr(4,nn)
304 zj =xsphr(5,nn)
305 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
306 IF(dd<dmin)THEN
307 ippv=jnod
308 dmin=dd
309 ENDIF
310 ENDIF
311 ENDIF
312 ENDDO
313C-------
314 ifpres=isphio(6,impose)
315 IF(ifpres/=0)THEN
316 npf = (npc(ifpres+1)-npc(ifpres))/2
317 ii = npc(ifpres)
318 IF (tt<=pld(ii)) THEN
319 pentp=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
320 pinfini =pld(ii+1)+pentp*(tt-pld(ii))
321 ELSEIF (tt>=pld(ii+2*(npf-1))) THEN
322 jj=ii+2*(npf-1)
323 pentp=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
324 pinfini =pld(jj+1)+max(-pld(jj+1),pentp*(tt-pld(jj)))
325 ELSE
326 DO ipt=1,npf-1
327 IF (pld(ii)<=tt.AND.tt<=pld(ii+2)) THEN
328 pentp=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
329 pinfini =pld(ii+1)+pentp*(tt-pld(ii))
330 GOTO 240
331 ENDIF
332 ii=ii+2
333 ENDDO
334 240 CONTINUE
335 ENDIF
336 pinfini=pinfini*vsphio(ivad+1)
337 ELSE
338 pinfini=vsphio(ivad+1)
339 ENDIF
340C------
341 vx=v(1,inod)
342 vy=v(2,inod)
343 vz=v(3,inod)
344 vn=vnormal(1,n)*vx+vnormal(2,n)*vy+vnormal(3,n)*vz
345
346 IF(vn>=zero)THEN
347 rhoi=spbuf(2,n)
348 sspi=wa(8,n)
349C
350 iprt =ipartsp(n)
351 imat =ipart(1,iprt)
352 pmin= pm(37,imat)
353C
354 vo=sphveln(1,n)
355 po=sphveln(2,n)
356 pinfini=pinfini-half*rhoi*vn*vn
357 lc =vsphio(ivad+2)
358 alp =half*sspi/max(em30,lc)*dt1
359 pn=po+(one - alp)*rhoi*sspi*(vn-vo)+alp*(pinfini-po)
360 pn=max(pn,pmin)
361C
362 ng=mod(kxsp(2,n),ngroup+1)
363 CALL initbuf(iparg ,ng ,
364 2 mtn ,nel ,nft ,kad ,ity ,
365 3 npt ,jale ,ismstr ,jeul ,jtur ,
366 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
367 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
368 6 irep ,iint ,igtyp ,israt ,isrot ,
369 7 icsen ,isorth ,isorthg ,ifailure,jsms )
370!
371 DO i=1,6
372 kk(i) = nel*(i-1)
373 ENDDO
374!
375 gbuf => elbuf_tab(ng)%GBUF
376 k=n-nft
377!
378 px=-( gbuf%SIG(kk(1)+k)
379 . +gbuf%SIG(kk(2)+k)
380 . +gbuf%SIG(kk(3)+k))*third
381
382 gbuf%SIG(kk(1)+k)=gbuf%SIG(kk(1)+k)+px-pn
383 gbuf%SIG(kk(2)+k)=gbuf%SIG(kk(2)+k)+px-pn
384 gbuf%SIG(kk(3)+k)=gbuf%SIG(kk(3)+k)+px-pn
385 wa(1,n)=gbuf%SIG(kk(1)+k)
386 wa(2,n)=gbuf%SIG(kk(2)+k)
387 wa(3,n)=gbuf%SIG(kk(3)+k)
388C
389 dvol=spbuf(12,n)/max(em20,spbuf(2,n))
390 . -spbuf(12,n)/max(em20,wa(10,n))
391 wfextt=wfextt+(pn-px)*dvol
392 ELSE
393 ng=mod(kxsp(2,n),ngroup+1)
394 CALL initbuf(iparg ,ng ,
395 2 mtn ,nel ,nft ,kad ,ity ,
396 3 npt ,jale ,ismstr ,jeul ,jtur ,
397 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
398 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
399 6 irep ,iint ,igtyp ,israt ,isrot ,
400 7 icsen ,isorth ,isorthg ,ifailure,jsms )
401!
402 DO i=1,6
403 kk(i) = nel*(i-1)
404 ENDDO
405!
406 gbuf => elbuf_tab(ng)%GBUF
407 k=n-nft
408!
409 pn=-( gbuf%SIG(kk(1)+k)
410 . +gbuf%SIG(kk(2)+k)
411 . +gbuf%SIG(kk(3)+k))*third
412 ENDIF
413 sphveln(1,n)= vn
414 sphveln(2,n)= pn
415 ENDIF
416 ENDIF
417 ENDDO
418C-------------------------------------------
419!$OMP ATOMIC
420 wfext=wfext+wfextt
421C-------------------------------------------
422 RETURN
423 END
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
sponfprs
subroutine sponfprs(x, v, a, ms, spbuf, itab, kxsp, ixsp, nod2sp, isphio, vsphio, npc, pld, pm, iparg, elbuf_tab, ipart, ipartsp, waspact, vnormal, wa, sphveln, war, wfext)
Definition sponfprs.F:39