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

Go to the source code of this file.

Functions/Subroutines

subroutine scfint_reg (nloc_dmg, var_reg, nel, off, vol, nc1, nc2, nc3, nc4, nc5, nc6, nc7, nc8, 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

◆ scfint_reg()

subroutine scfint_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) nc1,
integer, dimension(nel) nc2,
integer, dimension(nel) nc3,
integer, dimension(nel) nc4,
integer, dimension(nel) nc5,
integer, dimension(nel) nc6,
integer, dimension(nel) nc7,
integer, dimension(nel) nc8,
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 )
Parameters
[in,out]nloc_dmgNon-local data structure
[in]nelNumber of elements in the group
nc1Node 1 of the elements
nc2Node 2 of the elements
nc3Node 3 of the elements
nc4Node 4 of the elements
nc5Node 5 of the elements
nc6Node 6 of the elements
nc7Node 7 of the elements
nc8Node 8 of the elements
[in]imatMaterial internal Id
[in]itaskNumber of the thread
[in]nftAddress of the first element of the group
[in]nlayNumber of layers in the elements
[in,out]bufnltsNon-local specific buffer for thickshell thickness

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