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