OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
s4fint_reg.F File Reference
#include "implicit_f.inc"
#include "parit_c.inc"
#include "scr02_c.inc"
#include "scr18_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine s4fint_reg (nloc_dmg, var_reg, nel, off, vol, nc1, nc2, nc3, nc4, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, imat, itask, dt2t, vol0, nft)

Function/Subroutine Documentation

◆ s4fint_reg()

subroutine s4fint_reg ( type(nlocal_str_), target nloc_dmg,
intent(in) var_reg,
integer nel,
intent(in) off,
intent(in) vol,
integer, dimension(nel) nc1,
integer, dimension(nel) nc2,
integer, dimension(nel) nc3,
integer, dimension(nel) nc4,
intent(in) px1,
intent(in) px2,
intent(in) px3,
intent(in) px4,
intent(in) py1,
intent(in) py2,
intent(in) py3,
intent(in) py4,
intent(in) pz1,
intent(in) pz2,
intent(in) pz3,
intent(in) pz4,
integer imat,
integer itask,
intent(inout) dt2t,
intent(in) vol0,
integer, intent(in) nft )

Definition at line 30 of file s4fint_reg.F.

38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE nlocal_reg_mod
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "parit_c.inc"
50#include "scr02_c.inc"
51#include "scr18_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER, INTENT(IN) :: NFT
56 INTEGER :: NEL, IMAT, ITASK
57 INTEGER, DIMENSION(NEL) :: NC1,NC2,NC3,NC4
58 my_real, INTENT(INOUT) ::
59 . dt2t
60 my_real, DIMENSION(NEL), INTENT(IN) ::
61 . vol,off,var_reg,vol0,
62 . px1,px2,px3,px4,py1,py2,py3,py4,pz1,pz2,pz3,pz4
63 TYPE(NLOCAL_STR_), TARGET :: NLOC_DMG
64C-----------------------------------------------
65C L o c a l V a r i a b l e s
66C-----------------------------------------------
67 INTEGER I,II,K,NNOD,N1,N2,N3,N4,L_NLOC
68 my_real
69 . l2,ntn,ntn_unl,ntn_vnl,xi,ntvar,a,
70 . b1,b2,b3,b4,zeta,sspnl,dtnl,le_max,maxstif
71 my_real, DIMENSION(NEL) ::
72 . f1,f2,f3,f4,lc
73 my_real, DIMENSION(:) ,ALLOCATABLE ::
74 . btb11,btb12,btb13,btb14,btb22,btb23,btb24,
75 . btb33,btb34,btb44,sti1,sti2,sti3,sti4
76 INTEGER, DIMENSION(:), ALLOCATABLE ::
77 . POS1,POS2,POS3,POS4
78 my_real, POINTER, DIMENSION(:) ::
79 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
80c-----------------------------------------------------------------------
81c VAR_REG : variable a regulariser (local cumulated plastic strain)
82c NTVAR = NT * VAR_REG
83C=======================================================================
84 l2 = nloc_dmg%LEN(imat)**2
85 xi = nloc_dmg%DAMP(imat)
86 nnod = nloc_dmg%NNOD
87 l_nloc = nloc_dmg%L_NLOC
88 zeta = nloc_dmg%DENS(imat)
89 sspnl = nloc_dmg%SSPNL(imat)
90 le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
91 ntn = four*four
92 lc(1:nel) = zero
93 ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
94 . btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),pos1(nel),
95 . pos2(nel),pos3(nel),pos4(nel))
96 ! For nodal timestep
97 IF (nodadt > 0) THEN
98 ! Non-local nodal stifness
99 ALLOCATE(sti1(nel),sti2(nel),sti3(nel),sti4(nel))
100 ! Non-local mass
101 mass => nloc_dmg%MASS(1:l_nloc)
102 ! Initial non-local mass
103 mass0 => nloc_dmg%MASS0(1:l_nloc)
104 ENDIF
105 vnl => nloc_dmg%VNL(1:l_nloc)
106 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
107 unl => nloc_dmg%UNL(1:l_nloc)
108c
109 !-----------------------------------------------------------------------
110 ! Computation of the element volume and the BtB matrix product
111 !-----------------------------------------------------------------------
112 ! Loop over elements
113# include "vectorize.inc"
114 DO i=1,nel
115c
116 ! Recovering the nodes of the tetra element
117 n1 = nloc_dmg%IDXI(nc1(i))
118 n2 = nloc_dmg%IDXI(nc2(i))
119 n3 = nloc_dmg%IDXI(nc3(i))
120 n4 = nloc_dmg%IDXI(nc4(i))
121c
122 ! Recovering the positions of the first d.o.fs of each nodes
123 pos1(i) = nloc_dmg%POSI(n1)
124 pos2(i) = nloc_dmg%POSI(n2)
125 pos3(i) = nloc_dmg%POSI(n3)
126 pos4(i) = nloc_dmg%POSI(n4)
127c
128 ! Computation of the product BtxB
129 btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
130 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
131 btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
132 btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
133 btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
134 btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
135 btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
136 btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
137 btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
138 btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
139c
140 ENDDO
141c
142 !-----------------------------------------------------------------------
143 ! Computation of non-local forces
144 !-----------------------------------------------------------------------
145 ! Loop over elements
146# include "vectorize.inc"
147 DO i = 1, nel
148c
149 ! If the element is not broken, normal computation
150 IF (off(i)/=zero) THEN
151c
152 ! Computing the product NtN*UNL
153 ntn_unl = (unl(pos1(i)) + unl(pos2(i)) + unl(pos3(i)) + unl(pos4(i))) / ntn
154c
155 ! Computing the product XDAMP*NtN*VNL
156 ntn_vnl = (vnl(pos1(i)) + vnl(pos2(i)) + vnl(pos3(i)) + vnl(pos4(i))) / ntn
157 IF (nodadt > 0) THEN
158 ntn_vnl = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
159 . sqrt(mass(pos2(i))/mass0(pos2(i))),
160 . sqrt(mass(pos3(i))/mass0(pos3(i))),
161 . sqrt(mass(pos4(i))/mass0(pos4(i))))*ntn_vnl
162 ENDIF
163c
164 ! Computation of the product LEN**2 * BtxB
165 b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i))
166 . + btb13(i)*unl(pos3(i)) + btb14(i)*unl(pos4(i)))
167c
168 b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i))
169 . + btb23(i)*unl(pos3(i)) + btb24(i)*unl(pos4(i)))
170c
171 b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i))
172 . + btb33(i)*unl(pos3(i)) + btb34(i)*unl(pos4(i)))
173c
174 b4 = l2 * vol(i) * ( btb14(i)*unl(pos1(i)) + btb24(i)*unl(pos2(i))
175 . + btb34(i)*unl(pos3(i)) + btb44(i)*unl(pos4(i)))
176c
177 ! Multiplication by the volume of the element
178 ntn_unl = ntn_unl * vol(i)
179 ntn_vnl = ntn_vnl * xi * vol(i)
180c
181 ! Introducing the internal variable to be regularized
182 ntvar = var_reg(i)*fourth* vol(i)
183c
184 ! Computing the elementary non-local forces
185 a = ntn_unl + ntn_vnl - ntvar
186 f1(i) = a + b1
187 f2(i) = a + b2
188 f3(i) = a + b3
189 f4(i) = a + b4
190c
191 ! Computing nodal equivalent stiffness
192 IF (nodadt > 0) THEN
193 sti1(i) = (abs(l2*btb11(i) + one/ntn) + abs(l2*btb12(i) + one/ntn) + abs(l2*btb13(i) + one/ntn) +
194 . abs(l2*btb14(i) + one/ntn))* vol(i)
195 sti2(i) = (abs(l2*btb12(i) + one/ntn) + abs(l2*btb22(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) +
196 . abs(l2*btb24(i) + one/ntn))* vol(i)
197 sti3(i) = (abs(l2*btb13(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) + abs(l2*btb33(i) + one/ntn) +
198 . abs(l2*btb34(i) + one/ntn))* vol(i)
199 sti4(i) = (abs(l2*btb14(i) + one/ntn) + abs(l2*btb24(i) + one/ntn) + abs(l2*btb34(i) + one/ntn) +
200 . abs(l2*btb44(i) + one/ntn))* vol(i)
201 ENDIF
202c
203 ! If the element is broken, computation of absorbing forces
204 ELSE
205c
206 ! Initial element characteristic length
207 lc(i) = vol0(i)**third
208c
209 IF (nodadt > 0) THEN
210 ! Non-local absorbing forces
211 f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*zeta*sspnl*half*
212 . (vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
213 f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*zeta*sspnl*half*
214 . (vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
215 f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*zeta*sspnl*half*
216 . (vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
217 f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*zeta*sspnl*half*
218 . (vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)
219 ! Computing nodal equivalent stiffness
220 sti1(i) = em20
221 sti2(i) = em20
222 sti3(i) = em20
223 sti4(i) = em20
224 ELSE
225 ! Non-local absorbing forces
226 f1(i) = zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
227 f2(i) = zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
228 f3(i) = zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
229 f4(i) = zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)
230 ENDIF
231 ENDIF
232 ENDDO
233c-----------------------------------------------------------------------
234c Assemblage
235c-----------------------------------------------------------------------
236 ! If PARITH/OFF
237 IF (iparit == 0) THEN
238 fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
239 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
240 DO i=1,nel
241 ! Assembling the forces in the classic way
242 fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
243 fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
244 fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
245 fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
246 IF (nodadt > 0) THEN
247 ! Spectral radius of stiffness matrix
248 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
249 ! Computing nodal stiffness
250 stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
251 stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
252 stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
253 stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
254 ENDIF
255 ENDDO
256c
257 ! If PARITH/ON
258 ELSE
259c
260 DO i=1,nel
261 ii = i + nft
262c
263 ! Spectral radius of stiffness matrix
264 IF (nodadt > 0) THEN
265 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
266 ENDIF
267c
268 k = nloc_dmg%IADS(1,ii)
269 nloc_dmg%FSKY(k,1) = -f1(i)
270 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
271c
272 k = nloc_dmg%IADS(3,ii)
273 nloc_dmg%FSKY(k,1) = -f2(i)
274 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
275c
276 k = nloc_dmg%IADS(6,ii)
277 nloc_dmg%FSKY(k,1) = -f3(i)
278 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
279c
280 k = nloc_dmg%IADS(5,ii)
281 nloc_dmg%FSKY(k,1) = -f4(i)
282 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
283c
284 ENDDO
285 ENDIF
286c
287 !-----------------------------------------------------------------------
288 ! Computing non-local timestep
289 !-----------------------------------------------------------------------
290 IF (nodadt == 0) THEN
291 DO i = 1,nel
292 ! If the element is not broken, normal computation
293 IF (off(i)/=zero) THEN
294 ! Non-local critical time-step in the plane
295 dtnl = (two*(min(vol(i)**third,le_max))*sqrt(three*zeta))/
296 . sqrt(twelve*l2 + (min(vol(i)**third,le_max))**2)
297 ! Retaining the minimal value
298 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl)
299 ENDIF
300 ENDDO
301 ENDIF
302c
303c-----------
304 ! Deallocation of tables
305 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
306 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
307 IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
308 IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
309 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
310 IF (ALLOCATED(btb23)) DEALLOCATE(btb23)
311 IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
312 IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
313 IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
314 IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
315 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
316 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
317 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
318 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
319 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
320 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
321 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
322 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
323c
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
nlocal_reg_mod
Definition nlocal_reg_mod.F:119