OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
clusterf.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!|| clusterf ../engine/source/output/cluster/clusterf.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- uses -----------------------------------------------------
28!|| cluster_mod ../engine/share/modules/cluster_mod.F
29!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
30!|| h3d_mod ../engine/share/modules/h3d_mod.F
31!||====================================================================
32 SUBROUTINE clusterf(CLUSTER ,ELBUF_TAB,X ,A ,AR ,
33 . SKEW ,IXS ,IPARG ,FCLUSTER,MCLUSTER ,
34 . H3D_DATA ,GEO )
35C-----------------------------------------------
36C M o d u l e s
37C-----------------------------------------------
38 USE elbufdef_mod
39 USE cluster_mod
40 USE h3d_mod
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "param_c.inc"
49#include "units_c.inc"
50#include "comlock.inc"
51#include "com01_c.inc"
52#include "com04_c.inc"
53#include "com08_c.inc"
54#include "scr14_c.inc"
55#include "scr17_c.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER IXS(NIXS,*),IPARG(NPARG,*)
60 my_real ,DIMENSION(3,*) :: X,A,AR,FCLUSTER,MCLUSTER
61 my_real ,DIMENSION(LSKEW,*) :: skew
62 TYPE (CLUSTER_) ,DIMENSION(NCLUSTER) :: CLUSTER
63 TYPE (ELBUF_STRUCT_),DIMENSION(NGROUP) :: ELBUF_TAB
64 TYPE (H3D_DATABASE) :: H3D_DATA
65C-----------------------------------------------
66C L o c a l V a r i a b l e s
67C-----------------------------------------------
68 INTEGER I,J,K,IL,IEL,NG,NFT,NNOD,ISKN,N,N1,N2,N3,N4,NINDX,IFAIL,IPID
69 INTEGER CLUSTERNOD(NCLUSTER),LCLUSTER(NCLUSTER),LCL(NCLUSTER),
70 . iad(ncluster)
71 INTEGER INDX(NCLUSTER)
72 my_real, DIMENSION(NPROPG,*) :: GEO
73 my_real, DIMENSION(3) :: fbot,ftop,mbot,mtop,m1,xg,x1,x2
74 my_real, DIMENSION(3,NCLUSTER) :: vn,vx,vy
75 my_real :: fn,ft,mr,mb,dmg,xm,ym,zm,dx1,dy1,dz1,dx2,dy2,dz2,
76 . fx,fy,fz,momx,momy,momz,norm,critf,critm,drx,dry,drz,
77 . sx,sy,sz,tx,ty,tz
78 my_real, DIMENSION(NCLUSTER) :: tthick
79C=======================================================================
80 nindx = 0
81 tthick(1:ncluster) = zero
82c
83 DO i = 1, ncluster
84 IF (cluster(i)%OFF == 0) cycle
85 nnod = cluster(i)%NNOD
86 iskn = cluster(i)%SKEW
87 ifail= cluster(i)%IFAIL
88c------
89c center of moments on lower face
90c------
91 x1(1:3) = zero
92 DO j = 1,nnod
93 n1 = cluster(i)%NOD1(j)
94 x1(1) = x1(1) + x(1,n1)
95 x1(2) = x1(2) + x(2,n1)
96 x1(3) = x1(3) + x(3,n1)
97 END DO
98 xm = x1(1) / nnod
99 ym = x1(2) / nnod
100 zm = x1(3) / nnod
101c------
102c local skew
103c------
104 IF (ifail > 0 .and. iskn == 0) THEN ! local skew fixed on bottom surf
105c
106c calculate normal direction of the local skew
107 vn(1,i) = zero
108 vn(2,i) = zero
109 vn(3,i) = zero
110c
111 IF (cluster(i)%TYPE == 1) THEN ! cohesive 3D elements
112 DO j = 1,cluster(i)%NEL
113 ng = cluster(i)%NG(j)
114 iel = cluster(i)%ELEM(j)
115 nft = iparg(3,ng)
116 ipid = ixs(10,nft+iel)
117 n1 = ixs(2,nft+iel)
118 n2 = ixs(3,nft+iel)
119 n3 = ixs(4,nft+iel)
120 n4 = ixs(5,nft+iel)
121 tthick(i) = geo(41,ipid)
122 sx = x(1,n3) - x(1,n1)
123 sy = x(2,n3) - x(2,n1)
124 sz = x(3,n3) - x(3,n1)
125 tx = x(1,n4) - x(1,n2)
126 ty = x(2,n4) - x(2,n2)
127 tz = x(3,n4) - x(3,n2)
128 vn(1,i) = vn(1,i) + sy*tz - sz*ty
129 vn(2,i) = vn(2,i) + sz*tx - sx*tz
130 vn(3,i) = vn(3,i) + sx*ty - sy*tx
131 END DO
132c
133 ELSE ! spring cluster
134 n1 = cluster(i)%NOD1(nnod)
135 n2 = cluster(i)%NOD1(1)
136 sx = xm - x(1,n1)
137 sy = ym - x(2,n1)
138 sz = zm - x(3,n1)
139 tx = xm - x(1,n2)
140 ty = ym - x(2,n2)
141 tz = zm - x(3,n2)
142 vn(1,i) = vn(1,i) + sy*tz - sz*ty
143 vn(2,i) = vn(2,i) + sz*tx - sx*tz
144 vn(3,i) = vn(3,i) + sx*ty - sy*tx
145 DO j = 1,nnod-1
146 n1 = cluster(i)%NOD1(j)
147 n2 = cluster(i)%NOD1(j+1)
148 sx = xm - x(1,n1)
149 sy = ym - x(2,n1)
150 sz = zm - x(3,n1)
151 tx = xm - x(1,n2)
152 ty = ym - x(2,n2)
153 tz = zm - x(3,n2)
154 vn(1,i) = vn(1,i) + sy*tz - sz*ty
155 vn(2,i) = vn(2,i) + sz*tx - sx*tz
156 vn(3,i) = vn(3,i) + sx*ty - sy*tx
157 END DO
158 END IF ! cluster type
159c
160 norm = one / sqrt(vn(1,i)**2 + vn(2,i)**2 + vn(3,i)**2)
161 vn(1,i) = vn(1,i)*norm
162 vn(2,i) = vn(2,i)*norm
163 vn(3,i) = vn(3,i)*norm
164c
165c calculate X and Y directions of the local skew
166c
167 n1 = cluster(i)%NOD1(1)
168 n2 = cluster(i)%NOD1(2)
169 vx(1,i) = x(1,n1) - xm
170 vx(2,i) = x(2,n1) - ym
171 vx(3,i) = x(3,n1) - zm
172 vy(1,i) = vn(2,i)*vx(3,i) - vn(3,i)*vx(2,i)
173 vy(2,i) = vn(3,i)*vx(1,i) - vn(1,i)*vx(3,i)
174 vy(3,i) = vn(1,i)*vx(2,i) - vn(2,i)*vx(1,i)
175 norm = one / sqrt(vy(1,i)**2 + vy(2,i)**2 + vy(3,i)**2)
176 vy(1,i) = vy(1,i)*norm
177 vy(2,i) = vy(2,i)*norm
178 vy(3,i) = vy(3,i)*norm
179 vx(1,i) = vy(2,i)*vn(3,i) - vy(3,i)*vn(2,i)
180 vx(2,i) = vy(3,i)*vn(1,i) - vy(1,i)*vn(3,i)
181 vx(3,i) = vy(1,i)*vn(2,i) - vy(2,i)*vn(1,i)
182 norm = one / sqrt(vx(1,i)**2 + vx(2,i)**2 + vx(3,i)**2)
183 vx(1,i) = vx(1,i)*norm
184 vx(2,i) = vx(2,i)*norm
185 vx(3,i) = vx(3,i)*norm
186c
187 ENDIF ! IFAIL > 0 .and. ISKN == 0
188c------
189c Forces
190c------
191 fbot = zero
192 ftop = zero
193 DO j = 1,nnod
194 n1 = cluster(i)%NOD1(j)
195 n2 = cluster(i)%NOD2(j)
196 fbot(1) = fbot(1) + a(1,n1)
197 fbot(2) = fbot(2) + a(2,n1)
198 fbot(3) = fbot(3) + a(3,n1)
199 ftop(1) = ftop(1) + a(1,n2)
200 ftop(2) = ftop(2) + a(2,n2)
201 ftop(3) = ftop(3) + a(3,n2)
202 END DO
203c------
204c Moments
205c------
206 mbot = zero
207 mtop = zero
208c
209 IF (cluster(i)%TYPE == 1 .and. iskn == 0 .and. tthick(i) > zero) THEN
210 DO j = 1,nnod
211 n1 = cluster(i)%NOD1(j)
212 n2 = cluster(i)%NOD2(j)
213
214 drx = x(1,n2) - xm
215 dry = x(2,n2) - ym
216 drz = sign(tthick(i), x(3,n2) - zm)
217 mtop(1) = mtop(1) + dry*a(3,n2) - drz*a(2,n2)
218 mtop(2) = mtop(2) + drz*a(1,n2) - drx*a(3,n2)
219 mtop(3) = mtop(3) + drx*a(2,n2) - dry*a(1,n2)
220c
221 drx = x(1,n1) - xm
222 dry = x(2,n1) - ym
223 mbot(1) = mbot(1) + dry*a(3,n1)
224 mbot(2) = mbot(2) - drx*a(3,n1)
225 mbot(3) = mbot(3) + drx*a(2,n1) - dry*a(1,n1)
226 END DO ! NNOD
227 ELSE
228 DO j = 1,nnod
229 n1 = cluster(i)%NOD1(j)
230 n2 = cluster(i)%NOD2(j)
231
232 drx = x(1,n2) - xm
233 dry = x(2,n2) - ym
234 drz = x(3,n2) - zm
235 mtop(1) = mtop(1) + dry*a(3,n2) - drz*a(2,n2)
236 mtop(2) = mtop(2) + drz*a(1,n2) - drx*a(3,n2)
237 mtop(3) = mtop(3) + drx*a(2,n2) - dry*a(1,n2)
238c
239 drx = x(1,n1) - xm
240 dry = x(2,n1) - ym
241 drz = x(3,n1) - zm
242 mbot(1) = mbot(1) + dry*a(3,n1) - drz*a(2,n1)
243 mbot(2) = mbot(2) + drz*a(1,n1) - drx*a(3,n1)
244 mbot(3) = mbot(3) + drx*a(2,n1) - dry*a(1,n1)
245 END DO ! NNOD
246 END IF
247c------
248 IF (cluster(i)%TYPE == 1) THEN ! Brick cluster
249 fx = (ftop(1) - fbot(1))*half
250 fy = (ftop(2) - fbot(2))*half
251 fz = (ftop(3) - fbot(3))*half
252 momx = (mtop(1) - mbot(1))*half
253 momy = (mtop(2) - mbot(2))*half
254 momz = (mtop(3) - mbot(3))*half
255 ELSE ! Spring cluster
256 fx = ftop(1)
257 fy = ftop(2)
258 fz = ftop(3)
259 momx = mtop(1)
260 momy = mtop(2)
261 momz = mtop(3)
262 DO j = 1,nnod
263 n1 = cluster(i)%NOD1(j)
264 n2 = cluster(i)%NOD2(j)
265 momx = momx + ar(1,n2)
266 momy = momy + ar(2,n2)
267 momz = momz + ar(3,n2)
268 END DO
269 ENDIF
270c
271 cluster(i)%FOR(1) = fx
272 cluster(i)%FOR(2) = fy
273 cluster(i)%FOR(3) = fz
274 cluster(i)%MOM(1) = momx
275 cluster(i)%MOM(2) = momy
276 cluster(i)%MOM(3) = momz
277c------
278 ENDDO ! I = 1, NCLUSTER
279c
280c------------------------------
281c Cluster failure
282c---------------------------------
283 nindx = 0
284 indx = 0
285c
286 DO i = 1, ncluster
287
288 IF (cluster(i)%OFF == 0) THEN
289 cluster(i)%FOR(1) = zero
290 cluster(i)%FOR(2) = zero
291 cluster(i)%FOR(3) = zero
292 cluster(i)%MOM(1) = zero
293 cluster(i)%MOM(2) = zero
294 cluster(i)%MOM(3) = zero
295 END IF
296
297 nnod = cluster(i)%NNOD
298 iskn = cluster(i)%SKEW
299 ifail= cluster(i)%IFAIL
300
301c---
302c IF (CLUSTER(I)%TYPE == 1) THEN ! check local failure
303cc IF (CLUSTER(I)%IFAIL >= 10) THEN
304cc check local element failure : /FAIL
305cc IFAIL = IFAIl - 10
306cc
307c NOFF = 0
308c DO J = 1,CLUSTER(I)%NEL
309c NG = CLUSTER(I)%NG(J)
310c IEL = CLUSTER(I)%ELEM(J)
311c IF (ELBUF_TAB(NG)%GBUF%OFF(IEL) == ZERO) NOFF = NOFF + 1
312c END DO
313c IF (NOFF == CLUSTER(I)%NEL) THEN
314c CLUSTER(I)%OFF = 0
315c NINDX = NINDX+1
316c INDX(NINDX) = I
317c IDEL7NOK = 1
318c ENDIF
319c ENDIF
320c ENDIF
321c---
322 IF (ifail > 0) THEN
323 IF (iskn > 0) THEN
324 fbot(1) = cluster(i)%FOR(1)*skew(1,iskn) +
325 . cluster(i)%FOR(2)*skew(2,iskn) +
326 . cluster(i)%FOR(3)*skew(3,iskn)
327 fbot(2) = cluster(i)%FOR(1)*skew(4,iskn) +
328 . cluster(i)%FOR(2)*skew(5,iskn) +
329 . cluster(i)%FOR(3)*skew(6,iskn)
330 fbot(3) = cluster(i)%FOR(1)*skew(7,iskn) +
331 . cluster(i)%FOR(2)*skew(8,iskn) +
332 . cluster(i)%FOR(3)*skew(9,iskn)
333 m1(1) = cluster(i)%MOM(1)*skew(1,iskn) +
334 . cluster(i)%MOM(2)*skew(2,iskn) +
335 . cluster(i)%MOM(3)*skew(3,iskn)
336 m1(2) = cluster(i)%MOM(1)*skew(4,iskn) +
337 . cluster(i)%MOM(2)*skew(5,iskn) +
338 . cluster(i)%MOM(3)*skew(6,iskn)
339 m1(3) = cluster(i)%MOM(1)*skew(7,iskn) +
340 . cluster(i)%MOM(2)*skew(8,iskn) +
341 . cluster(i)%MOM(3)*skew(9,iskn)
342 ELSE
343 fbot(1) = cluster(i)%FOR(1)*vx(1,i) +
344 . cluster(i)%FOR(2)*vx(2,i) +
345 . cluster(i)%FOR(3)*vx(3,i)
346 fbot(2) = cluster(i)%FOR(1)*vy(1,i) +
347 . cluster(i)%FOR(2)*vy(2,i) +
348 . cluster(i)%FOR(3)*vy(3,i)
349 fbot(3) = cluster(i)%FOR(1)*vn(1,i) +
350 . cluster(i)%FOR(2)*vn(2,i) +
351 . cluster(i)%FOR(3)*vn(3,i)
352 m1(1) = cluster(i)%MOM(1)*vx(1,i) +
353 . cluster(i)%MOM(2)*vx(2,i) +
354 . cluster(i)%MOM(3)*vx(3,i)
355 m1(2) = cluster(i)%MOM(1)*vy(1,i) +
356 . cluster(i)%MOM(2)*vy(2,i) +
357 . cluster(i)%MOM(3)*vy(3,i)
358 m1(3) = cluster(i)%MOM(1)*vn(1,i) +
359 . cluster(i)%MOM(2)*vn(2,i) +
360 . cluster(i)%MOM(3)*vn(3,i)
361
362
363 ENDIF ! IF (ISKN > 0) THEN
364 fn = abs(fbot(3))
365 ft = sqrt(fbot(1)*fbot(1) + fbot(2)*fbot(2))
366 mr = abs(m1(3))
367 mb = sqrt(m1(1)*m1(1) + m1(2)*m1(2))
368 ENDIF ! IF (IFAIL > 0) THEN
369
370c---------------------------
371 IF (ifail == 1) THEN
372c Monodirectional + one direction
373 critf = max(fn/cluster(i)%FMAX(1),ft/cluster(i)%FMAX(2))
374 critm = max(mr/cluster(i)%MMAX(1),mb/cluster(i)%MMAX(2))
375 dmg = max(critf,critm)
376
377 ELSEIF (ifail == 2) THEN
378c Monodirectional + all directions
379 dmg = fourth*(min(one+em10, fn/cluster(i)%FMAX(1)) +
380 . min(one+em10, ft/cluster(i)%FMAX(2)) +
381 . min(one+em10, mr/cluster(i)%MMAX(1)) +
382 . min(one+em10, mb/cluster(i)%MMAX(2)))
383
384 ELSEIF (ifail == 3) THEN
385c Multidirectional failure
386 dmg =
387 . cluster(i)%AX(1)*(fn/cluster(i)%FMAX(1))**cluster(i)%NX(1)
388 . + cluster(i)%AX(2)*(ft/cluster(i)%FMAX(2))**cluster(i)%NX(2)
389 . + cluster(i)%AX(3)*(mr/cluster(i)%MMAX(1))**cluster(i)%NX(3)
390 . + cluster(i)%AX(4)*(mb/cluster(i)%MMAX(2))**cluster(i)%NX(4)
391 ELSE ! no fail
392 dmg = zero
393 ENDIF
394 cluster(i)%FAIL = dmg
395c---------------------------
396 IF (dmg > one) THEN
397
398 nindx = nindx+1
399 indx(nindx) = i
400 idel7nok = 1
401 cluster(i)%OFF = 0
402 cluster(i)%FOR(1) = zero
403 cluster(i)%FOR(2) = zero
404 cluster(i)%FOR(3) = zero
405 cluster(i)%MOM(1) = zero
406 cluster(i)%MOM(2) = zero
407 cluster(i)%MOM(3) = zero
408 IF (cluster(i)%TYPE == 1) THEN
409 DO j = 1,cluster(i)%NEL
410 ng = cluster(i)%NG(j)
411 iel = cluster(i)%ELEM(j)
412 elbuf_tab(ng)%GBUF%OFF(iel) = zero
413 END DO
414 ELSE ! spring cluster
415 ! set OFF flag to zero for all elements
416 ENDIF
417 ENDIF
418c
419 ENDDO ! I = 1, NCLUSTER
420c---------------------------
421 IF (anim_v(19) + h3d_data%N_VECT_CLUST_FORCE > 0) THEN
422 DO i = 1, ncluster
423 nnod = cluster(i)%NNOD
424 DO j = 1,nnod
425 n = cluster(i)%NOD1(j)
426 fcluster(1,n) = cluster(i)%FOR(1)
427 fcluster(2,n) = cluster(i)%FOR(2)
428 fcluster(3,n) = cluster(i)%FOR(3)
429 n = cluster(i)%NOD2(j)
430 fcluster(1,n) = cluster(i)%FOR(1)
431 fcluster(2,n) = cluster(i)%FOR(2)
432 fcluster(3,n) = cluster(i)%FOR(3)
433 ENDDO
434 ENDDO ! I = 1, NCLUSTER
435 ENDIF
436 IF (anim_v(20) + h3d_data%N_VECT_CLUST_MOM > 0) THEN
437 DO i = 1, ncluster
438 nnod = cluster(i)%NNOD
439 DO j = 1,nnod
440 n = cluster(i)%NOD1(j)
441 mcluster(1,n) = cluster(i)%MOM(1)
442 mcluster(2,n) = cluster(i)%MOM(2)
443 mcluster(3,n) = cluster(i)%MOM(3)
444 n = cluster(i)%NOD2(j)
445 mcluster(1,n) = cluster(i)%MOM(1)
446 mcluster(2,n) = cluster(i)%MOM(2)
447 mcluster(3,n) = cluster(i)%MOM(3)
448 ENDDO
449 ENDDO ! I = 1, NCLUSTER
450 ENDIF
451C-----------------------------------------------
452 IF (nindx > 0) THEN
453 DO j=1,nindx
454#include "lockon.inc"
455 WRITE(iout ,1000) cluster(indx(j))%ID
456 WRITE(istdo,1100) cluster(indx(j))%ID,tt
457#include "lockoff.inc"
458 END DO
459 ENDIF
460C-----------------------------------------------
461 1000 FORMAT(5x,'DELETE ELEMENT CLUSTER,ID=',i10)
462 1100 FORMAT(5x,'DELETE ELEMENT CLUSTER,ID=',i10,', AT TIME ',1pe16.9)
463C-----------
464 RETURN
465 END
clusterf
subroutine clusterf(cluster, elbuf_tab, x, a, ar, skew, ixs, iparg, fcluster, mcluster, h3d_data, geo)
Definition clusterf.F:35
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
cluster_mod
Definition cluster_mod.F:40
h3d_mod
Definition h3d_mod.F:308