OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbafint_reg.F File Reference
#include "implicit_f.inc"
#include "parit_c.inc"
#include "com20_c.inc"
#include "com08_c.inc"
#include "scr02_c.inc"
#include "scr18_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ cbafint_reg()

subroutine cbafint_reg ( type(nlocal_str_), target nloc_dmg,
intent(inout) var_reg,
intent(in) thk,
integer nel,
intent(in) off,
intent(in) area,
integer, dimension(nel) nc1,
integer, dimension(nel) nc2,
integer, dimension(nel) nc3,
integer, dimension(nel) nc4,
type(buf_nloc_), target bufnl,
integer imat,
integer nddl,
integer itask,
integer ng,
integer jft,
integer jlt,
intent(in) x13,
intent(in) y13,
intent(in) x24,
intent(in) y24,
intent(inout) dt2t,
intent(in) thk0,
intent(in) area0,
integer, intent(in) nft )

Definition at line 31 of file cbafint_reg.F.

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