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 . dx, dy, dz, 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
80 ! Coefficient for non-local stability to take into account damping
81 my_real, PARAMETER :: cdamp = 0.7d0
82c-----------------------------------------------------------------------
83c VAR_REG : variable a regulariser (local cumulated plastic strain)
84c NTVAR = NT * VAR_REG
85C=======================================================================
86 l2 = nloc_dmg%LEN(imat)**2
87 xi = nloc_dmg%DAMP(imat)
88 nnod = nloc_dmg%NNOD
89 l_nloc = nloc_dmg%L_NLOC
90 zeta = nloc_dmg%DENS(imat)
91 sspnl = nloc_dmg%SSPNL(imat)
92 le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
93 ntn = four*four
94 lc(1:nel) = zero
95 ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
96 . btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),pos1(nel),
97 . pos2(nel),pos3(nel),pos4(nel))
98 ! For nodal timestep
99 IF (nodadt > 0) THEN
100 ! Non-local nodal stifness
101 ALLOCATE(sti1(nel),sti2(nel),sti3(nel),sti4(nel))
102 ! Non-local mass
103 mass => nloc_dmg%MASS(1:l_nloc)
104 ! initial non-local mass
105 mass0 => nloc_dmg%MASS0(1:l_nloc)
106 ENDIF
107 vnl => nloc_dmg%VNL(1:l_nloc)
108 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
109 unl => nloc_dmg%UNL(1:l_nloc)
110c
111 !-----------------------------------------------------------------------
112 ! Computation of the element volume and the BtB matrix product
113 !-----------------------------------------------------------------------
114 ! Loop over elements
115# include "vectorize.inc"
116 DO i=1,nel
117c
118 ! Recovering the nodes of the tetra element
119 n1 = nloc_dmg%IDXI(nc1(i))
120 n2 = nloc_dmg%IDXI(nc2(i))
121 n3 = nloc_dmg%IDXI(nc3(i))
122 n4 = nloc_dmg%IDXI(nc4(i))
123c
124 ! Recovering the positions of the first d.o.fs of each nodes
125 pos1(i) = nloc_dmg%POSI(n1)
126 pos2(i) = nloc_dmg%POSI(n2)
127 pos3(i) = nloc_dmg%POSI(n3)
128 pos4(i) = nloc_dmg%POSI(n4)
129c
130 ! Computation of the product BtxB
131 btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
132 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
133 btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
134 btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
135 btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
136 btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
137 btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
138 btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
139 btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
140 btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
141c
142 ENDDO
143c
144 !-----------------------------------------------------------------------
145 ! Computation of non-local forces
146 !-----------------------------------------------------------------------
147 ! Loop over elements
148# include "vectorize.inc"
149 DO i = 1, nel
150c
151 ! If the element is not broken, normal computation
152 IF (off(i)/=zero) THEN
153c
154 ! Computing the product NtN*UNL
155 ntn_unl = (unl(pos1(i)) + unl(pos2(i)) + unl(pos3(i)) + unl(pos4(i))) / ntn
156c
157 ! Computing the product XDAMP*NtN*VNL
158 ntn_vnl = (vnl(pos1(i)) + vnl(pos2(i)) + vnl(pos3(i)) + vnl(pos4(i))) / ntn
159 IF (nodadt > 0) THEN
160 ntn_vnl = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
161 . sqrt(mass(pos2(i))/mass0(pos2(i))),
162 . sqrt(mass(pos3(i))/mass0(pos3(i))),
163 . sqrt(mass(pos4(i))/mass0(pos4(i))))*ntn_vnl
164 ENDIF
165c
166 ! Computation of the product LEN**2 * BtxB
167 b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i))
168 . + btb13(i)*unl(pos3(i)) + btb14(i)*unl(pos4(i)))
169c
170 b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i))
171 . + btb23(i)*unl(pos3(i)) + btb24(i)*unl(pos4(i)))
172c
173 b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i))
174 . + btb33(i)*unl(pos3(i)) + btb34(i)*unl(pos4(i)))
175c
176 b4 = l2 * vol(i) * ( btb14(i)*unl(pos1(i)) + btb24(i)*unl(pos2(i))
177 . + btb34(i)*unl(pos3(i)) + btb44(i)*unl(pos4(i)))
178c
179 ! Multiplication by the volume of the element
180 ntn_unl = ntn_unl * vol(i)
181 ntn_vnl = ntn_vnl * xi * vol(i)
182c
183 ! Introducing the internal variable to be regularized
184 ntvar = var_reg(i)*fourth* vol(i)
185c
186 ! Computing the elementary non-local forces
187 a = ntn_unl + ntn_vnl - ntvar
188 f1(i) = a + b1
189 f2(i) = a + b2
190 f3(i) = a + b3
191 f4(i) = a + b4
192c
193 ! Computing nodal equivalent stiffness
194 IF (nodadt > 0) THEN
195 sti1(i) = (abs(l2*btb11(i) + one/ntn) + abs(l2*btb12(i) + one/ntn) + abs(l2*btb13(i) + one/ntn) +
196 . abs(l2*btb14(i) + one/ntn))* vol(i)
197 sti2(i) = (abs(l2*btb12(i) + one/ntn) + abs(l2*btb22(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) +
198 . abs(l2*btb24(i) + one/ntn))* vol(i)
199 sti3(i) = (abs(l2*btb13(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) + abs(l2*btb33(i) + one/ntn) +
200 . abs(l2*btb34(i) + one/ntn))* vol(i)
201 sti4(i) = (abs(l2*btb14(i) + one/ntn) + abs(l2*btb24(i) + one/ntn) + abs(l2*btb34(i) + one/ntn) +
202 . abs(l2*btb44(i) + one/ntn))* vol(i)
203 ENDIF
204c
205 ! If the element is broken, computation of absorbing forces
206 ELSE
207c
208 ! Initial element characteristic length
209 lc(i) = vol0(i)**third
210c
211 IF (nodadt > 0) THEN
212 ! Non-local absorbing forces
213 f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*zeta*sspnl*half*
214 . (vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
215 f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*zeta*sspnl*half*
216 . (vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
217 f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*zeta*sspnl*half*
218 . (vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
219 f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*zeta*sspnl*half*
220 . (vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)
221 ! Computing nodal equivalent stiffness
222 sti1(i) = em20
223 sti2(i) = em20
224 sti3(i) = em20
225 sti4(i) = em20
226 ELSE
227 ! Non-local absorbing forces
228 f1(i) = zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
229 f2(i) = zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
230 f3(i) = zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
231 f4(i) = zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)
232 ENDIF
233 ENDIF
234 ENDDO
235c-----------------------------------------------------------------------
236c Assemblage
237c-----------------------------------------------------------------------
238 ! If PARITH/OFF
239 IF (iparit == 0) THEN
240 fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
241 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
242 DO i=1,nel
243 ! Assembling the forces in the classic way
244 fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
245 fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
246 fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
247 fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
248 IF (nodadt > 0) THEN
249 ! Spectral radius of stiffness matrix
250 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
251 ! computing nodal stiffness
252 stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
253 stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
254 stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
255 stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
256 ENDIF
257 ENDDO
258c
259 ! If PARITH/ON
260 ELSE
261c
262 DO i=1,nel
263 ii = i + nft
264c
265 ! Spectral radius of stiffness matrix
266 IF (nodadt > 0) THEN
267 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
268 ENDIF
269c
270 k = nloc_dmg%IADS(1,ii)
271 nloc_dmg%FSKY(k,1) = -f1(i)
272 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
273c
274 k = nloc_dmg%IADS(3,ii)
275 nloc_dmg%FSKY(k,1) = -f2(i)
276 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
277c
278 k = nloc_dmg%IADS(6,ii)
279 nloc_dmg%FSKY(k,1) = -f3(i)
280 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
281c
282 k = nloc_dmg%IADS(5,ii)
283 nloc_dmg%FSKY(k,1) = -f4(i)
284 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
285c
286 ENDDO
287 ENDIF
288c
289 !-----------------------------------------------------------------------
290 ! Computing non-local timestep
291 !-----------------------------------------------------------------------
292 IF (nodadt == 0) THEN
293 DO i = 1,nel
294 ! If the element is not broken, normal computation
295 IF (off(i)/=zero) THEN
296 ! Non-local critical time-step in the plane
297 dtnl = (two*(min(vol(i)**third,le_max))*sqrt(three*zeta))/
298 . sqrt(twelve*l2 + (min(vol(i)**third,le_max))**2)
299 ! Retaining the minimal value
300 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl)
301 ENDIF
302 ENDDO
303 ENDIF
304c
305c-----------
306 ! Deallocation of tables
307 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
308 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
309 IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
310 IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
311 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
312 IF (ALLOCATED(btb23)) DEALLOCATE(btb23)
313 IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
314 IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
315 IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
316 IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
317 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
318 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
319 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
320 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
321 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
322 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
323 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
324 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
325c
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:117