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