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

Go to the source code of this file.

Functions/Subroutines

subroutine s6cfint_reg (nloc_dmg, var_reg, nel, off, vol, nc1, nc2, nc3, nc4, nc5, nc6, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, imat, itask, dt2t, vol0, nft, nlay, ws, as, area, bufnlts)

Function/Subroutine Documentation

◆ s6cfint_reg()

subroutine s6cfint_reg ( type(nlocal_str_), intent(inout), target nloc_dmg,
intent(inout) var_reg,
integer, intent(in) nel,
intent(in) off,
intent(in) vol,
integer, dimension(nel), intent(in) nc1,
integer, dimension(nel), intent(in) nc2,
integer, dimension(nel), intent(in) nc3,
integer, dimension(nel), intent(in) nc4,
integer, dimension(nel), intent(in) nc5,
integer, dimension(nel), intent(in) nc6,
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, intent(in) imat,
integer, intent(in) itask,
intent(inout) dt2t,
intent(in) vol0,
integer, intent(in) nft,
integer, intent(in) nlay,
intent(in) ws,
intent(in) as,
intent(in) area,
type(buf_nlocts_), intent(inout), target bufnlts )

Definition at line 31 of file s6cfint_reg.F.

41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE nlocal_reg_mod
45 USE elbufdef_mod
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C G l o b a l P a r a m e t e r s
52C-----------------------------------------------
53#include "parit_c.inc"
54#include "scr02_c.inc"
55#include "scr18_c.inc"
56#include "com08_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER, INTENT(IN) ::
61 . NFT,NLAY,NEL,IMAT,ITASK
62 INTEGER, INTENT(IN), DIMENSION(NEL) ::
63 . NC1,NC2,NC3,NC4,NC5,NC6
64 my_real, INTENT(INOUT) ::
65 . dt2t
66 my_real, DIMENSION(9,9), INTENT(IN) ::
67 . ws,as
68 my_real, DIMENSION(NEL,NLAY), INTENT(INOUT) ::
69 . var_reg
70 my_real, DIMENSION(NEL), INTENT(IN) ::
71 . vol,off,vol0,px1,px2,px3,px4,area,
72 . py1,py2,py3,py4,pz1,pz2,pz3,pz4
73 TYPE(NLOCAL_STR_), INTENT(INOUT), TARGET :: NLOC_DMG
74 TYPE(BUF_NLOCTS_), INTENT(INOUT), TARGET :: BUFNLTS
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I,II,J,K,NNOD,N1,N2,N3,N4,N5,N6,
79 . L_NLOC,NDDL,NDNOD
80 my_real
81 . l2,ntn,ntn_unl,ntn_vnl,xi,ntvar,a,dtnl,le_max,
82 . b1,b2,b3,b4,b5,b6,zeta,sspnl,maxstif,
83 . bth1,bth2,nth1,nth2,dt2p,dtnod,k1,k2,k12,
84 . dtnl_th
85 my_real, DIMENSION(NEL,NLAY) ::
86 . f1,f2,f3,f4,f5,f6
87 my_real, DIMENSION(NEL) ::
88 . lc,thk,lthk,pxx1,pxx2,pxx3,
89 . pxx4,pxx5,pxx6,pyy1,pyy2,pyy3,
90 . pyy4,pyy5,pyy6,pzz1,pzz2,pzz3,
91 . pzz4,pzz5,pzz6
92 my_real, DIMENSION(:) ,ALLOCATABLE ::
93 . btb11,btb12,btb13,btb14,btb15,btb16,
94 . btb22,btb23,btb24,btb25,btb26,btb33,
95 . btb34,btb35,btb36,btb44,btb45,btb46,
96 . btb55,btb56,btb66
97 my_real, DIMENSION(:,:) ,ALLOCATABLE ::
98 . sti1,sti2,sti3,sti4,sti5,sti6
99 INTEGER, DIMENSION(:), ALLOCATABLE ::
100 . POS1,POS2,POS3,POS4,POS5,POS6
101 my_real, POINTER, DIMENSION(:) ::
102 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
103 my_real, POINTER, DIMENSION(:,:) ::
104 . massth,unlth,vnlth,fnlth
105 my_real, DIMENSION(:,:), ALLOCATABLE ::
106 . stifnlth,dtn
107 my_real, target :: nothing(1)
108 ! Safety coefficient for non-local stability vs mechanical stability
109 ! (it has been slightly increased vs nloc_dmg_init.F)
110 my_real, PARAMETER :: csta = 40.0d0
111 ! Coefficient for non-local stability to take into account damping
112 my_real, PARAMETER :: cdamp = 0.7d0
113 my_real
114 . zs(10,9)
115 ! Position of nodes in the thickshell thickness
116 DATA zs /
117 1 0. ,0. ,0. ,
118 1 0. ,0. ,0. ,
119 1 0. ,0. ,0. ,
120 1 0. ,
121 2 -1. ,0. ,1. ,
122 2 0. ,0. ,0. ,
123 2 0. ,0. ,0. ,
124 2 0. ,
125 3 -1. ,-.549193338482966,0.549193338482966,
126 3 1. ,0. ,0. ,
127 3 0. ,0. ,0. ,
128 3 0. ,
129 4 -1. ,-.600558677589454,0. ,
130 4 0.600558677589454,1. ,0. ,
131 4 0. ,0. ,0. ,
132 4 0. ,
133 5 -1. ,-.812359691877328,-.264578928334038,
134 5 0.264578928334038,0.812359691877328,1. ,
135 5 0. ,0. ,0. ,
136 5 0. ,
137 6 -1. ,-.796839450334708,-.449914286274731,
138 6 0. ,0.449914286274731,0.796839450334708,
139 6 1. ,0. ,0. ,
140 6 0. ,
141 7 -1. ,-.898215824685518,-.584846546513270,
142 7 -.226843756241524,0.226843756241524,0.584846546513270,
143 7 0.898215824685518,1. ,0. ,
144 7 0. ,
145 8 -1. ,-.878478166955581,-.661099443664978,
146 8 -.354483526205989,0. ,0.354483526205989,
147 8 0.661099443664978,0.878478166955581,1. ,
148 8 0. ,
149 9 -1. ,-.936320479015252,-.735741735638020,
150 9 -.491001129763160,-.157505717044458,0.157505717044458,
151 9 0.491001129763160,0.735741735638020,0.936320479015252,
152 9 1. /
153C=======================================================================
154 stifnl => nothing
155 l2 = nloc_dmg%LEN(imat)**2
156 xi = nloc_dmg%DAMP(imat)
157 nnod = nloc_dmg%NNOD
158 l_nloc = nloc_dmg%L_NLOC
159 zeta = nloc_dmg%DENS(imat)
160 sspnl = nloc_dmg%SSPNL(imat)
161 le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
162 ntn = six*six
163 lc(1:nel) = zero
164 ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb15(nel),
165 . btb16(nel),btb22(nel),btb23(nel),btb24(nel),btb25(nel),
166 . btb26(nel),btb33(nel),btb34(nel),btb35(nel),btb36(nel),
167 . btb44(nel),btb45(nel),btb46(nel),btb55(nel),btb56(nel),
168 . btb66(nel),pos1(nel) ,pos2(nel) ,pos3(nel) ,pos4(nel) ,
169 . pos5(nel) ,pos6(nel) )
170 ! For nodal timestep
171 IF (nodadt > 0) THEN
172 ! Non-local nodal stifness
173 ALLOCATE(sti1(nel,nlay),sti2(nel,nlay),sti3(nel,nlay),
174 . sti4(nel,nlay),sti5(nel,nlay),sti6(nel,nlay))
175 ! Non-local mass
176 mass => nloc_dmg%MASS(1:l_nloc)
177 ! Initial non-local mass
178 mass0 => nloc_dmg%MASS0(1:l_nloc)
179 ELSE
180 NULLIFY(mass)
181 NULLIFY(mass0)
182 ALLOCATE(sti1(1,1),sti2(1,1),sti3(1,1),
183 . sti4(1,1),sti5(1,1),sti6(1,1))
184 ENDIF
185 vnl => nloc_dmg%VNL(1:l_nloc)
186 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
187 unl => nloc_dmg%UNL(1:l_nloc)
188c
189 !-----------------------------------------------------------------------
190 ! Computation of the element volume and the BtB matrix product
191 !-----------------------------------------------------------------------
192 ! Loop over elements
193# include "vectorize.inc"
194 DO i=1,nel
195c
196 ! Recovering the nodes of the brick element
197 n1 = nloc_dmg%IDXI(nc1(i))
198 n2 = nloc_dmg%IDXI(nc2(i))
199 n3 = nloc_dmg%IDXI(nc3(i))
200 n4 = nloc_dmg%IDXI(nc4(i))
201 n5 = nloc_dmg%IDXI(nc5(i))
202 n6 = nloc_dmg%IDXI(nc6(i))
203c
204 ! Recovering the positions of the first d.o.fs of each nodes
205 pos1(i) = nloc_dmg%POSI(n1)
206 pos2(i) = nloc_dmg%POSI(n2)
207 pos3(i) = nloc_dmg%POSI(n3)
208 pos4(i) = nloc_dmg%POSI(n4)
209 pos5(i) = nloc_dmg%POSI(n5)
210 pos6(i) = nloc_dmg%POSI(n6)
211c
212 ! Computation of derivatives of shape functions
213 pxx1(i) = px1(i)-px4(i)
214 pyy1(i) = py1(i)-py4(i)
215 pzz1(i) = pz1(i)-pz4(i)
216c
217 pxx2(i) = px2(i)-px4(i)
218 pyy2(i) = py2(i)-py4(i)
219 pzz2(i) = pz2(i)-pz4(i)
220c
221 pxx3(i) = px3(i)-px4(i)
222 pyy3(i) = py3(i)-py4(i)
223 pzz3(i) = pz3(i)-pz4(i)
224c
225 pxx4(i) = px1(i)+px4(i)
226 pyy4(i) = py1(i)+py4(i)
227 pzz4(i) = pz1(i)+pz4(i)
228c
229 pxx5(i) = px2(i)+px4(i)
230 pyy5(i) = py2(i)+py4(i)
231 pzz5(i) = pz2(i)+pz4(i)
232c
233 pxx6(i) = px3(i)+px4(i)
234 pyy6(i) = py3(i)+py4(i)
235 pzz6(i) = pz3(i)+pz4(i)
236c
237 ! Computation of the product BtxB
238 btb11(i) = pxx1(i)**2 + pyy1(i)**2 + pzz1(i)**2
239 btb12(i) = pxx1(i)*pxx2(i) + pyy1(i)*pyy2(i) + pzz1(i)*pzz2(i)
240 btb13(i) = pxx1(i)*pxx3(i) + pyy1(i)*pyy3(i) + pzz1(i)*pzz3(i)
241 btb14(i) = pxx1(i)*pxx4(i) + pyy1(i)*pyy4(i) + pzz1(i)*pzz4(i)
242 btb15(i) = pxx1(i)*pxx5(i) + pyy1(i)*pyy5(i) + pzz1(i)*pzz5(i)
243 btb16(i) = pxx1(i)*pxx6(i) + pyy1(i)*pyy6(i) + pzz1(i)*pzz6(i)
244c
245 btb22(i) = pxx2(i)**2 + pyy2(i)**2 + pzz2(i)**2
246 btb23(i) = pxx2(i)*pxx3(i) + pyy2(i)*pyy3(i) + pzz2(i)*pzz3(i)
247 btb24(i) = pxx2(i)*pxx4(i) + pyy2(i)*pyy4(i) + pzz2(i)*pzz4(i)
248 btb25(i) = pxx2(i)*pxx5(i) + pyy2(i)*pyy5(i) + pzz2(i)*pzz5(i)
249 btb26(i) = pxx2(i)*pxx6(i) + pyy2(i)*pyy6(i) + pzz2(i)*pzz6(i)
250c
251 btb33(i) = pxx3(i)**2 + pyy3(i)**2 + pzz3(i)**2
252 btb34(i) = pxx3(i)*pxx4(i) + pyy3(i)*pyy4(i) + pzz3(i)*pzz4(i)
253 btb35(i) = pxx3(i)*pxx5(i) + pyy3(i)*pyy5(i) + pzz3(i)*pzz5(i)
254 btb36(i) = pxx3(i)*pxx6(i) + pyy3(i)*pyy6(i) + pzz3(i)*pzz6(i)
255c
256 btb44(i) = pxx4(i)**2 + pyy4(i)**2 + pzz4(i)**2
257 btb45(i) = pxx4(i)*pxx5(i) + pyy4(i)*pyy5(i) + pzz4(i)*pzz5(i)
258 btb46(i) = pxx4(i)*pxx6(i) + pyy4(i)*pyy6(i) + pzz4(i)*pzz6(i)
259c
260 btb55(i) = pxx5(i)**2 + pyy5(i)**2 + pzz5(i)**2
261 btb56(i) = pxx5(i)*pxx6(i) + pyy5(i)*pyy6(i) + pzz5(i)*pzz6(i)
262c
263 btb66(i) = pxx6(i)**2 + pyy6(i)**2 + pzz6(i)**2
264c
265 ENDDO
266c
267 !-----------------------------------------------------------------------
268 ! Pre-treatment non-local regularization in the thickshell thickness
269 !-----------------------------------------------------------------------
270 IF ((l2>zero).AND.(nlay > 1)) THEN
271c
272 ! Compute thickshell thickness
273 DO i = 1,nel
274 thk(i) = vol(i)/area(i)
275 lthk(i) = (zs(nlay+1,nlay)-zs(nlay,nlay))*thk(i)*half
276 ENDDO
277c
278 ! Allocation of the velocities predictor
279 nddl = nlay
280 IF (nodadt > 0) THEN
281 ALLOCATE(stifnlth(nel,nddl+1))
282 ALLOCATE(dtn(nel,nddl+1))
283 ELSE
284 ALLOCATE(dtn(1,1))
285 ALLOCATE(stifnlth(1,1))
286 dtn(1,1) = ep20
287 stifnlth(1,1) = ep20
288 ENDIF
289 ndnod = nddl+1
290c
291 ! Pointing the non-local values in the thickness of the corresponding element
292 massth => bufnlts%MASSTH(1:nel,1:ndnod)
293 unlth => bufnlts%UNLTH(1:nel ,1:ndnod)
294 vnlth => bufnlts%VNLTH(1:nel ,1:ndnod)
295 fnlth => bufnlts%FNLTH(1:nel ,1:ndnod)
296c
297 DO k = 1,ndnod
298 DO i = 1,nel
299 ! Resetting non-local forces
300 fnlth(i,k) = zero
301 ! Resetting non-local nodal stiffness
302 IF (nodadt > 0) THEN
303 stifnlth(i,k) = em20
304 ENDIF
305 ENDDO
306 ENDDO
307c
308 ! Computation of non-local forces in the shell thickness
309 DO k = 1, nddl
310c
311 ! Computation of shape functions value
312 nth1 = (as(k,nddl) - zs(k+1,nddl)) /
313 . (zs(k,nddl) - zs(k+1,nddl))
314 nth2 = (as(k,nddl) - zs(k,nddl)) /
315 . (zs(k+1,nddl) - zs(k,nddl))
316c
317 ! Loop over elements
318 DO i = 1,nel
319c
320 ! Computation of B-matrix values
321 bth1 = (one/(zs(k,nddl) - zs(k+1,nddl)))*(two/thk(i))
322 bth2 = (one/(zs(k+1,nddl) - zs(k,nddl)))*(two/thk(i))
323c
324 ! Computation of the non-local K matrix
325 k1 = l2*(bth1**2) + nth1**2
326 k12 = l2*(bth1*bth2)+ (nth1*nth2)
327 k2 = l2*(bth2**2) + nth2**2
328c
329 ! Computation of the non-local forces
330 fnlth(i,k) = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1)
331 . + xi*((nth1**2)*vnlth(i,k)
332 . + (nth1*nth2)*vnlth(i,k+1))
333 . - (nth1*var_reg(i,k)))*half*ws(k,nddl)*vol(i)
334 fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i,k) + k2*unlth(i,k+1)
335 . + xi*(nth1*nth2*vnlth(i,k)
336 . + (nth2**2)*vnlth(i,k+1))
337 . - nth2*var_reg(i,k))*half*ws(k,nddl)*vol(i)
338c
339 ! Computation of non-local nodal stiffness
340 IF (nodadt > 0) THEN
341 stifnlth(i,k) = stifnlth(i,k) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)
342 stifnlth(i,k+1) = stifnlth(i,k+1) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)
343 ENDIF
344c
345 ENDDO
346 ENDDO
347c
348 ! Updating non-local mass with /DT/NODA
349 IF (nodadt > 0) THEN
350C
351 ! Initial computation of the nodal timestep
352 dtnod = ep20
353 DO k = 1,ndnod
354 DO i = 1,nel
355 dtn(i,k) = dtfac1(11)*cdamp*sqrt(two*massth(i,k)/max(stifnlth(i,k),em20))
356 dtnod = min(dtn(i,k),dtnod)
357 ENDDO
358 ENDDO
359C
360 ! /DT/NODA/CSTX - Constant timestep with added mass
361 IF ((idtmin(11)==3).OR.(idtmin(11)==4).OR.(idtmin(11)==8)) THEN
362 ! Added mass computation if necessary
363 IF (dtnod < dtmin1(11)*(sqrt(csta))) THEN
364 DO k = 1,ndnod
365 DO i = 1,nel
366 IF (dtn(i,k) < dtmin1(11)) THEN
367 dt2p = dtmin1(11)/(dtfac1(11)*cdamp)
368 massth(i,k) = max(massth(i,k),csta*half*stifnlth(i,k)*dt2p*dt2p*onep00001)
369 ENDIF
370 ENDDO
371 ENDDO
372 ENDIF
373 dtnod = dtmin1(11)*(sqrt(csta))
374 ENDIF
375C
376 ! Classical nodal timestep check
377 IF (dtnod < dt2t) THEN
378 dt2t = min(dt2t,dtnod)
379 ENDIF
380 ENDIF
381c
382 DO k = 1,ndnod
383 DO i = 1,nel
384 ! Updating the non-local in-thickness velocities
385 vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt12
386 ENDDO
387 ENDDO
388c
389 DO k = 1,ndnod
390 DO i = 1,nel
391 ! Computing the non-local in-thickness cumulated values
392 unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt1
393 ENDDO
394 ENDDO
395c
396 ! Transfert at the integration point
397 DO k = 1, nddl
398 !Computation of shape functions value
399 nth1 = (as(k,nddl) - zs(k+1,nddl))/
400 . (zs(k,nddl) - zs(k+1,nddl))
401 nth2 = (as(k,nddl) - zs(k,nddl))/
402 . (zs(k+1,nddl) - zs(k,nddl))
403 ! Loop over elements
404 DO i = 1,nel
405 !Integration points non-local variables
406 var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i,k+1)
407 ENDDO
408 ENDDO
409 ENDIF
410c
411 !-----------------------------------------------------------------------
412 ! Computation of non-local forces
413 !-----------------------------------------------------------------------
414 ! Loop over elements
415 DO k = 1,nlay
416c
417 ! Loop over elements
418# include "vectorize.inc"
419 DO i = 1, nel
420c
421 ! If the element is not broken, normal computation
422 IF (off(i) /= zero) THEN
423c
424 ! Computing the product NtN*UNL
425 ntn_unl = (unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) + unl(pos4(i)+k-1)
426 . + unl(pos5(i)+k-1) + unl(pos6(i)+k-1)) / ntn
427c
428 ! Computing the product XDAMP*NtN*VNL
429 ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) + vnl(pos4(i)+k-1)
430 . + vnl(pos5(i)+k-1) + vnl(pos6(i)+k-1)) / ntn
431 IF (nodadt > 0) THEN
432 ntn_vnl = min(sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1)),
433 . sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1)),
434 . sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1)),
435 . sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1)),
436 . sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1)),
437 . sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1)))*ntn_vnl
438 ENDIF
439c
440 ! Computation of the product LEN**2 * BtxB
441 b1 = l2 * vol(i) * ws(k,nlay) *half * ( btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
442 . + btb13(i)*unl(pos3(i)+k-1) + btb14(i)*unl(pos4(i)+k-1) + btb15(i)*unl(pos5(i)+k-1)
443 . + btb16(i)*unl(pos6(i)+k-1) )
444c
445 b2 = l2 * vol(i) * ws(k,nlay) *half * ( btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
446 . + btb23(i)*unl(pos3(i)+k-1) + btb24(i)*unl(pos4(i)+k-1) + btb25(i)*unl(pos5(i)+k-1)
447 . + btb26(i)*unl(pos6(i)+k-1) )
448c
449 b3 = l2 * vol(i) * ws(k,nlay) *half * ( btb13(i)*unl(pos1(i)+k-1) + btb23(i)*unl(pos2(i)+k-1)
450 . + btb33(i)*unl(pos3(i)+k-1) + btb34(i)*unl(pos4(i)+k-1) + btb35(i)*unl(pos5(i)+k-1)
451 . + btb36(i)*unl(pos6(i)+k-1) )
452c
453 b4 = l2 * vol(i) * ws(k,nlay) *half * ( btb14(i)*unl(pos1(i)+k-1) + btb24(i)*unl(pos2(i)+k-1)
454 . + btb34(i)*unl(pos3(i)+k-1) + btb44(i)*unl(pos4(i)+k-1) + btb45(i)*unl(pos5(i)+k-1)
455 . + btb46(i)*unl(pos6(i)+k-1) )
456c
457 b5 = l2 * vol(i) * ws(k,nlay) *half * ( btb15(i)*unl(pos1(i)+k-1) + btb25(i)*unl(pos2(i)+k-1)
458 . + btb35(i)*unl(pos3(i)+k-1) + btb45(i)*unl(pos4(i)+k-1) + btb55(i)*unl(pos5(i)+k-1)
459 . + btb56(i)*unl(pos6(i)+k-1) )
460c
461 b6 = l2 * vol(i) * ws(k,nlay) *half * ( btb16(i)*unl(pos1(i)+k-1) + btb26(i)*unl(pos2(i)+k-1)
462 . + btb36(i)*unl(pos3(i)+k-1) + btb46(i)*unl(pos4(i)+k-1) + btb56(i)*unl(pos5(i)+k-1)
463 . + btb66(i)*unl(pos6(i)+k-1) )
464c
465 ! Multiplication by the volume of the element
466 ntn_unl = ntn_unl * vol(i) * ws(k,nlay) * half
467 ntn_vnl = ntn_vnl * xi * vol(i) * ws(k,nlay) * half
468c
469 ! Introducing the internal variable to be regularized
470 ntvar = var_reg(i,k)*one_over_6* vol(i) * ws(k,nlay) * half
471c
472 ! Computing the elementary non-local forces
473 a = ntn_unl + ntn_vnl - ntvar
474 f1(i,k) = a + b1
475 f2(i,k) = a + b2
476 f3(i,k) = a + b3
477 f4(i,k) = a + b4
478 f5(i,k) = a + b5
479 f6(i,k) = a + b6
480c
481 ! Computing nodal equivalent stiffness
482 IF (nodadt > 0) THEN
483 sti1(i,k) = (abs(l2*btb11(i) + one/ntn) + abs(l2*btb12(i) + one/ntn) + abs(l2*btb13(i) + one/ntn) +
484 . abs(l2*btb14(i) + one/ntn) + abs(l2*btb15(i) + one/ntn) + abs(l2*btb16(i) + one/ntn))
485 . *vol(i)*ws(k,nlay)*half
486 sti2(i,k) = (abs(l2*btb12(i) + one/ntn) + abs(l2*btb22(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) +
487 . abs(l2*btb24(i) + one/ntn) + abs(l2*btb25(i) + one/ntn) + abs(l2*btb26(i) + one/ntn))
488 . *vol(i)*ws(k,nlay)*half
489 sti3(i,k) = (abs(l2*btb13(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) + abs(l2*btb33(i) + one/ntn) +
490 . abs(l2*btb34(i) + one/ntn) + abs(l2*btb35(i) + one/ntn) + abs(l2*btb36(i) + one/ntn))
491 . *vol(i)*ws(k,nlay)*half
492 sti4(i,k) = (abs(l2*btb14(i) + one/ntn) + abs(l2*btb24(i) + one/ntn) + abs(l2*btb34(i) + one/ntn) +
493 . abs(l2*btb44(i) + one/ntn) + abs(l2*btb45(i) + one/ntn) + abs(l2*btb46(i) + one/ntn))
494 . *vol(i)*ws(k,nlay)*half
495 sti5(i,k) = (abs(l2*btb15(i) + one/ntn) + abs(l2*btb25(i) + one/ntn) + abs(l2*btb35(i) + one/ntn) +
496 . abs(l2*btb45(i) + one/ntn) + abs(l2*btb55(i) + one/ntn) + abs(l2*btb56(i) + one/ntn))
497 . *vol(i)*ws(k,nlay)*half
498 sti6(i,k) = (abs(l2*btb16(i) + one/ntn) + abs(l2*btb26(i) + one/ntn) + abs(l2*btb36(i) + one/ntn) +
499 . abs(l2*btb46(i) + one/ntn) + abs(l2*btb56(i) + one/ntn) + abs(l2*btb66(i) + one/ntn))
500 . *vol(i)*ws(k,nlay)*half
501 ENDIF
502c
503 ! If the element is broken, the non-local wave is absorbed
504 ELSE
505c
506 ! Initial element characteristic length
507 lc(i) = (vol0(i)*ws(k,nlay)*half)**third
508c
509 IF (nodadt > 0) THEN
510
511 ! Non-local absorbing forces
512 f1(i,k) = sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*half*
513 . (vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
514 f2(i,k) = sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*half*
515 . (vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
516 f3(i,k) = sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*zeta*sspnl*half*
517 . (vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
518 f4(i,k) = sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*half*
519 . (vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
520 f5(i,k) = sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1))*zeta*sspnl*half*
521 . (vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
522 f6(i,k) = sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1))*zeta*sspnl*half*
523 . (vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)
524 ! Computing nodal equivalent stiffness
525 sti1(i,k) = em20
526 sti2(i,k) = em20
527 sti3(i,k) = em20
528 sti4(i,k) = em20
529 sti5(i,k) = em20
530 sti6(i,k) = em20
531 ELSE
532 ! Non-local absorbing forces
533 f1(i,k) = zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
534 f2(i,k) = zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
535 f3(i,k) = zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
536 f4(i,k) = zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
537 f5(i,k) = zeta*sspnl*half*(vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
538 f6(i,k) = zeta*sspnl*half*(vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)
539 ENDIF
540 ENDIF
541 ENDDO
542 ENDDO
543c
544 !-----------------------------------------------------------------------
545 ! Assembling of the non-local forces
546 !-----------------------------------------------------------------------
547c
548 ! If PARITH/OFF
549 IF (iparit == 0) THEN
550 ! Recovering non-local internal forces
551 fnl => nloc_dmg%FNL(1:l_nloc,itask+1) ! Non-local forces
552 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
553 ! Loop over elements
554 DO i=1,nel
555 ! Loop over non-local degrees of freedom
556# include "vectorize.inc"
557 DO k=1,nlay
558 ! Assembling the forces in the classic way
559 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
560 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
561 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
562 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
563 fnl(pos5(i)+k-1) = fnl(pos5(i)+k-1) - f5(i,k)
564 fnl(pos6(i)+k-1) = fnl(pos6(i)+k-1) - f6(i,k)
565 IF (nodadt > 0) THEN
566 ! Spectral radius of stiffness matrix
567 maxstif = max(sti1(i,k),sti2(i,k),sti3(i,k),
568 . sti4(i,k),sti5(i,k),sti6(i,k))
569 ! Computing nodal stiffness
570 stifnl(pos1(i)+k-1) = stifnl(pos1(i)+k-1) + maxstif
571 stifnl(pos2(i)+k-1) = stifnl(pos2(i)+k-1) + maxstif
572 stifnl(pos3(i)+k-1) = stifnl(pos3(i)+k-1) + maxstif
573 stifnl(pos4(i)+k-1) = stifnl(pos4(i)+k-1) + maxstif
574 stifnl(pos5(i)+k-1) = stifnl(pos5(i)+k-1) + maxstif
575 stifnl(pos6(i)+k-1) = stifnl(pos6(i)+k-1) + maxstif
576 ENDIF
577 ENDDO
578 ENDDO
579c
580 ! If PARITH/ON
581 ELSE
582 ! Loop over additional d.o.fs
583 DO j = 1,nlay
584c
585 ! Loop over elements
586 DO i=1,nel
587 ii = i + nft
588c
589 ! Spectral radius of stiffness matrix
590 IF (nodadt > 0) THEN
591 maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j),
592 . sti5(i,j),sti6(i,j))
593 ENDIF
594c
595 k = nloc_dmg%IADS(1,ii)
596 nloc_dmg%FSKY(k,j) = -f1(i,j)
597 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
598c
599 k = nloc_dmg%IADS(2,ii)
600 nloc_dmg%FSKY(k,j) = -f2(i,j)
601 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
602c
603 k = nloc_dmg%IADS(3,ii)
604 nloc_dmg%FSKY(k,j) = -f3(i,j)
605 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
606c
607 k = nloc_dmg%IADS(5,ii)
608 nloc_dmg%FSKY(k,j) = -f4(i,j)
609 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
610c
611 k = nloc_dmg%IADS(6,ii)
612 nloc_dmg%FSKY(k,j) = -f5(i,j)
613 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
614c
615 k = nloc_dmg%IADS(7,ii)
616 nloc_dmg%FSKY(k,j) = -f6(i,j)
617 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
618c
619 ENDDO
620 ENDDO
621 ENDIF
622c
623 !-----------------------------------------------------------------------
624 ! Computing non-local timestep
625 !-----------------------------------------------------------------------
626 IF (nodadt == 0) THEN
627 DO i = 1,nel
628 ! If the element is not broken, normal computation
629 IF (off(i)/=zero) THEN
630 ! Non-local critical time-step in the plane
631 dtnl = (two*(min((vol(i))**third,le_max))*sqrt(three*zeta))/
632 . sqrt(twelve*l2 + (min((vol(i))**third,le_max))**2)
633 ! Non-local critical time-step in the thickness
634 IF (nlay > 1) THEN
635 dtnl_th = (two*(min(lthk(i),le_max))*sqrt(three*zeta))/
636 . sqrt(twelve*l2 + (min(lthk(i),le_max))**2)
637 ELSE
638 dtnl_th = ep20
639 ENDIF
640 ! Retaining the minimal value
641 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl_th,dtfac1(1)*cdamp*dtnl)
642 ENDIF
643 ENDDO
644 ENDIF
645c
646 ! Deallocation of tables
647 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
648 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
649 IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
650 IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
651 IF (ALLOCATED(btb15)) DEALLOCATE(btb15)
652 IF (ALLOCATED(btb16)) DEALLOCATE(btb16)
653 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
654 IF (ALLOCATED(btb23)) DEALLOCATE(btb23)
655 IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
656 IF (ALLOCATED(btb25)) DEALLOCATE(btb25)
657 IF (ALLOCATED(btb26)) DEALLOCATE(btb26)
658 IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
659 IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
660 IF (ALLOCATED(btb35)) DEALLOCATE(btb35)
661 IF (ALLOCATED(btb36)) DEALLOCATE(btb36)
662 IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
663 IF (ALLOCATED(btb45)) DEALLOCATE(btb45)
664 IF (ALLOCATED(btb46)) DEALLOCATE(btb46)
665 IF (ALLOCATED(btb55)) DEALLOCATE(btb55)
666 IF (ALLOCATED(btb56)) DEALLOCATE(btb56)
667 IF (ALLOCATED(btb66)) DEALLOCATE(btb66)
668 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
669 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
670 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
671 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
672 IF (ALLOCATED(pos5)) DEALLOCATE(pos5)
673 IF (ALLOCATED(pos6)) DEALLOCATE(pos6)
674 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
675 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
676 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
677 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
678 IF (ALLOCATED(sti5)) DEALLOCATE(sti5)
679 IF (ALLOCATED(sti6)) DEALLOCATE(sti6)
680 IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
681c-----------
my_real
#define my_real
Definition cppsort.cpp:32
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
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