OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25neigh.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "param_c.inc"
#include "scr03_c.inc"
#include "scr17_c.inc"
#include "com01_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i25neigh (nrtm, nsn, nsv, irect, irtlm, mvoisin, evoisin, mseglo, msegtyp, itab, x, id, titr, igeo, nadmsr, admsr, adskyn, iadnor, nrtm_sh, iedge, nedge, ledge, lbound, edg_cos, nisub, lisub, addsubm, lisubm, inflg_subm, nisube, addsube, lisube, inflg_sube, noint, nmn, msr, nom_opt, ilev, mbinflg, ebinflg, ielem_m, idel_solid)
subroutine i25ini_gap_n (nrtm, irect, ixs, geo, ixc, ixtg, ixt, ixp, ipart, ipartc, iparttg, thk, thk_part, gap_n, gap_m, nmn, msr, gapn_m, gapmax_m, gapscale, igeo, msegtyp, gapmsav, ithk25)
subroutine i25neigh_seg_en (n1, ied1, n2, ied2, mvoisin, ne, ice, iself, i1, i2, irect, x, ie, nrtm, msegtyp, irr)
subroutine i25neigh_seg_e (n1, ied1, n2, ied2, n, ic, iself, i1, i2, irect, nrtm, msegtyp, mvoisin)
subroutine i25neigh_removeallbut1 (n, ic, iself, irect, x, i1, i2, irr)
subroutine i25neigh_seg_opp (ei, ej, irect, x, iop)
subroutine i25neigh_msg_err (i1, i2, itab, irr)

Function/Subroutine Documentation

◆ i25ini_gap_n()

subroutine i25ini_gap_n ( integer nrtm,
integer, dimension(4,*) irect,
integer, dimension(nixs,*) ixs,
geo,
integer, dimension(nixc,*) ixc,
integer, dimension(nixtg,*) ixtg,
integer, dimension(nixt,*) ixt,
integer, dimension(nixp,*) ixp,
integer, dimension(lipart1,*) ipart,
integer, dimension(*) ipartc,
integer, dimension(*) iparttg,
thk,
thk_part,
gap_n,
gap_m,
integer nmn,
integer, dimension(*) msr,
gapn_m,
gapmax_m,
gapscale,
integer, dimension(npropgi,*) igeo,
integer, dimension(*) msegtyp,
gapmsav,
integer ithk25 )

Definition at line 1132 of file i25neigh.F.

1138C-----------------------------------------------
1139C I m p l i c i t T y p e s
1140C-----------------------------------------------
1141#include "implicit_f.inc"
1142C-----------------------------------------------
1143C C o m m o n B l o c k s
1144C-----------------------------------------------
1145#include "com01_c.inc"
1146#include "com04_c.inc"
1147#include "param_c.inc"
1148#include "scr17_c.inc"
1149C-----------------------------------------------
1150C D u m m y A r g u m e n t s
1151C-----------------------------------------------
1152 INTEGER NRTM,IRECT(4,*), IXS(NIXS,*), IXC(NIXC,*),
1153 . IXTG(NIXTG,*), IXT(NIXT,*), IXP(NIXP,*),
1154 . IPART(LIPART1,*), IPARTC(*), IPARTTG(*),
1155 . MSR(*),NMN,IGEO(NPROPGI,*),MSEGTYP(*), ITHK25
1156C REAL
1157 my_real
1158 . geo(npropg,*), thk(*),thk_part(*),gap_n(4,*),gap_m(*),
1159 . gapn_m(*),gapmax_m, gapscale,gapmsav(*)
1160C-----------------------------------------------
1161C L o c a l V a r i a b l e s
1162C-----------------------------------------------
1163 INTEGER I,J,K,IW,I1,I2,I3,MG,M,IP,IGTYP
1164 my_real
1165 . dx
1166 my_real, DIMENSION(:), ALLOCATABLE :: wa
1167C init
1168C
1169 ALLOCATE(wa(numnod))
1170 DO i=1,numnod
1171 wa(i)=zero
1172 END DO
1173C
1174C------------------------------------
1175C GAP NOEUDS IN WA
1176C------------------------------------
1177 DO i=1,numelc
1178 mg=ixc(6,i)
1179 igtyp = igeo(11,mg)
1180 ip = ipartc(i)
1181 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
1182 dx=half*thk_part(ip)
1183 ELSEIF ( thk(i) /= zero .AND. iintthick == 0) THEN
1184 dx=half*thk(i)
1185 ELSEIF(igtyp == 17) THEN
1186 dx=half*thk(i)
1187 ELSE
1188 dx=half*geo(1,mg)
1189 ENDIF
1190 wa(ixc(2,i))=max(wa(ixc(2,i)),dx)
1191 wa(ixc(3,i))=max(wa(ixc(3,i)),dx)
1192 wa(ixc(4,i))=max(wa(ixc(4,i)),dx)
1193 wa(ixc(5,i))=max(wa(ixc(5,i)),dx)
1194 ENDDO
1195C
1196 DO i=1,numeltg
1197 mg=ixtg(5,i)
1198 igtyp = igeo(11,mg)
1199 ip = iparttg(i)
1200 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
1201 dx=half*thk_part(ip)
1202 ELSEIF ( thk(numelc+i) /= zero .AND. iintthick == 0) THEN
1203 dx=half*thk(numelc+i)
1204 ELSEIF(igtyp == 17) THEN
1205 dx=half*thk(numelc+i)
1206 ELSE
1207 dx=half*geo(1,mg)
1208 ENDIF
1209 wa(ixtg(2,i))=max(wa(ixtg(2,i)),dx)
1210 wa(ixtg(3,i))=max(wa(ixtg(3,i)),dx)
1211 wa(ixtg(4,i))=max(wa(ixtg(4,i)),dx)
1212 ENDDO
1213C
1214 DO i=1,numnod
1215 wa(i)=gapscale*wa(i)
1216 END DO
1217C--------exclude lines in main surfaces
1218C DO I=1,NUMELT
1219C MG=IXT(4,I)
1220C DX=HALF*SQRT(GEO(1,MG))
1221C WA(IXT(2,I))=MAX(WA(IXT(2,I)),DX)
1222C WA(IXT(3,I))=MAX(WA(IXT(3,I)),DX)
1223C ENDDO
1224C
1225C DO I=1,NUMELP
1226C MG=IXP(5,I)
1227C DX=HALF*SQRT(GEO(1,MG))
1228C WA(IXP(2,I))=MAX(WA(IXP(2,I)),DX)
1229C WA(IXP(3,I))=MAX(WA(IXP(3,I)),DX)
1230C ENDDO
1231C------------------------------------
1232C INI GAP_N (4), GAP_M is modified taking into account nodal gap for sorting
1233C------------------------------------
1234C -----due to the fact that if surf_M does not contain the defining w/ shell
1235C-------> should not take into account GAP_shell
1236 DO i=1,nrtm
1237 IF (msegtyp(i)==0) THEN
1238 DO j=1,4
1239 m=irect(j,i)
1240 wa(m) = zero
1241 END DO
1242 END IF !(MSEGTYP(I)==0) THEN
1243 END DO ! nrtm
1244C
1245 DO i=1,nmn
1246 m = msr(i)
1247 wa(m) = min(wa(m),gapmax_m)
1248 gapn_m(i) = wa(m)
1249 END DO
1250C
1251 DO i=1,nrtm
1252 gap_m(i) = zero
1253 DO j=1,4
1254 m=irect(j,i)
1255 gap_n(j,i)=wa(m)
1256 gap_m(i) = max(gap_m(i),wa(m))
1257 END DO
1258 END DO ! nrtm
1259C
1260 IF(ithk25==1) THEN
1261 DO i=1,nrtm
1262 gapmsav(i) = gap_m(i)
1263 ENDDO
1264 ENDIF
1265 DEALLOCATE(wa)
1266C
1267 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ i25neigh()

subroutine i25neigh ( integer nrtm,
integer nsn,
integer, dimension(*) nsv,
integer, dimension(4,nrtm) irect,
integer, dimension(4,nsn) irtlm,
integer, dimension(4,nrtm) mvoisin,
integer, dimension(4,nrtm) evoisin,
integer, dimension(nrtm) mseglo,
integer, dimension(nrtm) msegtyp,
integer, dimension(*) itab,
x,
integer id,
character(len=nchartitle) titr,
integer, dimension(npropgi,*) igeo,
integer nadmsr,
integer, dimension(4,nrtm) admsr,
integer, dimension(4*nrtm+1) adskyn,
integer, dimension(4,nrtm) iadnor,
integer nrtm_sh,
integer iedge,
integer nedge,
integer, dimension(nledge,*) ledge,
integer, dimension(*) lbound,
edg_cos,
integer nisub,
integer, dimension(*) lisub,
integer, dimension(*) addsubm,
integer, dimension(*) lisubm,
integer, dimension(*) inflg_subm,
integer nisube,
integer, dimension(*) addsube,
integer, dimension(*) lisube,
integer, dimension(*) inflg_sube,
integer noint,
integer nmn,
integer, dimension(*) msr,
integer, dimension(lnopt1,*) nom_opt,
integer ilev,
integer, dimension(*) mbinflg,
integer, dimension(*) ebinflg,
integer, dimension(2,nrtm), intent(in) ielem_m,
integer, intent(in) idel_solid )

