OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sfint_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!|| sfint_reg ../engine/source/elements/solid/solide/sfint_reg.F
25!||--- called by ------------------------------------------------------
26!|| sforc3 ../engine/source/elements/solid/solide/sforc3.F
27!|| szforc3 ../engine/source/elements/solid/solidez/szforc3.F
28!||--- uses -----------------------------------------------------
29!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
30!|| nlocal_reg_mod ../common_source/modules/nlocal_reg_mod.F
31!||====================================================================
32 SUBROUTINE sfint_reg(
33 1 NLOC_DMG ,VAR_REG ,NEL ,OFF ,
34 2 VOL ,NLOCS ,VOL0 ,
35 3 PX1 ,PX2 ,PX3 ,PX4 ,
36 4 PY1 ,PY2 ,PY3 ,PY4 ,
37 5 PZ1 ,PZ2 ,PZ3 ,PZ4 ,
38 6 IMAT ,ITASK ,DT2T ,NFT )
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE nlocal_reg_mod
43 USE elbufdef_mod
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "parit_c.inc"
52#include "scr02_c.inc"
53#include "scr18_c.inc"
54C-----------------------------------------------
55C D u m m y A r g u m e n t s
56C-----------------------------------------------
57 INTEGER, INTENT(IN) :: NFT !< Address of the first element of the group
58 INTEGER, INTENT(IN) :: NEL !< Number of elements in the group
59 INTEGER, INTENT(IN) :: IMAT !< Material internal Id
60 INTEGER, INTENT(IN) :: ITASK !< Number of the thread
61 my_real, INTENT(INOUT) :: DT2T !< Element critical timestep
62 my_real, DIMENSION(NEL), INTENT(IN) :: VOL !< Current volume of the element
63 my_real, DIMENSION(NEL), INTENT(IN) :: off !< Flag for element deletion
64 my_real, DIMENSION(NEL), INTENT(IN) :: var_reg !< Variable to regularize
65 my_real, DIMENSION(NEL), INTENT(IN) :: vol0 !< Initial volume of the element
66 my_real, DIMENSION(NEL), INTENT(IN) :: px1 !< Derivative of the first node shape function along x
67 my_real, DIMENSION(NEL), INTENT(IN) :: px2 !< Derivative of the second node shape function along x
68 my_real, DIMENSION(NEL), INTENT(IN) :: px3 !< Derivative of the third node shape function along x
69 my_real, DIMENSION(NEL), INTENT(IN) :: px4 !< Derivative of the fourth node shape function along x
70 my_real, DIMENSION(NEL), INTENT(IN) :: py1 !< Derivative of the first node shape function along y
71 my_real, DIMENSION(NEL), INTENT(IN) :: py2 !< Derivative of the second node shape function along y
72 my_real, DIMENSION(NEL), INTENT(IN) :: py3 !< Derivative of the third node shape function along y
73 my_real, DIMENSION(NEL), INTENT(IN) :: py4 !< Derivative of the fourth node shape function along y
74 my_real, DIMENSION(NEL), INTENT(IN) :: pz1 !< Derivative of the first node shape function along z
75 my_real, DIMENSION(NEL), INTENT(IN) :: pz2 !< Derivative of the second node shape function along z
76 my_real, DIMENSION(NEL), INTENT(IN) :: pz3 !< Derivative of the third node shape function along z
77 my_real, DIMENSION(NEL), INTENT(IN) :: pz4 !< Derivative of the fourth node shape function along z
78 TYPE(nlocal_str_), INTENT(INOUT), TARGET :: NLOC_DMG !< Non-local data structure
79 TYPE(buf_nlocs_), INTENT(IN) :: NLOCS !< Non-local effective nodes data structure
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER I,II,K,L,N1,N2,N3,N4,N5,N6,N7,N8,L_NLOC,
84 . nindx6,indx6(nel),nindx8,indx8(nel)
85 INTEGER, DIMENSION(:), ALLOCATABLE ::
86 . POS1,POS2,POS3,POS4,POS5,POS6,POS7,POS8
87 my_real
88 . L2,NTN_UNL,NTN_VNL,XI,NTVAR,A,
89 . b1,b2,b3,b4,b5,b6,b7,b8,zeta,sspnl,dtnl,le_max,
90 . maxstif,ntn6,ntn8
91 my_real, DIMENSION(:) ,ALLOCATABLE ::
92 . btb11,btb12,btb13,btb14,btb22,btb23,btb24,
93 . btb33,btb34,btb44,sti1,sti2,sti3,sti4,sti5,
94 . sti6,sti7,sti8,f1,f2,f3,f4,f5,f6,f7,f8,lc
95 my_real, POINTER, DIMENSION(:) ::
96 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
97 ! Coefficient for non-local stability to take into account damping
98 my_real, target :: nothing(1)
99 my_real, PARAMETER :: cdamp = 0.7d0
100c-----------------------------------------------------------------------
101c VAR_REG : variable a regulariser (local cumulated plastic strain)
102c NTVAR = NT * VAR_REG
103C=======================================================================
104c
105 nothing = zero
106 vnl => nothing
107 fnl => nothing
108 unl => nothing
109 stifnl => nothing
110 mass => nothing
111 mass0 => nothing
112 vnl0 => nothing
113
114 ! recover non-local parameters
115 l2 = nloc_dmg%LEN(imat)**2 !< Squared internal length
116 xi = nloc_dmg%DAMP(imat) !< Damping parameter
117 l_nloc = nloc_dmg%L_NLOC !< Size of non-local tables
118 zeta = nloc_dmg%DENS(imat) !< Density parameter
119 sspnl = nloc_dmg%SSPNL(imat) !< Non-local "Sound speed"
120 le_max = nloc_dmg%LE_MAX(imat) !< Maximal length of convergence
121c
122 ! Allocation of variable table
123 ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
124 . btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),pos1(nel),
125 . pos2(nel),pos3(nel),pos4(nel),pos5(nel),pos6(nel),pos7(nel),pos8(nel),
126 . f1(nel),f2(nel),f3(nel),f4(nel),f5(nel),f6(nel),f7(nel),f8(nel),lc(nel))
127c
128 ! Degenerated geometry index tables
129 nindx6 = 0
130 indx6(1:nel) = 0
131 ntn6 = six*six
132 nindx8 = 0
133 indx8(1:nel) = 0
134 ntn8 = eight*eight
135c
136 ! Local variables initialization
137 lc(1:nel) = zero
138 ! For nodal timestep
139 IF (nodadt > 0) THEN
140 ! Non-local nodal stifness
141 ALLOCATE(sti1(nel),sti2(nel),sti3(nel),sti4(nel),sti5(nel),sti6(nel),
142 . sti7(nel),sti8(nel))
143 ! Non-local mass
144 mass => nloc_dmg%MASS(1:l_nloc)
145 ! Initial non-local mass
146 mass0 => nloc_dmg%MASS0(1:l_nloc)
147 ENDIF
148 ! Current timestep non-local velocities
149 vnl => nloc_dmg%VNL(1:l_nloc)
150 ! Previous timestep non-local velocities
151 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
152 ! Non-local displacements
153 unl => nloc_dmg%UNL(1:l_nloc)
154c
155 !--------------------------------------------------------------------------------
156 ! Computation of the position of the non-local d.o.fs and the BtB matrix product
157 !--------------------------------------------------------------------------------
158 ! Loop over elements
159 DO i=1,nel
160 ! Computation of the product BtxB
161 btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
162 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
163 btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
164 btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
165 btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
166 btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
167 btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
168 btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
169 btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
170 btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
171 ENDDO
172c
173 ! Loop over elements
174 DO i=1,nel
175c
176 ! Degenerated brick elements
177 IF (nlocs%NL_ISOLNOD(i) == 6) THEN
178c
179 ! Indexing degenerated penta bricks
180 nindx6 = nindx6 + 1
181 indx6(nindx6) = i
182c
183 ! Recovering the nodes of the brick element
184 n1 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(1,i))
185 n2 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(2,i))
186 n3 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(3,i))
187 n4 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(4,i))
188 n5 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(5,i))
189 n6 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(6,i))
190c
191 ! Recovering the positions of the first d.o.fs of each nodes
192 pos1(i) = nloc_dmg%POSI(n1)
193 pos2(i) = nloc_dmg%POSI(n2)
194 pos3(i) = nloc_dmg%POSI(n3)
195 pos4(i) = nloc_dmg%POSI(n4)
196 pos5(i) = nloc_dmg%POSI(n5)
197 pos6(i) = nloc_dmg%POSI(n6)
198c
199 ! Classic hexa brick elements
200 ELSEIF (nlocs%NL_ISOLNOD(i) == 8) THEN
201c
202 ! Indexing classic hexa bricks
203 nindx8 = nindx8 + 1
204 indx8(nindx8) = i
205c
206 ! Recovering the nodes of the brick element
207 n1 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(1,i))
208 n2 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(2,i))
209 n3 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(3,i))
210 n4 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(4,i))
211 n5 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(5,i))
212 n6 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(6,i))
213 n7 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(7,i))
214 n8 = nloc_dmg%IDXI(nlocs%NL_SOLNOD(8,i))
215c
216 ! Recovering the positions of the first d.o.fs of each nodes
217 pos1(i) = nloc_dmg%POSI(n1)
218 pos2(i) = nloc_dmg%POSI(n2)
219 pos3(i) = nloc_dmg%POSI(n3)
220 pos4(i) = nloc_dmg%POSI(n4)
221 pos5(i) = nloc_dmg%POSI(n5)
222 pos6(i) = nloc_dmg%POSI(n6)
223 pos7(i) = nloc_dmg%POSI(n7)
224 pos8(i) = nloc_dmg%POSI(n8)
225c
226 ENDIF
227c
228 ENDDO
229c
230 !-----------------------------------------------------------------------
231 ! Computation of non-local forces
232 !-----------------------------------------------------------------------
233 ! Loop over classic brick elements
234#include "vectorize.inc"
235 DO ii = 1, nindx8
236c
237 ! Number of the element
238 i = indx8(ii)
239c
240 ! If the element is not broken, normal computation
241 IF (off(i) /= zero) THEN
242c
243 ! Computing the product NtN*UNL
244 ntn_unl = (unl(pos1(i)) + unl(pos2(i)) + unl(pos3(i)) + unl(pos4(i))
245 . + unl(pos5(i)) + unl(pos6(i)) + unl(pos7(i)) + unl(pos8(i))) / ntn8
246c
247 ! Computing the product NtN*VNL
248 ntn_vnl = (vnl(pos1(i)) + vnl(pos2(i)) + vnl(pos3(i)) + vnl(pos4(i))
249 . + vnl(pos5(i)) + vnl(pos6(i)) + vnl(pos7(i)) + vnl(pos8(i))) / ntn8
250 IF (nodadt > 0) THEN
251 ntn_vnl = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
252 . sqrt(mass(pos2(i))/mass0(pos2(i))),
253 . sqrt(mass(pos3(i))/mass0(pos3(i))),
254 . sqrt(mass(pos4(i))/mass0(pos4(i))),
255 . sqrt(mass(pos5(i))/mass0(pos5(i))),
256 . sqrt(mass(pos6(i))/mass0(pos6(i))),
257 . sqrt(mass(pos7(i))/mass0(pos7(i))),
258 . sqrt(mass(pos8(i))/mass0(pos8(i))))*ntn_vnl
259 ENDIF
260c
261 ! Computation of the product LEN**2 * BtxB
262 b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i))
263 . + btb13(i)*unl(pos3(i)) + btb14(i)*unl(pos4(i)) - btb13(i)*unl(pos5(i))
264 . - btb14(i)*unl(pos6(i)) - btb11(i)*unl(pos7(i)) - btb12(i)*unl(pos8(i)))
265c
266 b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i))
267 . + btb23(i)*unl(pos3(i)) + btb24(i)*unl(pos4(i)) - btb23(i)*unl(pos5(i))
268 . - btb24(i)*unl(pos6(i)) - btb12(i)*unl(pos7(i)) - btb22(i)*unl(pos8(i)))
269c
270 b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i))
271 . + btb33(i)*unl(pos3(i)) + btb34(i)*unl(pos4(i)) - btb33(i)*unl(pos5(i))
272 . - btb34(i)*unl(pos6(i)) - btb13(i)*unl(pos7(i)) - btb23(i)*unl(pos8(i)))
273c
274 b4 = l2 * vol(i) * ( btb14(i)*unl(pos1(i)) + btb24(i)*unl(pos2(i))
275 . + btb34(i)*unl(pos3(i)) + btb44(i)*unl(pos4(i)) - btb34(i)*unl(pos5(i))
276 . - btb44(i)*unl(pos6(i)) - btb14(i)*unl(pos7(i)) - btb24(i)*unl(pos8(i)))
277c
278 b5 = l2 * vol(i) * ( -btb13(i)*unl(pos1(i)) - btb23(i)*unl(pos2(i))
279 . - btb33(i)*unl(pos3(i)) - btb34(i)*unl(pos4(i)) + btb33(i)*unl(pos5(i))
280 . + btb34(i)*unl(pos6(i)) + btb13(i)*unl(pos7(i)) + btb23(i)*unl(pos8(i)))
281c
282 b6 = l2 * vol(i) * ( -btb14(i)*unl(pos1(i)) - btb24(i)*unl(pos2(i))
283 . - btb34(i)*unl(pos3(i)) - btb44(i)*unl(pos4(i)) + btb34(i)*unl(pos5(i))
284 . + btb44(i)*unl(pos6(i)) + btb14(i)*unl(pos7(i)) + btb24(i)*unl(pos8(i)))
285c
286 b7 = l2 * vol(i) * ( -btb11(i)*unl(pos1(i)) - btb12(i)*unl(pos2(i))
287 . - btb13(i)*unl(pos3(i)) - btb14(i)*unl(pos4(i)) + btb13(i)*unl(pos5(i))
288 . + btb14(i)*unl(pos6(i)) + btb11(i)*unl(pos7(i)) + btb12(i)*unl(pos8(i)))
289c
290 b8 = l2 * vol(i) * ( -btb12(i)*unl(pos1(i)) - btb22(i)*unl(pos2(i))
291 . - btb23(i)*unl(pos3(i)) - btb24(i)*unl(pos4(i)) + btb23(i)*unl(pos5(i))
292 . + btb24(i)*unl(pos6(i)) + btb12(i)*unl(pos7(i)) + btb22(i)*unl(pos8(i)))
293c
294 ! Multiplication by the volume of the element (and damping parameter XI)
295 ntn_unl = ntn_unl * vol(i)
296 ntn_vnl = ntn_vnl * xi * vol(i)
297c
298 ! Introducing the internal variable to be regularized
299 ntvar = var_reg(i)*one_over_8* vol(i)
300c
301 ! Computing the elementary non-local forces
302 a = ntn_unl + ntn_vnl - ntvar
303 f1(i) = a + b1
304 f2(i) = a + b2
305 f3(i) = a + b3
306 f4(i) = a + b4
307 f5(i) = a + b5
308 f6(i) = a + b6
309 f7(i) = a + b7
310 f8(i) = a + b8
311c
312 ! Computing nodal equivalent stiffness
313 IF (nodadt > 0) THEN
314 sti1(i) = (abs(l2*btb11(i) + one/ntn8) + abs(l2*btb12(i) + one/ntn8) + abs(l2*btb13(i) + one/ntn8) +
315 . abs(l2*btb14(i) + one/ntn8) + abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb14(i) + one/ntn8) +
316 . abs(-l2*btb11(i) + one/ntn8) + abs(-l2*btb12(i) + one/ntn8))*vol(i)
317 sti2(i) = (abs(l2*btb12(i) + one/ntn8) + abs(l2*btb22(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) +
318 . abs(l2*btb24(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8) +
319 . abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb22(i) + one/ntn8))*vol(i)
320 sti3(i) = (abs(l2*btb13(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) + abs(l2*btb33(i) + one/ntn8) +
321 . abs(l2*btb34(i) + one/ntn8) + abs(-l2*btb33(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) +
322 . abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8))*vol(i)
323 sti4(i) = (abs(l2*btb14(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) +
324 . abs(l2*btb44(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) + abs(-l2*btb44(i) + one/ntn8) +
325 . abs(-l2*btb14(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8))*vol(i)
326 sti5(i) = (abs(-l2*btb13(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) + abs(-l2*btb33(i) + one/ntn8) +
327 . abs(-l2*btb34(i) + one/ntn8) + abs(l2*btb33(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) +
328 . abs(l2*btb13(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8))*vol(i)
329 sti6(i) = (abs(-l2*btb14(i) + one/ntn8) + abs(-l2*btb24(i) + one/ntn8) + abs(-l2*btb34(i) + one/ntn8) +
330 . abs(-l2*btb44(i) + one/ntn8) + abs(l2*btb34(i) + one/ntn8) + abs(l2*btb44(i) + one/ntn8) +
331 . abs(l2*btb14(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8))*vol(i)
332 sti7(i) = (abs(-l2*btb11(i) + one/ntn8) + abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb13(i) + one/ntn8) +
333 . abs(-l2*btb14(i) + one/ntn8) + abs(l2*btb13(i) + one/ntn8) + abs(l2*btb14(i) + one/ntn8) +
334 . abs(l2*btb11(i) + one/ntn8) + abs(l2*btb12(i) + one/ntn8))*vol(i)
335 sti8(i) = (abs(-l2*btb12(i) + one/ntn8) + abs(-l2*btb22(i) + one/ntn8) + abs(-l2*btb23(i) + one/ntn8) +
336 . abs(-l2*btb24(i) + one/ntn8) + abs(l2*btb23(i) + one/ntn8) + abs(l2*btb24(i) + one/ntn8) +
337 . abs(l2*btb12(i) + one/ntn8) + abs(l2*btb22(i) + one/ntn8))*vol(i)
338 ENDIF
339c
340 ! If the element is broken, the non-local wave is absorbed
341 ELSE
342c
343 ! Initial element characteristic length
344 lc(i) = vol0(i)**third
345c
346 IF (nodadt > 0) THEN
347
348 ! Non-local absorbing forces
349 f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*zeta*sspnl*half*
350 . (vnl(pos1(i))+vnl0(pos1(i)))*(three/four)*(lc(i)**2)
351 f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*zeta*sspnl*half*
352 . (vnl(pos2(i))+vnl0(pos2(i)))*(three/four)*(lc(i)**2)
353 f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*zeta*sspnl*half*
354 . (vnl(pos3(i))+vnl0(pos3(i)))*(three/four)*(lc(i)**2)
355 f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*zeta*sspnl*half*
356 . (vnl(pos4(i))+vnl0(pos4(i)))*(three/four)*(lc(i)**2)
357 f5(i) = sqrt(mass(pos5(i))/mass0(pos5(i)))*zeta*sspnl*half*
358 . (vnl(pos5(i))+vnl0(pos5(i)))*(three/four)*(lc(i)**2)
359 f6(i) = sqrt(mass(pos6(i))/mass0(pos6(i)))*zeta*sspnl*half*
360 . (vnl(pos6(i))+vnl0(pos6(i)))*(three/four)*(lc(i)**2)
361 f7(i) = sqrt(mass(pos7(i))/mass0(pos7(i)))*zeta*sspnl*half*
362 . (vnl(pos7(i))+vnl0(pos7(i)))*(three/four)*(lc(i)**2)
363 f8(i) = sqrt(mass(pos8(i))/mass0(pos8(i)))*zeta*sspnl*half*
364 . (vnl(pos8(i))+vnl0(pos8(i)))*(three/four)*(lc(i)**2)
365 ! Computing nodal equivalent stiffness
366 sti1(i) = em20
367 sti2(i) = em20
368 sti3(i) = em20
369 sti4(i) = em20
370 sti5(i) = em20
371 sti6(i) = em20
372 sti7(i) = em20
373 sti8(i) = em20
374 ELSE
375 ! Non-local absorbing forces
376 f1(i) = zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(three/four)*(lc(i)**2)
377 f2(i) = zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(three/four)*(lc(i)**2)
378 f3(i) = zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(three/four)*(lc(i)**2)
379 f4(i) = zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(three/four)*(lc(i)**2)
380 f5(i) = zeta*sspnl*half*(vnl(pos5(i))+vnl0(pos5(i)))*(three/four)*(lc(i)**2)
381 f6(i) = zeta*sspnl*half*(vnl(pos6(i))+vnl0(pos6(i)))*(three/four)*(lc(i)**2)
382 f7(i) = zeta*sspnl*half*(vnl(pos7(i))+vnl0(pos7(i)))*(three/four)*(lc(i)**2)
383 f8(i) = zeta*sspnl*half*(vnl(pos8(i))+vnl0(pos8(i)))*(three/four)*(lc(i)**2)
384 ENDIF
385 ENDIF
386 ENDDO
387c
388 ! Loop over penta degenerated bricks
389#include "vectorize.inc"
390 DO ii = 1, nindx6
391c
392 ! Number of the element
393 i = indx6(ii)
394c
395 ! If the element is not broken, normal computation
396 IF (off(i) /= zero) THEN
397c
398 ! Computing the product NtN*UNL
399 ntn_unl = (unl(pos1(i)) + unl(pos2(i)) + unl(pos3(i)) +
400 . unl(pos4(i)) + unl(pos5(i)) + unl(pos6(i))) / ntn6
401c
402 ! Computing the product NtN*VNL
403 ntn_vnl = (vnl(pos1(i)) + vnl(pos2(i)) + vnl(pos3(i)) +
404 . vnl(pos4(i)) + vnl(pos5(i)) + vnl(pos6(i))) / ntn6
405 IF (nodadt > 0) THEN
406 ntn_vnl = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
407 . sqrt(mass(pos2(i))/mass0(pos2(i))),
408 . sqrt(mass(pos3(i))/mass0(pos3(i))),
409 . sqrt(mass(pos4(i))/mass0(pos4(i))),
410 . sqrt(mass(pos5(i))/mass0(pos5(i))),
411 . sqrt(mass(pos6(i))/mass0(pos6(i))))*ntn_vnl
412 ENDIF
413c
414 ! Computation of the product LEN**2 * BtxB
415 b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i))
416 . + btb13(i)*unl(pos3(i)) - btb13(i)*unl(pos4(i))
417 . - btb12(i)*unl(pos5(i)) - btb11(i)*unl(pos6(i)) )
418c
419 b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i))
420 . + btb23(i)*unl(pos3(i)) - btb23(i)*unl(pos4(i))
421 . - btb22(i)*unl(pos5(i)) - btb12(i)*unl(pos6(i)) )
422c
423 b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i))
424 . + btb33(i)*unl(pos3(i)) - btb33(i)*unl(pos4(i))
425 . - btb23(i)*unl(pos5(i)) - btb13(i)*unl(pos6(i)) )
426c
427 b4 = l2 * vol(i) * (-btb13(i)*unl(pos1(i)) - btb23(i)*unl(pos2(i))
428 . - btb33(i)*unl(pos3(i)) + btb33(i)*unl(pos4(i))
429 . + btb23(i)*unl(pos5(i)) + btb13(i)*unl(pos6(i)) )
430c
431 b5 = l2 * vol(i) * (-btb12(i)*unl(pos1(i)) - btb22(i)*unl(pos2(i))
432 . - btb23(i)*unl(pos3(i)) + btb23(i)*unl(pos4(i))
433 . + btb22(i)*unl(pos5(i)) + btb12(i)*unl(pos6(i)) )
434c
435 b6 = l2 * vol(i) * (-btb11(i)*unl(pos1(i)) - btb12(i)*unl(pos2(i))
436 . - btb13(i)*unl(pos3(i)) + btb13(i)*unl(pos4(i))
437 . + btb12(i)*unl(pos5(i)) + btb11(i)*unl(pos6(i)) )
438c
439 ! Multiplication by the volume of the element (and damping parameter XI)
440 ntn_unl = ntn_unl * vol(i)
441 ntn_vnl = ntn_vnl * xi * vol(i)
442c
443 ! Introducing the internal variable to be regularized
444 ntvar = var_reg(i)*one_over_6* vol(i)
445c
446 ! Computing the elementary non-local forces
447 a = ntn_unl + ntn_vnl - ntvar
448 f1(i) = a + b1
449 f2(i) = a + b2
450 f3(i) = a + b3
451 f4(i) = a + b4
452 f5(i) = a + b5
453 f6(i) = a + b6
454c
455 ! Computing nodal equivalent stiffness
456 IF (nodadt > 0) THEN
457 sti1(i) = (abs(l2*btb11(i) + one/ntn6) + abs(l2*btb12(i) + one/ntn6) +
458 . abs(l2*btb13(i) + one/ntn6) + abs(-l2*btb13(i) + one/ntn6) +
459 . abs(-l2*btb12(i) + one/ntn6) + abs(-l2*btb11(i) + one/ntn6))*vol(i)
460 sti2(i) = (abs(l2*btb12(i) + one/ntn6) + abs(l2*btb22(i) + one/ntn6) +
461 . abs(l2*btb23(i) + one/ntn6) + abs(-l2*btb23(i) + one/ntn6) +
462 . abs(-l2*btb22(i) + one/ntn6) + abs(-l2*btb12(i) + one/ntn6))*vol(i)
463 sti3(i) = (abs(l2*btb13(i) + one/ntn6) + abs(l2*btb23(i) + one/ntn6) +
464 . abs(l2*btb33(i) + one/ntn6) + abs(-l2*btb33(i) + one/ntn6) +
465 . abs(-l2*btb23(i) + one/ntn6) + abs(-l2*btb13(i) + one/ntn6))*vol(i)
466 sti4(i) = (abs(-l2*btb13(i) + one/ntn6) + abs(-l2*btb23(i) + one/ntn6) +
467 . abs(-l2*btb33(i) + one/ntn6) + abs(l2*btb33(i) + one/ntn6) +
468 . abs(l2*btb23(i) + one/ntn6) + abs(l2*btb13(i) + one/ntn6))*vol(i)
469 sti5(i) = (abs(-l2*btb12(i) + one/ntn6) + abs(-l2*btb22(i) + one/ntn6) +
470 . abs(-l2*btb23(i) + one/ntn6) + abs(l2*btb23(i) + one/ntn6) +
471 . abs(l2*btb22(i) + one/ntn6) + abs(l2*btb12(i) + one/ntn6))*vol(i)
472 sti6(i) = (abs(-l2*btb11(i) + one/ntn6) + abs(-l2*btb12(i) + one/ntn6) +
473 . abs(-l2*btb13(i) + one/ntn6) + abs(l2*btb13(i) + one/ntn6) +
474 . abs(l2*btb12(i) + one/ntn6) + abs(l2*btb11(i) + one/ntn6))*vol(i)
475 ENDIF
476c
477 ELSE
478c
479 ! Initial element characteristic length
480 lc(i) = vol0(i)**third
481c
482 IF (nodadt > 0) THEN
483
484 ! Non-local absorbing forces
485 f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*zeta*sspnl*half*
486 . (vnl(pos1(i))+vnl0(pos1(i)))*(two*third)*(lc(i)**2)
487 f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*zeta*sspnl*half*
488 . (vnl(pos2(i))+vnl0(pos2(i)))*(two*third)*(lc(i)**2)
489 f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*zeta*sspnl*half*
490 . (vnl(pos3(i))+vnl0(pos3(i)))*(two*third)*(lc(i)**2)
491 f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*zeta*sspnl*half*
492 . (vnl(pos4(i))+vnl0(pos4(i)))*(two*third)*(lc(i)**2)
493 f5(i) = sqrt(mass(pos5(i))/mass0(pos5(i)))*zeta*sspnl*half*
494 . (vnl(pos5(i))+vnl0(pos5(i)))*(two*third)*(lc(i)**2)
495 f6(i) = sqrt(mass(pos6(i))/mass0(pos6(i)))*zeta*sspnl*half*
496 . (vnl(pos6(i))+vnl0(pos6(i)))*(two*third)*(lc(i)**2)
497 ! Computing nodal equivalent stiffness
498 sti1(i) = em20
499 sti2(i) = em20
500 sti3(i) = em20
501 sti4(i) = em20
502 sti5(i) = em20
503 sti6(i) = em20
504 ELSE
505 ! Non-local absorbing forces
506 f1(i) = zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(two*third)*(lc(i)**2)
507 f2(i) = zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(two*third)*(lc(i)**2)
508 f3(i) = zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(two*third)*(lc(i)**2)
509 f4(i) = zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(two*third)*(lc(i)**2)
510 f5(i) = zeta*sspnl*half*(vnl(pos5(i))+vnl0(pos5(i)))*(two*third)*(lc(i)**2)
511 f6(i) = zeta*sspnl*half*(vnl(pos6(i))+vnl0(pos6(i)))*(two*third)*(lc(i)**2)
512 ENDIF
513 ENDIF
514 ENDDO
515c
516c-----------------------------------------------------------------------
517c Assemblage
518c-----------------------------------------------------------------------
519 ! If PARITH/OFF
520 IF (iparit == 0) THEN
521 fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
522 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
523c
524 ! Loop over classic hexa brick elements
525#include "vectorize.inc"
526 DO ii=1,nindx8
527 i = indx8(ii)
528 ! Assembling the forces in the classic way
529 fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
530 fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
531 fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
532 fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
533 fnl(pos5(i)) = fnl(pos5(i)) - f5(i)
534 fnl(pos6(i)) = fnl(pos6(i)) - f6(i)
535 fnl(pos7(i)) = fnl(pos7(i)) - f7(i)
536 fnl(pos8(i)) = fnl(pos8(i)) - f8(i)
537 IF (nodadt > 0) THEN
538 ! Spectral radius of stiffness matrix
539 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i),sti5(i),sti6(i),sti7(i),sti8(i))
540 ! Computing nodal stiffness
541 stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
542 stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
543 stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
544 stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
545 stifnl(pos5(i)) = stifnl(pos5(i)) + maxstif
546 stifnl(pos6(i)) = stifnl(pos6(i)) + maxstif
547 stifnl(pos7(i)) = stifnl(pos7(i)) + maxstif
548 stifnl(pos8(i)) = stifnl(pos8(i)) + maxstif
549 ENDIF
550 ENDDO
551c
552 ! Loop over penta degenerated bricks
553#include "vectorize.inc"
554 DO ii=1,nindx6
555 i = indx6(ii)
556 ! Assembling the forces in the classic way
557 fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
558 fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
559 fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
560 fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
561 fnl(pos5(i)) = fnl(pos5(i)) - f5(i)
562 fnl(pos6(i)) = fnl(pos6(i)) - f6(i)
563 IF (nodadt > 0) THEN
564 ! Spectral radius of stiffness matrix
565 maxstif = max(sti1(i),sti2(i),sti3(i),
566 . sti4(i),sti5(i),sti6(i))
567 ! Computing nodal stiffness
568 stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
569 stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
570 stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
571 stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
572 stifnl(pos5(i)) = stifnl(pos5(i)) + maxstif
573 stifnl(pos6(i)) = stifnl(pos6(i)) + maxstif
574 ENDIF
575 ENDDO
576c
577 ! If PARITH/ON
578 ELSE
579c
580 ! Loop over classic hexa bricks
581 DO l=1,nindx8
582 i = indx8(l)
583 ii = i + nft
584c
585 ! Spectral radius of stiffness matrix
586 IF (nodadt > 0) THEN
587 maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i),sti5(i),sti6(i),sti7(i),sti8(i))
588 ENDIF
589c
590 k = nloc_dmg%IADS(1,ii)
591 nloc_dmg%FSKY(k,1) = -f1(i)
592 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
593c
594 k = nloc_dmg%IADS(2,ii)
595 nloc_dmg%FSKY(k,1) = -f2(i)
596 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
597c
598 k = nloc_dmg%IADS(3,ii)
599 nloc_dmg%FSKY(k,1) = -f3(i)
600 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
601c
602 k = nloc_dmg%IADS(4,ii)
603 nloc_dmg%FSKY(k,1) = -f4(i)
604 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
605c
606 k = nloc_dmg%IADS(5,ii)
607 nloc_dmg%FSKY(k,1) = -f5(i)
608 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
609c
610 k = nloc_dmg%IADS(6,ii)
611 nloc_dmg%FSKY(k,1) = -f6(i)
612 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
613c
614 k = nloc_dmg%IADS(7,ii)
615 nloc_dmg%FSKY(k,1) = -f7(i)
616 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
617c
618 k = nloc_dmg%IADS(8,ii)
619 nloc_dmg%FSKY(k,1) = -f8(i)
620 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
621c
622 ENDDO
623c
624 ! Loop over penta degenerated bricks
625 DO l = 1,nindx6
626 i = indx6(l)
627 ii = i + nft
628c
629 ! Spectral radius of stiffness matrix
630 IF (nodadt > 0) THEN
631 maxstif = max(sti1(i),sti2(i),sti3(i),
632 . sti4(i),sti5(i),sti6(i))
633 ENDIF
634c
635 k = nloc_dmg%IADS(1,ii)
636 nloc_dmg%FSKY(k,1) = -f1(i)
637 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
638c
639 k = nloc_dmg%IADS(2,ii)
640 nloc_dmg%FSKY(k,1) = -f2(i)
641 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
642c
643 k = nloc_dmg%IADS(3,ii)
644 nloc_dmg%FSKY(k,1) = -f3(i)
645 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
646c
647 k = nloc_dmg%IADS(5,ii)
648 nloc_dmg%FSKY(k,1) = -f4(i)
649 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
650c
651 k = nloc_dmg%IADS(6,ii)
652 nloc_dmg%FSKY(k,1) = -f5(i)
653 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
654c
655 k = nloc_dmg%IADS(7,ii)
656 nloc_dmg%FSKY(k,1) = -f6(i)
657 IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
658c
659 ENDDO
660 ENDIF
661c
662 !-----------------------------------------------------------------------
663 ! Computing non-local timestep
664 !-----------------------------------------------------------------------
665 IF (nodadt == 0) THEN
666 DO i = 1,nel
667 ! If the element is not broken, normal computation
668 IF (off(i)/=zero) THEN
669 ! Non-local critical time-step in the plane
670 dtnl = (two*(min(vol(i)**third,le_max))*sqrt(three*zeta))/
671 . sqrt(twelve*l2 + (min(vol(i)**third,le_max))**2)
672 ! Retaining the minimal value
673 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl)
674 ENDIF
675 ENDDO
676 ENDIF
677c
678 ! Deallocation of tables
679 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
680 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
681 IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
682 IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
683 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
684 IF (ALLOCATED(btb23)) DEALLOCATE(btb23)
685 IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
686 IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
687 IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
688 IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
689 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
690 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
691 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
692 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
693 IF (ALLOCATED(pos5)) DEALLOCATE(pos5)
694 IF (ALLOCATED(pos6)) DEALLOCATE(pos6)
695 IF (ALLOCATED(pos7)) DEALLOCATE(pos7)
696 IF (ALLOCATED(pos8)) DEALLOCATE(pos8)
697 IF (ALLOCATED(f1)) DEALLOCATE(f1)
698 IF (ALLOCATED(f2)) DEALLOCATE(f2)
699 IF (ALLOCATED(f3)) DEALLOCATE(f3)
700 IF (ALLOCATED(f4)) DEALLOCATE(f4)
701 IF (ALLOCATED(f5)) DEALLOCATE(f5)
702 IF (ALLOCATED(f6)) DEALLOCATE(f6)
703 IF (ALLOCATED(f7)) DEALLOCATE(f7)
704 IF (ALLOCATED(f8)) DEALLOCATE(f8)
705 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
706 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
707 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
708 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
709 IF (ALLOCATED(sti5)) DEALLOCATE(sti5)
710 IF (ALLOCATED(sti6)) DEALLOCATE(sti6)
711 IF (ALLOCATED(sti7)) DEALLOCATE(sti7)
712 IF (ALLOCATED(sti8)) DEALLOCATE(sti8)
713 IF (ALLOCATED(lc)) DEALLOCATE(lc)
714c-----------
715 END
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
sfint_reg
subroutine sfint_reg(nloc_dmg, var_reg, nel, off, vol, nlocs, vol0, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, imat, itask, dt2t, nft)
Definition sfint_reg.F:39
nlocal_reg_mod::nlocal_str_
Definition nlocal_reg_mod.F:124