OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
chvis3.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!|| chvis3 ../engine/source/elements/shell/coque/chvis3.F
25!||--- called by ------------------------------------------------------
26!|| cforc3 ../engine/source/elements/shell/coque/cforc3.F
27!|| cforc3_crk ../engine/source/elements/xfem/cforc3_crk.F
28!||====================================================================
29 SUBROUTINE chvis3(
30 2 JFT ,JLT ,THK ,HOUR,OFF ,PX1 ,PX2 ,PY1 ,PY2 ,
31 3 IXC ,DT1C,SSP ,RHO ,STI ,VX1 ,VX2 ,VX3 ,VX4 ,VY1 ,
32 4 VY2 ,VY3 ,VY4 ,VZ1 ,VZ2 ,VZ3 ,VZ4 ,AREA,THK0,VHX ,
33 5 VHY ,SHF ,Z2 ,EANI,STIR,VISCMX,G ,A11 ,
34 6 H1 ,H2 ,H3 ,YM ,NU ,THK02 ,ALPE ,H11 ,
35 7 H12 ,H13 ,H21 ,H22 ,H23 ,H31 ,H32 ,H33 ,
36 8 B11 ,B12 ,B13 ,B14 ,B21 ,B22 ,B23 ,B24 ,
37 9 RX1 ,RX2 ,RX3 ,RX4 ,RY1 ,RY2 ,RY3 ,RY4 ,
38 A IPARTC,PARTSAV,
39 B IHBE ,NFT ,ISMSTR,KFTS ,
40 C SRH1 ,SRH2 ,SRH3 ,IGTYP ,
41 D IGMAT ,A11R,NEL)
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46#include "comlock.inc"
47C-----------------------------------------------
48C G l o b a l P a r a m e t e r s
49C-----------------------------------------------
50#include "mvsiz_p.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "com04_c.inc"
55#include "com06_c.inc"
56#include "param_c.inc"
57#include "remesh_c.inc"
58#include "scr02_c.inc"
59#include "scr06_c.inc"
60#include "scr14_c.inc"
61#include "scr16_c.inc"
62#include "scr17_c.inc"
63#include "sms_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER IXC(NIXC,*),IPARTC(*), JFT, JLT,
68 . IHBE ,NFT ,ISMSTR,KFTS,IGTYP,
69 . IGMAT,NEL
70C REAL
71 my_real
72 . THK(*), HOUR(NEL,5), OFF(*),PARTSAV(NPSAV,*),
73 . PX1(*), PX2(*), PY1(*), PY2(*),DT1C(*),EANI(*),
74 . SSP(MVSIZ), RHO(MVSIZ),STI(MVSIZ),STIR(MVSIZ),
75 . VX1(MVSIZ), VX2(MVSIZ), VX3(MVSIZ), VX4(MVSIZ),
76 . VY1(MVSIZ), VY2(MVSIZ), VY3(MVSIZ), VY4(MVSIZ),
77 . VZ1(MVSIZ), VZ2(MVSIZ), VZ3(MVSIZ), VZ4(MVSIZ),
78 . AREA(MVSIZ),THK0(MVSIZ),VHX(MVSIZ), VHY(MVSIZ),
79 . shf(mvsiz),z2(mvsiz),viscmx(mvsiz),g(mvsiz),
80 . h1(mvsiz), h2(mvsiz), h3(mvsiz), a11(mvsiz),
81 . b11(mvsiz), b12(mvsiz), b13(mvsiz), b14(mvsiz),
82 . b21(mvsiz), b22(mvsiz), b23(mvsiz), b24(mvsiz),
83 . h11(mvsiz), h12(mvsiz), h13(mvsiz), h14(mvsiz),
84 . h21(mvsiz), h22(mvsiz), h23(mvsiz), h24(mvsiz),
85 . h31(mvsiz), h32(mvsiz), h33(mvsiz), h34(mvsiz),
86 . rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz),
87 . ry1(mvsiz), ry2(mvsiz), ry3(mvsiz), ry4(mvsiz),
88 . thk02(mvsiz),ym(mvsiz), nu(mvsiz), alpe(mvsiz),
89 . srh1(*) ,srh2(*) ,srh3(*),a11r(*)
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I,MX,J, II, IC, JST(MVSIZ+1)
94C REAL
95 my_real
96 . H1L(MVSIZ), H2L(MVSIZ), H3L(MVSIZ), H1Q(MVSIZ), H2Q(MVSIZ),
97 . H3Q(MVSIZ), HG1(MVSIZ), HG2(MVSIZ),
98 . HH1(MVSIZ), HH2(MVSIZ), HH3(MVSIZ),
99 . A1(MVSIZ), A2(MVSIZ), A3(MVSIZ), A4(MVSIZ),
100 . A5(MVSIZ), A6(MVSIZ), A7(MVSIZ), A8(MVSIZ),
101 . B1(MVSIZ), B2(MVSIZ), EHOU(MVSIZ),
102 . HOUR1(MVSIZ), HOUR2(MVSIZ), HOUR3(MVSIZ),
103 . GAMA1(MVSIZ), GAMA2(MVSIZ), GAMA3(MVSIZ), GAMA4(MVSIZ)
104 my_real
105 . px1v, px2v, py1v, py2v, ssphvl, hvish1, hvish2,zz2,
106 . shfpr3, fac, ehourt,
107 . r0,r1,t2a,tsa,td,vv,mm,hour1p,hour2p,hour3p,
108 . px1vp,py1vp,px2vp,py2vp,hour1a,hour2a,hour3a,tsaphi,
109 . tdphi,shfpr3d,hvish1d,hvish2d,r0d,r1d,hh3d,
110 . sr2d2,srshfpr3, inv9, inv12,scale(mvsiz),iz
111C=======================================================================
112 ehourt = zero
113 sr2d2 = sqrt(two)* half
114 inv9 = one_over_9
115 inv12 = one_over_12
116
117!! IGTYP = IGEO(11,MAT(1))
118C
119 IF(ismstr/=1.AND.ismstr/=11.AND.ihbe>=1)THEN
120 DO i=jft,jlt
121 px1v = px1(i)*vhx(i)
122 px2v = px2(i)*vhx(i)
123 py1v = py1(i)*vhy(i)
124 py2v = py2(i)*vhy(i)
125 gama1(i)= off(i)*( one - px1v-py1v)
126 gama3(i)= off(i)*( one + px1v+py1v)
127 gama2(i)= off(i)*(-one - px2v-py2v)
128 gama4(i)= off(i)*(-one + px2v+py2v)
129 ENDDO
130 ELSE
131 DO i=jft,jlt
132 gama1(i)= off(i)
133 gama3(i)= off(i)
134 gama2(i)= -off(i)
135 gama4(i)= -off(i)
136 ENDDO
137 ENDIF
138C----------------------------------------
139 DO i=jft,jlt
140 shfpr3 = shf(i)/(three*(one + nu(i)))
141 hvish1 = hvisc*h1(i)
142 hvish2 = hvisc*h2(i)
143 r0 = fourth*rho(i)
144 r1 = r0*hundred
145 r0 = r0*hvlin
146 a1(i) = r1*hvish1
147 a2(i) = r0*sr2d2*srh1(i)
148 srshfpr3 = sqrt(shfpr3)
149 a3(i) = r1*hvish2*srshfpr3
150 a4(i) = r0*sr2d2*srh2(i)*srshfpr3
151 hh3(i) = helas*h3(i)
152 a5(i) = hh3(i)*r1*zep072169
153 hh3(i) = sr2d2*srh3(i)
154 a6(i) = hh3(i)*r0*zep072169
155 r0 = fourth*ym(i)*helas
156 a7(i) = h1(i)*r0
157 a8(i) = h2(i)*r0*shfpr3
158 ENDDO
159C-------------------------------------
160 DO i=jft,jlt
161 t2a = thk02(i)*area(i)
162 tsa = sqrt(t2a)
163 h1q(i) = a1(i)*tsa
164 h1l(i) = a2(i)*ssp(i)*tsa
165 h2q(i) = a3(i)*thk02(i)
166 h2l(i) = a4(i)*ssp(i)*thk02(i)
167 h3q(i) = a5(i)*t2a
168 h3l(i) = a6(i)*ssp(i)*t2a
169 td = thk0(i)*dt1c(i)
170 hh1(i) = a7(i)*td
171 b1(i) = px1(i)*px1(i)+py1(i)*py1(i)
172 b2(i) = px2(i)*px2(i)+py2(i)*py2(i)
173 hh2(i) = a8(i)*thk02(i)*td/(b1(i)+b2(i))
174 ENDDO
175C-----------------------------------
176C TRIANGLES
177C-----------------------------------
178 DO i=jft,jlt
179 IF(ixc(4,i)==ixc(5,i))THEN
180 h1q(i)=zero
181 h1l(i)=zero
182 h2q(i)=zero
183 h2l(i)=zero
184 h3q(i)=zero
185 h3l(i)=zero
186 hh1(i)=zero
187 hh2(i)=zero
188 END IF
189 END DO
190C-----------------------------------
191C STIFFNESS - DT
192C-----------------------------------
193 DO i=jft,jlt
194 scale(i) = zero
195 ENDDO
196 IF(nodadt/=0.OR.idt1sh==1.OR.idtmins==2)THEN
197 DO i=jft,jlt
198 scale(i) = max(gama1(i)*gama1(i),gama2(i)*gama2(i),
199 . gama3(i)*gama3(i),gama4(i)*gama4(i)) *
200 . dt1c(i)*max(hh1(i)+h1l(i),hh2(i)+h2l(i),h3l(i))
201 . /max(dt1c(i)*dt1c(i),em20)
202 sti(i) = sti(i) + scale(i)
203 ENDDO
204 !! IF(IDTMINS==2 )IGMAT = 1
205C
206 IF(igtyp == 52 .OR.
207 . ((igtyp == 11 .OR. igtyp == 17 .OR. igtyp == 51)
208 . .AND. igmat > 0 )) THEN
209 IF(nadmesh==0)THEN
210 DO i=jft,jlt
211 IF (off(i)==zero) THEN
212 sti(i) = zero
213 stir(i) = zero
214 ELSE
215 vv = viscmx(i) * viscmx(i) * alpe(i)
216 fac = max(b1(i),b2(i))/(area(i) * vv)
217 sti(i) = sti(i) + fac * thk0(i) * a11(i)
218 stir(i) = fac * a11r(i)*one_over_12*thk0(i)**3 +
219 . fac * a11(i)*thk0(i)*area(i)*inv9 +
220 . fac*scale(i)*(one_over_12*thk0(i)**2 + area(i)*inv9)
221 ENDIF
222 ENDDO
223 ELSE
224 DO i=jft,jlt
225 IF (off(i)==zero) THEN
226 sti(i) = zero
227 stir(i) = zero
228 ELSE
229 vv = viscmx(i) * viscmx(i) * alpe(i)
230 fac = max(b1(i),b2(i))/(area(i) * vv)
231 sti(i) = sti(i) + fac * thk0(i) * a11(i)
232 stir(i) = fac * a11r(i)*one_over_12*thk0(i)**3 +
233 . fac*scale(i)*one_over_12*thk0(i)**2
234 ENDIF
235 ENDDO
236 END IF
237 ELSE
238 IF(nadmesh==0)THEN
239 DO i=jft,jlt
240 IF (off(i)==zero) THEN
241 sti(i) = zero
242 stir(i) = zero
243 ELSE
244 vv = viscmx(i) * viscmx(i) * alpe(i)
245 sti(i) = sti(i) + max(b1(i),b2(i))
246 . * thk0(i) * a11(i) / (area(i) * vv)
247 stir(i) = sti(i) * (thk02(i)*inv12 + area(i)*inv9)
248c .. + 0.5 * SHF(I) * AREA(I) * G(I)/A11(I))
249 ENDIF
250 ENDDO
251 ELSE
252 DO i=jft,jlt
253 IF (off(i)==zero) THEN
254 sti(i) = zero
255 stir(i) = zero
256 ELSE
257 vv = viscmx(i) * viscmx(i) * alpe(i)
258 sti(i) = sti(i) + max(b1(i),b2(i))
259 . * thk0(i) * a11(i) / (area(i) * vv)
260 stir(i) = sti(i) * thk02(i)*inv12
261 ENDIF
262 ENDDO
263 END IF
264 endif! IGTYP
265 ENDIF
266C-----------------------------------
267C ANTI-HOURGLASS MEMBRANE FORCES
268C-----------------------------------
269 IF(ismstr==1.OR.ismstr==11.OR.ihbe<1)THEN
270 DO i=jft,jlt
271 hg1(i)=(vx1(i)-vx2(i)+vx3(i)-vx4(i))*off(i)
272 hg2(i)=(vy1(i)-vy2(i)+vy3(i)-vy4(i))*off(i)
273 ENDDO
274 ELSE
275 DO i=jft,jlt
276 hg1(i)=vx1(i)*gama1(i)+vx2(i)*gama2(i)
277 . +vx3(i)*gama3(i)+vx4(i)*gama4(i)
278 hg2(i)=vy1(i)*gama1(i)+vy2(i)*gama2(i)
279 . +vy3(i)*gama3(i)+vy4(i)*gama4(i)
280 ENDDO
281 ENDIF
282 DO i=jft,jlt
283 hour(i,1)=hour(i,1)+hg1(i)*hh1(i)
284 hour(i,2)=hour(i,2)+hg2(i)*hh1(i)
285 hour1a =hour(i,1)+hg1(i)*(h1l(i)+h1q(i)*abs(hg1(i)))
286 h11(i)=hour1a*gama1(i)
287 h12(i)=hour1a*gama2(i)
288 h13(i)=hour1a*gama3(i)
289C H14(I)=HOUR1A*GAMA4(I)
290 hour2a =hour(i,2)+hg2(i)*(h1l(i)+h1q(i)*abs(hg2(i)))
291 h21(i)=hour2a*gama1(i)
292 h22(i)=hour2a*gama2(i)
293 h23(i)=hour2a*gama3(i)
294C H24(I)=HOUR2A*GAMA4(I)
295 ehou(i) = hour1a*hg1(i) + hour2a*hg2(i)
296 ENDDO
297C
298C----------------------------------
299C ANTI-HOURGLASS BENDING FORCES
300C----------------------------------
301 IF(ismstr==1.OR.ismstr==11.OR.ihbe<1)THEN
302 DO i=jft,jlt
303 hg1(i)=(vz1(i)-vz2(i)+vz3(i)-vz4(i))*off(i)
304 ENDDO
305 ELSE
306 DO i=jft,jlt
307 hg1(i)=vz1(i)*gama1(i)+vz2(i)*gama2(i)
308 . +vz3(i)*gama3(i)+vz4(i)*gama4(i)
309 ENDDO
310 ENDIF
311 DO i=jft,jlt
312 hour(i,3)=hour(i,3)+hg1(i)*hh2(i)
313 hour3a =hour(i,3)+hg1(i)*(h2l(i)+h2q(i)*abs(hg1(i)))
314 h31(i)=hour3a*gama1(i)
315 h32(i)=hour3a*gama2(i)
316 h33(i)=hour3a*gama3(i)
317C H34(I)=HOUR3(I)*GAMA4(I)
318 ehou(i) = ehou(i) + hour3a*hg1(i)
319 ENDDO
320C
321C------------------------
322C ANTI-HOURGLASS MOMENTS
323C------------------------
324 DO i=jft,jlt
325 hg1(i)=rx1(i)-rx2(i)+rx3(i)-rx4(i)
326 hg2(i)=ry1(i)-ry2(i)+ry3(i)-ry4(i)
327 ENDDO
328C
329 DO i=jft,jlt
330 hour(i,4)=hg1(i)*(h3l(i)+h3q(i)*abs(hg1(i)))
331 hour(i,5)=hg2(i)*(h3l(i)+h3q(i)*abs(hg2(i)))
332 ehou(i) = ehou(i) +
333 . hour(i,4)*hg1(i) + hour(i,5)*hg2(i)
334 ehou(i) = dt1c(i) * ehou(i) * off(i)
335C EHOURT = EHOURT + EHOU(I)
336 b11(i)= hour(i,4)*off(i)
337 b12(i)=-hour(i,4)*off(i)
338 b13(i)= hour(i,4)*off(i)
339 b14(i)=-hour(i,4)*off(i)
340 b21(i)= hour(i,5)*off(i)
341 b22(i)=-hour(i,5)*off(i)
342 b23(i)= hour(i,5)*off(i)
343 b24(i)=-hour(i,5)*off(i)
344 ENDDO
345C separate loop that does not vectorize w/o simd directive
346 DO i=jft,jlt
347 ehourt = ehourt + ehou(i)
348 ENDDO
349C
350 ic=1
351 jst(ic)=jft
352 DO j=jft+1,jlt
353 IF (ipartc(j)/=ipartc(j-1)) THEN
354 ic=ic+1
355 jst(ic)=j
356 ENDIF
357 ENDDO
358
359 jst(ic+1)=jlt+1
360 IF(ic==1) THEN
361 mx = ipartc(jft)
362 DO i=jft,jlt
363 partsav(8,mx)=partsav(8,mx) + ehou(i)
364 ENDDO
365
366 ELSEIF(ic==2.AND.kfts>0)THEN
367 mx = ipartc(jft)
368 DO i=jft, kfts-1
369 partsav(8,mx)=partsav(8,mx) + ehou(i)
370 ENDDO
371 mx = ipartc(jlt)
372 DO i=kfts,jlt
373 partsav(8,mx)=partsav(8,mx) + ehou(i)
374 ENDDO
375
376 ELSE
377C
378 DO ii=1,ic
379 mx=ipartc(jst(ii))
380 IF (jst(ii+1)-jst(ii)>15) THEN
381 DO j=jst(ii),jst(ii+1)-1
382 partsav(8,mx)=partsav(8,mx)+ehou(j)
383 ENDDO
384 ELSE
385 DO j=jst(ii),jst(ii+1)-1
386 partsav(8,mx)=partsav(8,mx)+ehou(j)
387 ENDDO
388 ENDIF
389 ENDDO
390 ENDIF
391C
392!$OMP ATOMIC
393 ehour = ehour + ehourt
394C
395 DO i=jft,jlt
396 eani(nft+numels+i) = eani(nft+numels+i) + ehou(i)
397 ENDDO
398C
399 RETURN
400 END
chvis3
subroutine chvis3(jft, jlt, thk, hour, off, px1, px2, py1, py2, ixc, dt1c, ssp, rho, sti, vx1, vx2, vx3, vx4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, area, thk0, vhx, vhy, shf, z2, eani, stir, viscmx, g, a11, h1, h2, h3, ym, nu, thk02, alpe, h11, h12, h13, h21, h22, h23, h31, h32, h33, b11, b12, b13, b14, b21, b22, b23, b24, rx1, rx2, rx3, rx4, ry1, ry2, ry3, ry4, ipartc, partsav, ihbe, nft, ismstr, kfts, srh1, srh2, srh3, igtyp, igmat, a11r, nel)
Definition chvis3.F:42
max
#define max(a, b)
Definition macros.h:21