Definition at line 36 of file i25neigh.F.

45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE my_alloc_mod
49 USE message_mod
50 USE names_and_titles_mod , ONLY : nchartitle
51C-----------------------------------------------
52C I m p l i c i t T y p e s
53C-----------------------------------------------
54#include "implicit_f.inc"
55C-----------------------------------------------
56C C o m m o n B l o c k s
57C-----------------------------------------------
58#include "com04_c.inc"
59#include "units_c.inc"
60#include "param_c.inc"
61#include "scr03_c.inc"
62#include "scr17_c.inc"
63C-----------------------------------------------
64C D u m m y A r g u m e n t s
65C-----------------------------------------------
66 INTEGER NRTM,NSN,NMN,NSV(*),IRECT(4,NRTM),MVOISIN(4,NRTM),EVOISIN(4,NRTM),
67 . MSEGLO(NRTM),IRTLM(4,NSN),MSEGTYP(NRTM),ITAB(*),
68 . IGEO(NPROPGI,*), NADMSR, ADMSR(4,NRTM), ADSKYN(4*NRTM+1),
69 . IADNOR(4,NRTM), NRTM_SH, IEDGE, NEDGE, LEDGE(NLEDGE,*),
70 . LBOUND(*),
71 . NISUB, LISUB(*), ADDSUBM(*), LISUBM(*), INFLG_SUBM(*),
72 . NISUBE, ADDSUBE(*), LISUBE(*), INFLG_SUBE(*), MSR(*),
73 . ILEV, MBINFLG(*), EBINFLG(*)
74 INTEGER NOINT, NOM_OPT(LNOPT1,*)
75 INTEGER ID
76 my_real
77 . x(3,*), edg_cos
78 CHARACTER(LEN=NCHARTITLE) :: TITR
79 INTEGER , INTENT(IN) :: IDEL_SOLID ! solid erosion
80 INTEGER , INTENT(IN) :: IELEM_M(2,NRTM) ! elements connected to the segment (especially solid)
81C-----------------------------------------------
82C L o c a l V a r i a b l e s
83C-----------------------------------------------
84 INTEGER I,J,K,L,IW,I1,I2,I3,I4,M,N,NMAX,E_MAX,E_ID,N_EI,
85 1 NE0,NRTM0,ISH,NVERTX, IVERTX, JVERTX,
86 2 J1,J2,J3,J4,K1,K2,L1,L2,KPERM1(4),KPERM2(4),IRR,
87 3 NFT,JLT,MI,II(4),JJ(4),IEDG,JEDG,IOK,
88 4 IBASE, KBASE, KOLD, FIN, NE, EJ, NEDGE_TMP,
89 5 SOL_EDGE, SH_EDGE, ISTORE,
90 6 LL, KK, CUR, NEXT, JSUB, KSUB, IMS1, IMS2, IMS3, IMS4,
91 7 NISUBN, INFLG, MAXADD, STAT
92 INTEGER IX1, IX2, IX3, IX4,
93 , JX1, JX2, JX3, JX4,
94 . NA, NB, EA, EB,JM,MJ,DN_EI
95 INTEGER, DIMENSION(:),ALLOCATABLE :: MVOI
96 INTEGER, DIMENSION(:,:),ALLOCATABLE :: EIDNOD
97 my_real
98 . xa1,xa2,xa3,xa4,
99 . ya1,ya2,ya3,ya4,
100 . za1,za2,za3,za4,
101 . xb1,xb2,xb3,xb4,
102 . yb1,yb2,yb3,yb4,
103 . zb1,zb2,zb3,zb4,
104 . x01,y01,z01,x02,y02,z02,
105 . xna, yna, zna, xnb, ynb, znb, aaa, bbb, ang
106 INTEGER WORK(70000)
107 INTEGER, DIMENSION(:,:),ALLOCATABLE :: CLEF,LEDGE_TMP1,LEDGE_TMP2
108 INTEGER, DIMENSION(:),ALLOCATABLE :: INDEX, ITAG, NTAG, ROOT
109 INTEGER, DIMENSION(:), ALLOCATABLE :: MLOC,KAD,ADDSUBN,ADDSUBN_TMP,
110 . LISUBN,INFLG_SUBN,IXSUB
111 INTEGER, DIMENSION(:,:,:),ALLOCATABLE :: MNEIGH_SOLID
112C-----------------------------------------------
113C E x t e r n a l F u n c t i o n s
114C-----------------------------------------------
115 INTEGER BITSET, BITGET
116 EXTERNAL bitset, bitget
117C-----------------------------------------------
118 DATA kperm1/1,3,5,7/
119 DATA kperm2/2,4,6,8/
120C-----------------------------------------------
121 CALL my_alloc(itag,numnod)
122 CALL my_alloc(ntag,4*nrtm)
123 CALL my_alloc(root,4*nrtm)
124C
125 DO i=1,nrtm
126 mseglo(i)=i
127 ENDDO
128C-------------------------------------------------------
129C Neighbhoors researching
130C shell segs have been duplicated w/ inverse order
131C for the moment all antisymmetry surface will be stored at the end
132C for Solids : searching of neighbhoors is only done for external segments
133C-----------------------------------------------
134C
135 DO i=1,numnod
136 itag(i)=0
137 ENDDO
138 DO i=1,nrtm
139 IF(ielem_m(2,i) == 0) THEN
140
141 DO j=1,3
142 m=irect(j,i)
143 itag(m)=itag(m)+1
144 END DO
145 IF (irect(4,i)/=irect(3,i))THEN
146 m= irect(4,i)
147 itag(m)=itag(m)+1
148 END IF
149 ENDIF
150 END DO
151C-----MSEGTYP (> NRTM for i=NRTM0+1,NRTM0+NRTM_SH) -> IM2SH = MSEGTYP-NRTM---------
152C-----------max number of connected segment per node
153 nmax=0
154 DO i=1,numnod
155 nmax=max(nmax,itag(i))
156 itag(i)=0
157 ENDDO
158 ALLOCATE(mvoi(nmax+10),eidnod(nmax,numnod))
159 mvoi(1:nmax+10)=0
160 eidnod(1:nmax,1:numnod)=0
161C------------ini- E_ids of each node
162 DO i=1,nrtm
163 IF(ielem_m(2,i) == 0) THEN
164
165 DO j=1,3
166 m=irect(j,i)
167 itag(m)=itag(m)+1
168 eidnod(itag(m),m)=i
169 END DO
170 IF (irect(4,i)/=irect(3,i)) THEN
171 m= irect(4,i)
172 itag(m)=itag(m)+1
173 eidnod(itag(m),m)=i
174 END IF
175 ENDIF
176 END DO
177C------------MVOISIN-(seg number) ---
178 e_max=4
179 DO i=1,nrtm
180 DO j=1,e_max
181 mvoisin(j,i)=0
182 END DO
183 END DO
184C
185 DO i=1,nrtm
186 IF(ielem_m(2,i) == 0) THEN
187 n_ei=0
188C----seg 1-2------
189 i1 =irect(1,i)
190 i2 =irect(2,i)
191 CALL i25neigh_seg_en(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),mvoisin,
192 1 n_ei,mvoi ,i ,i1 ,i2 ,irect,
193 2 x ,mvoisin(1,i),nrtm,msegtyp,irr )
194 IF (irr >0) CALL i25neigh_msg_err(i1,i2,itab,irr)
195C----seg 2-3------
196 i1 =irect(2,i)
197 i2 =irect(3,i)
198 CALL i25neigh_seg_en(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),mvoisin,
199 1 n_ei,mvoi ,i ,i1 ,i2 ,irect,
200 2 x ,mvoisin(1,i),nrtm,msegtyp,irr )
201 IF (irr >0) CALL i25neigh_msg_err(i1,i2,itab,irr)
202C----seg 3-4------
203 i1 =irect(3,i)
204 i2 =irect(4,i)
205 CALL i25neigh_seg_en(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),mvoisin,
206 1 n_ei,mvoi ,i ,i1 ,i2 ,irect,
207 2 x ,mvoisin(1,i),nrtm,msegtyp,irr )
208 IF (irr >0) CALL i25neigh_msg_err(i1,i2,itab,irr)
209C----seg 1-4------
210 i1 =irect(4,i)
211 i2 =irect(1,i)
212 CALL i25neigh_seg_en(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),mvoisin,
213 1 n_ei,mvoi ,i ,i1 ,i2 ,irect,
214 2 x ,mvoisin(1,i),nrtm,msegtyp,irr )
215 IF (irr >0) CALL i25neigh_msg_err(i1,i2,itab,irr)
216C
217 mvoi(1:n_ei)=0 ! Reset MVOI
218C
219 ENDIF
220 END DO !I=1,NRTM
221C
222C--------------------------------------------
223 DO i=1,nrtm
224 DO j=1,4
225 l = mvoisin(j,i)
226 IF(l/=0) THEN
227 DO k=1,4
228 IF(mvoisin(k,l)==i) GOTO 120
229 END DO
230 WRITE(istdo,'(A,/,10I10)')
231 . 'i25inisu_nei - internal error : a segment is not neighboring its neighbor segments...',
232 . i,(itab(irect(m,i)),m=1,4),
233 . l,(itab(irect(m,l)),m=1,4)
234 120 CONTINUE
235 END IF
236 END DO
237 END DO !I=1,NRTM
238C--------------------------------------------
239 DEALLOCATE(mvoi,eidnod)
240C--------------------------------------------
241C------------------------------------------------------
242C Solid erosion :
243C if internal segments of solids are defined : need to build ADMSR for all NRTM
244C Store all segments connected to the same edge for solid segments
245C-----------------------------------------------------
246 IF(idel_solid > 0) THEN
247 DO i=1,numnod
248 itag(i)=0
249 ENDDO
250 DO i=1,nrtm
251 IF(ielem_m(1,i) <= numels) THEN
252
253 DO j=1,3
254 m=irect(j,i)
255 itag(m)=itag(m)+1
256 END DO
257 IF (irect(4,i)/=irect(3,i))THEN
258 m= irect(4,i)
259 itag(m)=itag(m)+1
260 END IF
261 ENDIF
262 END DO
263C
264 nmax=0
265 DO i=1,numnod
266 nmax=max(nmax,itag(i))
267 itag(i)=0
268 ENDDO
269 ALLOCATE(eidnod(nmax,numnod))
270 eidnod(1:nmax,1:numnod)=0
271 ALLOCATE(mvoi(nmax+10))
272 mvoi(1:nmax+10)=0
273C
274 ALLOCATE(mneigh_solid(nmax,4,nrtm))
275 mneigh_solid(1:nmax,1:4,1:nrtm) = 0
276 DO i=1,nrtm
277 IF(ielem_m(1,i) <= numels) THEN
278
279 DO j=1,3
280 m=irect(j,i)
281 itag(m)=itag(m)+1
282 eidnod(itag(m),m)=i
283 END DO
284 IF (irect(4,i)/=irect(3,i)) THEN
285 m= irect(4,i)
286 itag(m)=itag(m)+1
287 eidnod(itag(m),m)=i
288 END IF
289 ENDIF
290 END DO
291C
292C-------------neighbour segments-----------
293C
294 DO i=1,nrtm
295 IF(ielem_m(1,i) <= numels) THEN
296 n_ei=0
297C----seg 1-2------
298 i1 =irect(1,i)
299 i2 =irect(2,i)
300 CALL i25neigh_seg_e(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),n_ei,mvoi,i,
301 . i1,i2,irect,nrtm,msegtyp ,mvoisin)
302 mneigh_solid(1:n_ei,1,i) = mvoi(1:n_ei)
303 ne0 =n_ei
304C----seg 2-3------
305 i1 =irect(2,i)
306 i2 =irect(3,i)
307 CALL i25neigh_seg_e(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),n_ei,mvoi,i,
308 . i1,i2,irect,nrtm,msegtyp ,mvoisin)
309 dn_ei = n_ei - ne0
310 mneigh_solid(1:dn_ei,2,i) = mvoi(ne0+1:n_ei)
311 ne0 =n_ei
312C----seg 3-4------
313 i1 =irect(3,i)
314 i2 =irect(4,i)
315 CALL i25neigh_seg_e(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),n_ei,mvoi,i,
316 . i1,i2,irect,nrtm,msegtyp ,mvoisin)
317 dn_ei = n_ei - ne0
318 mneigh_solid(1:dn_ei,3,i) = mvoi(ne0+1:n_ei)
319 ne0 =n_ei
320C----seg 1-4------
321 i1 =irect(4,i)
322 i2 =irect(1,i)
323 CALL i25neigh_seg_e(itag(i1),eidnod(1,i1),itag(i2),eidnod(1,i2),n_ei,mvoi,i,
324 . i1,i2,irect,nrtm,msegtyp ,mvoisin)
325 dn_ei = n_ei - ne0
326 mneigh_solid(1:dn_ei,4,i) = mvoi(ne0+1:n_ei)
327C
328 mvoi(1:n_ei)=0 ! Reset MVOI
329 ENDIF
330 ENDDO
331
332 DEALLOCATE(mvoi,eidnod)
333
334 ENDIF
335 DEALLOCATE(itag)
336
337C--------------------------------------------
338
339
340 nvertx=0
341 DO i=1,nrtm
342 DO j=1,3
343 nvertx =nvertx+1
344 admsr(j,i) =nvertx
345 root(nvertx)=nvertx
346 END DO
347 IF(irect(4,i)/=irect(3,i))THEN
348 nvertx =nvertx+1
349 admsr(4,i) =nvertx
350 root(nvertx)=nvertx
351 ELSE
352 admsr(4,i)=admsr(3,i)
353 END IF
354 END DO
355C
356 DO i=1,nrtm
357C
358 ii(1:4)=irect(1:4,i)
359C
360 DO iedg=1,4
361C
362 IF(ii(4)==ii(3).AND.iedg==3) cycle
363 i1=iedg
364 i2=mod(iedg,4)+1
365C
366 mi=mvoisin(iedg,i)
367 IF(mi/=0)THEN
368C
369 jj(1:4)=irect(1:4,mi)
370C
371 iok=0
372 DO jedg=1,4
373C
374 IF(jj(4)==jj(3).AND.jedg==3) cycle
375C
376 j1=jedg
377 j2=mod(jedg,4)+1
378 IF(ii(i2)==jj(j1).AND.ii(i1)==jj(j2))THEN
379
380 ivertx=min(root(admsr(i1,i)),root(admsr(j2,mi)))
381 jvertx=max(root(admsr(i1,i)),root(admsr(j2,mi)))
382 IF(jvertx/=ivertx) root(jvertx)=ivertx
383
384 ivertx=min(root(admsr(i2,i)),root(admsr(j1,mi)))
385 jvertx=max(root(admsr(i2,i)),root(admsr(j1,mi)))
386 IF(jvertx/=ivertx) root(jvertx)=ivertx
387
388 iok=1
389
390 evoisin(iedg,i)=jedg
391
392 ELSEIF(ii(i1)==jj(j1).AND.ii(i2)==jj(j2))THEN
393 print *,'i25inisu_nei - internal error : non-consistent neighboring segment'
394
395 END IF
396 END DO
397
398 IF(iok==0)
399 . WRITE(istdo,*) 'i25inisu_nei - internal error : no common edge w/neighboring segment',
400 . itab(ii(1)),itab(ii(2)),itab(ii(3)),itab(ii(4)),
401 . itab(jj(1)),itab(jj(2)),itab(jj(3)),itab(jj(4))
402 END IF ! IF(MI/=0)THEN
403C
404 IF(ielem_m(1,i) <= numels.AND.idel_solid > 0) THEN
405 DO k=1,nmax
406 mj = mneigh_solid(k,iedg,i)
407 IF(mj/=0.AND.mj/=mi) THEN
408C
409 jj(1:4)=irect(1:4,mj)
410C
411 iok=0
412 DO jedg=1,4
413C
414 IF(jj(4)==jj(3).AND.jedg==3) cycle
415C
416 j1=jedg
417 j2=mod(jedg,4)+1
418 IF(ii(i2)==jj(j1).AND.ii(i1)==jj(j2))THEN
419
420 ivertx=min(root(admsr(i1,i)),root(admsr(j2,mj)))
421 jvertx=max(root(admsr(i1,i)),root(admsr(j2,mj)))
422 IF(jvertx/=ivertx) root(jvertx)=ivertx
423
424 ivertx=min(root(admsr(i2,i)),root(admsr(j1,mj)))
425 jvertx=max(root(admsr(i2,i)),root(admsr(j1,mj)))
426 IF(jvertx/=ivertx) root(jvertx)=ivertx
427
428 iok=1
429
430c ELSEIF(II(I1)==JJ(J1).AND.II(I2)==JJ(J2))THEN
431c print *,'i25inisu_nei - internal error : non-consistent neighboring segment'
432
433 END IF
434 END DO
435
436c IF(IOK==0)
437c . WRITE(ISTDO,*) 'i25inisu_nei - internal error : no common edge w/neighboring segment',
438c . itab(ii(1)),itab(ii(2)),itab(ii(3)),itab(ii(4)),
439c . itab(jj(1)),itab(jj(2)),itab(jj(3)),itab(jj(4))
440c END IF ! IF(MJ/=0.AND.MJ/=MI)THEN
441 ENDIF
442 ENDDO
443 ENDIF
444 END DO
445C
446 END DO
447C
448 IF(idel_solid > 0) DEALLOCATE(mneigh_solid)
449C
450 DO i=1,nvertx
451 j=i
452 DO WHILE(root(j) < j)
453 j=root(j)
454 END DO
455 root(i)=j
456 END DO
457C
458 nadmsr=0
459 ntag(1:nvertx)=0
460 DO i=1,nvertx
461 j=root(i)
462 IF(ntag(j)==0)THEN
463 nadmsr =nadmsr+1
464 ntag(j)=nadmsr
465 END IF
466 END DO
467C
468 DO i=1,nrtm
469 admsr(1,i)=root(admsr(1,i))
470 admsr(1,i)=ntag(admsr(1,i))
471 admsr(2,i)=root(admsr(2,i))
472 admsr(2,i)=ntag(admsr(2,i))
473 admsr(3,i)=root(admsr(3,i))
474 admsr(3,i)=ntag(admsr(3,i))
475 admsr(4,i)=root(admsr(4,i))
476 admsr(4,i)=ntag(admsr(4,i))
477 END DO
478C-------------------------------------------------
479C Compute addresses for Parith/ON assembling of the normals
480C-------------------------------------------------
481 ntag(1:4*nrtm)=0
482 DO i=1,nrtm
483 i1 = abs(admsr(1,i))
484 i2 = abs(admsr(2,i))
485 i3 = abs(admsr(3,i))
486 i4 = abs(admsr(4,i))
487 ntag(i1) = ntag(i1) + 1
488 ntag(i2) = ntag(i2) + 1
489 ntag(i3) = ntag(i3) + 1
490 IF(i4/=i3) ntag(i4) = ntag(i4) + 1
491 END DO
492
493 adskyn(1) = 1
494 DO n=1,nadmsr
495 adskyn(n+1) = adskyn(n)+ntag(n)
496 END DO
497
498 DO n=1,nrtm
499 i1 = abs(admsr(1,n))
500 i2 = abs(admsr(2,n))
501 i3 = abs(admsr(3,n))
502 i4 = abs(admsr(4,n))
503 IF(irect(3,n)/=irect(4,n))THEN
504 iadnor(1,n)=adskyn(i1)
505 adskyn(i1) = adskyn(i1)+1
506 iadnor(2,n)=adskyn(i2)
507 adskyn(i2) = adskyn(i2)+1
508 iadnor(3,n)=adskyn(i3)
509 adskyn(i3) = adskyn(i3)+1
510 iadnor(4,n)=adskyn(i4)
511 adskyn(i4) = adskyn(i4)+1
512 ELSE
513 iadnor(1,n)=adskyn(i1)
514 adskyn(i1) = adskyn(i1)+1
515 iadnor(2,n)=adskyn(i2)
516 adskyn(i2) = adskyn(i2)+1
517 iadnor(3,n)=adskyn(i3)
518 adskyn(i3) = adskyn(i3)+1
519 END IF
520 END DO
521
522C
523C Reset ADSKYN
524 adskyn(1) = 1
525 DO n=1,nadmsr
526 adskyn(n+1) = adskyn(n)+ntag(n)
527 END DO
528C
529 DEALLOCATE(ntag,root)
530C
531 nrtm0=nrtm-nrtm_sh
532C
533 nedge =0
534 lbound(1:nadmsr)=0
535C-----------------------------------------------------------------------
536C LEDGE(1:NLEDGE)
537C 1: "left" segment number <=> index dans IRECT(1:NRTM)
538C 2: local number of the edge within "left" segment
539C 3: "right" segment number
540C 4: local number of the edge within "right" segment ! Squeezing T, X.. shapes wrt rupture
541C 5: I1 = numero local du noeud 1 dans NSV(1:NSN)
542C 6: I2 = numero local du noeud 1 dans NSV(1:NSN)
543C 7: Type of the edge :
544C < 0 <=> l'arete est uniquement main, pas second
545C +/-1 arete de solide
546C +/-2 arete de coque
547C for engine:
548C 8: Global ID (id in starter)
549C 9: Weight = 1 => secondary secnd handled by current domain
550C 10: orientation flag left seg
551C 11: first node of left seg
552C 12: second node left seg
553C 13: orientation flag right seg
554C 14: first node of right seg
555C 15: second node right seg
556
557C-----------------------------------------------------------------------
558 ALLOCATE(ledge_tmp1(nledge,4*nrtm))
559C
560 DO i=1,nrtm
561 IF(ielem_m(2,i) == 0) THEN
562 IF(i <= nrtm0)THEN
563 ibase=i
564 ELSE
565 ibase=-msegtyp(i)
566 ! IF(IBASE > 0) THEN
567 IF(ibase > nrtm)ibase=ibase-nrtm
568 ! END IF
569 END IF
570
571 DO j=1,4
572 i1=irect(j,i)
573 i2=irect(mod(j,4)+1,i)
574C
575 k=mvoisin(j,i)
576C
577 IF(k <= nrtm0)THEN
578 kbase=k ! y-compris K=0
579 ELSE
580 kbase=-msegtyp(k)
581 ! IF(KBASE > 0) THEN
582 IF(kbase > nrtm)kbase=kbase-nrtm
583 ! END IF
584 END IF
585C
586 IF(kbase < ibase)THEN
587
588 IF(.NOT.(i1==i2.AND.j==3))THEN
589C
590 nedge=nedge+1
591 ledge_tmp1(1,nedge)=i
592 ledge_tmp1(2,nedge)=j
593 ledge_tmp1(3,nedge)=k
594 ledge_tmp1(4,nedge)=0
595 IF(itab(i1) < itab(i2))THEN
596 ledge_tmp1(5,nedge)=i1
597 ledge_tmp1(6,nedge)=i2
598 ELSE
599 ledge_tmp1(5,nedge)=i2
600 ledge_tmp1(6,nedge)=i1
601 END IF
602C
603 IF(k/=0)THEN
604 IF(msegtyp(i)==0.AND.msegtyp(k)==0)THEN ! Solid only
605 ledge_tmp1(7,nedge)=1 ! arete solide-solide
606 ELSE
607 ledge_tmp1(7,nedge)=2 ! arete coque-coque ou coque-solide
608 END IF
609 ELSE ! Bord
610 ledge_tmp1(7,nedge)=2
611 END IF
612C
613 IF(k/=0)THEN
614 DO l=1,4
615 k1=irect(l,k)
616 k2=irect(mod(l,4)+1,k)
617 IF(.NOT.(k1==k2.AND.l==3).AND.((k1==i1.AND.k2==i2).OR.(k2==i1.AND.k1==i2))) THEN
618 ledge_tmp1(4,nedge)=l
619 END IF
620 END DO
621 IF(ledge_tmp1(4,nedge)==0)THEN
622 WRITE(istdo,'(A)')
623 . 'i25inisu_nei - internal error : could not find the edge on neighboring element'
624 END IF
625 END IF
626C
627 END IF
628
629
630 END IF
631C
632 IF(k==0)THEN
633 IF(.NOT.(i1==i2.AND.j==3))THEN
634C
635 lbound(admsr(j,i))=1
636 lbound(admsr(mod(j,4)+1,i))=1
637C
638 END IF
639 END IF
640 END DO
641 ENDIF
642 END DO
643C-----------------------------------------------------------------------
644C
645C Only 1 edge upon upper and lower skin is retained
646C Except for T, X... connections <=> keeping all edges between vertices (cf ADMSR)
647C
648C x-------x-------x
649C | Only 1 Edge
650C x-------x-------x
651C
652C
653C x
654C |
655C edge | edge
656C x-------x-------x
657C | T shape <=> 3 edges, X shape <=> 4 edges, ...
658C x-------x-------x
659C edge
660C
661C-----------------------------------------------------------------------
662C
663C For T, X ... shapes
664C If 1 element fails, free edges may appear while it is not really the case
665C x x
666C | |
667C edge | edge free edge | edge
668C x-------x-------x x-------x
669C | ====> |
670C x-------x-------x x-------x
671C edge free edge
672C It is considered as it is not a real problem
673C
674C-----------------------------------------------------------------------
675 ALLOCATE(clef(4,nedge),ledge_tmp2(nledge,nedge),index(2*nedge))
676 DO i=1,nedge
677
678 i1 = ledge_tmp1(5,i)
679 i2 = ledge_tmp1(6,i)
680C
681 ne =ledge_tmp1(1,i)
682 ej =ledge_tmp1(2,i)
683
684 IF(ne <= nrtm0)THEN
685 ibase=ne
686 ELSE
687 ibase=-msegtyp(ne)
688 IF(ibase > nrtm)ibase=ibase-nrtm
689 END IF
690C
691 ne =ledge_tmp1(3,i)
692 ej =ledge_tmp1(4,i)
693
694 IF(ne <= nrtm0)THEN
695 kbase=ne ! y-compris NE=0
696 ELSE
697 kbase=-msegtyp(ne)
698 IF(kbase > nrtm)kbase=kbase-nrtm
699 END IF
700
701 index(i) = i
702
703 clef(1,i) = i1
704 clef(2,i) = i2
705 clef(3,i) = ibase
706 clef(4,i) = kbase
707
708 END DO
709C
710 CALL my_orders(0,work,clef,index,nedge,4)
711C
712 nedge_tmp = nedge
713 ledge_tmp2(1:nledge,1:nedge)=ledge_tmp1(1:nledge,1:nedge)
714C
715 nedge = 0
716C
717 i=1
718 DO WHILE(i <= nedge_tmp)
719 nedge = nedge+1
720 kold = index(i)
721 ledge_tmp1(1:nledge,nedge)=ledge_tmp2(1:nledge,kold)
722
723 fin = 0
724 DO WHILE (fin == 0 .AND. i < nedge_tmp)
725 i = i+1
726 k = index(i)
727 IF(clef(1,k)/=clef(1,kold).OR.
728 . clef(2,k)/=clef(2,kold).OR.
729 . clef(3,k)/=clef(3,kold).OR.
730 . clef(4,k)/=clef(4,kold))THEN
731 fin=1
732 END IF
733 kold=k
734 END DO
735
736 IF(i==nedge_tmp .AND. fin==0) i=i+1 ! Exit
737
738 END DO
739
740C
741 nedge_tmp = nedge
742 nedge = 0
743
744 sol_edge =iedge/10 ! solids
745 sh_edge =iedge-10*sol_edge ! shells
746C
747 DO i=1,nedge_tmp
748C
749C IEDG= 1 : Edge to edge contact is activated using the external border edges from surf_ID1 and surf_ID2
750C & Edge to edge contact is not activated for edges supported by solids only
751C IEDG= 2 : Edge to edge contact is activated using all edges supported by shell elements from surf_ID1 and surf_ID2
752C & Edge to edge contact is not activated for edges supported by solids only
753C IEDG= 10 : Edge to edge contact is activated using sharp edges between contact solid segments
754C & Edge to edge contact is not activated for edges supported by shells only
755C IEDG= 20: Edge to edge contact is activated using all edges between contact solid segments
756C & Edge to edge contact is not activated for edges supported by shells only
757C And then the combinations 11, 12, 21 and 22 (let, 1st digit for solids, 2nd for shells)
758C
759 istore=0
760 IF(ledge_tmp1(7,i)==1)THEN ! edge between contact solid segments
761 IF(sol_edge==1.OR.sol_edge==3)THEN
762C
763 na =ledge_tmp1(1,i)
764 ea =ledge_tmp1(2,i)
765 nb =ledge_tmp1(3,i)
766 eb =ledge_tmp1(4,i)
767
768 IF(na==0 .OR. nb==0) THEN
769 print *,' internal error - i25neigh'
770 END IF
771
772 ix1 =irect(1,na)
773 ix2 =irect(2,na)
774 ix3 =irect(3,na)
775 ix4 =irect(4,na)
776C
777 xa1=x(1,ix1)
778 ya1=x(2,ix1)
779 za1=x(3,ix1)
780 xa2=x(1,ix2)
781 ya2=x(2,ix2)
782 za2=x(3,ix2)
783 xa3=x(1,ix3)
784 ya3=x(2,ix3)
785 za3=x(3,ix3)
786 xa4=x(1,ix4)
787 ya4=x(2,ix4)
788 za4=x(3,ix4)
789C
790 x01 = xa3 - xa1
791 y01 = ya3 - ya1
792 z01 = za3 - za1
793C
794 x02 = xa4 - xa2
795 y02 = ya4 - ya2
796 z02 = za4 - za2
797C
798 xna = y01*z02 - z01*y02
799 yna = z01*x02 - x01*z02
800 zna = x01*y02 - y01*x02
801C
802 aaa=one/max(em30,sqrt(xna*xna+yna*yna+zna*zna))
803 xna = -xna*aaa
804 yna = -yna*aaa
805 zna = -zna*aaa
806C
807 jx1 =irect(1,nb)
808 jx2 =irect(2,nb)
809 jx3 =irect(3,nb)
810 jx4 =irect(4,nb)
811C
812 xb1=x(1,jx1)
813 yb1=x(2,jx1)
814 zb1=x(3,jx1)
815 xb2=x(1,jx2)
816 yb2=x(2,jx2)
817 zb2=x(3,jx2)
818 xb3=x(1,jx3)
819 yb3=x(2,jx3)
820 zb3=x(3,jx3)
821 xb4=x(1,jx4)
822 yb4=x(2,jx4)
823 zb4=x(3,jx4)
824C
825 x01 = xb3 - xb1
826 y01 = yb3 - yb1
827 z01 = zb3 - zb1
828C
829 x02 = xb4 - xb2
830 y02 = yb4 - yb2
831 z02 = zb4 - zb2
832C
833 xnb = y01*z02 - z01*y02
834 ynb = z01*x02 - x01*z02
835 znb = x01*y02 - y01*x02
836C
837 bbb=one/max(em30,sqrt(xnb*xnb+ynb*ynb+znb*znb))
838 xnb = -xnb*bbb
839 ynb = -ynb*bbb
840 znb = -znb*bbb
841C
842 ang = xna*xnb+yna*ynb+zna*znb
843 IF (ang < edg_cos) THEN
844 istore=1 ! arete vive
845 ELSEIF(sol_edge==1)THEN
846 ledge_tmp1(7,i)=-1
847 istore=1 ! all edges on main side
848 END IF
849 ELSEIF(sol_edge==2)THEN
850 istore=1
851 END IF
852 ELSEIF(sh_edge/=0)THEN
853 istore=1 ! all shell edges are stored even if IEDG=1 (because of sorting & rupture)
854 END IF
855
856
857 IF(istore==1)THEN
858 nedge = nedge+1
859 ledge(1:nledge,nedge)=ledge_tmp1(1:nledge,i)
860 END IF
861
862 END DO
863C
864 DEALLOCATE(clef,index,ledge_tmp1,ledge_tmp2)
865C
866c DO I=1,NRTM
867c print *,'I,MSEGTYP(I)=',I,MSEGTYP(I),itab(irect(1:4,i))
868c print *,'MVOISIN(1,I)=',(MVOISIN(J,I),J=1,4)
869c END DO
870c DO I=1,NEDGE
871c print *,I,LEDGE(1,I),LEDGE(2,I),LEDGE(3,I),LEDGE(4,I),ITAB(LEDGE(5,I)),ITAB(LEDGE(6,I)),LEDGE(7,I)
872c END DO
873C
874C-----------------------------------------
875C Flag d'appartenance des aretes a S1/S2
876C-----------------------------------------
877 IF(nedge/=0 .AND. ilev==2)THEN
878 ebinflg(1:nedge)=0
879 DO i=1,nedge
880 ne =ledge(1,i)
881 inflg=mbinflg(ne)
882 ims1=bitget(inflg,0)
883 ims2=bitget(inflg,1)
884 IF(ims1/=0) ebinflg(i)=bitset(ebinflg(i),0)
885 IF(ims2/=0) ebinflg(i)=bitset(ebinflg(i),1)
886 ne=ledge(3,i)
887 IF(ne/=0)THEN
888 inflg=mbinflg(ne)
889 ims1=bitget(inflg,0)
890 ims2=bitget(inflg,1)
891 IF(ims1/=0) ebinflg(i)=bitset(ebinflg(i),0)
892 IF(ims2/=0) ebinflg(i)=bitset(ebinflg(i),1)
893 END IF
894 END DO
895 END IF
896C
897 IF(nedge/=0.AND.nisub/=0)THEN
898
899 ALLOCATE (mloc(numnod),kad(max(nmn,nedge)),stat=stat)
900 mloc(1:numnod)=0
901 DO i=1,nmn
902 n = msr(i)
903 mloc(n) = i
904 END DO
905
906 ALLOCATE(addsubn_tmp(nmn+1), addsubn(nmn+1))
907C-----------------------------------------
908C Calcul de ADDSUBN, LISUBN, INFLG_SUBN pour les noeuds mains
909C-----------------------------------------
910 addsubn_tmp(1:nmn+1)=0
911 DO i=1,nrtm
912 DO j=1,4
913 IF(.NOT.(j==3.AND.irect(3,i)==irect(4,i)))THEN
914 n = mloc(irect(j,i))
915 addsubn_tmp(n)=addsubn_tmp(n)+addsubm(i+1)-addsubm(i)
916 END IF
917 END DO
918 END DO
919C
920 cur=1
921 DO n=1,nmn
922 next = cur+addsubn_tmp(n)
923 addsubn_tmp(n) = cur
924 cur = next
925 END DO
926 addsubn_tmp(1+nmn)=cur
927C
928 nisubn=addsubn_tmp(1+nmn)-1
929 ALLOCATE (lisubn(nisubn),inflg_subn(nisubn),stat=stat)
930 inflg_subn(1:nisubn)=0
931C
932C utilise KAD(1:NMN) pour remplir LISUBN, INFLG_SUBN
933 DO n=1,nmn
934 kad(n)=addsubn_tmp(n)
935 END DO
936 DO i=1,nrtm
937 DO j=1,4
938 IF(.NOT.(j==3.AND.irect(3,i)==irect(4,i)))THEN
939 n = mloc(irect(j,i))
940 DO kk=addsubm(i),addsubm(i+1)-1
941 lisubn(kad(n)) =lisubm(kk)
942 inflg_subn(kad(n))=inflg_subm(kk)
943 kad(n)=kad(n)+1
944 END DO
945 END IF
946 END DO
947 END DO
948C-----------------------------------------
949C Compactage de LISUBN ET Combinaison des flags INFLG_SUBN
950C-----------------------------------------
951 maxadd = 0
952 DO n=1,nmn
953 maxadd=max(maxadd,addsubn_tmp(n+1)-addsubn_tmp(n))
954 END DO
955 ALLOCATE(ixsub(2*maxadd),stat=stat)
956C
957 DO n=1,nmn
958 kad(n)=addsubn_tmp(n)
959 END DO
960C
961 DO n=1,nmn
962 DO ll = 1,addsubn_tmp(n+1)-addsubn_tmp(n)
963 ixsub(ll) = ll
964 ENDDO
965C
966C 8.1356: Subscript #1 of the array LISUBN has value 12641 which is greater than the upper bound of 12640
967 IF(addsubn_tmp(n+1)-addsubn_tmp(n) > 1 .AND. addsubn_tmp(n) <=nisubn ) THEN
968 CALL my_orders(0,work,lisubn(addsubn_tmp(n)),ixsub,addsubn_tmp(n+1)-addsubn_tmp(n),1)
969 ENDIF
970C
971 cur =0
972 DO ll=addsubn_tmp(n),addsubn_tmp(n+1)-1
973 kk = addsubn_tmp(n)-1+ixsub(ll-addsubn_tmp(n)+1)
974C
975 IF(lisubn(kk)/=cur)THEN
976C
977C Combines INFLG
978 inflg=inflg_subn(kk)
979 ims1=bitget(inflg,0)
980 IF(ims1/=0) inflg_subn(kad(n))=
981 . bitset(inflg_subn(kad(n)),0)
982 ims2=bitget(inflg,1)
983 IF(ims2/=0) inflg_subn(kad(n))=
984 . bitset(inflg_subn(kad(n)),1)
985C
986 cur = lisubn(kk)
987 lisubn(kad(n))=cur
988 kad(n)=kad(n)+1
989 ELSE
990C
991C Combines INFLG
992 inflg=inflg_subn(kk)
993 ims1=bitget(inflg,0)
994 IF(ims1/=0) inflg_subn(kad(n)-1)=
995 . bitset(inflg_subn(kad(n)-1),0)
996 ims2=bitget(inflg,1)
997 IF(ims2/=0) inflg_subn(kad(n)-1)=
998 . bitset(inflg_subn(kad(n)-1),1)
999 END IF
1000 END DO
1001C
1002 addsubn(n)=kad(n)-addsubn_tmp(n)
1003 END DO
1004C
1005 cur=1
1006 DO n=1,nmn
1007 next = cur+addsubn(n)
1008 addsubn(n) = cur
1009 cur = next
1010 END DO
1011 addsubn(1+nmn)=cur
1012C
1013 DO n=1,nmn
1014 DO kk=addsubn(n),addsubn(n+1)-1
1015 lisubn(kk) =lisubn(addsubn_tmp(n)+kk-addsubn(n))
1016 inflg_subn(kk)=inflg_subn(addsubn_tmp(n)+kk-addsubn(n))
1017 END DO
1018 END DO
1019C-----------------------------------------
1020C Calcul de ADDSUBE, LISUBE, INFLG_SUBE a partir de ADDSUBN, LISUBN, INFLG_SUBN
1021C-----------------------------------------
1022 addsube(1:nedge+1) = 0
1023 inflg_sube(1:nisube)=0
1024C
1025 DO ne=1,nedge
1026 i1 =mloc(ledge(5,ne))
1027 i2 =mloc(ledge(6,ne))
1028
1029 ll =addsubn(i1)
1030 kk =addsubn(i2)
1031 DO WHILE(ll<addsubn(i1+1))
1032 jsub=lisubn(ll)
1033 DO WHILE(kk<addsubn(i2+1))
1034 ksub=lisubn(kk)
1035 IF(ksub==jsub)THEN
1036 addsube(ne)=addsube(ne)+1
1037 kk=kk+1
1038 ELSE IF(ksub<jsub)THEN
1039 kk=kk+1
1040 ELSE
1041 GO TO 100
1042 END IF
1043 END DO
1044 100 CONTINUE
1045 ll=ll+1
1046 END DO
1047 END DO
1048C
1049 cur=1
1050 DO ne=1,nedge
1051 next = cur+addsube(ne)
1052 addsube(ne)= cur
1053 cur = next
1054 END DO
1055 addsube(1+nedge)=cur
1056C
1057C utilise KAD(1:NEDGE) pour remplir LISUBE, INFLG_SUBE
1058 DO ne=1,nedge
1059 kad(ne)=addsube(ne)
1060 END DO
1061
1062 DO ne=1,nedge
1063 i1 =mloc(ledge(5,ne))
1064 i2 =mloc(ledge(6,ne))
1065
1066 ll =addsubn(i1)
1067 kk =addsubn(i2)
1068 DO WHILE(ll<addsubn(i1+1))
1069 jsub=lisubn(ll)
1070 ims1 = bitget(inflg_subn(ll),0)
1071 ims2 = bitget(inflg_subn(ll),1)
1072 DO WHILE(kk<addsubn(i2+1))
1073 ksub=lisubn(kk)
1074 ims3 = bitget(inflg_subn(kk),0)
1075 ims4 = bitget(inflg_subn(kk),1)
1076 IF(ksub==jsub)THEN
1077 lisube(kad(ne))=jsub
1078 IF(ims1==1.AND.ims3==1) ! edge belongs to S1 is I1 and I2 belong to S1
1079 . inflg_sube(kad(ne))=
1080 . bitset(inflg_sube(kad(ne)),0)
1081 IF(ims2==1.AND.ims4==1) ! edge belongs to S2 is I1 and I2 belong to S2
1082 . inflg_sube(kad(ne))=
1083 . bitset(inflg_sube(kad(ne)),1)
1084 kad(ne) =kad(ne)+1
1085 kk=kk+1
1086 ELSE IF(ksub<jsub)THEN
1087 kk=kk+1
1088 ELSE
1089 GO TO 200
1090 END IF
1091 END DO
1092 200 CONTINUE
1093 ll=ll+1
1094 END DO
1095 END DO
1096
1097 DEALLOCATE (mloc,kad,addsubn_tmp,addsubn,lisubn,inflg_subn,ixsub)
1098C-------------------------------------
1099 IF(ipri>=6) THEN
1100 WRITE(iout,1000)
1101 WRITE(iout,1010)noint
1102 WRITE(iout,'(10I10)')
1103 . (nom_opt(1,ninter+lisub(jsub)),jsub=1,nisub)
1104 WRITE(iout,1030)
1105 DO ne=1,nedge
1106 jsub=addsube(ne)
1107 n =addsube(ne+1)-addsube(ne)
1108 IF(n>0)THEN
1109 WRITE(iout,'(3I10)')ne,itab(ledge(5,ne)),itab(ledge(6,ne))
1110 WRITE(iout,'(30X,2I10)')
1111 . (lisube(jsub-1+k),inflg_sube(jsub-1+k),k=1,n)
1112 END IF
1113 END DO
1114 END IF
1115 END IF ! IF(NEDGE/=0.AND.NISUB/=0)THEN
1116C-------------------------------------
1117 1000 FORMAT( /1x,' STRUCTURE OF SUB-INTERFACES OUTPUT TO TH'/
1118 . 1x,' ----------------------------------------'// )
1119 1010 FORMAT( /1x,' INTERFACE ID . . . . . . . . . . . . . .',i10/,
1120 . ' -> LIST OF SUB-INTERFACES IDS : ')
1121 1030 FORMAT(/,' EDGE NODE 1 NODE 2'/
1122 . ' NUMBER ID ID'/
1123 . ' -> LIST OF SUB-INTERFACES (LOCAL NUMBERS IN INTERFACE)'/)
1124C-------------------------------------
1125 RETURN
bitget
integer function bitget(i, n)
Definition bitget.F:37
bitset
integer function bitset(i, n)
Definition bitget.F:66
i25neigh_msg_err
subroutine i25neigh_msg_err(i1, i2, itab, irr)
Definition i25neigh.F:1704
i25neigh_seg_e
subroutine i25neigh_seg_e(n1, ied1, n2, ied2, n, ic, iself, i1, i2, irect, nrtm, msegtyp, mvoisin)
Definition i25neigh.F:1384
i25neigh_seg_en
subroutine i25neigh_seg_en(n1, ied1, n2, ied2, mvoisin, ne, ice, iself, i1, i2, irect, x, ie, nrtm, msegtyp, irr)
Definition i25neigh.F:1281
my_orders
void my_orders(int *mode, int *iwork, int *data, int *index, int *n, int *irecl)
Definition my_orders.c:82
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003

◆ i25neigh_msg_err()

subroutine i25neigh_msg_err ( integer i1,
integer i2,
integer, dimension(*) itab,
integer irr )

Definition at line 1703 of file i25neigh.F.

1704 USE message_mod
1705C-----------------------------------------------
1706C I m p l i c i t T y p e s
1707C-----------------------------------------------
1708#include "implicit_f.inc"
1709C-----------------------------------------------
1710C C o m m o n B l o c k s
1711C-----------------------------------------------
1712#include "scr03_c.inc"
1713C-----------------------------------------------
1714C D u m m y A r g u m e n t s
1715C-----------------------------------------------
1716 INTEGER I1,I2,ITAB(*),IRR
1717C-----------------------------------------------
1718C L o c a l V a r i a b l e s
1719C-----------------------------------------------
1720 INTEGER I
1721C----Warning,ERROR out----------------
1722 IF(ipri==0) RETURN
1723#ifndef HYPERMESH_LIB
1724 IF (irr ==11) THEN
1725C-----multi-neibour but no valid one
1726 CALL ancmsg(msgid=1245,
1727 . msgtype=msgwarning,
1728 . anmode=aninfo_blind_2,
1729 . i2=itab(i1),i3=itab(i2))
1730c write(iout,*) '***Warning: No valid commun Seg with line:',
1731c + Itab(I1),Itab(I2)
1732 ELSEIF (irr ==12) THEN
1733C-----multi-neibour but no valid one
1734 CALL ancmsg(msgid=1246,
1735 . msgtype=msgerror,
1736 . anmode=aninfo_blind_2,
1737 . i2=itab(i1),i3=itab(i2))
1738c write(iout,*) '***ERROR: Too many commun Seg with line:',
1739c + Itab(I1),Itab(I2)
1740c CALL ARRET(2)
1741 END IF
1742#endif
1743C----6---------------------------------------------------------------7---------8
1744 RETURN
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889

◆ i25neigh_removeallbut1()

subroutine i25neigh_removeallbut1 ( integer n,
integer, dimension(*) ic,
integer iself,
integer, dimension(4,*) irect,
x,
integer i1,
integer i2,
integer irr )

Definition at line 1478 of file i25neigh.F.

1479C----6---------------------------------------------------------------7---------8
1480C I m p l i c i t T y p e s
1481C-----------------------------------------------
1482#include "implicit_f.inc"
1483C-----------------------------------------------
1484C C o m m o n B l o c k s
1485C-----------------------------------------------
1486C-----------------------------------------------------------------
1487C D u m m y A r g u m e n t s
1488C-----------------------------------------------
1489 INTEGER N,IC(*),ISELF,IRECT(4,*),I1,I2,IE,IRR
1490C REAL
1491 my_real
1492 . x(3,*)
1493C-----------------------------------------------
1494C L o c a l V a r i a b l e s
1495C-----------------------------------------------
1496 INTEGER I,J,J1,J2,INV,ipr
1497 my_real
1498 . s,yj(3,n),yi(3),y0(3),yjni(n),angle(n),smin,norm,
1499 . nxi,nyi,nzi,x12,y12,z12,nxj,nyj,nzj,smax
1500C-----------------------------------------------
1501C --------YI = N_iself ^ 12
1502 CALL norma4n(nxi,nyi,nzi,norm,irect(1,iself) ,x )
1503 x12= x(1,i2)-x(1,i1)
1504 y12= x(2,i2)-x(2,i1)
1505 z12= x(3,i2)-x(3,i1)
1506 yi(1)=nyi*z12-nzi*y12
1507 yi(2)=nzi*x12-nxi*z12
1508 yi(3)=nxi*y12-nyi*x12
1509 CALL normv3(yi,norm)
1510 j=0
1511 DO i=1,n
1512 ie=ic(i)
1513 CALL norma4n(nxj,nyj,nzj,norm,irect(1,ie) ,x )
1514C----YJ = N_ie ^ 21
1515 yj(1,i)=-nyj*z12+nzj*y12
1516 yj(2,i)=-nzj*x12+nxj*z12
1517 yj(3,i)=-nxj*y12+nyj*x12
1518 CALL normv3(yj(1,i),norm)
1519 yjni(i)=nxi*yj(1,i)+nyi*yj(2,i)+nzi*yj(3,i)
1520 IF (yjni(i)>=zero) j=j+1
1521 angle(i)=yi(1)*yj(1,i)+yi(2)*yj(2,i)+yi(3)*yj(3,i)
1522 END DO
1523C
1524 smax=-onep01
1525 j1=0
1526C--------groupe YJ*Ni>=0 :concave keep angle (max_cos) only
1527 DO i=1,n
1528 IF (yjni(i)<zero) cycle
1529 IF (angle(i)>=-one) THEN
1530 IF(smax < angle(i)) THEN
1531 smax=angle(i)
1532 j1=i
1533 END IF
1534 END IF !(ANGLE(I)>=-ONE) THEN
1535 END DO
1536C------angle >180------
1537 IF(j1==0.AND.j >0) THEN
1538 smin=ep10
1539 DO i=1,n
1540 IF (yjni(i)<zero) cycle
1541 IF (smin > angle(i)) THEN
1542 smin=angle(i)
1543 j1=i
1544 END IF
1545 END DO
1546 END IF
1547C--------same side
1548 IF (j==n) THEN
1549C--------only groupe YJ*Ni<0(convex) and no valid one before
1550 ELSEIF(j==0.OR.j1==0) THEN
1551C------angle >180- first-----
1552 smax=-onep01
1553 DO i=1,n
1554 IF (yjni(i)>=zero) cycle
1555 IF(angle(i)< -one .AND.smax < angle(i)) THEN
1556 smax=angle(i)
1557 j1=i
1558 END IF
1559 END DO
1560C ------------------
1561 IF (j1==0) then
1562 smin=ep10
1563 DO i=1,n
1564 IF (yjni(i)>=zero) cycle
1565C--------groupe YJ*Ni<0 :convex keep angle (min_cos) only
1566 IF(angle(i)>= -one .AND. smin > angle(i)) THEN
1567 smin=angle(i)
1568 j1=i
1569 END IF
1570 END DO
1571 END IF !(J1==0) then
1572 END IF !(J==N) then
1573C ------still no valid one-----------
1574 IF (j1==0) then
1575c print *,'***warning** No valid Neighbour Segs of',N,' candidats'
1576 irr = 11
1577 ic(1)=0
1578 ELSE
1579 ic(1)=ic(j1)
1580 END IF
1581 DO i=2,n
1582 ic(i)=0
1583 END DO
1584C----6---------------------------------------------------------------7---------8
1585 RETURN
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
normv3
subroutine normv3(v, norm)
Definition i24tools.F:129
norma4n
subroutine norma4n(n1, n2, n3, area, irect, x)
Definition norma1.F:132

◆ i25neigh_seg_e()

subroutine i25neigh_seg_e ( integer n1,
integer, dimension(*) ied1,
integer n2,
integer, dimension(*) ied2,
integer n,
integer, dimension(*) ic,
integer iself,
integer i1,
integer i2,
integer, dimension(4,*) irect,
integer nrtm,
integer, dimension(*) msegtyp,
integer, dimension(4,*) mvoisin )

Definition at line 1382 of file i25neigh.F.

1384C----6---------------------------------------------------------------7---------8
1385C I m p l i c i t T y p e s
1386C-----------------------------------------------
1387#include "implicit_f.inc"
1388C-----------------------------------------------
1389C C o m m o n B l o c k s
1390C-----------------------------------------------
1391C-----------------------------------------------------------------
1392C D u m m y A r g u m e n t s
1393C-----------------------------------------------
1394 INTEGER N1,IED1(*),N2,IED2(*),N,IC(*),ISELF,
1395 . I1,I2,IRECT(4,*),NRTM,MSEGTYP(*),MVOISIN(4,*)
1396C-----------------------------------------------
1397c FUNCTION: find neighbour segment and which share the same nodes I1,I2
1398c ---neighbour node array will be built after taking into account of treatment w/ IC
1399c (remains only one segment at the end)
1400c
1401c Note:
1402c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
1403c
1404c TYPE NAME FUNCTION
1405c I N1,IED1(N1) - number and neighbour segment id list of node I1
1406c I N2,IED2(N2) - number and neighbour segment id list of node I2
1407c O N,IC(N) - Number and neighbour segment id list of segment id ISELF
1408c I ISELF,I1,I2 - input segment id ISELF w/ nodes I1,I2 (I25NEIGH_un nodes)
1409c I IRECT(4,*) - connectivity of segment id *
1410C-----------------------------------------------
1411C External function
1412C-----------------------------------------------
1413 LOGICAL INTAB,SAME_SEG
1414 EXTERNAL intab,same_seg
1415C-----------------------------------------------
1416C L o c a l V a r i a b l e s
1417C-----------------------------------------------
1418 INTEGER I,J,NEW,K,M,LING,NE,NN,JJ,J1,J2
1419 DATA ling/0/
1420C----add I25NEIGH_un ID--at end--------------------------
1421 IF (i1==i2) THEN
1422C----add 0 in IC as convention ----------------------------
1423 CALL add_id(n,ic,ling)
1424 ELSE
1425 ne=n
1426 DO j=1,n2
1427 new=ied2(j)
1428 IF (new==iself) cycle
1429 DO jj=1,4
1430 IF(i2==irect(jj,new))THEN
1431 j2 = jj
1432 EXIT
1433 END IF
1434 END DO
1435 IF (intab(n1,ied1,new)) THEN
1436 j1 = 0
1437 DO jj=1,4
1438 IF(i1==irect(jj,new))THEN
1439 j1 = jj
1440 EXIT
1441 END IF
1442 END DO
1443 if(j1==0) then
1444 print *,'I25NEIGH_SEG_E - internal error',i1,irect(1:4,new)
1445 end if
1446 IF( (j1-j2==1.OR.
1447 . (irect(3,new)/=irect(4,new).AND.j1-j2==-3).OR.
1448 . (irect(3,new)==irect(4,new).AND.j1-j2==-2))
1449 . .AND. .NOT.same_seg(irect(1,iself),irect(1,new))
1450C
1451C A coating shell can not be neighboring a non coating shell (cf bumper)
1452C
1453 . .AND. .NOT. ((iabs(msegtyp(iself)) > nrtm .AND. iabs(msegtyp(new)) <= nrtm).OR.
1454 . (iabs(msegtyp(iself)) <= nrtm .AND. iabs(msegtyp(new)) > nrtm))) THEN
1455C Triangles ::
1456 IF(irect(3,new)==irect(4,new).AND.j2==3.AND.j1==1) j2=4
1457C
1458C Does not consider NEW if NEW is already neighbor of another segment.
1459 IF(mvoisin(j2,new)==0) CALL add_id(n,ic,new)
1460C
1461 end if
1462 END IF
1463 END DO
1464C----add 0 for IC if find nothing -> consisting w/ ICN------------------------
1465 IF (ne==n) CALL add_id(n,ic,ling)
1466 END IF !(I1==I2) THEN
1467C----6---------------------------------------------------------------7---------8
1468 RETURN
same_seg
logical function same_seg(irect1, irect2)
Definition i24tools.F:329
intab
logical function intab(nic, ic, n)
Definition i24tools.F:95
add_id
subroutine add_id(n, ic, id)
Definition i24tools.F:30

◆ i25neigh_seg_en()

subroutine i25neigh_seg_en ( integer n1,
integer, dimension(*) ied1,
integer n2,
integer, dimension(*) ied2,
integer, dimension(4,*) mvoisin,
integer ne,
integer, dimension(*) ice,
integer iself,
integer i1,
integer i2,
integer, dimension(4,*) irect,
x,
integer, dimension(*) ie,
integer nrtm,
integer, dimension(*) msegtyp,
integer irr )

Definition at line 1278 of file i25neigh.F.

1281C----6---------------------------------------------------------------7---------8
1282C I m p l i c i t T y p e s
1283C-----------------------------------------------
1284#include "implicit_f.inc"
1285C-----------------------------------------------
1286C C o m m o n B l o c k s
1287C-----------------------------------------------
1288C-----------------------------------------------------------------
1289C D u m m y A r g u m e n t s
1290C-----------------------------------------------
1291 INTEGER N1,IED1(*),N2,IED2(*),NE,ICE(*),ISELF,IRR,
1292 . I1,I2,IRECT(4,*),IE(*), NRTM, MSEGTYP(*), MVOISIN(4,*)
1293 my_real
1294 . x(3,*)
1295C-----------------------------------------------
1296c FUNCTION: find neighbour segment and nodes which share the same nodes I1,I2
1297c
1298c Note:
1299c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
1300c
1301c TYPE NAME FUNCTION
1302c I N1,IED1(N1) - number and neighbour segment id list of node I1
1303c I N2,IED2(N2) - number and neighbour segment id list of node I2
1304c O NE,ICE(NE) - Number and neighbour segment id list of segment id ISELF
1305c I ISELF,I1,I2 - input segment id ISELF w/ nodes I1,I2 (I25NEIGH_un nodes)
1306c I IRECT(4,*) - connectivity of segment id *
1307c I X(3,*) - node coordinates
1308c O IE(NE) - final (reduced) neighbour segment,array
1309C-----------------------------------------------
1310C L o c a l V a r i a b l e s
1311C-----------------------------------------------
1312 INTEGER I,J,J1,J2,JJ,NEW,K,M,NE0,DNE,IOP
1313C---------------------
1314 irr = 0
1315 ne0=ne
1316C
1317C MVOISIN already filled due to symetrization ::
1318 IF(ie(ne0+1)/=0) THEN
1319 ne=ne0+1
1320 RETURN
1321 END IF
1322C-------------neighbour segments-----------
1323 CALL i25neigh_seg_e(n1,ied1,n2,ied2,ne,ice,iself,
1324 . i1,i2,irect,nrtm,msegtyp ,mvoisin)
1325C-------------treatment of multi-neighbours (T shape,shell) segments => reduce to 1----------
1326 dne = ne-ne0
1327 IF (dne > 1) THEN
1328 CALL i25neigh_removeallbut1(dne,ice(ne0+1),iself,irect,x ,i1,i2,irr)
1329 ne=ne0+1
1330 END IF
1331 IF (ice(ne)>0 )THEN
1332C
1333C shell symmetric segment will be removed here !
1334 CALL i25neigh_seg_opp(iself,ice(ne),irect,x ,iop)
1335 IF (iop > 0 ) ice(ne) = 0
1336 END IF !(ICE(NE)>0 )THEN
1337C-------------after convention--------------
1338 IF (ne-ne0 > 1) THEN
1339C print *,'!!!error (report to developer)!!!',(NE-NE0)
1340 irr=12
1341 ne=ne0+1
1342 END IF
1343
1344 DO i=1+ne0,ne
1345 ie(i)=ice(i)
1346 END DO
1347
1348 IF(ne-ne0 == 1 .AND. ie(ne)/=0)THEN
1349C
1350C Enforce symmetry ::
1351 DO jj=1,4
1352 IF(i2==irect(jj,ie(ne)))THEN
1353 j2 = jj
1354 EXIT
1355 END IF
1356 END DO
1357 DO jj=1,4
1358 IF(i1==irect(jj,ie(ne)))THEN
1359 j1 = jj
1360 EXIT
1361 END IF
1362 END DO
1363C Triangles ::
1364 IF(irect(3,ie(ne))==irect(4,ie(ne)).AND.j2==3.AND.j1==1) j2=4
1365C
1366 mvoisin(j2,ie(ne))=iself
1367C
1368 END IF
1369C----6---------------------------------------------------------------7---------8
1370 RETURN
i25neigh_removeallbut1
subroutine i25neigh_removeallbut1(n, ic, iself, irect, x, i1, i2, irr)
Definition i25neigh.F:1479
i25neigh_seg_opp
subroutine i25neigh_seg_opp(ei, ej, irect, x, iop)
Definition i25neigh.F:1595

◆ i25neigh_seg_opp()

subroutine i25neigh_seg_opp ( integer ei,
integer ej,
integer, dimension(4,*) irect,
x,
integer iop )

Definition at line 1594 of file i25neigh.F.

1595C----6---------------------------------------------------------------7---------8
1596C I m p l i c i t T y p e s
1597C-----------------------------------------------
1598#include "implicit_f.inc"
1599C-----------------------------------------------
1600C C o m m o n B l o c k s
1601C-----------------------------------------------
1602C-----------------------------------------------------------------
1603C D u m m y A r g u m e n t s
1604C-----------------------------------------------
1605 INTEGER EI,EJ,IRECT(4,*),IOP
1606C REAL
1607 my_real
1608 . x(3,*)
1609C----if the normals of two segments are opposite or near opposite
1610C----if two common segs, will also be eliminated
1611C-----------------------------------------------
1612C External function
1613C-----------------------------------------------
1614 LOGICAL INTAB
1615 EXTERNAL intab
1616C-----------------------------------------------
1617C L o c a l V a r i a b l e s
1618C-----------------------------------------------
1619 INTEGER I,J,NN,IMIN,JMIN,IRMIN,JRMIN,IRI(4),IRJ(4)
1620C-----------------------------------------------
1621 irmin=max(irect(1,ei),irect(2,ei),irect(3,ei),irect(4,ei))
1622 IF(irect(3,ei)/=irect(4,ei))THEN
1623 DO i=1,4
1624 IF(irect(i,ei) <= irmin)THEN
1625 irmin=irect(i,ei)
1626 imin =i
1627 END IF
1628 END DO
1629 iri(1)=irect(imin,ei)
1630 iri(2)=irect(mod(imin ,4)+1,ei)
1631 iri(3)=irect(mod(imin+1,4)+1,ei)
1632 iri(4)=irect(mod(imin+2,4)+1,ei)
1633 ELSE
1634 DO i=1,3
1635 IF(irect(i,ei) <= irmin)THEN
1636 irmin=irect(i,ei)
1637 imin =i
1638 END IF
1639 END DO
1640 iri(1)=irect(imin,ei)
1641 iri(2)=irect(mod(imin ,3)+1,ei)
1642 iri(3)=irect(mod(imin+1,3)+1,ei)
1643 iri(4)=iri(3)
1644 END IF
1645C-----------------------------------------------
1646 jrmin=max(irect(1,ej),irect(2,ej),irect(3,ej),irect(4,ej))
1647 IF(irect(3,ej)/=irect(4,ej))THEN
1648 DO j=1,4
1649 IF(irect(j,ej) <= jrmin)THEN
1650 jrmin=irect(j,ej)
1651 jmin =j
1652 END IF
1653 END DO
1654 irj(1)=irect(jmin,ej)
1655 irj(2)=irect(mod(jmin ,4)+1,ej)
1656 irj(3)=irect(mod(jmin+1,4)+1,ej)
1657 irj(4)=irect(mod(jmin+2,4)+1,ej)
1658 ELSE
1659 DO j=1,3
1660 IF(irect(j,ej) <= jrmin)THEN
1661 jrmin=irect(j,ej)
1662 jmin =j
1663 END IF
1664 END DO
1665 irj(1)=irect(jmin,ej)
1666 irj(2)=irect(mod(jmin ,3)+1,ej)
1667 irj(3)=irect(mod(jmin+1,3)+1,ej)
1668 irj(4)=irj(3)
1669 END IF
1670C-----------------------------------------------
1671 iop=0
1672 IF(irj(3)/=irj(4) .AND. iri(3)/=iri(4))THEN
1673 IF(iri(1)==irj(1))THEN
1674 IF(iri(2)==irj(4))THEN
1675 IF(iri(3)==irj(3))THEN
1676 IF(iri(4)==irj(2))THEN
1677 iop=1
1678 END IF
1679 END IF
1680 END IF
1681 END IF
1682 ELSEIF(irj(3)==irj(4) .AND. iri(3)==iri(4))THEN
1683 IF(iri(1)==irj(1))THEN
1684 IF(iri(2)==irj(3))THEN
1685 IF(iri(3)==irj(2))THEN
1686 iop=1
1687 END IF
1688 END IF
1689 END IF
1690 END IF
1691C----6---------------------------------------------------------------7---------8
1692 RETURN