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