OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sfor_n2s4.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!|| sfor_n2s4 ../engine/source/elements/solid/solide/sfor_n2s4.F
25!||--- called by ------------------------------------------------------
26!|| s6for_distor ../engine/source/elements/thickshell/solide6c/s6for_distor.F90
27!|| s8for_distor ../engine/source/elements/solid/solide/s8for_distor.F
28!||--- calls -----------------------------------------------------
29!|| ssort_n4 ../engine/source/elements/solid/solide/ssort_n4.f
30!||====================================================================
31 SUBROUTINE sfor_n2s4( XI, YI, ZI, STIF,
32 . X1, X2, X3, X4,
33 . Y1, Y2, Y3, Y4,
34 . Z1, Z2, Z3, Z4,
35 . VX1, VX2, VX3, VX4,
36 . VY1, VY2, VY3, VY4,
37 . VZ1, VZ2, VZ3, VZ4,
38 . FOR_T1, FOR_T2, FOR_T3, FOR_T4,
39 . FORC_N, LL , IFCTL, IFC1 ,
40 . PENMIN, PENREF, MARGE, FQMAX,
41 . STIF0, NEL , VC,E_DISTOR,
42 . DT1 )
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 "mvsiz_p.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER, INTENT (IN) :: NEL
55 INTEGER, INTENT (OUT) :: IFCTL
56 INTEGER, DIMENSION(MVSIZ), INTENT (IN) :: IFC1
57 my_real, INTENT (IN) :: FQMAX,DT1
58 my_real, DIMENSION(MVSIZ), INTENT (IN) ::
59 . X1, X2, X3, X4,
60 . Y1, Y2, Y3, Y4,
61 . Z1, Z2, Z3, Z4,
62 . VX1, VX2, VX3, VX4,
63 . VY1, VY2, VY3, VY4,
64 . VZ1, VZ2, VZ3, VZ4,
65 . xi, yi, zi, penmin,
66 . penref, ll, stif0, marge
67 my_real, DIMENSION(MVSIZ), INTENT (INOUT) :: stif
68 my_real, DIMENSION(MVSIZ,3), INTENT (INOUT) :: forc_n,
69 . for_t1, for_t2, for_t3, for_t4
70 my_real, DIMENSION(MVSIZ,3), INTENT (IN) :: vc
71 my_real, DIMENSION(NEL), INTENT (INOUT) :: e_distor
72C-----------------------------------------------
73C C o m m o n B l o c k s
74C-----------------------------------------------
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I,J,IFCTL1,ITG(MVSIZ),IFC2(MVSIZ),ie
79C 12
80 my_real
81 . nx(mvsiz),ny(mvsiz),nz(mvsiz),
82 . xsc(mvsiz),ysc(mvsiz),zsc(mvsiz),pene(mvsiz),
83 . xa(mvsiz),ya(mvsiz),za(mvsiz),hj(mvsiz,4),
84 . xb(mvsiz),yb(mvsiz),zb(mvsiz),fn(mvsiz),
85 . xc(mvsiz),yc(mvsiz),zc(mvsiz),area(mvsiz),
86 . la(mvsiz),lb(mvsiz),lc(mvsiz),fkt,
87 . rx, ry, rz, sx, sy, sz,vxc,vyc,vzc,
88 . x42,y42, z42, x31, y31, z31,fx,fy,fz,
89 . sax,say,saz,sbx,sby,sbz,scx,scy,scz,
90 . trx,try,trz,tsx,tsy,tsz,ttx,tty,ttz,
91 . tr2,ts2,tt2,aaa,bbb,vr,vs,vt,nnx,nny,nnz,
92 . xab,xbc,xca,yab,ybc,yca,zab,zbc,zca,
93 . xia, xib, xic, yia, yib, yic,
94 . zia, zib, zic, h0,norm,s2,fac,
95 . f_q,f_c,kts,zerom,tx,ty,tz,pene3,pene4,pendr,lj,
96 . dx,dy,dz,dn
97C----------------------------
98 ifc2(1:nel) = ifc1(1:nel)
99 CALL ssort_n4(xi, yi, zi , marge,
100 . x1, x2, x3, x4,
101 . y1, y2, y3, y4,
102 . z1, z2, z3, z4,
103 . ifc2, stif0, nel)
104 ifctl = 0
105 ifctl1 = 0
106C-------diff in velocity as 1er sorting
107 DO i=1,nel
108 IF (ifc2(i)==0) cycle
109 xsc(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
110 ysc(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
111 zsc(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
112 ifctl1=1
113 END DO
114C
115 IF (ifctl1==1) THEN
116 zerom = -two*em03
117 itg(1:nel) = 1
118 pene(1:nel) = zero
119 DO i=1,nel
120 IF (ifc2(i)==0) cycle
121 rx =x2(i)+x3(i)-x1(i)-x4(i)
122 ry =y2(i)+y3(i)-y1(i)-y4(i)
123 rz =z2(i)+z3(i)-z1(i)-z4(i)
124 sx =x3(i)+x4(i)-x1(i)-x2(i)
125 sy =y3(i)+y4(i)-y1(i)-y2(i)
126 sz =z3(i)+z4(i)-z1(i)-z2(i)
127 nx(i)=ry*sz - rz*sy
128 ny(i)=rz*sx - rx*sz
129 nz(i)=rx*sy - ry*sx
130 area(i) = nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i)
131 norm=one/max(em20,sqrt(area(i)))
132 nx(i)=nx(i)*norm
133 ny(i)=ny(i)*norm
134 nz(i)=nz(i)*norm
135 bbb =(xsc(i)-xi(i))*nx(i) +
136 . (ysc(i)-yi(i))*ny(i) +
137 . (zsc(i)-zi(i))*nz(i) -penmin(i)
138 pene(i) = max(zero,-bbb)
139 IF (four*area(i)<penmin(i)*ll(i))
140 . pene(i)=min(pene(i),em01*penmin(i))
141 IF (area(i)<em20) pene(i) =zero
142 ENDDO
143! 4-------3
144! |\ /|
145! | \ 3 / |
146! | \ /2 |
147! | 5 |
148! | 4/ \ |
149! | / 1 \ |
150! |/ \|
151! 1-------2
152 DO i=1,nel
153 IF(pene(i) == zero ) cycle
154 IF(ifc2(i)>=3) THEN ! degenerated quad
155 SELECT CASE (ifc2(i))
156 CASE (3)
157 xb(i) = x1(i)
158 yb(i) = y1(i)
159 zb(i) = z1(i)
160 xc(i) = x2(i)
161 yc(i) = y2(i)
162 zc(i) = z2(i)
163 xa(i) = x3(i)
164 ya(i) = y3(i)
165 za(i) = z3(i)
166 CASE (4)
167 xb(i) = x2(i)
168 yb(i) = y2(i)
169 zb(i) = z2(i)
170 xc(i) = x3(i)
171 yc(i) = y3(i)
172 zc(i) = z3(i)
173 xa(i) = x4(i)
174 ya(i) = y4(i)
175 za(i) = z4(i)
176 CASE (5)
177 xb(i) = x4(i)
178 yb(i) = y4(i)
179 zb(i) = z4(i)
180 xc(i) = x1(i)
181 yc(i) = y1(i)
182 zc(i) = z1(i)
183 xa(i) = x2(i)
184 ya(i) = y2(i)
185 za(i) = z2(i)
186 CASE (6)
187 xb(i) = x3(i)
188 yb(i) = y3(i)
189 zb(i) = z3(i)
190 xc(i) = x4(i)
191 yc(i) = y4(i)
192 zc(i) = z4(i)
193 xa(i) = x1(i)
194 ya(i) = y1(i)
195 za(i) = z1(i)
196 END SELECT
197! pene(i) = em01*pene(i)
198 ELSE
199 xa(i) = xsc(i)
200 ya(i) = ysc(i)
201 za(i) = zsc(i)
202C-------- if other sub-tria recompute normal
203 bbb =(x3(i)-xi(i))*nx(i) +
204 . (y3(i)-yi(i))*ny(i) +
205 . (z3(i)-zi(i))*nz(i) -penmin(i)
206 pene3 = max(zero,-bbb)
207 bbb =(x4(i)-xi(i))*nx(i) +
208 . (y4(i)-yi(i))*ny(i) +
209 . (z4(i)-zi(i))*nz(i) -penmin(i)
210 pene4 = max(zero,-bbb)
211 IF (pene3>pene(i).AND.pene4>pene(i)) THEN
212 itg(i) = 3
213 xb(i) = x3(i)
214 yb(i) = y3(i)
215 zb(i) = z3(i)
216 xc(i) = x4(i)
217 yc(i) = y4(i)
218 zc(i) = z4(i)
219 ELSEIF (pene3>pene(i)) THEN
220 itg(i) = 2
221 xb(i) = x2(i)
222 yb(i) = y2(i)
223 zb(i) = z2(i)
224 xc(i) = x3(i)
225 yc(i) = y3(i)
226 zc(i) = z3(i)
227 ELSEIF (pene4>pene(i)) THEN
228 itg(i) = 4
229 xb(i) = x4(i)
230 yb(i) = y4(i)
231 zb(i) = z4(i)
232 xc(i) = x1(i)
233 yc(i) = y1(i)
234 zc(i) = z1(i)
235 ELSE
236 xb(i) = x1(i)
237 yb(i) = y1(i)
238 zb(i) = z1(i)
239 xc(i) = x2(i)
240 yc(i) = y2(i)
241 zc(i) = z2(i)
242 END IF
243 END IF ! IFC2(I)>=3
244 ENDDO
245 DO i=1,nel
246 IF(pene(i) == zero.OR.itg(i)==1) cycle
247 rx =xb(i)-xa(i)
248 ry =yb(i)-ya(i)
249 rz =zb(i)-za(i)
250 sx =xc(i)-xa(i)
251 sy =yc(i)-ya(i)
252 sz =zc(i)-za(i)
253 nx(i)=ry*sz - rz*sy
254 ny(i)=rz*sx - rx*sz
255 nz(i)=rx*sy - ry*sx
256 area(i) = nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i)
257 norm=one/max(em20,sqrt(area(i)))
258 nx(i)=nx(i)*norm
259 ny(i)=ny(i)*norm
260 nz(i)=nz(i)*norm
261 bbb = (x3(i)-xi(i))*nx(i) +
262 . (y3(i)-yi(i))*ny(i) +
263 . (z3(i)-zi(i))*nz(i) -penmin(i)
264 pene(i) = max(zero,-bbb)
265 IF (four*area(i)<penmin(i)*ll(i))
266 . pene(i)=min(pene(i),em01*penmin(i))
267 ENDDO
268 DO i=1,nel
269 IF(pene(i) == zero) cycle
270 xab = xb(i)-xa(i)
271 yab = yb(i)-ya(i)
272 zab = zb(i)-za(i)
273 xbc = xc(i)-xb(i)
274 ybc = yc(i)-yb(i)
275 zbc = zc(i)-zb(i)
276 xca = xa(i)-xc(i)
277 yca = ya(i)-yc(i)
278 zca = za(i)-zc(i)
279
280 xia = xa(i)-xi(i)
281 yia = ya(i)-yi(i)
282 zia = za(i)-zi(i)
283 xib = xb(i)-xi(i)
284 yib = yb(i)-yi(i)
285 zib = zb(i)-zi(i)
286 xic = xc(i)-xi(i)
287 yic = yc(i)-yi(i)
288 zic = zc(i)-zi(i)
289 sx = - yab*zca + zab*yca
290 sy = - zab*xca + xab*zca
291 sz = - xab*yca + yab*xca
292 s2 = sx*sx+sy*sy+sz*sz
293 sax = yib*zic - zib*yic
294 say = zib*xic - xib*zic
295 saz = xib*yic - yib*xic
296 la(i) = (sx*sax+sy*say+sz*saz)/s2
297 sbx = yic*zia - zic*yia
298 sby = zic*xia - xic*zia
299 sbz = xic*yia - yic*xia
300 lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
301 lc(i) = one - la(i) - lb(i)
302 lj = min(la(i),lb(i),lc(i))
303 IF (lj<zerom) pene(i)=min(pene(i),penmin(i))
304 IF(la(i)<zero)THEN
305 IF(lb(i)<zero)THEN
306 la(i) = zero
307 lb(i) = zero
308 lc(i) = one
309 ELSEIF(lc(i)<zero)THEN
310 lc(i) = zero
311 la(i) = zero
312 lb(i) = one
313 ELSE
314 la(i) = zero
315 aaa = lb(i) + lc(i)
316 lb(i) = lb(i)/aaa
317 lc(i) = lc(i)/aaa
318 ENDIF
319 ELSEIF(lb(i)<zero)THEN
320 IF(lc(i)<zero)THEN
321 lb(i) = zero
322 lc(i) = zero
323 la(i) = one
324 ELSE
325 lb(i) = zero
326 aaa = lc(i) + la(i)
327 lc(i) = lc(i)/aaa
328 la(i) = la(i)/aaa
329 ENDIF
330 ELSEIF(lc(i)<zero)THEN
331 lc(i) = zero
332 aaa = la(i) + lb(i)
333 la(i) = la(i)/aaa
334 lb(i) = lb(i)/aaa
335 ENDIF
336 ENDDO
337 DO i=1,nel
338 IF(pene(i) == zero) cycle
339 IF(ifc2(i)>=3) THEN ! degenerated quad
340 SELECT CASE (ifc2(i))
341 CASE (3)
342 hj(i,1) = lb(i)
343 hj(i,2) = lc(i)
344 hj(i,3) = la(i)
345 hj(i,4) = zero
346 CASE (4)
347 hj(i,2) = lb(i)
348 hj(i,3) = lc(i)
349 hj(i,4) = la(i)
350 hj(i,1) = zero
351 CASE (5)
352 hj(i,4) = lb(i)
353 hj(i,1) = lc(i)
354 hj(i,2) = la(i)
355 hj(i,3) = zero
356 CASE (6)
357 hj(i,3) = lb(i)
358 hj(i,4) = lc(i)
359 hj(i,1) = la(i)
360 hj(i,2) = zero
361 END SELECT
362 ELSE
363 h0 = fourth * la(i)
364 SELECT CASE (itg(i))
365 CASE (1)
366 hj(i,1) = h0 + lb(i)
367 hj(i,2) = h0 + lc(i)
368 hj(i,3) = h0
369 hj(i,4) = h0
370 CASE (2)
371 hj(i,2) = h0 + lb(i)
372 hj(i,3) = h0 + lc(i)
373 hj(i,4) = h0
374 hj(i,1) = h0
375 CASE (3)
376 hj(i,3) = h0 + lb(i)
377 hj(i,4) = h0 + lc(i)
378 hj(i,1) = h0
379 hj(i,2) = h0
380 CASE (4)
381 hj(i,4) = h0 + lb(i)
382 hj(i,1) = h0 + lc(i)
383 hj(i,2) = h0
384 hj(i,3) = h0
385 END SELECT
386 END IF
387 ENDDO
388 f_q = ep02
389 DO i=1,nel
390 IF(pene(i) == zero) cycle
391 pendr = (pene(i)/penref(i))**2
392 fac = min(fqmax,f_q*pendr)
393 fn(i) = (fac+one)*stif0(i)*pene(i)
394 fkt = one+three*fac
395 stif(i) =max(stif(i),fkt*stif0(i))
396 ENDDO
397 DO i=1,nel
398 IF(pene(i) == zero) cycle
399 dx = vc(i,1) - hj(i,1)*vx1(i) - hj(i,2)*vx2(i)
400 . - hj(i,3)*vx3(i) - hj(i,4)*vx4(i)
401 dy = vc(i,2) - hj(i,1)*vy1(i) - hj(i,2)*vy2(i)
402 . - hj(i,3)*vy3(i) - hj(i,4)*vy4(i)
403 dz = vc(i,3) - hj(i,1)*vz1(i) - hj(i,2)*vz2(i)
404 . - hj(i,3)*vz3(i) - hj(i,4)*vz4(i)
405 dn = (nx(i)*dx + ny(i)*dy + nz(i)*dz)*dt1
406 e_distor(i) = e_distor(i) - fn(i)*dn
407 ENDDO
408 DO i=1,nel
409 IF (pene(i) ==zero) cycle
410 fx = nx(i)*fn(i)
411 fy = ny(i)*fn(i)
412 fz = nz(i)*fn(i)
413 forc_n(i,1) = forc_n(i,1) - fx
414 forc_n(i,2) = forc_n(i,2) - fy
415 forc_n(i,3) = forc_n(i,3) - fz
416 for_t1(i,1) = for_t1(i,1) + fx*hj(i,1)
417 for_t1(i,2) = for_t1(i,2) + fy*hj(i,1)
418 for_t1(i,3) = for_t1(i,3) + fz*hj(i,1)
419 for_t2(i,1) = for_t2(i,1) + fx*hj(i,2)
420 for_t2(i,2) = for_t2(i,2) + fy*hj(i,2)
421 for_t2(i,3) = for_t2(i,3) + fz*hj(i,2)
422 for_t3(i,1) = for_t3(i,1) + fx*hj(i,3)
423 for_t3(i,2) = for_t3(i,2) + fy*hj(i,3)
424 for_t3(i,3) = for_t3(i,3) + fz*hj(i,3)
425 for_t4(i,1) = for_t4(i,1) + fx*hj(i,4)
426 for_t4(i,2) = for_t4(i,2) + fy*hj(i,4)
427 for_t4(i,3) = for_t4(i,3) + fz*hj(i,4)
428 ENDDO
429 END IF ! (IFCTL1==1) THEN
430c-----------
431 RETURN
432 END
433
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
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
sfor_n2s4
subroutine sfor_n2s4(xi, yi, zi, stif, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, vx1, vx2, vx3, vx4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, for_t1, for_t2, for_t3, for_t4, forc_n, ll, ifctl, ifc1, penmin, penref, marge, fqmax, stif0, nel, vc, e_distor, dt1)
Definition sfor_n2s4.F:43
ssort_n4
subroutine ssort_n4(xi, yi, zi, marge, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, ifc1, stif, nel)
Definition ssort_n4.F:33