OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
scfint_reg.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!|| scfint_reg ../engine/source/elements/thickshell/solidec/scfint_reg.F
25!||--- called by ------------------------------------------------------
26!|| scforc3 ../engine/source/elements/thickshell/solidec/scforc3.F
27!||--- uses -----------------------------------------------------
28!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
29!|| nlocal_reg_mod ../common_source/modules/nlocal_reg_mod.f
30!||====================================================================
31 SUBROUTINE scfint_reg(
32 1 NLOC_DMG,VAR_REG ,NEL ,OFF ,
33 2 VOL ,NLOCS ,AREA ,
34 3 PX1 ,PX2 ,PX3 ,PX4 ,
35 4 PY1 ,PY2 ,PY3 ,PY4 ,
36 5 PZ1 ,PZ2 ,PZ3 ,PZ4 ,
37 6 IMAT ,ITASK ,DT2T ,VOL0 ,
38 7 NFT ,NLAY ,WS ,AS ,
39 8 BUFNLTS )
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE nlocal_reg_mod
44 USE elbufdef_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "parit_c.inc"
53#include "scr02_c.inc"
54#include "scr18_c.inc"
55#include "com08_c.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER, INTENT(IN) :: NFT !< Address of the first element of the group
60 INTEGER, INTENT(IN) :: NLAY !< Number of layers in the elements
61 INTEGER, INTENT(IN) :: NEL !< Number of elements in the group
62 INTEGER, INTENT(IN) :: IMAT !< Material internal Id
63 INTEGER, INTENT(IN) :: ITASK !< Number of the thread
64 my_real, INTENT(INOUT) :: DT2T !< Element critical timestep
65 my_real, DIMENSION(9,9), INTENT(IN) :: WS !< Weight of integration points through thickness
66 my_real, DIMENSION(9,9), INTENT(IN) :: AS !< Position of integration points through thickness
67 my_real, DIMENSION(NEL,NLAY), INTENT(INOUT) :: var_reg !< Variable to regularize
68 my_real, DIMENSION(NEL), INTENT(IN) :: vol !< Current volume of the element
69 my_real, DIMENSION(NEL), INTENT(IN) :: off !< Flag for element deletion
70 my_real, DIMENSION(NEL), INTENT(IN) :: vol0 !< Initial volume of the element
71 my_real, DIMENSION(NEL), INTENT(IN) :: area !< Area of the thickshell element
72 my_real, DIMENSION(NEL), INTENT(IN) :: px1 !< Derivative of the first node shape function along x
73 my_real, DIMENSION(NEL), INTENT(IN) :: px2 !< Derivative of the second node shape function along x
74 my_real, DIMENSION(NEL), INTENT(IN) :: px3 !< Derivative of the third node shape function along x
75 my_real, DIMENSION(NEL), INTENT(IN) :: px4 !< Derivative of the fourth node shape function along x
76 my_real, DIMENSION(NEL), INTENT(IN) :: py1 !< Derivative of the first node shape function along y
77 my_real, DIMENSION(NEL), INTENT(IN) :: py2 !< Derivative of the second node shape function along y
78 my_real, DIMENSION(NEL), INTENT(IN) :: py3 !< Derivative of the third node shape function along y
79 my_real, DIMENSION(NEL), INTENT(IN) :: py4 !< Derivative of the fourth node shape function along y
80 my_real, DIMENSION(NEL), INTENT(IN) :: pz1 !< Derivative of the first node shape function along z
81 my_real, DIMENSION(NEL), INTENT(IN) :: pz2 !< Derivative of the second node shape function along z
82 my_real, DIMENSION(NEL), INTENT(IN) :: pz3 !< Derivative of the third node shape function along z
83 my_real, DIMENSION(NEL), INTENT(IN) :: pz4 !< Derivative of the fourth node shape function along z
84 TYPE(nlocal_str_), INTENT(INOUT), TARGET :: NLOC_DMG !< Non-local data structure
85 TYPE(buf_nlocts_), INTENT(INOUT), TARGET :: BUFNLTS !< Non-local specific buffer for thickshell thickness
86 TYPE(buf_nlocs_) , INTENT(IN) :: NLOCS !< Non-local effective nodes data structure
87C-----------------------------------------------
88C L o c a l V a r i a b l e s
89C-----------------------------------------------
90 INTEGER I,II,J,K,L,N1,N2,N3,N4,N5,N6,N7,N8,
91 . l_nloc,nddl,ndnod,nindx6,indx6(nel),
92 . nindx8,indx8(nel)
93 INTEGER :: NODA_DT
94 INTEGER, DIMENSION(:), ALLOCATABLE ::
95 . POS1,POS2,POS3,POS4,POS5,POS6,POS7,POS8
96 my_real
97 . L2,NTN_UNL,NTN_VNL,XI,NTVAR,A,DTNL,LE_MAX,
98 . b1,b2,b3,b4,b5,b6,b7,b8,zeta,sspnl,maxstif,
99 . bth1,bth2,nth1,nth2,dt2p,dtnod,k1,k2,k12,
100 . dtnl_th,ntn6,ntn8
101 my_real, DIMENSION(:) ,ALLOCATABLE ::
102 . btb11,btb12,btb13,btb14,btb22,btb23,btb24,
103 . btb33,btb34,btb44,lc,thk,lthk
104 my_real, DIMENSION(:,:) ,ALLOCATABLE ::
105 . f1,f2,f3,f4,f5,f6,f7,f8,stifnlth,dtn,
106 . sti1,sti2,sti3,sti4,sti5,sti6,sti7,sti8
107 my_real, POINTER, DIMENSION(:) ::
108 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
109 my_real, POINTER, DIMENSION(:,:) ::
110 . massth,unlth,vnlth,fnlth
111 ! Safety coefficient for non-local stability vs mechanical stability
112 my_real, PARAMETER :: csta = 40.0d0
113 ! Coefficient for non-local stability to take into account damping
114 my_real, PARAMETER :: cdamp = 0.7d0
115 my_real
116 . zs(10,9)
117 ! Position of nodes in the thickshell thickness
118 DATA zs /
119 1 0. ,0. ,0. ,
120 1 0. ,0. ,0. ,
121 1 0. ,0. ,0. ,
122 1 0. ,
123 2 -1. ,0. ,1. ,
124 2 0. ,0. ,0. ,
125 2 0. ,0. ,0. ,
126 2 0. ,
127 3 -1. ,-.549193338482966,0.549193338482966,
128 3 1. ,0. ,0. ,
129 3 0. ,0. ,0. ,
130 3 0. ,
131 4 -1. ,-.600558677589454,0. ,
132 4 0.600558677589454,1. ,0. ,
133 4 0. ,0. ,0. ,
134 4 0. ,
135 5 -1. ,-.812359691877328,-.264578928334038,
136 5 0.264578928334038,0.812359691877328,1. ,
137 5 0. ,0. ,0. ,
138 5 0. ,
139 6 -1. ,-.796839450334708,-.449914286274731,
140 6 0. ,0.449914286274731,0.796839450334708,
141 6 1. ,0. ,0. ,
142 6 0. ,
143 7 -1. ,-.898215824685518,-.584846546513270,
144 7 -.226843756241524,0.226843756241524,0.584846546513270,
145 7 0.898215824685518,1. ,0. ,
146 7 0. ,
147 8 -1. ,-.878478166955581,-.661099443664978,
148 8 -.354483526205989,0. ,0.354483526205989,
149 8 0.661099443664978,0.878478166955581,1. ,
150 8 0. ,
151 9 -1. ,-.936320479015252,-.735741735638020,
152 9 -.491001129763160,-.157505717044458,0.157505717044458,
153 9 0.491001129763160,0.735741735638020,0.936320479015252,
154 9 1. /
155C=======================================================================
156 noda_dt = nodadt
157c
158 ! Recover non-local parameters
159 l2 = nloc_dmg%LEN(imat)**2 !< Squared internal length
160 xi = nloc_dmg%DAMP(imat) !< Damping parameter
161 l_nloc = nloc_dmg%L_NLOC !< Size of non-local tables
162 zeta = nloc_dmg%DENS(imat) !< Density parameter
163 sspnl = nloc_dmg%SSPNL(imat) !< Non-local "Sound speed"
164 le_max = nloc_dmg%LE_MAX(imat) !< Maximal length of convergence
165c
166 ! Allocation of variable table
167 ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
168 . btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),
169 . f1(nel,nlay),f2(nel,nlay),f3(nel,nlay),f4(nel,nlay),
170 . f5(nel,nlay),f6(nel,nlay),f7(nel,nlay),f8(nel,nlay),
171 . pos1(nel),pos2(nel),pos3(nel),pos4(nel),pos5(nel),
172 . pos6(nel),pos7(nel),pos8(nel),lc(nel),thk(nel),lthk(nel))
173c
174 ! Degenerated geometry index tables
175 nindx6 = 0
176 indx6(1:nel) = 0
177 ntn6 = six*six
178 nindx8 = 0
179 indx8(1:nel) = 0
180 ntn8 = eight*eight
181c
182 ! Local variables initialization
183 lc(1:nel) = zero
184 ! For nodal timestep
185 IF (noda_dt > 0) THEN
186 ! Non-local nodal stifness
187 ALLOCATE(sti1(nel,nlay),sti2(nel,nlay),sti3(nel,nlay),sti4(nel,nlay),
188 . sti5(nel,nlay),sti6(nel,nlay),sti7(nel,nlay),sti8(nel,nlay))
189 ! Non-local mass
190 mass => nloc_dmg%MASS(1:l_nloc)
191 ! Initial non-local mass
192 mass0 => nloc_dmg%MASS0(1:l_nloc)
193 ELSE
194 NULLIFY(mass,mass0)
195 ALLOCATE(sti1(1,1),sti2(1,1),sti3(1,1),sti4(1,1),
196 . sti5(1,1),sti6(1,1),sti7(1,1),sti8(1,1))
197 ENDIF
198 ! Current timestep non-local velocities
199 vnl => nloc_dmg%VNL(1:l_nloc)
200 ! Previous timestep non-local velocities
201 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
202 ! Non-local displacements
203 unl => nloc_dmg%UNL(1:l_nloc)
204c
205 !--------------------------------------------------------------------------------
206 ! Computation of the position of the non-local d.o.fs and the BtB matrix product
207 !--------------------------------------------------------------------------------
208 ! Loop over elements
209 DO i=1,nel
210 ! Computation of the product BtxB
211 btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
212 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
213 btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
214 btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
215 btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
216 btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
217 btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
218 btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
219 btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
220 btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
221 ENDDO
222c
223 ! Loop over elements
224 DO i=1,nel
225c
226 ! Degenerated brick elements
227 IF (nlocs%NL_ISOLNOD(i) == 6) THEN
228c
229 ! Indexing degenerated penta bricks
230 nindx6 = nindx6 + 1
231 indx6(nindx6) = i
232c
233 ! Recovering the nodes of the brick element
234 n1 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(1,i))
235 n2 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(2,i))
236 n3 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(3,i))
237 n4 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(4,i))
238 n5 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(5,i))
239 n6 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(6,i))
240c
241 ! Recovering the positions of the first d.o.fs of each nodes
242 pos1(i) = nloc_dmg%POSI(n1)
243 pos2(i) = nloc_dmg%POSI(n2)
244 pos3(i) = nloc_dmg%POSI(n3)
245 pos4(i) = nloc_dmg%POSI(n4)
246 pos5(i) = nloc_dmg%POSI(n5)
247 pos6(i) = nloc_dmg%POSI(n6)
248c
249 ! Classic hexa brick elements
250 ELSEIF (nlocs%NL_ISOLNOD(i) == 8) THEN
251c
252 ! Indexing classic hexa bricks
253 nindx8 = nindx8 + 1
254 indx8(nindx8) = i
255c
256 ! Recovering the nodes of the brick element
257 n1 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(1,i))
258 n2 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(2,i))
259 n3 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(3,i))
260 n4 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(4,i))
261 n5 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(5,i))
262 n6 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(6,i))
263 n7 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(7,i))
264 n8 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(8,i))
265c
266 ! Recovering the positions of the first d.o.fs of each nodes
267 pos1(i) = nloc_dmg%POSI(n1)
268 pos2(i) = nloc_dmg%POSI(n2)
269 pos3(i) = nloc_dmg%POSI(n3)
270 pos4(i) = nloc_dmg%POSI(n4)
271 pos5(i) = nloc_dmg%POSI(n5)
272 pos6(i) = nloc_dmg%POSI(n6)
273 pos7(i) = nloc_dmg%POSI(n7)
274 pos8(i) = nloc_dmg%POSI(n8)
275c
276 ENDIF
277 ENDDO
278c
279 !-----------------------------------------------------------------------
280 ! Pre-treatment non-local regularization in the thickshell thickness
281 !-----------------------------------------------------------------------
282 IF ((l2>zero).AND.(nlay > 1)) THEN
283c
284 ! Compute thickshell thickness
285 DO i = 1,nel
286 thk(i) = vol(i)/area(i)
287 lthk(i) = (zs(nlay+1,nlay)-zs(nlay,nlay))*thk(i)*half
288 ENDDO
289c
290 ! Allocation of the velocities predictor
291 nddl = nlay
292 IF (noda_dt > 0) THEN
293 ALLOCATE(stifnlth(nel,nddl+1))
294 ALLOCATE(dtn(nel,nddl+1))
295 ELSE
296 ALLOCATE(dtn(1,1))
297 dtn = ep20
298 ALLOCATE(stifnlth(1,1))
299 stifnlth = ep20
300 ENDIF
301 ndnod = nddl+1
302c
303 ! Pointing the non-local values in the thickness of the corresponding element
304 massth => bufnlts%MASSTH(1:nel,1:ndnod)
305 unlth => bufnlts%UNLTH(1:nel ,1:ndnod)
306 vnlth => bufnlts%VNLTH(1:nel ,1:ndnod)
307 fnlth => bufnlts%FNLTH(1:nel ,1:ndnod)
308c
309 DO k = 1,ndnod
310 DO i = 1,nel
311 ! Resetting non-local forces
312 fnlth(i,k) = zero
313 ! Resetting non-local nodal stiffness
314 IF (noda_dt > 0) THEN
315 stifnlth(i,k) = em20
316 ENDIF
317 ENDDO
318 ENDDO
319c
320 ! Computation of non-local forces in the shell thickness
321 DO k = 1, nddl
322c
323 ! Computation of shape functions value
324 nth1 = (as(k,nddl) - zs(k+1,nddl)) /
325 . (zs(k,nddl) - zs(k+1,nddl))
326 nth2 = (as(k,nddl) - zs(k,nddl)) /
327 . (zs(k+1,nddl) - zs(k,nddl))
328c
329 ! Loop over elements
330 DO i = 1,nel
331c
332 ! Computation of B-matrix values
333 bth1 = (one/(zs(k,nddl) - zs(k+1,nddl)))*(two/thk(i))
334 bth2 = (one/(zs(k+1,nddl) - zs(k,nddl)))*(two/thk(i))
335c
336 ! Computation of the non-local K matrix
337 k1 = l2*(bth1**2) + nth1**2
338 k12 = l2*(bth1*bth2)+ (nth1*nth2)
339 k2 = l2*(bth2**2) + nth2**2
340c
341 ! Computation of the non-local forces
342 fnlth(i,k) = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1)
343 . + xi*((nth1**2)*vnlth(i,k)
344 . + (nth1*nth2)*vnlth(i,k+1))
345 . - (nth1*var_reg(i,k)))*half*ws(k,nddl)*vol(i)
346 fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i,k) + k2*unlth(i,k+1)
347 . + xi*(nth1*nth2*vnlth(i,k)
348 . + (nth2**2)*vnlth(i,k+1))
349 . - nth2*var_reg(i,k))*half*ws(k,nddl)*vol(i)
350c
351 ! Computation of non-local nodal stiffness
352 IF (noda_dt > 0) THEN
353 stifnlth(i,k) = stifnlth(i,k) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)
354 stifnlth(i,k+1) = stifnlth(i,k+1) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)
355 ENDIF
356c
357 ENDDO
358 ENDDO
359c
360 ! Updating non-local mass with /DT/NODA
361 IF (noda_dt > 0) THEN
362C
363 ! Initial computation of the nodal timestep
364 dtnod = ep20
365 DO k = 1,ndnod
366 DO i = 1,nel
367 dtn(i,k) = dtfac1(11)*cdamp*sqrt(two*massth(i,k)/max(stifnlth(i,k),em20))
368 dtnod = min(dtn(i,k),dtnod)
369 ENDDO
370 ENDDO
371C
372 ! /DT/NODA/CSTX - Constant timestep with added mass
373 IF ((idtmin(11)==3).OR.(idtmin(11)==4).OR.(idtmin(11)==8)) THEN
374 ! Added mass computation if necessary
375 IF (dtnod < dtmin1(11)*(sqrt(csta))) THEN
376 DO k = 1,ndnod
377 DO i = 1,nel
378 IF (dtn(i,k) < dtmin1(11)) THEN
379 dt2p = dtmin1(11)/(dtfac1(11)*cdamp)
380 massth(i,k) = max(massth(i,k),csta*half*stifnlth(i,k)*dt2p*dt2p*onep00001)
381 ENDIF
382 ENDDO
383 ENDDO
384 ENDIF
385 dtnod = dtmin1(11)*(sqrt(csta))
386 ENDIF
387C
388 ! Classical nodal timestep check
389 IF (dtnod < dt2t) THEN
390 dt2t = min(dt2t,dtnod)
391 ENDIF
392 ENDIF
393c
394 DO k = 1,ndnod
395 DO i = 1,nel
396 ! Updating the non-local in-thickness velocities
397 vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt12
398 ENDDO
399 ENDDO
400c
401 DO k = 1,ndnod
402 DO i = 1,nel
403 ! Computing the non-local in-thickness cumulated values
404 unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt1
405 ENDDO
406 ENDDO
407c
408 ! Transfert at the integration point
409 DO k = 1, nddl
410 !Computation of shape functions value
411 nth1 = (as(k,nddl) - zs(k+1,nddl))/
412 . (zs(k,nddl) - zs(k+1,nddl))
413 nth2 = (as(k,nddl) - zs(k,nddl))/
414 . (zs(k+1,nddl) - zs(k,nddl))
415 ! Loop over elements
416 DO i = 1,nel
417 !Integration points non-local variables
418 var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i,k+1)
419 ENDDO
420 ENDDO
421 ENDIF
422c
423 !-----------------------------------------------------------------------
424 ! Computation of non-local forces
425 !-----------------------------------------------------------------------
426 ! Loop over layers
427 DO k = 1,nlay
428c
429 ! Loop over classic hexa bricks
430#include "vectorize.inc"
431 DO ii = 1, nindx8
432c
433 ! Number of the element
434 i = indx8(ii)
435c
436 ! If the element is not broken, normal computation
437 IF (off(i) /= zero) THEN
438c
439 ! Computing the product NtN*UNL
440 ntn_unl = (unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) + unl(pos4(i)+k-1)
441 . + unl(pos5(i)+k-1) + unl(pos6(i)+k-1) + unl(pos7(i)+k-1) + unl(pos8(i)+k-1)) / ntn8
442c
443 ! Computing the product NtN*VNL
444 ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) + vnl(pos4(i)+k-1)
445 . + vnl(pos5(i)+k-1) + vnl(pos6(i)+k-1) + vnl(pos7(i)+k-1) + vnl(pos8(i)+k-1)) / ntn8
446 IF (noda_dt > 0) THEN
447 ntn_vnl = min(sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1)),
448 . sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1)),
449 . sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1)),
450 . sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1)),
451 . sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1)),
452 . sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1)),
453 . sqrt(mass(pos7(i)+k-1)/mass0(pos7(i)+k-1)),
454 . sqrt(mass(pos8(i)+k-1)/mass0(pos8(i)+k-1)))*ntn_vnl
455 ENDIF
456c
457 ! Computation of the product LEN**2 * BtxB
458 b1 = l2 * vol(i) * ws(k,nlay) *half * ( btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
459 . + btb13(i)*unl(pos3(i)+k-1) + btb14(i)*unl(pos4(i)+k-1) - btb13(i)*unl(pos5(i)+k-1)
460 . - btb14(i)*unl(pos6(i)+k-1) - btb11(i)*unl(pos7(i)+k-1) - btb12(i)*unl(pos8(i)+k-1))
461c
462 b2 = l2 * vol(i) * ws(k,nlay) *half * ( btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
463 . + btb23(i)*unl(pos3(i)+k-1) + btb24(i)*unl(pos4(i)+k-1) - btb23(i)*unl(pos5(i)+k-1)
464 . - btb24(i)*unl(pos6(i)+k-1) - btb12(i)*unl(pos7(i)+k-1) - btb22(i)*unl(pos8(i)+k-1))
465c
466 b3 = l2 * vol(i) * ws(k,nlay) *half * ( btb13(i)*unl(pos1(i)+k-1) + btb23(i)*unl(pos2(i)+k-1)
467 . + btb33(i)*unl(pos3(i)+k-1) + btb34(i)*unl(pos4(i)+k-1) - btb33(i)*unl(pos5(i)+k-1)
468 . - btb34(i)*unl(pos6(i)+k-1) - btb13(i)*unl(pos7(i)+k-1) - btb23(i)*unl(pos8(i)+k-1))
469c
470 b4 = l2 * vol(i) * ws(k,nlay) *half * ( btb14(i)*unl(pos1(i)+k-1) + btb24(i)*unl(pos2(i)+k-1)
471 . + btb34(i)*unl(pos3(i)+k-1) + btb44(i)*unl(pos4(i)+k-1) - btb34(i)*unl(pos5(i)+k-1)
472 . - btb44(i)*unl(pos6(i)+k-1) - btb14(i)*unl(pos7(i)+k-1) - btb24(i)*unl(pos8(i)+k-1))
473c
474 b5 = l2 * vol(i) * ws(k,nlay) *half * ( -btb13(i)*unl(pos1(i)+k-1)- btb23(i)*unl(pos2(i)+k-1)
475 . - btb33(i)*unl(pos3(i)+k-1) - btb34(i)*unl(pos4(i)+k-1) + btb33(i)*unl(pos5(i)+k-1)
476 . + btb34(i)*unl(pos6(i)+k-1) + btb13(i)*unl(pos7(i)+k-1) + btb23(i)*unl(pos8(i)+k-1))
477c
478 b6 = l2 * vol(i) * ws(k,nlay) *half * ( -btb14(i)*unl(pos1(i)+k-1)- btb24(i)*unl(pos2(i)+k-1)
479 . - btb34(i)*unl(pos3(i)+k-1) - btb44(i)*unl(pos4(i)+k-1) + btb34(i)*unl(pos5(i)+k-1)
480 . + btb44(i)*unl(pos6(i)+k-1) + btb14(i)*unl(pos7(i)+k-1) + btb24(i)*unl(pos8(i)+k-1))
481c
482 b7 = l2 * vol(i) * ws(k,nlay) *half * ( -btb11(i)*unl(pos1(i)+k-1)- btb12(i)*unl(pos2(i)+k-1)
483 . - btb13(i)*unl(pos3(i)+k-1) - btb14(i)*unl(pos4(i)+k-1) + btb13(i)*unl(pos5(i)+k-1)
484 . + btb14(i)*unl(pos6(i)+k-1) + btb11(i)*unl(pos7(i)+k-1) + btb12(i)*unl(pos8(i)+k-1))
485c
486 b8 = l2 * vol(i) * ws(k,nlay) *half * ( -btb12(i)*unl(pos1(i)+k-1)- btb22(i)*unl(pos2(i)+k-1)
487 . - btb23(i)*unl(pos3(i)+k-1) - btb24(i)*unl(pos4(i)+k-1) + btb23(i)*unl(pos5(i)+k-1)
488 . + btb24(i)*unl(pos6(i)+k-1) + btb12(i)*unl(pos7(i)+k-1) + btb22(i)*unl(pos8(i)+k-1))
489c
490 ! Multiplication by the volume of the element (and damping parameter XI)
491 ntn_unl = ntn_unl * vol(i) * ws(k,nlay) * half
492 ntn_vnl = ntn_vnl * xi * vol(i) * ws(k,nlay) * half
493c
494 ! Introducing the internal variable to be regularized
495 ntvar = var_reg(i,k)*one_over_8* vol(i) * ws(k,nlay) * half
496c
497 ! Computing the elementary non-local forces
498 a = ntn_unl + ntn_vnl - ntvar
499 f1(i,k) = a + b1
500 f2(i,k) = a + b2
501 f3(i,k) = a + b3
502 f4(i,k) = a + b4
503 f5(i,k) = a + b5
504 f6(i,k) = a + b6
505 f7(i,k) = a + b7
506 f8(i,k) = a + b8
507c
508 ! Computing nodal equivalent stiffness
509 IF (noda_dt > 0) THEN
510 sti1(i,k) = (abs(l2*btb11(i) + one/ntn8) + abs(l2*btb12(i) + one/ntn8) + abs(l2*btb13(i) + one/ntn8) +
511 . abs(l2*btb14(i) + one/ntn8) + abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb14(i) + one/ntn8) +
512 . abs(-l2*btb11(i) + one/ntn8) + abs(-l2*btb12(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
513 sti2(i,k) = (abs(l2*btb12(i) + one/ntn8) + abs(l2*btb22(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) +
514 . abs(l2*btb24(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8) +
515 . abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb22(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
516 sti3(i,k) = (abs(l2*btb13(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) + abs(l2*btb33(i) + one/ntn8) +
517 . abs(l2*btb34(i) + one/ntn8) + abs(-l2*btb33(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) +
518 . abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
519 sti4(i,k) = (abs(l2*btb14(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) +
520 . abs(l2*btb44(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) + abs(-l2*btb44(i) + one/ntn8) +
521 . abs(-l2*btb14(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
522 sti5(i,k) = (abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) + abs(-l2*btb33(i) + one/ntn8) +
523 . abs(-l2*btb34(i) + one/ntn8) + abs(l2*btb33(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) +
524 . abs(l2*btb13(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
525 sti6(i,k) = (abs(-l2*btb14(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) +
526 . abs(-l2*btb44(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) + abs(l2*btb44(i) + one/ntn8) +
527 . abs(l2*btb14(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
528 sti7(i,k) = (abs(-l2*btb11(i) + one/ntn8) + abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb13(i) + one/ntn8) +
529 . abs(-l2*btb14(i) + one/ntn8) + abs(l2*btb13(i) + one/ntn8) + abs(l2*btb14(i) + one/ntn8) +
530 . abs(l2*btb11(i) + one/ntn8) + abs(l2*btb12(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
531 sti8(i,k) = (abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb22(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) +
532 . abs(-l2*btb24(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8) +
533 . abs(l2*btb12(i) + one/ntn8) + abs(l2*btb22(i) + one/ntn8))*vol(i)*ws(k,nlay)*half
534 ENDIF
535c
536 ! If the element is broken, the non-local wave is absorbed
537 ELSE
538c
539 ! Initial element characteristic length
540 lc(i) = (vol0(i)*ws(k,nlay)*half)**third
541c
542 IF (noda_dt > 0) THEN
543 ! Non-local absorbing forces
544 f1(i,k) = sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*half*
545 . (vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(three/four)*(lc(i)**2)
546 f2(i,k) = sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*half*
547 . (vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(three/four)*(lc(i)**2)
548 f3(i,k) = sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*zeta*sspnl*half*
549 . (vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(three/four)*(lc(i)**2)
550 f4(i,k) = sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*half*
551 . (vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(three/four)*(lc(i)**2)
552 f5(i,k) = sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1))*zeta*sspnl*half*
553 . (vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(three/four)*(lc(i)**2)
554 f6(i,k) = sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1))*zeta*sspnl*half*
555 . (vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(three/four)*(lc(i)**2)
556 f7(i,k) = sqrt(mass(pos7(i)+k-1)/mass0(pos7(i)+k-1))*zeta*sspnl*half*
557 . (vnl(pos7(i)+k-1)+vnl0(pos7(i)+k-1))*(three/four)*(lc(i)**2)
558 f8(i,k) = sqrt(mass(pos8(i)+k-1)/mass0(pos8(i)+k-1))*zeta*sspnl*half*
559 . (vnl(pos8(i)+k-1)+vnl0(pos8(i)+k-1))*(three/four)*(lc(i)**2)
560 ! Computing nodal equivalent stiffness
561 sti1(i,k) = em20
562 sti2(i,k) = em20
563 sti3(i,k) = em20
564 sti4(i,k) = em20
565 sti5(i,k) = em20
566 sti6(i,k) = em20
567 sti7(i,k) = em20
568 sti8(i,k) = em20
569 ELSE
570 ! Non-local absorbing forces
571 f1(i,k) = zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(three/four)*(lc(i)**2)
572 f2(i,k) = zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(three/four)*(lc(i)**2)
573 f3(i,k) = zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(three/four)*(lc(i)**2)
574 f4(i,k) = zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(three/four)*(lc(i)**2)
575 f5(i,k) = zeta*sspnl*half*(vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(three/four)*(lc(i)**2)
576 f6(i,k) = zeta*sspnl*half*(vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(three/four)*(lc(i)**2)
577 f7(i,k) = zeta*sspnl*half*(vnl(pos7(i)+k-1)+vnl0(pos7(i)+k-1))*(three/four)*(lc(i)**2)
578 f8(i,k) = zeta*sspnl*half*(vnl(pos8(i)+k-1)+vnl0(pos8(i)+k-1))*(three/four)*(lc(i)**2)
579 ENDIF
580 ENDIF
581 ENDDO
582c
583 ! Loop over penta degenerated bricks
584#include "vectorize.inc"
585 DO ii = 1, nindx6
586c
587 ! Number of the element
588 i = indx6(ii)
589c
590 ! If the element is not broken, normal computation
591 IF (off(i) /= zero) THEN
592c
593 ! Computing the product NtN*UNL
594 ntn_unl = (unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) +
595 . unl(pos4(i)+k-1) + unl(pos5(i)+k-1) + unl(pos6(i)+k-1)) / ntn6
596c
597 ! Computing the product NtN*VNL
598 ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) +
599 . vnl(pos4(i)+k-1) + vnl(pos5(i)+k-1) + vnl(pos6(i)+k-1)) / ntn6
600 IF (noda_dt > 0) THEN
601 ntn_vnl = (sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*vnl(pos1(i)+k-1) +
602 . sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*vnl(pos2(i)+k-1) +
603 . sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*vnl(pos3(i)+k-1) +
604 . sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*vnl(pos4(i)+k-1) +
605 . sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1))*vnl(pos5(i)+k-1) +
606 . sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1))*vnl(pos6(i)+k-1)) / ntn6
607 ENDIF
608c
609 ! Computation of the product LEN**2 * BtxB
610 b1 = l2 * vol(i) * ws(k,nlay) *half * ( btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
611 . + btb13(i)*unl(pos3(i)+k-1) - btb13(i)*unl(pos4(i)+k-1) - btb12(i)*unl(pos5(i)+k-1)
612 . - btb11(i)*unl(pos6(i)+k-1) )
613c
614 b2 = l2 * vol(i) * ws(k,nlay) *half * ( btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
615 . + btb23(i)*unl(pos3(i)+k-1) - btb23(i)*unl(pos4(i)+k-1) - btb22(i)*unl(pos5(i)+k-1)
616 . - btb12(i)*unl(pos6(i)+k-1) )
617c
618 b3 = l2 * vol(i) * ws(k,nlay) *half * ( btb13(i)*unl(pos1(i)+k-1) + btb23(i)*unl(pos2(i)+k-1)
619 . + btb33(i)*unl(pos3(i)+k-1) - btb33(i)*unl(pos4(i)+k-1) - btb23(i)*unl(pos5(i)+k-1)
620 . - btb13(i)*unl(pos6(i)+k-1) )
621c
622 b4 = l2 * vol(i) * ws(k,nlay) *half * (-btb13(i)*unl(pos1(i)+k-1) - btb23(i)*unl(pos2(i)+k-1)
623 . - btb33(i)*unl(pos3(i)+k-1) + btb33(i)*unl(pos4(i)+k-1) + btb23(i)*unl(pos5(i)+k-1)
624 . + btb13(i)*unl(pos6(i)+k-1) )
625c
626 b5 = l2 * vol(i) * ws(k,nlay) *half * ( -btb12(i)*unl(pos1(i)+k-1)- btb22(i)*unl(pos2(i)+k-1)
627 . - btb23(i)*unl(pos3(i)+k-1) + btb23(i)*unl(pos4(i)+k-1) + btb22(i)*unl(pos5(i)+k-1)
628 . + btb12(i)*unl(pos6(i)+k-1) )
629c
630 b6 = l2 * vol(i) * ws(k,nlay) *half * ( -btb11(i)*unl(pos1(i)+k-1)- btb12(i)*unl(pos2(i)+k-1)
631 . - btb13(i)*unl(pos3(i)+k-1) + btb13(i)*unl(pos4(i)+k-1) + btb12(i)*unl(pos5(i)+k-1)
632 . + btb11(i)*unl(pos6(i)+k-1) )
633c
634 ! Multiplication by the volume of the element (and the damping parameter XI)
635 ntn_unl = ntn_unl * vol(i) * ws(k,nlay) * half
636 ntn_vnl = ntn_vnl * xi * vol(i) * ws(k,nlay) * half
637c
638 ! Introducing the internal variable to be regularized
639 ntvar = var_reg(i,k)*one_over_6* vol(i) * ws(k,nlay) * half
640c
641 ! Computing the elementary non-local forces
642 a = ntn_unl + ntn_vnl - ntvar
643 f1(i,k) = a + b1
644 f2(i,k) = a + b2
645 f3(i,k) = a + b3
646 f4(i,k) = a + b4
647 f5(i,k) = a + b5
648 f6(i,k) = a + b6
649c
650 ! Computing nodal equivalent stiffness
651 IF (noda_dt > 0) THEN
652 sti1(i,k) = (abs(l2*btb11(i) + one/ntn6) + abs(l2*btb12(i) + one/ntn6) +
653 . abs(l2*btb13(i) + one/ntn6) + abs(-l2*btb13(i) + one/ntn6) +
654 . abs(-l2*btb12(i) + one/ntn6) + abs(-l2*btb11(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
655 sti2(i,k) = (abs(l2*btb12(i) + one/ntn6) + abs(l2*btb22(i) + one/ntn6) +
656 . abs(l2*btb23(i) + one/ntn6) + abs(-l2*btb23(i) + one/ntn6) +
657 . abs(-l2*btb22(i) + one/ntn6) + abs(-l2*btb12(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
658 sti3(i,k) = (abs(l2*btb13(i) + one/ntn6) + abs(l2*btb23(i) + one/ntn6) +
659 . abs(l2*btb33(i) + one/ntn6) + abs(-l2*btb33(i) + one/ntn6) +
660 . abs(-l2*btb23(i) + one/ntn6) + abs(-l2*btb13(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
661 sti4(i,k) = (abs(-l2*btb13(i) + one/ntn6) + abs(-l2*btb23(i) + one/ntn6) +
662 . abs(-l2*btb33(i) + one/ntn6) + abs(l2*btb33(i) + one/ntn6) +
663 . abs(l2*btb23(i) + one/ntn6) + abs(l2*btb13(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
664 sti5(i,k) = (abs(-l2*btb12(i) + one/ntn6) + abs(-l2*btb22(i) + one/ntn6) +
665 . abs(-l2*btb23(i) + one/ntn6) + abs(l2*btb23(i) + one/ntn6) +
666 . abs(l2*btb22(i) + one/ntn6) + abs(l2*btb12(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
667 sti6(i,k) = (abs(-l2*btb11(i) + one/ntn6) + abs(-l2*btb12(i) + one/ntn6) +
668 . abs(-l2*btb13(i) + one/ntn6) + abs(l2*btb13(i) + one/ntn6) +
669 . abs(l2*btb12(i) + one/ntn6) + abs(l2*btb11(i) + one/ntn6))*vol(i)*ws(k,nlay)*half
670 ENDIF
671c
672 ELSE
673c
674 ! Initial element characteristic length
675 lc(i) = (vol0(i)*ws(k,nlay)*half)**third
676c
677 IF (noda_dt > 0) THEN
678 ! Non-local absorbing forces
679 f1(i,k) = sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*half*
680 . (vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
681 f2(i,k) = sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*half*
682 . (vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
683 f3(i,k) = sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*zeta*sspnl*half*
684 . (vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
685 f4(i,k) = sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*half*
686 . (vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
687 f5(i,k) = sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1))*zeta*sspnl*half*
688 . (vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
689 f6(i,k) = sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1))*zeta*sspnl*half*
690 . (vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)
691 ! Computing nodal equivalent stiffness
692 sti1(i,k) = em20
693 sti2(i,k) = em20
694 sti3(i,k) = em20
695 sti4(i,k) = em20
696 sti5(i,k) = em20
697 sti6(i,k) = em20
698 ELSE
699 ! Non-local absorbing forces
700 f1(i,k) = zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
701 f2(i,k) = zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
702 f3(i,k) = zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
703 f4(i,k) = zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
704 f5(i,k) = zeta*sspnl*half*(vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
705 f6(i,k) = zeta*sspnl*half*(vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)
706 ENDIF
707 ENDIF
708 ENDDO
709 ENDDO
710c
711 !-----------------------------------------------------------------------
712 ! Assembling of the non-local forces
713 !-----------------------------------------------------------------------
714c
715 ! If PARITH/OFF
716 IF (iparit == 0) THEN
717 ! Recovering non-local internal forces
718 fnl => nloc_dmg%FNL(1:l_nloc,itask+1) ! Non-local forces
719c IF (NODA_DT > 0) STIFNL => NLOC_DMG%STIFNL(1:L_NLOC,ITASK+1) ! Non-local equivalent nodal stiffness
720c
721 ! Loop over non-local degrees of freedom
722 DO k=1,nlay
723 ! Loop over classic hexa brick elements
724#include "vectorize.inc"
725 DO ii=1,nindx8
726 i = indx8(ii)
727 ! Assembling the forces in the classic way
728 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
729 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
730 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
731 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
732 fnl(pos5(i)+k-1) = fnl(pos5(i)+k-1) - f5(i,k)
733 fnl(pos6(i)+k-1) = fnl(pos6(i)+k-1) - f6(i,k)
734 fnl(pos7(i)+k-1) = fnl(pos7(i)+k-1) - f7(i,k)
735 fnl(pos8(i)+k-1) = fnl(pos8(i)+k-1) - f8(i,k)
736 IF (noda_dt > 0) THEN
737 ! Spectral radius of stiffness matrix
738 maxstif = max(sti1(i,k),sti2(i,k),sti3(i,k),sti4(i,k),
739 . sti5(i,k),sti6(i,k),sti7(i,k),sti8(i,k))
740 ! Computing nodal stiffness
741 nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) + maxstif
742 nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) + maxstif
743 nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) + maxstif
744 nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) + maxstif
745 nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) + maxstif
746 nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) + maxstif
747 nloc_dmg%STIFNL(pos7(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos7(i)+k-1,itask+1) + maxstif
748 nloc_dmg%STIFNL(pos8(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos8(i)+k-1,itask+1) + maxstif
749 ENDIF
750 ENDDO
751c
752 ! Loop over classic hexa brick elements
753#include "vectorize.inc"
754 DO ii=1,nindx6
755 i = indx6(ii)
756 ! Assembling the forces in the classic way
757 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
758 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
759 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
760 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
761 fnl(pos5(i)+k-1) = fnl(pos5(i)+k-1) - f5(i,k)
762 fnl(pos6(i)+k-1) = fnl(pos6(i)+k-1) - f6(i,k)
763 IF (noda_dt > 0) THEN
764 ! Spectral radius of stiffness matrix
765 maxstif = max(sti1(i,k),sti2(i,k),sti3(i,k),
766 . sti4(i,k),sti5(i,k),sti6(i,k))
767 ! Computing nodal stiffness
768 nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) + maxstif
769 nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) + maxstif
770 nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) + maxstif
771 nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) + maxstif
772 nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) + maxstif
773 nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) + maxstif
774 ENDIF
775 ENDDO
776 ENDDO
777c
778 ! If PARITH/ON
779 ELSE
780 ! Loop over additional d.o.fs
781 DO j = 1,nlay
782c
783 ! Loop over classic hexa bricks
784 DO l = 1,nindx8
785 i = indx8(l)
786 ii = i + nft
787c
788 ! Spectral radius of stiffness matrix
789 IF (noda_dt > 0) THEN
790 maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j),
791 . sti5(i,j),sti6(i,j),sti7(i,j),sti8(i,j))
792 ENDIF
793c
794 k = nloc_dmg%IADS(1,ii)
795 nloc_dmg%FSKY(k,j) = -f1(i,j)
796 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
797c
798 k = nloc_dmg%IADS(2,ii)
799 nloc_dmg%FSKY(k,j) = -f2(i,j)
800 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
801c
802 k = nloc_dmg%IADS(3,ii)
803 nloc_dmg%FSKY(k,j) = -f3(i,j)
804 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
805c
806 k = nloc_dmg%IADS(4,ii)
807 nloc_dmg%FSKY(k,j) = -f4(i,j)
808 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
809c
810 k = nloc_dmg%IADS(5,ii)
811 nloc_dmg%FSKY(k,j) = -f5(i,j)
812 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
813c
814 k = nloc_dmg%IADS(6,ii)
815 nloc_dmg%FSKY(k,j) = -f6(i,j)
816 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
817c
818 k = nloc_dmg%IADS(7,ii)
819 nloc_dmg%FSKY(k,j) = -f7(i,j)
820 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
821c
822 k = nloc_dmg%IADS(8,ii)
823 nloc_dmg%FSKY(k,j) = -f8(i,j)
824 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
825c
826 ENDDO
827c
828 ! Loop over penta degenerated bricks
829 DO l = 1,nindx6
830 i = indx6(l)
831 ii = i + nft
832c
833 ! Spectral radius of stiffness matrix
834 IF (noda_dt > 0) THEN
835 maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),
836 . sti4(i,j),sti5(i,j),sti6(i,j))
837 ENDIF
838c
839 k = nloc_dmg%IADS(1,ii)
840 nloc_dmg%FSKY(k,j) = -f1(i,j)
841 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
842c
843 k = nloc_dmg%IADS(2,ii)
844 nloc_dmg%FSKY(k,j) = -f2(i,j)
845 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
846c
847 k = nloc_dmg%IADS(3,ii)
848 nloc_dmg%FSKY(k,j) = -f3(i,j)
849 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
850c
851 k = nloc_dmg%IADS(5,ii)
852 nloc_dmg%FSKY(k,j) = -f4(i,j)
853 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
854c
855 k = nloc_dmg%IADS(6,ii)
856 nloc_dmg%FSKY(k,j) = -f5(i,j)
857 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
858c
859 k = nloc_dmg%IADS(7,ii)
860 nloc_dmg%FSKY(k,j) = -f6(i,j)
861 IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
862c
863 ENDDO
864 ENDDO
865 ENDIF
866c
867 !-----------------------------------------------------------------------
868 ! Computing non-local timestep
869 !-----------------------------------------------------------------------
870 IF (noda_dt == 0) THEN
871 DO i = 1,nel
872 ! If the element is not broken, normal computation
873 IF (off(i)/=zero) THEN
874 ! Non-local critical time-step in the plane
875 dtnl = (two*(min((vol(i))**third,le_max))*sqrt(three*zeta))/
876 . sqrt(twelve*l2 + (min((vol(i))**third,le_max))**2)
877 ! Non-local critical time-step in the thickness
878 IF ((l2>zero).AND.(nlay > 1)) THEN
879 dtnl_th = (two*(min(lthk(i),le_max))*sqrt(three*zeta))/
880 . sqrt(twelve*l2 + (min(lthk(i),le_max))**2)
881 ELSE
882 dtnl_th = ep20
883 ENDIF
884 ! Keeping the minimal value
885 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl_th,dtfac1(1)*cdamp*dtnl)
886 ENDIF
887 ENDDO
888 ENDIF
889c
890 ! Deallocation of tables
891 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
892 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
893 IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
894 IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
895 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
896 IF (ALLOCATED(btb23)) DEALLOCATE(btb23)
897 IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
898 IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
899 IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
900 IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
901 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
902 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
903 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
904 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
905 IF (ALLOCATED(pos5)) DEALLOCATE(pos5)
906 IF (ALLOCATED(pos6)) DEALLOCATE(pos6)
907 IF (ALLOCATED(pos7)) DEALLOCATE(pos7)
908 IF (ALLOCATED(pos8)) DEALLOCATE(pos8)
909 IF (ALLOCATED(f1)) DEALLOCATE(f1)
910 IF (ALLOCATED(f2)) DEALLOCATE(f2)
911 IF (ALLOCATED(f3)) DEALLOCATE(f3)
912 IF (ALLOCATED(f4)) DEALLOCATE(f4)
913 IF (ALLOCATED(f5)) DEALLOCATE(f5)
914 IF (ALLOCATED(f6)) DEALLOCATE(f6)
915 IF (ALLOCATED(f7)) DEALLOCATE(f7)
916 IF (ALLOCATED(f8)) DEALLOCATE(f8)
917 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
918 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
919 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
920 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
921 IF (ALLOCATED(sti5)) DEALLOCATE(sti5)
922 IF (ALLOCATED(sti6)) DEALLOCATE(sti6)
923 IF (ALLOCATED(sti7)) DEALLOCATE(sti7)
924 IF (ALLOCATED(sti8)) DEALLOCATE(sti8)
925 IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
926 IF (ALLOCATED(dtn)) DEALLOCATE(dtn)
927 IF (ALLOCATED(lc)) DEALLOCATE(lc)
928 IF (ALLOCATED(thk)) DEALLOCATE(thk)
929 IF (ALLOCATED(lthk)) DEALLOCATE(lthk)
930c-----------
931 END
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
scfint_reg
subroutine scfint_reg(nloc_dmg, var_reg, nel, off, vol, nlocs, area, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, imat, itask, dt2t, vol0, nft, nlay, ws, as, bufnlts)
Definition scfint_reg.F:40
nlocal_reg_mod::nlocal_str_
Definition nlocal_reg_mod.F:124