OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cfint_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!|| cfint_reg ../engine/source/elements/shell/coque/cfint_reg.F
25!||--- called by ------------------------------------------------------
26!|| cforc3 ../engine/source/elements/shell/coque/cforc3.F
27!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.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 cfint_reg(
33 1 NLOC_DMG,VAR_REG, THK, NEL,
34 2 OFF, AREA, NC1, NC2,
35 3 NC3, NC4, PX1, PY1,
36 4 PX2, PY2, BUFNL, IMAT,
37 5 NDDL, ITASK, DT2T, LE,
38 6 THK0, AREA0, NFT)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
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 "com20_c.inc"
53#include "com08_c.inc"
54#include "scr02_c.inc"
55#include "scr18_c.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER, INTENT(IN) :: NFT
60 INTEGER :: NEL,IMAT,NDDL,ITASK
61 INTEGER, DIMENSION(NEL) :: NC1,NC2,NC3,NC4
62 my_real, DIMENSION(NEL,NDDL), INTENT(INOUT)::
63 . VAR_REG
64 my_real, DIMENSION(NEL), INTENT(IN) ::
65 . area,off,px1,py1,px2,py2,thk,le,thk0,area0
66 my_real, INTENT(INOUT) ::
67 . dt2t
68 TYPE(nlocal_str_),TARGET :: NLOC_DMG
69 TYPE(BUF_NLOC_) ,TARGET :: BUFNL
70C-----------------------------------------------
71C L o c a l V a r i a b l e s
72C-----------------------------------------------
73 INTEGER I,II,K,N1,N2,N3,N4,L_NLOC,J,NDNOD
74 my_real
75 . l2,ntn,ntn_unl,ntn_vnl,xi,ntvar,a,zeta,
76 . b1,b2,b3,b4,sspnl,nth1,nth2,bth1,bth2,k1,k2,k12,
77 . dtnl_th,dtnl,le_max,maxstif,dtnod,dt2p
78 my_real, DIMENSION(:,:), ALLOCATABLE ::
79 . f1,f2,f3,f4,sti1,sti2,sti3,sti4
80 my_real, DIMENSION(:) ,ALLOCATABLE ::
81 . btb11,btb12,btb22,vol
82 INTEGER, DIMENSION(:), ALLOCATABLE ::
83 . POS1,POS2,POS3,POS4
84 my_real, POINTER, DIMENSION(:) ::
85 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
86 my_real, POINTER, DIMENSION(:,:) ::
87 . massth,unlth,vnlth,fnlth
88 my_real, DIMENSION(:,:), ALLOCATABLE ::
89 . stifnlth,dtn
90c-----------------------------------------------------------------------
91c VAR_REG : variable a regulariser (local cumulated plastic strain)
92c NTVAR = NT * VAR_REG
93C=======================================================================
94 ! Recovering non-local parameters
95 l2 = nloc_dmg%LEN(imat)**2 ! Non-local internal length ** 2
96 xi = nloc_dmg%DAMP(imat) ! Non-local damping parameter
97 zeta = nloc_dmg%DENS(imat) ! Non-local density
98 sspnl = nloc_dmg%SSPNL(imat) ! Non-local sound speed
99 l_nloc = nloc_dmg%L_NLOC ! Length of non-local tables
100 le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
101 ! Allocation of elementary forces vectors
102 ALLOCATE(f1(nel,nddl),f2(nel,nddl),f3(nel,nddl),f4(nel,nddl))
103 ! Only for nodal timestep
104 IF (nodadt > 0) THEN
105 ! Nodal stiffness
106 ALLOCATE(sti1(nel,nddl),sti2(nel,nddl),sti3(nel,nddl),sti4(nel,nddl))
107 ! Non-local mass
108 mass => nloc_dmg%MASS(1:l_nloc)
109 ! Initial non-local mass
110 mass0 => nloc_dmg%MASS0(1:l_nloc)
111 ENDIF
112 ALLOCATE(btb11(nel),btb12(nel),btb22(nel),vol(nel),
113 . pos1(nel),pos2(nel),pos3(nel),pos4(nel))
114 ! Recovering non-local data
115 vnl => nloc_dmg%VNL(1:l_nloc) ! Non-local variable velocities
116 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc) ! Non-local variable velocities
117 unl => nloc_dmg%UNL(1:l_nloc) ! Non-local cumulated variable
118 ! Constant for quad elements
119 ntn = four*four
120c
121 !-----------------------------------------------------------------------
122 ! Computation of the element volume and the BtB matrix product
123 !-----------------------------------------------------------------------
124 ! Loop over elements
125# include "vectorize.inc"
126 DO i=1,nel
127c
128 ! Recovering the nodes of the triangle element
129 n1 = nloc_dmg%IDXI(nc1(i))
130 n2 = nloc_dmg%IDXI(nc2(i))
131 n3 = nloc_dmg%IDXI(nc3(i))
132 n4 = nloc_dmg%IDXI(nc4(i))
133c
134 ! Recovering the positions of the first d.o.fs of each nodes
135 pos1(i) = nloc_dmg%POSI(n1)
136 pos2(i) = nloc_dmg%POSI(n2)
137 pos3(i) = nloc_dmg%POSI(n3)
138 pos4(i) = nloc_dmg%POSI(n4)
139c
140 ! Computation of the element volume
141 vol(i) = area(i)*thk(i)
142c
143 ! Computation of the product BtxB
144 btb11(i) = px1(i)**2 + py1(i)**2
145 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i)
146 btb22(i) = px2(i)**2 + py2(i)**2
147c
148 ENDDO
149c
150 !-----------------------------------------------------------------------
151 ! Pre-treatment non-local regularization in the thickness
152 !-----------------------------------------------------------------------
153 ! Only if NDDL > 1
154 IF ((nddl > 1).AND.(l2>zero)) THEN
155c
156 ! Allocation of the velocities predictor
157 IF (nddl > 2) THEN
158 IF (nodadt > 0) THEN
159 ALLOCATE(stifnlth(nel,nddl+1))
160 ALLOCATE(dtn(nel,nddl+1))
161 ENDIF
162 ndnod = nddl+1
163 ELSE
164 IF (nodadt > 0) THEN
165 ALLOCATE(stifnlth(nel,nddl))
166 ALLOCATE(dtn(nel,nddl))
167 ENDIF
168 ndnod = nddl
169 ENDIF
170c
171 ! Pointing the non-local values in the thickness of the corresponding element
172 massth => bufnl%MASSTH(1:nel,1:ndnod)
173 unlth => bufnl%UNLTH(1:nel,1:ndnod)
174 vnlth => bufnl%VNLTH(1:nel,1:ndnod)
175 fnlth => bufnl%FNLTH(1:nel,1:ndnod)
176c
177 DO k = 1,ndnod
178 DO i = 1,nel
179 ! Resetting non-local forces
180 fnlth(i,k) = zero
181 ! Resetting non-local nodal stiffness
182 IF (nodadt > 0) THEN
183 stifnlth(i,k) = em20
184 ENDIF
185 ENDDO
186 ENDDO
187c
188 ! Computation of non-local forces in the shell thickness
189 DO k = 1, nddl
190c
191 ! Computation of shape functions value
192 IF ((nddl==2).AND.(k==2)) THEN
193 nth1 = (z0(k,nddl) - zth(k,nddl)) / (zth(k-1,nddl) - zth(k,nddl))
194 nth2 = (z0(k,nddl) - zth(k-1,nddl)) / (zth(k,nddl) - zth(k-1,nddl))
195 ELSE
196 nth1 = (z0(k,nddl) - zth(k+1,nddl)) / (zth(k,nddl) - zth(k+1,nddl))
197 nth2 = (z0(k,nddl) - zth(k,nddl)) / (zth(k+1,nddl) - zth(k,nddl))
198 ENDIF
199c
200 ! Loop over elements
201 DO i = 1,nel
202 ! Computation of B-matrix values
203 IF ((nddl==2).AND.(k==2)) THEN
204 bth1 = (one/(zth(k-1,nddl) - zth(k,nddl)))*(one/thk(i))
205 bth2 = (one/(zth(k,nddl) - zth(k-1,nddl)))*(one/thk(i))
206 ELSE
207 bth1 = (one/(zth(k,nddl) - zth(k+1,nddl)))*(one/thk(i))
208 bth2 = (one/(zth(k+1,nddl) - zth(k,nddl)))*(one/thk(i))
209 ENDIF
210c
211 ! Computation of the non-local K matrix
212 k1 = l2*(bth1**2) + nth1**2
213 k12 = l2*(bth1*bth2)+ (nth1*nth2)
214 k2 = l2*(bth2**2) + nth2**2
215c
216 ! Computation of the non-local forces
217 IF ((nddl==2).AND.(k==2)) THEN
218 fnlth(i,k-1) = fnlth(i,k-1) + (k1*unlth(i,k-1) + k12*unlth(i,k)
219 . + xi*((nth1**2)*vnlth(i,k-1)
220 . + (nth1*nth2)*vnlth(i,k))
221 . - (nth1*var_reg(i,k)))*vol(i)*wf(k,nddl)
222 fnlth(i,k) = fnlth(i,k) + (k12*unlth(i,k-1) + k2*unlth(i,k)
223 . + xi*(nth1*nth2*vnlth(i,k-1)
224 . + (nth2**2)*vnlth(i,k))
225 . - nth2*var_reg(i,k))*vol(i)*wf(k,nddl)
226 ELSE
227 fnlth(i,k) = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1)
228 . + xi*((nth1**2)*vnlth(i,k)
229 . + (nth1*nth2)*vnlth(i,k+1))
230 . - (nth1*var_reg(i,k)))*vol(i)*wf(k,nddl)
231 fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i,k) + k2*unlth(i,k+1)
232 . + xi*(nth1*nth2*vnlth(i,k)
233 . + (nth2**2)*vnlth(i,k+1))
234 . - nth2*var_reg(i,k))*vol(i)*wf(k,nddl)
235 ENDIF
236c
237 ! Computation of non-local nodal stiffness
238 IF (nodadt > 0) THEN
239 IF ((nddl==2).AND.(k==2)) THEN
240 stifnlth(i,k-1) = stifnlth(i,k-1) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
241 stifnlth(i,k) = stifnlth(i,k) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
242 ELSE
243 stifnlth(i,k) = stifnlth(i,k) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
244 stifnlth(i,k+1) = stifnlth(i,k+1) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
245 ENDIF
246 ENDIF
247c
248 ENDDO
249 ENDDO
250c
251 ! Updating non-local mass with /DT/NODA
252 IF (nodadt > 0) THEN
253C
254 ! Initial computation of the nodal timestep
255 dtnod = ep20
256 DO k = 1,ndnod
257 DO i = 1,nel
258 dtn(i,k) = dtfac1(11)*cdamp*sqrt(two*massth(i,k)/max(stifnlth(i,k),em20))
259 dtnod = min(dtn(i,k),dtnod)
260 ENDDO
261 ENDDO
262C
263 ! /DT/NODA/CSTX - Constant timestep with added mass
264 IF ((idtmin(11)==3).OR.(idtmin(11)==4).OR.(idtmin(11)==8)) THEN
265 ! Added mass computation if necessary
266 IF (dtnod < dtmin1(11)*(sqrt(csta))) THEN
267 DO k = 1,ndnod
268 DO i = 1,nel
269 IF (dtn(i,k) < dtmin1(11)) THEN
270 dt2p = dtmin1(11)/(dtfac1(11)*cdamp)
271 massth(i,k) = max(massth(i,k),csta*half*stifnlth(i,k)*dt2p*dt2p*onep00001)
272 ENDIF
273 ENDDO
274 ENDDO
275 dtnod = dtmin1(11)*(sqrt(csta))
276 ENDIF
277 ENDIF
278C
279 ! Classical nodal timestep check
280 IF (dtnod < dt2t) THEN
281 dt2t = min(dt2t,dtnod)
282 ENDIF
283 ENDIF
284c
285 DO k = 1,ndnod
286 DO i = 1,nel
287 ! Updating the non-local in-thickness velocities
288 vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt12
289 ENDDO
290 ENDDO
291c
292 DO k = 1,ndnod
293 DO i = 1,nel
294 ! Computing the non-local in-thickness cumulated values
295 unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt1
296 ENDDO
297 ENDDO
298c
299 ! Transfert at the integration point
300 DO k = 1, nddl
301 !Computation of shape functions value
302 IF ((nddl==2).AND.(k==2)) THEN
303 nth1 = (z0(k,nddl) - zth(k,nddl))/(zth(k-1,nddl) - zth(k,nddl))
304 nth2 = (z0(k,nddl) - zth(k-1,nddl)) /(zth(k,nddl) - zth(k-1,nddl))
305 ELSE
306 nth1 = (z0(k,nddl) - zth(k+1,nddl))/(zth(k,nddl) - zth(k+1,nddl))
307 nth2 = (z0(k,nddl) - zth(k,nddl)) /(zth(k+1,nddl) - zth(k,nddl))
308 ENDIF
309 ! Loop over elements
310 DO i = 1,nel
311 !Integration points non-local variables
312 IF ((nddl==2).AND.(k==2)) THEN
313 var_reg(i,k) = nth1*unlth(i,k-1) + nth2*unlth(i,k)
314 ELSE
315 var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i,k+1)
316 ENDIF
317 ENDDO
318 ENDDO
319 ENDIF
320c
321 !-----------------------------------------------------------------------
322 ! Computation of the elementary non-local forces
323 !-----------------------------------------------------------------------
324 ! Loop over additional degrees of freedom
325 DO k = 1, nddl
326c
327 ! Loop over elements
328# include "vectorize.inc"
329 DO i = 1, nel
330c
331 ! If the element is not broken, normal computation
332 IF (off(i) /= zero) THEN
333 ! Computation of the product LEN**2 * BtxB
334 ! Warning: the derivatives of the shape function does not take into account the volume of the element.
335 ! That is why a factor (1/VOL**2) is added in B1, B2, B3 and B4.
336 b1 = (l2 / vol(i))*wf(k,nddl)*(btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
337 . - btb11(i)*unl(pos3(i)+k-1) - btb12(i)*unl(pos4(i)+k-1))
338c
339 b2 = (l2 / vol(i))*wf(k,nddl)*(btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
340 . - btb12(i)*unl(pos3(i)+k-1) - btb22(i)*unl(pos4(i)+k-1))
341c
342 b3 = (l2 / vol(i))*wf(k,nddl)*(-btb11(i)*unl(pos1(i)+k-1) - btb12(i)*unl(pos2(i)+k-1)
343 . + btb11(i)*unl(pos3(i)+k-1) + btb12(i)*unl(pos4(i)+k-1))
344c
345 b4 = (l2 / vol(i))*wf(k,nddl)*(-btb12(i)*unl(pos1(i)+k-1) - btb22(i)*unl(pos2(i)+k-1)
346 . + btb12(i)*unl(pos3(i)+k-1) + btb22(i)*unl(pos4(i)+k-1))
347c
348 ! Computing the product NtN*UNL
349 ntn_unl = (unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) + unl(pos4(i)+k-1))/ntn
350c
351 ! Computing the product XDAMP*NtN*VNL
352 ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) + vnl(pos4(i)+k-1))/ntn
353 IF (nodadt > 0) THEN
354 ntn_vnl = min(sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1)),
355 . sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1)),
356 . sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1)),
357 . sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1)))*ntn_vnl
358 ENDIF
359c
360 ! Multiplication by the volume of the element
361 ntn_unl = ntn_unl * vol(i) * wf(k,nddl)
362 ntn_vnl = ntn_vnl * xi * vol(i) * wf(k,nddl)
363c
364 ! Introducing the internal variable to be regularized
365 ntvar = var_reg(i,k)*fourth*vol(i)*wf(k,nddl)
366c
367 ! Computing the elementary non-local forces
368 a = ntn_unl + ntn_vnl - ntvar
369 f1(i,k) = a + b1
370 f2(i,k) = a + b2
371 f3(i,k) = a + b3
372 f4(i,k) = a + b4
373c
374 ! Computing nodal equivalent stiffness
375 IF (nodadt > 0) THEN
376 sti1(i,k) = wf(k,nddl)*(abs((l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
377 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
378 . abs(-(l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
379 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)))
380 sti2(i,k) = wf(k,nddl)*(abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
381 . abs((l2/vol(i))*btb22(i) + one/ntn*vol(i)) +
382 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
383 . abs(-(l2/vol(i))*btb22(i) + one/ntn*vol(i)))
384 sti3(i,k) = wf(k,nddl)*(abs(-(l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
385 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
386 . abs((l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
387 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)))
388 sti4(i,k) = wf(k,nddl)*(abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
389 . abs(-(l2/vol(i))*btb22(i) + one/ntn*vol(i)) +
390 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
391 . abs((l2/vol(i))*btb22(i) + one/ntn*vol(i)))
392 ENDIF
393c
394 ! If the element is broken, the non-local wave is absorbed
395 ELSE
396 IF (nodadt > 0) THEN
397 ! Non-local absorbing forces
398 f1(i,k) = wf(k,nddl)*sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*
399 . half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*sqrt(area0(i))*thk0(i)
400 f2(i,k) = wf(k,nddl)*sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*
401 . half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*sqrt(area0(i))*thk0(i)
402 f3(i,k) = wf(k,nddl)*sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*zeta*sspnl*
403 . half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*sqrt(area0(i))*thk0(i)
404 f4(i,k) = wf(k,nddl)*sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*
405 . half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*sqrt(area0(i))*thk0(i)
406 ! Computing nodal equivalent stiffness
407 sti1(i,k) = em20
408 sti2(i,k) = em20
409 sti3(i,k) = em20
410 sti4(i,k) = em20
411 ELSE
412 ! Non-local absorbing forces
413 f1(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*sqrt(area0(i))*thk0(i)
414 f2(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*sqrt(area0(i))*thk0(i)
415 f3(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*sqrt(area0(i))*thk0(i)
416 f4(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*sqrt(area0(i))*thk0(i)
417 ENDIF
418 ENDIF
419 ENDDO
420 ENDDO
421c
422 !-----------------------------------------------------------------------
423 ! Assembling of the non-local forces
424 !-----------------------------------------------------------------------
425c
426 ! If PARITH/OFF
427 IF (iparit == 0) THEN
428 ! Recovering non-local internal forces
429 fnl => nloc_dmg%FNL(1:l_nloc,itask+1) ! Non-local forces
430 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
431 ! Loop over elements
432 DO i=1,nel
433 ! Loop over non-local degrees of freedom
434# include "vectorize.inc"
435 DO k = 1,nddl
436 ! Assembling non-local forces
437 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
438 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
439 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
440 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
441 IF (nodadt > 0) THEN
442 ! Spectral radius of stiffness matrix
443 maxstif = max(sti1(i,k),sti2(i,k),sti3(i,k),sti4(i,k))
444 ! Computing nodal stiffness
445 stifnl(pos1(i)+k-1) = stifnl(pos1(i)+k-1) + maxstif
446 stifnl(pos2(i)+k-1) = stifnl(pos2(i)+k-1) + maxstif
447 stifnl(pos3(i)+k-1) = stifnl(pos3(i)+k-1) + maxstif
448 stifnl(pos4(i)+k-1) = stifnl(pos4(i)+k-1) + maxstif
449 ENDIF
450 ENDDO
451 ENDDO
452c
453 ! If PARITH/ON
454 ELSE
455 ! Loop over additional d.o.fs
456 DO j = 1,nddl
457c
458 ! Loop over elements
459 DO i=1,nel
460 ii = i + nft
461c
462 ! Spectral radius of stiffness matrix
463 IF (nodadt > 0) THEN
464 maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j))
465 ENDIF
466c
467 k = nloc_dmg%IADC(1,ii)
468 nloc_dmg%FSKY(k,j) = -f1(i,j)
469 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
470c
471 k = nloc_dmg%IADC(2,ii)
472 nloc_dmg%FSKY(k,j) = -f2(i,j)
473 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
474c
475 k = nloc_dmg%IADC(3,ii)
476 nloc_dmg%FSKY(k,j) = -f3(i,j)
477 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
478c
479 k = nloc_dmg%IADC(4,ii)
480 nloc_dmg%FSKY(k,j) = -f4(i,j)
481 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
482c
483 ENDDO
484c
485 ENDDO
486c
487 ENDIF
488c
489 !-----------------------------------------------------------------------
490 ! Computing non-local timestep
491 !-----------------------------------------------------------------------
492 IF (nodadt == 0) THEN
493 DO i = 1,nel
494 ! If the element is not broken, normal computation
495 IF (off(i) /= zero) THEN
496 ! Non-local critical time-step in the plane
497 dtnl = (two*(min(le(i),le_max))*sqrt(three*zeta))/
498 . sqrt(twelve*l2 + (min(le(i),le_max))**2)
499 ! Non-local critical time-step in the thickness
500 IF (nddl>1) THEN
501 IF (nddl > 2) THEN
502 dtnl_th = (two*(min(thk(i)/nddl,le_max))*sqrt(three*zeta))/
503 . sqrt(twelve*l2 + (min(thk(i)/nddl,le_max))**2)
504 ELSE
505 dtnl_th = (two*(min(thk(i),le_max))*sqrt(three*zeta))/
506 . sqrt(twelve*l2 + (min(thk(i),le_max))**2)
507 ENDIF
508 ELSE
509 dtnl_th = ep20
510 ENDIF
511 ! Retaining the minimal value
512 dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl_th,dtfac1(1)*cdamp*dtnl)
513 ENDIF
514 ENDDO
515 ENDIF
516c
517 ! Deallocation of tables
518 IF (ALLOCATED(f1)) DEALLOCATE(f1)
519 IF (ALLOCATED(f2)) DEALLOCATE(f2)
520 IF (ALLOCATED(f3)) DEALLOCATE(f3)
521 IF (ALLOCATED(f4)) DEALLOCATE(f4)
522 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
523 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
524 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
525 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
526 IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
527 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
528 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
529 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
530 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
531 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
532 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
533 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
534 IF (ALLOCATED(vol)) DEALLOCATE(vol)
535c
536 END
subroutine cfint_reg(nloc_dmg, var_reg, thk, nel, off, area, nc1, nc2, nc3, nc4, px1, py1, px2, py2, bufnl, imat, nddl, itask, dt2t, le, thk0, area0, nft)
Definition cfint_reg.F:39
subroutine area(d1, x, x2, y, y2, eint, stif0)
#define min(a, b)
Definition macros.h:20
#define max(a, b)
Definition macros.h:21