OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25pen3_edg.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!|| i25pen3_edg ../engine/source/interfaces/intsort/i25pen3_edg.F
25!||--- called by ------------------------------------------------------
26!|| i25sto_edg ../engine/source/interfaces/intsort/i25sto_edg.F
27!||--- uses -----------------------------------------------------
28!|| tri7box ../engine/share/modules/tri7box.F
29!||====================================================================
30 SUBROUTINE i25pen3_edg(JLT ,CAND_S ,CAND_M ,DRAD ,IGAP0 ,
31 . NEDGE ,LEDGE ,MARGE ,GAPE ,GAP_E_L ,
32 . IGAP ,X ,IRECT ,PENE ,ADMSR ,
33 . EDG_BISECTOR,VTX_BISECTOR,ITAB,
34 . XREM_EDGE,S1_XREM,S2_XREM,DGAPLOAD)
35
36C-----------------------------------------------
37C M o d u l e s
38C-----------------------------------------------
39C USE TRI25EBOX
40 USE tri7box
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "mvsiz_p.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "com04_c.inc"
53#include "param_c.inc"
54#include "i25edge_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER, INTENT(IN) :: S1_XREM,S2_XREM
59 my_real, INTENT(IN) :: XREM_EDGE(S1_XREM,S2_XREM)
60 INTEGER JLT, IGAP0, NEDGE, IGAP
61 INTEGER IRECT(4,*), CAND_S(*), CAND_M(*), LEDGE(NLEDGE,NEDGE), ADMSR(4,*), ITAB(NUMNOD)
62 my_real , INTENT(IN) :: dgapload ,drad
63 my_real
64 . marge,
65 . x(3,numnod), gape(*), gap_e_l(*), pene(mvsiz)
66 real*4 edg_bisector(3,4,*), vtx_bisector(3,2,*)
67C-----------------------------------------------
68C L o c a l V a r i a b l e s
69C-----------------------------------------------
70 INTEGER I, IG,N1,N2,M1,M2,NI,L,J, IE, JE, IL, JL
71 my_real
72 . XXS1(MVSIZ) ,XXS2(MVSIZ) ,XYS1(MVSIZ) ,XYS2(MVSIZ) ,XZS1(MVSIZ) ,XZS2(MVSIZ) ,
73 . XXM1(MVSIZ) ,XXM2(MVSIZ) ,XYM1(MVSIZ) ,XYM2(MVSIZ) ,XZM1(MVSIZ) ,XZM2(MVSIZ) ,
74 . xs12,ys12,zs12,xm12,ym12,zm12,xa,xb,
75 . xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
76 . xx,yy,zz,als,alm,det,gap2,
77 . xmax1,ymax1,zmax1,xmax2,ymax2,zmax2,
78 . xmin1,ymin1,zmin1,xmin2,ymin2,zmin2,gapv(mvsiz),pene_l(mvsiz),vdt(mvsiz),
79 . aaa, dx, dy, dz, dd, ex, ey ,ez, nni, ni2, invcos
80 my_real
81 . maxgapv
82C-----------------------------------------------
83CDIR$ NOFUSION
84 DO i=1,jlt
85 ie = cand_m(i)
86
87 IF(cand_s(i) <= nedge) THEN
88 je = cand_s(i)
89 IF(igap0 == 0) THEN
90 gapv(i)=two*gape(ie)+gape(je)
91 ELSE
92 gapv(i)=two*(gape(ie)+gape(je))
93 END IF
94 IF(igap == 3)
95 . gapv(i)=min(gap_e_l(ie)+gap_e_l(je),gapv(i)) ! under-estimated ...
96 gapv(i)=max(drad,gapv(i)+dgapload)+marge
97 ELSE
98 IF(igap0 == 0) THEN
99 gapv(i)=two*gape(ie)+xrem_edge(e_gap,cand_s(i)-nedge )
100 ELSE
101 gapv(i)=two*(gape(ie)+xrem_edge(e_gap,cand_s(i)-nedge ))
102 END IF
103 IF(igap == 3)
104 . gapv(i)=min(gap_e_l(ie)+xrem_edge(e_gapl,cand_s(i)-nedge),gapv(i))
105 gapv(i)=max(drad,gapv(i)+dgapload)+marge
106 ENDIF
107 ENDDO
108
109C--------------------------------------------------------
110C F = [A*X1+(1-A)*X2-B*X3-(1-B)*X4]^2 + [..Y..]^2 + [..Z..]^2
111C DF/DA = 0 = (X1-X2)(A(X1-X2)+X2-X4 +B(X4-X3))+...
112C DF/DA = 0 = A(X1-X2)^2 +X2-X4 + B(X1-X2)(X4-X3))+...
113C DF/DA = 0 = A[(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
114C + B[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
115C + (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
116C DF/DB = 0 = (X4-X3)(A(X1-X2)+X2-X4 +B(X4-X3))+...
117C DF/DB = 0 = B[(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
118C + A[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
119C + (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
120C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
121C XM2 = [(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
122C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
123C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
124C XB = (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
125C A XS2 + B XSM + XA = 0
126C A XSM + B XM2 + XB = 0
127C
128C A = -(XA + B XSM)/XS2
129C -(XA + B XSM)*XSM + B XM2*XS2 + XB*XS2 = 0
130C -B XSM*XSM + B XM2*XS2 + XB*XS2-XA*XSM = 0
131C B*(XM2*XS2 - XSM*XSM) = -XB*XS2+XA*XSM
132C B = (XA*XSM-XB*XS2) / (XM2*XS2 - XSM*XSM)
133C A = (XB*XSM-XA*XM2) / (XM2*XS2 - XSM*XSM)
134C--------------------------------------------------------
135CDIR$ NOFUSION
136 DO i=1,jlt
137 IF(cand_s(i)<=nedge) THEN
138 n1 = ledge(5,cand_s(i))
139 n2 = ledge(6,cand_s(i))
140
141 xxs1(i) = x(1,n1)
142 xys1(i) = x(2,n1)
143 xzs1(i) = x(3,n1)
144 xxs2(i) = x(1,n2)
145 xys2(i) = x(2,n2)
146 xzs2(i) = x(3,n2)
147 ELSE
148 ni = cand_s(i) - nedge
149 xxs1(i) = xrem_edge(e_x1,ni)
150 xys1(i) = xrem_edge(e_y1,ni)
151 xzs1(i) = xrem_edge(e_z1,ni)
152 xxs2(i) = xrem_edge(e_x2,ni)
153 xys2(i) = xrem_edge(e_y2,ni)
154 xzs2(i) = xrem_edge(e_z2,ni)
155 END IF
156
157 m1 = ledge(5,cand_m(i))
158 m2 = ledge(6,cand_m(i))
159
160 xxm1(i) = x(1,m1)
161 xym1(i) = x(2,m1)
162 xzm1(i) = x(3,m1)
163 xxm2(i) = x(1,m2)
164 xym2(i) = x(2,m2)
165 xzm2(i) = x(3,m2)
166 END DO
167C
168C--------------------------------------------------------
169C calcul d'un minorant de la distance
170C--------------------------------------------------------
171CDIR$ NOFUSION
172 DO i=1,jlt
173 xmax1 = max(xxs1(i),xxs2(i))
174 ymax1 = max(xys1(i),xys2(i))
175 zmax1 = max(xzs1(i),xzs2(i))
176 xmax2 = max(xxm1(i),xxm2(i))
177 ymax2 = max(xym1(i),xym2(i))
178 zmax2 = max(xzm1(i),xzm2(i))
179 xmin1 = min(xxs1(i),xxs2(i))
180 ymin1 = min(xys1(i),xys2(i))
181 zmin1 = min(xzs1(i),xzs2(i))
182 xmin2 = min(xxm1(i),xxm2(i))
183 ymin2 = min(xym1(i),xym2(i))
184 zmin2 = min(xzm1(i),xzm2(i))
185 dd = max(xmin1-xmax2,ymin1-ymax2,zmin1-zmax2,
186 . xmin2-xmax1,ymin2-ymax1,zmin2-zmax1)
187 IF(dd > gapv(i))THEN
188 pene(i) = zero
189 cycle
190 ENDIF
191
192c calcul de la distance^2
193
194 xm12 = xxm2(i)-xxm1(i)
195 ym12 = xym2(i)-xym1(i)
196 zm12 = xzm2(i)-xzm1(i)
197 xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
198
199 xs12 = xxs2(i)-xxs1(i)
200 ys12 = xys2(i)-xys1(i)
201 zs12 = xzs2(i)-xzs1(i)
202 xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
203 xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
204 xs2m2 = xxm2(i)-xxs2(i)
205 ys2m2 = xym2(i)-xys2(i)
206 zs2m2 = xzm2(i)-xzs2(i)
207
208 xa = xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
209 xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2
210 det = xm2*xs2 - xsm*xsm
211 det = max(em20,det)
212C
213 alm = (xa*xsm-xb*xs2) / det
214 xs2 = max(xs2,em20)
215 xm2 = max(xm2,em20)
216 alm=min(one,max(zero,alm))
217 als = -(xa + alm*xsm) / xs2
218 als=min(one,max(zero,als))
219 alm = -(xb + als*xsm) / xm2
220 alm=min(one,max(zero,alm))
221
222C PENE = GAP^2 - DIST^2 UTILISE POUR TESTER SI NON NUL
223
224 xx = als*xxs1(i) + (one-als)*xxs2(i)
225 . - alm*xxm1(i) - (one-alm)*xxm2(i)
226 yy = als*xys1(i) + (one-als)*xys2(i)
227 . - alm*xym1(i) - (one-alm)*xym2(i)
228 zz = als*xzs1(i) + (one-als)*xzs2(i)
229 . - alm*xzm1(i) - (one-alm)*xzm2(i)
230
231 gap2=gapv(i)*gapv(i)
232 pene(i) = max(zero,gap2- xx*xx - yy*yy - zz*zz)
233C
234 END DO
235C
236 RETURN
237 END
i25pen3_edg
subroutine i25pen3_edg(jlt, cand_s, cand_m, drad, igap0, nedge, ledge, marge, gape, gap_e_l, igap, x, irect, pene, admsr, edg_bisector, vtx_bisector, itab, xrem_edge, s1_xrem, s2_xrem, dgapload)
Definition i25pen3_edg.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21