OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
rgwalc.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!|| rgwalc ../engine/source/constraints/general/rwall/rgwalc.F
25!||--- called by ------------------------------------------------------
26!|| rgwal0 ../engine/source/constraints/general/rwall/rgwal0.F
27!|| rgwal0_imp ../engine/source/constraints/general/rwall/rgwal0.F
28!||--- calls -----------------------------------------------------
29!|| my_barrier ../engine/source/system/machine.F
30!|| sum_6_float ../engine/source/system/parit.F
31!||====================================================================
32 SUBROUTINE rgwalc(X ,A ,V ,RWL ,NSW ,
33 . NSN ,ITIED,MSR ,MS ,WEIGHT,
34 . ICONT,FRWL6,IMP_S,NT_RW,IDDL ,
35 . IKC ,NDOF ,NODNX_SMS, WEIGHT_MD,
36 . WFEXT, WFEXT_MD)
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41#include "comlock.inc"
42C-----------------------------------------------
43C C o m m o n B l o c k s
44C-----------------------------------------------
45#include "com06_c.inc"
46#include "com08_c.inc"
47#include "scr11_c.inc"
48#include "impl1_c.inc"
49#include "sms_c.inc"
50C-----------------------------------------------
51C D u m m y A r g u m e n t s
52C-----------------------------------------------
53 INTEGER NSN, ITIED, MSR,ICONT,IMP_S,NT_RW
54 INTEGER NSW(*),WEIGHT(*), IDDL(*),IKC(*),NDOF(*), NODNX_SMS(*), WEIGHT_MD(*)
55 my_real X(*), A(*), V(*), RWL(*), MS(*)
56 DOUBLE PRECISION FRWL6(7,6)
57 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT, WFEXT_MD
58C-----------------------------------------------
59C L o c a l V a r i a b l e s
60C-----------------------------------------------
61 INTEGER M3, M2, M1, I, N, N3, N2, N1,J,JJ,K,NINDEX,
62 . INDEX(NSN)
63C REAL
64 my_real
65 . ra2, xwl, ywl, zwl, vxw, vyw, vzw, wfxtn, fact,
66 . wfxt, vx, vy, vz, ux, uy, uz, xc, yc, zc, dd1, dd,
67 . dp, xt, yt, zt, xx, xn, yn, zn, dv, da, dvt, fnxn, fnyn, fnzn,
68 . fnxt, fnyt, fnzt, fndfn, ftdft, fric, fric2, fcoe,
69 . xwlo, ywlo, zwlo,
70 . f1(nsn), f2(nsn), f3(nsn), f4(nsn), f5(nsn), f6(nsn), f7(nsn),
71 . wfxt2, wfxtn2, wewe2
72 DOUBLE PRECISION
73 . FRWL6_L(7,6)
74C-----------------------------------------------
75 ICONT=0
76 CALL my_barrier
77! We need an OMP barrier here because each thread initializes
78! the variables : without barrier ICONT can be set to 1 by a
79! thread in the !$OMP DO loop whereas another (late) thread initializes
80! ICONT to 0
81! ICONT is a shared variable
82 ra2=(half*rwl(7))**2
83 IF(msr==0)THEN
84 xwlo=rwl(4)
85 ywlo=rwl(5)
86 zwlo=rwl(6)
87 xwl=rwl(4)
88 ywl=rwl(5)
89 zwl=rwl(6)
90 vxw=zero
91 vyw=zero
92 vzw=zero
93 ELSE
94 m3=3*msr
95 m2=m3-1
96 m1=m2-1
97C changement formulation : plus d'impasse sur contribution force
98 vxw=v(m1)+a(m1)*dt12
99 vyw=v(m2)+a(m2)*dt12
100 vzw=v(m3)+a(m3)*dt12
101 xwlo=x(m1)
102 ywlo=x(m2)
103 zwlo=x(m3)
104 xwl=x(m1)+vxw*dt2
105 ywl=x(m2)+vyw*dt2
106 zwl=x(m3)+vzw*dt2
107 ENDIF
108 wfxt =zero
109 wfxtn=zero
110 wfxt2 =zero
111 wfxtn2=zero
112 nindex=0
113C----
114 IF(idtmins==0.AND.idtmins_int==0)THEN
115!$OMP DO
116 DO i=1,nsn
117 n=nsw(i)
118 n3=3*n
119 n2=n3-1
120 n1=n2-1
121 vx=v(n1)+a(n1)*dt12
122 vy=v(n2)+a(n2)*dt12
123 vz=v(n3)+a(n3)*dt12
124 ux=x(n1)+vx*dt2
125 uy=x(n2)+vy*dt2
126 uz=x(n3)+vz*dt2
127 xc=ux-xwl
128 yc=uy-ywl
129 zc=uz-zwl
130 dd1=xc**2+yc**2+zc**2
131 dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
132 dp=dd1-dd**2
133 IF(dp <= ra2)THEN
134 icont=1
135 nindex = nindex+1
136 index(nindex) = i
137 END IF
138 END DO
139!$OMP END DO
140 ELSE
141!$OMP DO
142 DO i=1,nsn
143 n=nsw(i)
144 IF(nodnx_sms(n)/=0)cycle
145 n3=3*n
146 n2=n3-1
147 n1=n2-1
148 vx=v(n1)+a(n1)*dt12
149 vy=v(n2)+a(n2)*dt12
150 vz=v(n3)+a(n3)*dt12
151 ux=x(n1)+vx*dt2
152 uy=x(n2)+vy*dt2
153 uz=x(n3)+vz*dt2
154 xc=ux-xwl
155 yc=uy-ywl
156 zc=uz-zwl
157 dd1=xc**2+yc**2+zc**2
158 dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
159 dp=dd1-dd**2
160 IF(dp <= ra2)THEN
161 icont=1
162 nindex = nindex+1
163 index(nindex) = i
164 END IF
165 END DO
166!$OMP END DO
167 END IF
168C
169 fact=one/dt12
170 DO j = 1,nindex
171 i = index(j)
172 n=nsw(i)
173 n3=3*n
174 n2=n3-1
175 n1=n2-1
176 xc=x(n1)-xwlo
177 yc=x(n2)-ywlo
178 zc=x(n3)-zwlo
179 dd1=xc**2+yc**2+zc**2
180 dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
181 dp=dd1-dd**2
182 xt=dd*rwl(1)
183 yt=dd*rwl(2)
184 zt=dd*rwl(3)
185 xx=sqrt(dp)
186 xn=(xc-xt)/xx
187 yn=(yc-yt)/xx
188 zn=(zc-zt)/xx
189 dv=(v(n1)-vxw)*xn+(v(n2)-vyw)*yn+(v(n3)-vzw)*zn
190 da=a(n1)*xn+a(n2)*yn+a(n3)*zn
191 dvt=dv+da*dt12
192 fnxn=dvt*xn*ms(n)
193 fnyn=dvt*yn*ms(n)
194 fnzn=dvt*zn*ms(n)
195 wfxtn = wfxtn - weight_md(n)*fact*
196 . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)
197 wewe2 = (1-weight_md(n))*weight(n)
198 wfxtn2 = wfxtn2 - wewe2*fact*
199 . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)
200 f1(j) = fnxn*weight_md(n)
201 f2(j) = fnyn*weight_md(n)
202 f3(j) = fnzn*weight_md(n)
203 f4(j) = ms(n)*weight_md(n)
204 IF(itied/=0)THEN
205 fnxt=((v(n1)-vxw)+a(n1)*dt12)*ms(n)-fnxn
206 fnyt=((v(n2)-vyw)+a(n2)*dt12)*ms(n)-fnyn
207 fnzt=((v(n3)-vzw)+a(n3)*dt12)*ms(n)-fnzn
208 fndfn=fnxn**2+fnyn**2+fnzn**2
209 ftdft=fnxt**2+fnyt**2+fnzt**2
210 fric=rwl(13)
211 fric2=fric**2
212 IF(ftdft<=fric2*fndfn.OR.itied==1) THEN
213C POINT SECND TIED
214 a(n1)=zero
215 a(n2)=zero
216 a(n3)=zero
217 v(n1)=vxw
218 v(n2)=vyw
219 v(n3)=vzw
220 IF (imp_s==1)THEN
221 IF(ndof(n)>0) THEN
222
223 jj=iddl(n)+1
224 IF (ikc(jj)==0)ikc(jj)=3
225 IF (ikc(jj+1)==0)ikc(jj+1)=3
226 IF (ikc(jj+2)==0)ikc(jj+2)=3
227 ENDIF
228 ENDIF
229 ELSE
230C POINT SECND SLIDING
231 fcoe=fric*sqrt(fndfn/ftdft)
232 fnxt=fcoe*fnxt
233 fnyt=fcoe*fnyt
234 fnzt=fcoe*fnzt
235 a(n1)=a(n1)-(da*xn+fnxt/(dt12*ms(n)))
236 a(n2)=a(n2)-(da*yn+fnyt/(dt12*ms(n)))
237 a(n3)=a(n3)-(da*zn+fnzt/(dt12*ms(n)))
238 v(n1)=v(n1)-dv*xn
239 v(n2)=v(n2)-dv*yn
240 v(n3)=v(n3)-dv*zn
241 IF (imp_s==1) THEN
242 IF(ndof(n)>0) THEN
243 v(n1)=-dv
244 a(n1)=xn
245 a(n2)=yn
246 a(n3)=zn
247 jj=iddl(n)+1
248 IF (ikc(jj)==0)ikc(jj)=10
249 ENDIF
250 ENDIF
251 wfxt=wfxt-
252 . ((v(n1)-vxw)*fnxt+(v(n2)-vyw)*fnyt+(v(n3)-vzw)*fnzt)
253 . *weight_md(n)
254 wewe2 = (1-weight_md(n))*weight(n)
255 wfxt2=wfxt2-
256 . ((v(n1)-vxw)*fnxt+(v(n2)-vyw)*fnyt+(v(n3)-vzw)*fnzt)
257 . *wewe2
258 ENDIF
259 f5(j) = fnxt*weight_md(n)
260 f6(j) = fnyt*weight_md(n)
261 f7(j) = fnzt*weight_md(n)
262 ELSE
263 f5(j) = zero
264 f6(j) = zero
265 f7(j) = zero
266 a(n1)=a(n1)-da*xn
267 a(n2)=a(n2)-da*yn
268 a(n3)=a(n3)-da*zn
269 v(n1)=v(n1)-dv*xn
270 v(n2)=v(n2)-dv*yn
271 v(n3)=v(n3)-dv*zn
272 IF (imp_s==1) THEN
273 IF(ndof(n)>0) THEN
274 v(n1)=-dv
275 a(n1)=xn
276 a(n2)=yn
277 a(n3)=zn
278 jj=iddl(n)+1
279 IF (ikc(jj)==0)ikc(jj)=10
280 ENDIF
281 ENDIF
282 ENDIF
283 ENDDO
284C
285 IF(imconv==1)THEN
286 wfxt=wfxt+half*dt1*wfxtn
287 wfxt2=wfxt2+half*dt1*wfxtn2
288!$OMP ATOMIC
289 wfext=wfext+wfxt
290 wfext_md=wfext_md+wfxt2
291 DO k = 1, 6
292 frwl6_l(1,k) = zero
293 frwl6_l(2,k) = zero
294 frwl6_l(3,k) = zero
295 frwl6_l(4,k) = zero
296 frwl6_l(5,k) = zero
297 frwl6_l(6,k) = zero
298 frwl6_l(7,k) = zero
299 END DO
300 CALL sum_6_float(1, nindex, f1, frwl6_l(1,1), 7)
301 CALL sum_6_float(1, nindex, f2, frwl6_l(2,1), 7)
302 CALL sum_6_float(1, nindex, f3, frwl6_l(3,1), 7)
303 CALL sum_6_float(1, nindex, f4, frwl6_l(4,1), 7)
304 CALL sum_6_float(1, nindex, f5, frwl6_l(5,1), 7)
305 CALL sum_6_float(1, nindex, f6, frwl6_l(6,1), 7)
306 CALL sum_6_float(1, nindex, f7, frwl6_l(7,1), 7)
307
308#include "lockon.inc"
309 DO k = 1, 6
310 frwl6(1,k) = frwl6(1,k)+frwl6_l(1,k)
311 frwl6(2,k) = frwl6(2,k)+frwl6_l(2,k)
312 frwl6(3,k) = frwl6(3,k)+frwl6_l(3,k)
313 frwl6(4,k) = frwl6(4,k)+frwl6_l(4,k)
314 frwl6(5,k) = frwl6(5,k)+frwl6_l(5,k)
315 frwl6(6,k) = frwl6(6,k)+frwl6_l(6,k)
316 frwl6(7,k) = frwl6(7,k)+frwl6_l(7,k)
317 END DO
318#include "lockoff.inc"
319 ENDIF
320
321 IF (imp_s==1) THEN
322 DO j=1,nindex
323 i = index(j)
324 n=nsw(i)
325 IF (ndof(n)>0) THEN
326C to be uncommented the day it is parallel in implicit
327!$OMP ATOMIC
328 nt_rw = nt_rw + 1
329 END IF
330 ENDDO
331 ENDIF
332C
333 RETURN
334 END
sum_6_float
subroutine sum_6_float(jft, jlt, f, f6, n)
Definition parit.F:64
rgwalc
subroutine rgwalc(x, a, v, rwl, nsw, nsn, itied, msr, ms, weight, icont, frwl6, imp_s, nt_rw, iddl, ikc, ndof, nodnx_sms, weight_md, wfext, wfext_md)
Definition rgwalc.F:37
my_barrier
subroutine my_barrier
Definition machine.F:31