OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3e.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!|| i25dst3e ../engine/source/interfaces/int25/i25dst3e.F
25!||--- called by ------------------------------------------------------
26!|| i25mainf ../engine/source/interfaces/int25/i25mainf.F
27!||====================================================================
28 SUBROUTINE i25dst3e(
29 1 JLT ,CAND_S,CAND_M,H1S ,H2S ,
30 2 H1M ,H2M ,NX ,NY ,NZ ,
31 3 STIF ,N1 ,N2 ,M1 ,M2 ,
32 4 JLT_NEW,XXS1 ,XXS2 ,XYS1 ,XYS2 ,
33 5 XZS1 ,XZS2 ,XXM1 ,XXM2 ,XYM1 ,
34 6 XYM2 ,XZM1 ,XZM2 ,VXS1 ,VXS2 ,
35 7 VYS1 ,VYS2 ,VZS1 ,VZS2 ,VXM1 ,
36 8 VXM2 ,VYM1 ,VYM2 ,VZM1 ,VZM2 ,
37 9 MS1 ,MS2 ,MM1 ,MM2 ,IEDGE ,
38 B NSMS ,INDEX ,INTFRIC ,IPARTFRICSI,IPARTFRICMI ,
39 C GAPVE ,EX ,EY ,EZ ,FX ,
40 D FY ,FZ ,LEDGE ,IRECT ,CAND_P ,
41 E IAM ,JAM ,IBM ,JBM ,IAS ,
42 F JAS ,IBS ,JBS ,ITAB ,EDGE_ID,
43 G DGAPLOADPMAX)
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53#include "param_c.inc"
54#include "sms_c.inc"
55C debug:
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER JLT,JLT_NEW, IEDGE,INTFRIC
60 INTEGER CAND_S(MVSIZ),CAND_M(MVSIZ),NSMS(MVSIZ),INDEX(*),
61 . N1(MVSIZ),N2(MVSIZ),M1(MVSIZ),M2(MVSIZ),
62 . IAM(MVSIZ),JAM(MVSIZ),IBM(MVSIZ),JBM(MVSIZ),
63 . IAS(MVSIZ),JAS(MVSIZ),IBS(MVSIZ),JBS(MVSIZ),
64 . IPARTFRICSI(MVSIZ),IPARTFRICMI(MVSIZ),
65 . IRECT(4,*),LEDGE(NLEDGE,*),ITAB(*)
66 INTEGER :: EDGE_ID(2,MVSIZ)
67 my_real
68 . H1S(*),H2S(*),H1M(*),H2M(*),NX(*),NY(*),NZ(*),STIF(*),
69 . XXS1(*), XXS2(*), XYS1(*), XYS2(*), XZS1(*) , XZS2(*),
70 . XXM1(*), XXM2(*) , XYM1(*), XYM2(*), XZM1(*), XZM2(*),
71 . VXS1(*), VXS2(*), VYS1(*), VYS2(*), VZS1(*), VZS2(*) ,
72 . VXM1(*), VXM2(*), VYM1(*), VYM2(*), VZM1(*), VZM2(*),
73 . MS1(*) ,MS2(*) ,MM1(*) ,MM2(*),
74 . gapve(*), ex(*), ey(*), ez(*), fx(*), fy(*), fz(*),
75 . cand_p(*)
76 my_real , INTENT(IN) :: dgaploadpmax
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I, SOL_EDGE, SH_EDGE, K
81 my_real
82 . pene2(mvsiz),
83 . xs12,ys12,zs12,xm12,ym12,zm12,xa,xb,
84 . xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
85 . det,gap2,p1,p2,nn(mvsiz)
86 INTEGER NTRIA(3,4)
87 DATA NTRIA/1,2,4,2,4,1,0,0,0,4,1,2/
88C-----------------------------------------------
89 JLT_NEW = 0
90C--------------------------------------------------------
91C
92C--------------------------------------------------------
93C F = [A*X1+(1-A)*X2-B*X3-(1-B)*X4]^2 + [..Y..]^2 + [..Z..]^2
94C DF/DA = 0 = (X1-X2)(A(X1-X2)+X2-X4 +B(X4-X3))+...
95C DF/DA = 0 = A(X1-X2)^2 +X2-X4 + B(X1-X2)(X4-X3))+...
96C DF/DA = 0 = A[(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
97C + B[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
98C + (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
99C DF/DB = 0 = (X4-X3)(A(X1-X2)+X2-X4 +B(X4-X3))+...
100C DF/DB = 0 = B[(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
101C + A[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
102C + (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
103C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
104C XM2 = [(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
105C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
106C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
107C XB = (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
108C A XS2 + B XSM + XA = 0
109C A XSM + B XM2 + XB = 0
110C
111C A = -(XA + B XSM)/XS2
112C -(XA + B XSM)*XSM + B XM2*XS2 + XB*XS2 = 0
113C -B XSM*XSM + B XM2*XS2 + XB*XS2-XA*XSM = 0
114C B*(XM2*XS2 - XSM*XSM) = -XB*XS2+XA*XSM
115C B = (XA*XSM-XB*XS2) / (XM2*XS2 - XSM*XSM)
116C A = (XB*XSM-XA*XM2) / (XM2*XS2 - XSM*XSM)
117C
118C IF B<0 => B=0
119C
120C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
121C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
122C A = - XA /XS2
123C B = 0
124C
125C ELSEIF B>1 => B=1
126C
127C B = 1
128C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
129C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
130C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
131C A = -(XA + XSM)/XS2
132C
133C IF A<0 => A=0
134C
135C
136C ELSEIF A>1 => A=1
137C
138C
139 pene2(1:jlt)=ep20
140C
141 DO i=1,jlt
142
143 xm12 = xxm2(i)-xxm1(i)
144 ym12 = xym2(i)-xym1(i)
145 zm12 = xzm2(i)-xzm1(i)
146 xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
147
148 xs12 = xxs2(i)-xxs1(i)
149 ys12 = xys2(i)-xys1(i)
150 zs12 = xzs2(i)-xzs1(i)
151 xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
152 xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
153 xs2m2 = xxm2(i)-xxs2(i)
154 ys2m2 = xym2(i)-xys2(i)
155 zs2m2 = xzm2(i)-xzs2(i)
156
157 xa = xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
158 xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2
159 det = xm2*xs2 - xsm*xsm
160
161 det = max(em20,det)
162C
163 h1m(i) = (xa*xsm-xb*xs2) / det
164 xs2 = max(xs2,em20)
165 xm2 = max(xm2,em20)
166 h1m(i)=min(one,max(zero,h1m(i)))
167 h1s(i) = -(xa + h1m(i)*xsm) / xs2
168 h1s(i)=min(one,max(zero,h1s(i)))
169 h1m(i) = -(xb + h1s(i)*xsm) / xm2
170 h1m(i)=min(one,max(zero,h1m(i)))
171
172 h2s(i) = one -h1s(i)
173 h2m(i) = one -h1m(i)
174C !!!!!!!!!!!!!!!!!!!!!!!
175C PENE = GAP^2 - DIST^2 UTILISE POUR TESTER SI NON NUL
176C!!!!!!!!!!!!!!!!!!!!!!!!
177 nx(i) = h1s(i)*xxs1(i) + h2s(i)*xxs2(i)
178 . - h1m(i)*xxm1(i) - h2m(i)*xxm2(i)
179 ny(i) = h1s(i)*xys1(i) + h2s(i)*xys2(i)
180 . - h1m(i)*xym1(i) - h2m(i)*xym2(i)
181 nz(i) = h1s(i)*xzs1(i) + h2s(i)*xzs2(i)
182 . - h1m(i)*xzm1(i) - h2m(i)*xzm2(i)
183 nn(i)= nx(i)*nx(i) + ny(i)*ny(i) + nz(i)*nz(i)
184
185 gap2 = (gapve(i)+dgaploadpmax)*(gapve(i)+dgaploadpmax)
186 pene2(i) = gap2 - nn(i)
187 pene2(i) = max(zero,pene2(i))
188
189 ENDDO
190C
191 sol_edge =iedge/10 ! solids
192 sh_edge =iedge-10*sol_edge ! shells
193C
194 IF(sh_edge/=0)THEN
195 DO i=1,jlt
196C
197C Free edges, looking vs positive normal only
198C
199C / S
200C / x
201C M /
202C <------x Sector with Zero force
203C n(M) |
204C |
205C |
206
207 p1=zero
208 IF(ibm(i)==0)
209 . p1=-(nx(i)*ex(i)+ny(i)*ey(i)+nz(i)*ez(i)) ! (n(M),SM) > 45 degrees
210
211 p2=zero
212 IF(ibs(i)==0)
213 . p2= nx(i)*fx(i)+ny(i)*fy(i)+nz(i)*fz(i) ! (n(S),MS) > 45 degrees
214
215 nn(i)=sqrt(nn(i))
216 IF(p1 > em04*nn(i) .OR. p2 > em04*nn(i))pene2(i)=zero ! Tolerance EM04
217
218 ENDDO
219 END IF
220C
221C---------------------------------------
222C
223 DO i=1,jlt
224 IF(pene2(i)==zero) cand_p(index(i))=zero ! Reset CAND_P
225 ENDDO
226C
227 DO i=1,jlt
228 IF(pene2(i)/=zero.AND.stif(i)/=zero)THEN
229 jlt_new = jlt_new + 1
230
231 edge_id(1,jlt_new) = edge_id(1,i)
232 edge_id(2,jlt_new) = edge_id(2,i)
233
234 cand_s(jlt_new) = cand_s(i)
235 cand_m(jlt_new) = cand_m(i)
236 edge_id(1,jlt_new) = edge_id(1,i)
237 edge_id(2,jlt_new) = edge_id(2,i)
238 n1(jlt_new) = n1(i)
239 n2(jlt_new) = n2(i)
240 m1(jlt_new) = m1(i)
241 m2(jlt_new) = m2(i)
242 ibm(jlt_new) = ibm(i)
243 h1s(jlt_new) = h1s(i)
244 h2s(jlt_new) = h2s(i)
245 h1m(jlt_new) = h1m(i)
246 h2m(jlt_new) = h2m(i)
247 nx(jlt_new) = nx(i)
248 ny(jlt_new) = ny(i)
249 nz(jlt_new) = nz(i)
250 stif(jlt_new)= stif(i)
251 gapve(jlt_new)= gapve(i)
252
253 xxs1(jlt_new) = xxs1(i)
254 xys1(jlt_new) = xys1(i)
255 xzs1(jlt_new) = xzs1(i)
256 xxs2(jlt_new) = xxs2(i)
257 xys2(jlt_new) = xys2(i)
258 xzs2(jlt_new) = xzs2(i)
259
260 vxs1(jlt_new) = vxs1(i)
261 vys1(jlt_new) = vys1(i)
262 vzs1(jlt_new) = vzs1(i)
263 vxs2(jlt_new) = vxs2(i)
264 vys2(jlt_new) = vys2(i)
265 vzs2(jlt_new) = vzs2(i)
266 vxm1(jlt_new) = vxm1(i)
267 vym1(jlt_new) = vym1(i)
268 vzm1(jlt_new) = vzm1(i)
269 vxm2(jlt_new) = vxm2(i)
270 vym2(jlt_new) = vym2(i)
271 vzm2(jlt_new) = vzm2(i)
272 ms1(jlt_new) = ms1(i)
273 ms2(jlt_new) = ms2(i)
274 mm1(jlt_new) = mm1(i)
275 mm2(jlt_new) = mm2(i)
276 IF(idtmins == 2) nsms(jlt_new)= nsms(i)
277 IF(intfric /= 0) THEN
278 ipartfricsi(jlt_new)=ipartfricsi(i)
279 ipartfricmi(jlt_new)=ipartfricmi(i)
280 ENDIF
281 index(jlt_new) = index(i)
282 ENDIF
283 ENDDO
284 RETURN
285 END
i25dst3e
subroutine i25dst3e(jlt, cand_s, cand_m, h1s, h2s, h1m, h2m, nx, ny, nz, stif, n1, n2, m1, m2, jlt_new, xxs1, xxs2, xys1, xys2, xzs1, xzs2, xxm1, xxm2, xym1, xym2, xzm1, xzm2, vxs1, vxs2, vys1, vys2, vzs1, vzs2, vxm1, vxm2, vym1, vym2, vzm1, vzm2, ms1, ms2, mm1, mm2, iedge, nsms, index, intfric, ipartfricsi, ipartfricmi, gapve, ex, ey, ez, fx, fy, fz, ledge, irect, cand_p, iam, jam, ibm, jbm, ias, jas, ibs, jbs, itab, edge_id, dgaploadpmax)
Definition i25dst3e.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21