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