OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25pen3e.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!|| i25pen3e ../starter/source/interfaces/inter3d1/i25pen3e.F
25!||--- called by ------------------------------------------------------
26!|| inint3 ../starter/source/interfaces/inter3d1/inint3.F
27!||====================================================================
28 SUBROUTINE i25pen3e(
29 1 JLT ,IEDGE ,CAND_S,CAND_M,
30 2 N1 ,N2 ,M1 ,M2 ,
31 3 XXS1 ,XXS2 ,XYS1 ,XYS2 ,
32 4 XZS1 ,XZS2 ,XXM1 ,XXM2 ,XYM1 ,
33 5 XYM2 ,XZM1 ,XZM2 ,GAPVE ,PENE ,
34 6 EX ,EY ,EZ ,FX ,FY ,
35 7 FZ ,LEDGE ,IRECT ,X ,ITAB )
36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C G l o b a l P a r a m e t e r s
42C-----------------------------------------------
43#include "mvsiz_p.inc"
44#include "param_c.inc"
45C-----------------------------------------------
46C D u m m y A r g u m e n t s
47C-----------------------------------------------
48 INTEGER JLT, IEDGE
49 INTEGER CAND_S(*), CAND_M(*), IRECT(4,*), LEDGE(NLEDGE,*), ITAB(*),
50 . N1(MVSIZ),N2(MVSIZ),M1(MVSIZ),M2(MVSIZ)
51 my_real
52 . XXS1(*), XXS2(*), XYS1(*), XYS2(*), XZS1(*) , XZS2(*),
53 . XXM1(*), XXM2(*) , XYM1(*), XYM2(*), XZM1(*), XZM2(*),
54 . GAPVE(*), PENE(*), X(3,*),
55 . ex(mvsiz), ey(mvsiz), ez(mvsiz), fx(mvsiz), fy(mvsiz), fz(mvsiz)
56C-----------------------------------------------
57C L o c a l V a r i a b l e s
58C-----------------------------------------------
59 INTEGER I, IA, JA, IB, JB, SOL_EDGE, SH_EDGE, K, NJNDX, N4A, N4B,
60 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
61 . JX1(MVSIZ), JX2(MVSIZ), JX3(MVSIZ), JX4(MVSIZ),
62 . JNDX(MVSIZ), I4A(MVSIZ), I4B(MVSIZ)
63 my_real
64 . H1S(MVSIZ),H2S(MVSIZ),H1M(MVSIZ),H2M(MVSIZ),NX(MVSIZ),NY(MVSIZ),NZ(MVSIZ),
65 . XS12,YS12,ZS12,XM12,YM12,ZM12,XA,XB,
66 . xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
67 . xx,yy,zz,als,alm,det,aaa,bbb,nn(mvsiz),p1,p2
68 my_real
69 . xa0(mvsiz),xa1(mvsiz),xa2(mvsiz),xa3(mvsiz),xa4(mvsiz),
70 . ya0(mvsiz),ya1(mvsiz),ya2(mvsiz),ya3(mvsiz),ya4(mvsiz),
71 . za0(mvsiz),za1(mvsiz),za2(mvsiz),za3(mvsiz),za4(mvsiz),
72 . xb0(mvsiz),xb1(mvsiz),xb2(mvsiz),xb3(mvsiz),xb4(mvsiz),
73 . yb0(mvsiz),yb1(mvsiz),yb2(mvsiz),yb3(mvsiz),yb4(mvsiz),
74 . zb0(mvsiz),zb1(mvsiz),zb2(mvsiz),zb3(mvsiz),zb4(mvsiz)
75 my_real
76 . x0a(mvsiz,4),y0a(mvsiz,4),z0a(mvsiz,4),
77 . x0b(mvsiz,4),y0b(mvsiz,4),z0b(mvsiz,4),
78 . xna(mvsiz,4), yna(mvsiz,4), zna(mvsiz,4), xnb(mvsiz,4), ynb(mvsiz,4), znb(mvsiz,4),
79 . xs, ys, zs, xm, ym, zm, da, db, cnvx, da1, db1, da2, db2,
80 . rzero, run, rdix, rem30, rep30,
81 . alp,xxs,xys,xzs,
82 . xi0,yi0,zi0,xi1,yi1,zi1,xi2,yi2,zi2,
83 . sx1,sy1,sz1,sx2,sy2,sz2
84 INTEGER NTRIA(3,4)
85 DATA NTRIA/1,2,4,2,4,1,0,0,0,4,1,2/
86C-----------------------------------------------
87C F = [A*X1+(1-A)*X2-B*X3-(1-B)*X4]^2 + [..Y..]^2 + [..Z..]^2
88C DF/DA = 0 = (X1-X2)(A(X1-X2)+X2-X4 +B(X4-X3))+...
89C DF/DA = 0 = A(X1-X2)^2 +X2-X4 + B(X1-X2)(X4-X3))+...
90C DF/DA = 0 = A[(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
91C + B[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
92C + (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
93C DF/DB = 0 = (X4-X3)(A(X1-X2)+X2-X4 +B(X4-X3))+...
94C DF/DB = 0 = B[(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
95C + A[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
96C + (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
97C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
98C XM2 = [(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
99C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
100C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
101C XB = (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
102C A XS2 + B XSM + XA = 0
103C A XSM + B XM2 + XB = 0
104C
105C A = -(XA + B XSM)/XS2
106C -(XA + B XSM)*XSM + B XM2*XS2 + XB*XS2 = 0
107C -B XSM*XSM + B XM2*XS2 + XB*XS2-XA*XSM = 0
108C B*(XM2*XS2 - XSM*XSM) = -XB*XS2+XA*XSM
109C B = (XA*XSM-XB*XS2) / (XM2*XS2 - XSM*XSM)
110C A = (XB*XSM-XA*XM2) / (XM2*XS2 - XSM*XSM)
111C
112C IF B<0 => B=0
113C
114C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
115C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
116C A = - XA /XS2
117C B = 0
118C
119C ELSEIF B>1 => B=1
120C
121C B = 1
122C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
123C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
124C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
125C A = -(XA + XSM)/XS2
126C
127C IF A<0 => A=0
128C
129C
130C ELSEIF A>1 => A=1
131C
132C
133 pene(1:jlt)=ep20
134C
135 DO i=1,jlt
136
137 xm12 = xxm2(i)-xxm1(i)
138 ym12 = xym2(i)-xym1(i)
139 zm12 = xzm2(i)-xzm1(i)
140 xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
141
142 xs12 = xxs2(i)-xxs1(i)
143 ys12 = xys2(i)-xys1(i)
144 zs12 = xzs2(i)-xzs1(i)
145 xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
146 xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
147 xs2m2 = xxm2(i)-xxs2(i)
148 ys2m2 = xym2(i)-xys2(i)
149 zs2m2 = xzm2(i)-xzs2(i)
150
151 xa = xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
152 xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2
153 det = xm2*xs2 - xsm*xsm
154 det = max(em20,det)
155C
156 h1m(i) = (xa*xsm-xb*xs2) / det
157 xs2 = max(xs2,em20)
158 xm2 = max(xm2,em20)
159 h1m(i)=min(one,max(zero,h1m(i)))
160 h1s(i) = -(xa + h1m(i)*xsm) / xs2
161 h1s(i)=min(one,max(zero,h1s(i)))
162 h1m(i) = -(xb + h1s(i)*xsm) / xm2
163 h1m(i)=min(one,max(zero,h1m(i)))
164
165 h2s(i) = one -h1s(i)
166 h2m(i) = one -h1m(i)
167C !!!!!!!!!!!!!!!!!!!!!!!
168C PENE = GAP^2 - DIST^2 USED TO TEST IF NON-NULL
169C!!!!!!!!!!!!!!!!!!!!!!!!
170 nx(i) = h1s(i)*xxs1(i) + h2s(i)*xxs2(i)
171 . - h1m(i)*xxm1(i) - h2m(i)*xxm2(i)
172 ny(i) = h1s(i)*xys1(i) + h2s(i)*xys2(i)
173 . - h1m(i)*xym1(i) - h2m(i)*xym2(i)
174 nz(i) = h1s(i)*xzs1(i) + h2s(i)*xzs2(i)
175 . - h1m(i)*xzm1(i) - h2m(i)*xzm2(i)
176
177 nn(i) = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2)
178
179 pene(i) = max(zero,gapve(i) - nn(i))
180
181 ENDDO
182C
183 sol_edge =iedge/10 ! solids
184 sh_edge =iedge-10*sol_edge ! shells
185C
186 IF(sh_edge/=0)THEN
187 DO i=1,jlt
188
189C
190C Free edges, looking vs positive normal only
191C
192C / S
193C / x
194C M /
195C <------x Sector with Zero force
196C n(M) \
197C \
198C \
199
200 p1=zero
201 IF(ledge(3,cand_m(i))==0)
202 . p1=-(nx(i)*ex(i)+ny(i)*ey(i)+nz(i)*ez(i)) ! (n(M),SM) > 45 degrees
203
204 p2=zero
205 IF(ledge(3,cand_s(i))==0)
206 . p2= nx(i)*fx(i)+ny(i)*fy(i)+nz(i)*fz(i) ! (n(S),MS) > 45 degrees
207
208C NN(I)=SQRT(NN(I)) ! already done above
209 IF(p1 > em04*nn(i) .OR. p2 > em04*nn(i))pene(i)=zero ! Tolerance EM04
210
211 ENDDO
212 END IF
213C
214C---------------------------------------
215 RETURN
216 END
i25pen3e
subroutine i25pen3e(jlt, iedge, cand_s, cand_m, n1, n2, m1, m2, xxs1, xxs2, xys1, xys2, xzs1, xzs2, xxm1, xxm2, xym1, xym2, xzm1, xzm2, gapve, pene, ex, ey, ez, fx, fy, fz, ledge, irect, x, itab)
Definition i25pen3e.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21