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, igap, thk_m, thk_m_scale)
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,
integer, intent(in) igap,
intent(in) thk_m,
intent(in) thk_m_scale )

Definition at line 1133 of file i25neigh.F.

1139 use element_mod , only :nixs,nixc,nixtg,nixt,nixp
1140C-----------------------------------------------
1141C I m p l i c i t T y p e s
1142C-----------------------------------------------
1143#include "implicit_f.inc"
1144C-----------------------------------------------
1145C C o m m o n B l o c k s
1146C-----------------------------------------------
1147#include "com01_c.inc"
1148#include "com04_c.inc"
1149#include "param_c.inc"
1150#include "scr17_c.inc"
1151C-----------------------------------------------
1152C D u m m y A r g u m e n t s
1153C-----------------------------------------------
1154 INTEGER NRTM,IRECT(4,*), IXS(NIXS,*), IXC(NIXC,*),
1155 . IXTG(NIXTG,*), IXT(NIXT,*), IXP(NIXP,*),
1156 . IPART(LIPART1,*), IPARTC(*), IPARTTG(*),
1157 . MSR(*),NMN,IGEO(NPROPGI,*),MSEGTYP(*), ITHK25
1158C REAL
1159 my_real
1160 . geo(npropg,*), thk(*),thk_part(*),gap_n(4,*),gap_m(*),
1161 . gapn_m(*),gapmax_m, gapscale,gapmsav(*)
1162 INTEGER, INTENT(IN) :: IGAP
1163 my_real, INTENT(IN) :: thk_m ,thk_m_scale
1164C-----------------------------------------------
1165C L o c a l V a r i a b l e s
1166C-----------------------------------------------
1167 INTEGER I,J,K,IW,I1,I2,I3,MG,M,IP,IGTYP
1168 my_real
1169 . dx
1170 my_real, DIMENSION(:), ALLOCATABLE :: wa
1171C init
1172C
1173 ALLOCATE(wa(numnod))
1174 DO i=1,numnod
1175 wa(i)=zero
1176 END DO
1177C
1178C------------------------------------
1179C GAP NODES IN WA
1180C------------------------------------
1181 DO i=1,numelc
1182 mg=ixc(6,i)
1183 igtyp = igeo(11,mg)
1184 ip = ipartc(i)
1185 IF ( igap == 5.AND.thk_m /= zero) THEN
1186 dx=half*thk_m
1187 ELSEIF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
1188 dx=half*thk_part(ip)
1189 ELSEIF ( thk(i) /= zero .AND. iintthick == 0) THEN
1190 dx=half*thk(i)
1191 ELSEIF(igtyp == 17) THEN
1192 dx=half*thk(i)
1193 ELSE
1194 dx=half*geo(1,mg)
1195 ENDIF
1196 wa(ixc(2,i))=max(wa(ixc(2,i)),dx)
1197 wa(ixc(3,i))=max(wa(ixc(3,i)),dx)
1198 wa(ixc(4,i))=max(wa(ixc(4,i)),dx)
1199 wa(ixc(5,i))=max(wa(ixc(5,i)),dx)
1200 ENDDO
1201C
1202 DO i=1,numeltg
1203 mg=ixtg(5,i)
1204 igtyp = igeo(11,mg)
1205 ip = iparttg(i)
1206 IF ( igap == 5.AND.thk_m /= zero) THEN
1207 dx=half*thk_m
1208 ELSEIF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
1209 dx=half*thk_part(ip)
1210 ELSEIF ( thk(numelc+i) /= zero .AND. iintthick == 0) THEN
1211 dx=half*thk(numelc+i)
1212 ELSEIF(igtyp == 17) THEN
1213 dx=half*thk(numelc+i)
1214 ELSE
1215 dx=half*geo(1,mg)
1216 ENDIF
1217 wa(ixtg(2,i))=max(wa(ixtg(2,i)),dx)
1218 wa(ixtg(3,i))=max(wa(ixtg(3,i)),dx)
1219 wa(ixtg(4,i))=max(wa(ixtg(4,i)),dx)
1220 ENDDO
1221C
1222 IF(igap==5.AND.thk_m /= zero) THEN
1223 DO i=1,numnod
1224 wa(i)=thk_m_scale*wa(i)
1225 END DO
1226 ELSE
1227 DO i=1,numnod
1228 wa(i)=gapscale*wa(i)
1229 END DO
1230 ENDIF
1231C--------exclude lines in main surfaces
1232C DO I=1,NUMELT
1233C MG=IXT(4,I)
1234C DX=HALF*SQRT(GEO(1,MG))
1235C WA(IXT(2,I))=MAX(WA(IXT(2,I)),DX)
1236C WA(IXT(3,I))=MAX(WA(IXT(3,I)),DX)
1237C ENDDO
1238C
1239C DO I=1,NUMELP
1240C MG=IXP(5,I)
1241C DX=HALF*SQRT(GEO(1,MG))
1242C WA(IXP(2,I))=MAX(WA(IXP(2,I)),DX)
1243C WA(IXP(3,I))=MAX(WA(IXP(3,I)),DX)
1244C ENDDO
1245C------------------------------------
1246C INI GAP_N (4), GAP_M is modified taking into account nodal gap for sorting
1247C------------------------------------
1248C -----due to the fact that if surf_M does not contain the defining w/ shell
1249C-------> should not take into account GAP_shell
1250 DO i=1,nrtm
1251 IF (msegtyp(i)==0) THEN
1252 DO j=1,4
1253 m=irect(j,i)
1254 wa(m) = zero
1255 END DO
1256 END IF !(MSEGTYP(I)==0) THEN
1257 END DO ! nrtm
1258C
1259 DO i=1,nmn
1260 m = msr(i)
1261 wa(m) = min(wa(m),gapmax_m)
1262 gapn_m(i) = wa(m)
1263 END DO
1264C
1265 DO i=1,nrtm
1266 gap_m(i) = zero
1267 DO j=1,4
1268 m=irect(j,i)
1269 gap_n(j,i)=wa(m)
1270 gap_m(i) = max(gap_m(i),wa(m))
1271
1272 END DO
1273 END DO ! nrtm
1274C
1275 IF(ithk25==1) THEN
1276 DO i=1,nrtm
1277 gapmsav(i) = gap_m(i)
1278 ENDDO
1279 ENDIF
1280 DEALLOCATE(wa)
1281C
1282 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 in 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 = Local NUE NUMB 1 in NSV (1: NSN)
542C 6: I2 = local number of node 1 in NSV (1: nsn)
543C 7: Type of the edge :
544C <0 <=> the arete is only main, not second
545C +/-1 edge of solid
546C +/-2 edge of shell
547C for engine:
548C 8: Global ID (id in starter)
549C 9: Weight = 1 => secondary second 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 OR SOLIDE COOL
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 of belonging to the areas 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 Calculation of ADDSUBN, LISUBN, INFLG_SUBN for the main nodes
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 Use Kad (1: NMN) to fill 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 Compacting LISUBN AND Combining INFLG_SUBN flags
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 Calculation 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 Use Kad (1: Nedge) to fill 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:1744
i25neigh_seg_e
subroutine i25neigh_seg_e(n1, ied1, n2, ied2, n, ic, iself, i1, i2, irect, nrtm, msegtyp, mvoisin)
Definition i25neigh.F:1400
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:1296
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:1257
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1008

◆ i25neigh_msg_err()

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

Definition at line 1743 of file i25neigh.F.

1744 USE message_mod
1745C-----------------------------------------------
1746C I m p l i c i t T y p e s
1747C-----------------------------------------------
1748#include "implicit_f.inc"
1749C-----------------------------------------------
1750C C o m m o n B l o c k s
1751C-----------------------------------------------
1752#include "scr03_c.inc"
1753C-----------------------------------------------
1754C D u m m y A r g u m e n t s
1755C-----------------------------------------------
1756 INTEGER I1,I2,ITAB(*),IRR
1757C-----------------------------------------------
1758C L o c a l V a r i a b l e s
1759C-----------------------------------------------
1760 INTEGER I
1761C----Warning,ERROR out----------------
1762 IF(ipri==0) RETURN
1763#ifndef HYPERMESH_LIB
1764 IF (irr ==11) THEN
1765C-----multi-neibour but no valid one
1766 CALL ancmsg(msgid=1245,
1767 . msgtype=msgwarning,
1768 . anmode=aninfo_blind_2,
1769 . i2=itab(i1),i3=itab(i2))
1770c write(iout,*) '***Warning: No valid commun Seg with line:',
1771c + Itab(I1),Itab(I2)
1772 ELSEIF (irr ==12) THEN
1773C-----multi-neibour but no valid one
1774 CALL ancmsg(msgid=1246,
1775 . msgtype=msgerror,
1776 . anmode=aninfo_blind_2,
1777 . i2=itab(i1),i3=itab(i2))
1778c write(iout,*) '***ERROR: Too many commun Seg with line:',
1779c + Itab(I1),Itab(I2)
1780c CALL ARRET(2)
1781 END IF
1782#endif
1783C----6---------------------------------------------------------------7---------8
1784 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:895

◆ 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 1518 of file i25neigh.F.

1519C----6---------------------------------------------------------------7---------8
1520C I m p l i c i t T y p e s
1521C-----------------------------------------------
1522#include "implicit_f.inc"
1523C-----------------------------------------------
1524C C o m m o n B l o c k s
1525C-----------------------------------------------
1526C-----------------------------------------------------------------
1527C D u m m y A r g u m e n t s
1528C-----------------------------------------------
1529 INTEGER N,IC(*),ISELF,IRECT(4,*),I1,I2,IE,IRR
1530C REAL
1531 my_real
1532 . x(3,*)
1533C-----------------------------------------------
1534C L o c a l V a r i a b l e s
1535C-----------------------------------------------
1536 INTEGER I,J,J1,J2,INV
1537 my_real
1538 . s,yj(3,n),yi(3),y0(3),yjni(n),angle(n),smin,norm,
1539 . nxi,nyi,nzi,x12,y12,z12,nxj,nyj,nzj,smax
1540C-----------------------------------------------
1541C --------YI = N_iself ^ 12
1542 CALL norma4n(nxi,nyi,nzi,norm,irect(1,iself) ,x )
1543 x12= x(1,i2)-x(1,i1)
1544 y12= x(2,i2)-x(2,i1)
1545 z12= x(3,i2)-x(3,i1)
1546 yi(1)=nyi*z12-nzi*y12
1547 yi(2)=nzi*x12-nxi*z12
1548 yi(3)=nxi*y12-nyi*x12
1549 CALL normv3(yi,norm)
1550 j=0
1551 DO i=1,n
1552 ie=ic(i)
1553 CALL norma4n(nxj,nyj,nzj,norm,irect(1,ie) ,x )
1554C----YJ = N_ie ^ 21
1555 yj(1,i)=-nyj*z12+nzj*y12
1556 yj(2,i)=-nzj*x12+nxj*z12
1557 yj(3,i)=-nxj*y12+nyj*x12
1558 CALL normv3(yj(1,i),norm)
1559 yjni(i)=nxi*yj(1,i)+nyi*yj(2,i)+nzi*yj(3,i)
1560 IF (yjni(i)>=zero) j=j+1
1561 angle(i)=yi(1)*yj(1,i)+yi(2)*yj(2,i)+yi(3)*yj(3,i)
1562 END DO
1563C
1564 smax=-onep01
1565 j1=0
1566C--------groupe YJ*Ni>=0 :concave keep angle (max_cos) only
1567 DO i=1,n
1568 IF (yjni(i)<zero) cycle
1569 IF (angle(i)>=-one) THEN
1570 IF(smax < angle(i)) THEN
1571 smax=angle(i)
1572 j1=i
1573 END IF
1574 END IF !(ANGLE(I)>=-ONE) THEN
1575 END DO
1576C------angle >180------
1577 IF(j1==0.AND.j >0) THEN
1578 smin=ep10
1579 DO i=1,n
1580 IF (yjni(i)<zero) cycle
1581 IF (smin > angle(i)) THEN
1582 smin=angle(i)
1583 j1=i
1584 END IF
1585 END DO
1586 END IF
1587C--------same side
1588 IF (j==n) THEN
1589C--------only groupe YJ*Ni<0(convex) and no valid one before
1590 ELSEIF(j==0.OR.j1==0) THEN
1591C------angle >180- first-----
1592 smax=-onep01
1593 DO i=1,n
1594 IF (yjni(i)>=zero) cycle
1595 IF(angle(i)< -one .AND.smax < angle(i)) THEN
1596 smax=angle(i)
1597 j1=i
1598 END IF
1599 END DO
1600C ------------------
1601 IF (j1==0) then
1602 smin=ep10
1603 DO i=1,n
1604 IF (yjni(i)>=zero) cycle
1605C--------groupe YJ*Ni<0 :convex keep angle (min_cos) only
1606 IF(angle(i)>= -one .AND. smin > angle(i)) THEN
1607 smin=angle(i)
1608 j1=i
1609 END IF
1610 END DO
1611 END IF !(J1==0) then
1612 END IF !(J==N) then
1613C ------still no valid one-----------
1614 IF (j1==0) then
1615c print *,'***warning** No valid Neighbour Segs of',N,' CANDIDATES'
1616 irr = 11
1617 ic(1)=0
1618 ELSE
1619 ic(1)=ic(j1)
1620 END IF
1621 DO i=2,n
1622 ic(i)=0
1623 END DO
1624C----6---------------------------------------------------------------7---------8
1625 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 1398 of file i25neigh.F.

1400C----6---------------------------------------------------------------7---------8
1401C I m p l i c i t T y p e s
1402C-----------------------------------------------
1403#include "implicit_f.inc"
1404C-----------------------------------------------
1405C C o m m o n B l o c k s
1406C-----------------------------------------------
1407C-----------------------------------------------------------------
1408C D u m m y A r g u m e n t s
1409C-----------------------------------------------
1410 INTEGER N1,IED1(*),N2,IED2(*),N,IC(*),ISELF,
1411 . I1,I2,IRECT(4,*),NRTM,MSEGTYP(*),MVOISIN(4,*)
1412C-----------------------------------------------
1413c FUNCTION: find neighbour segment and which share the same nodes I1,I2
1414c ---neighbour node array will be built after taking into account of treatment w/ IC
1415c (remains only one segment at the end)
1416c
1417c Note:
1418c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
1419c
1420c TYPE NAME FUNCTION
1421c I N1,IED1(N1) - number and neighbour segment id list of node I1
1422c I N2,IED2(N2) - number and neighbour segment id list of node I2
1423c O N,IC(N) - Number and neighbour segment id list of segment id ISELF
1424c I ISELF,I1,I2 - input segment id ISELF w/ nodes I1,I2 (I25NEIGH_un nodes)
1425c I IRECT(4,*) - connectivity of segment id *
1426C-----------------------------------------------
1427C External function
1428C-----------------------------------------------
1429 LOGICAL INTAB,SAME_SEG
1430 EXTERNAL intab,same_seg
1431C-----------------------------------------------
1432C L o c a l V a r i a b l e s
1433C-----------------------------------------------
1434 INTEGER I,J,NEW,K,M,LING,NE,NN,JJ,J1,J2,NOK,KK
1435 DATA ling/0/
1436C----add I25NEIGH_un ID--at end--------------------------
1437 IF (i1==i2) THEN
1438C----add 0 in IC as convention ----------------------------
1439 CALL add_id(n,ic,ling)
1440 ELSE
1441 ne=n
1442 DO j=1,n2
1443 new=ied2(j)
1444 IF (new==iself) cycle
1445 DO jj=1,4
1446 IF(i2==irect(jj,new))THEN
1447 j2 = jj
1448 EXIT
1449 END IF
1450 END DO
1451 IF (intab(n1,ied1,new)) THEN
1452 j1 = 0
1453 DO jj=1,4
1454 IF(i1==irect(jj,new))THEN
1455 j1 = jj
1456 EXIT
1457 END IF
1458 END DO
1459 if(j1==0) then
1460 print *,'I25NEIGH_SEG_E - internal error',i1,irect(1:4,new)
1461 end if
1462 IF( (j1-j2==1.OR.
1463 . (irect(3,new)/=irect(4,new).AND.j1-j2==-3).OR.
1464 . (irect(3,new)==irect(4,new).AND.j1-j2==-2))
1465 . .AND. .NOT.same_seg(irect(1,iself),irect(1,new))
1466C
1467C A coating shell can not be neighboring a non coating shell (cf bumper)
1468C
1469 . .AND. .NOT. ((iabs(msegtyp(iself)) > nrtm .AND. iabs(msegtyp(new)) <= nrtm).OR.
1470 . (iabs(msegtyp(iself)) <= nrtm .AND. iabs(msegtyp(new)) > nrtm))) THEN
1471C Triangles ::
1472 IF(irect(3,new)==irect(4,new).AND.j2==3.AND.j1==1) j2=4
1473 nok = 0
1474
1475 IF(irect(3,iself)==irect(4,iself)) THEN
1476 DO jj=1,4
1477 IF(irect(jj,iself)/=i1.AND.irect(jj,iself)/=i2) THEN
1478 DO kk=1,4
1479 IF(kk/=j1.AND.kk/=j2) THEN
1480 IF (irect(jj,iself)==irect(kk,new)) nok = 1
1481 ENDIF
1482 ENDDO
1483 ENDIF
1484 ENDDO
1485 ENDIF
1486 IF(irect(3,new)==irect(4,new)) THEN
1487 DO jj=1,4
1488 IF(jj/=j1.AND.jj/=j2) THEN
1489 DO kk=1,4
1490 IF(irect(kk,iself)/=i1.AND.irect(kk,iself)/=i2) THEN
1491 IF (irect(jj,new)==irect(kk,iself)) nok = 1
1492 ENDIF
1493 ENDDO
1494 ENDIF
1495 ENDDO
1496 ENDIF
1497C
1498C Does not consider NEW if NEW is already neighbor of another segment.
1499 IF(mvoisin(j2,new)==0.AND.nok==0) CALL add_id(n,ic,new)
1500C
1501 end if
1502 END IF
1503 END DO
1504C----add 0 for IC if find nothing -> consisting w/ ICN------------------------
1505 IF (ne==n) CALL add_id(n,ic,ling)
1506 END IF !(I1==I2) THEN
1507C----6---------------------------------------------------------------7---------8
1508 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 1293 of file i25neigh.F.

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

◆ i25neigh_seg_opp()

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

Definition at line 1634 of file i25neigh.F.

1635C----6---------------------------------------------------------------7---------8
1636C I m p l i c i t T y p e s
1637C-----------------------------------------------
1638#include "implicit_f.inc"
1639C-----------------------------------------------
1640C C o m m o n B l o c k s
1641C-----------------------------------------------
1642C-----------------------------------------------------------------
1643C D u m m y A r g u m e n t s
1644C-----------------------------------------------
1645 INTEGER EI,EJ,IRECT(4,*),IOP
1646C REAL
1647 my_real
1648 . x(3,*)
1649C----if the normals of two segments are opposite or near opposite
1650C----if two common segs, will also be eliminated
1651C-----------------------------------------------
1652C External function
1653C-----------------------------------------------
1654 LOGICAL INTAB
1655 EXTERNAL intab
1656C-----------------------------------------------
1657C L o c a l V a r i a b l e s
1658C-----------------------------------------------
1659 INTEGER I,J,NN,IMIN,JMIN,IRMIN,JRMIN,IRI(4),IRJ(4)
1660C-----------------------------------------------
1661 irmin=max(irect(1,ei),irect(2,ei),irect(3,ei),irect(4,ei))
1662 IF(irect(3,ei)/=irect(4,ei))THEN
1663 DO i=1,4
1664 IF(irect(i,ei) <= irmin)THEN
1665 irmin=irect(i,ei)
1666 imin =i
1667 END IF
1668 END DO
1669 iri(1)=irect(imin,ei)
1670 iri(2)=irect(mod(imin ,4)+1,ei)
1671 iri(3)=irect(mod(imin+1,4)+1,ei)
1672 iri(4)=irect(mod(imin+2,4)+1,ei)
1673 ELSE
1674 DO i=1,3
1675 IF(irect(i,ei) <= irmin)THEN
1676 irmin=irect(i,ei)
1677 imin =i
1678 END IF
1679 END DO
1680 iri(1)=irect(imin,ei)
1681 iri(2)=irect(mod(imin ,3)+1,ei)
1682 iri(3)=irect(mod(imin+1,3)+1,ei)
1683 iri(4)=iri(3)
1684 END IF
1685C-----------------------------------------------
1686 jrmin=max(irect(1,ej),irect(2,ej),irect(3,ej),irect(4,ej))
1687 IF(irect(3,ej)/=irect(4,ej))THEN
1688 DO j=1,4
1689 IF(irect(j,ej) <= jrmin)THEN
1690 jrmin=irect(j,ej)
1691 jmin =j
1692 END IF
1693 END DO
1694 irj(1)=irect(jmin,ej)
1695 irj(2)=irect(mod(jmin ,4)+1,ej)
1696 irj(3)=irect(mod(jmin+1,4)+1,ej)
1697 irj(4)=irect(mod(jmin+2,4)+1,ej)
1698 ELSE
1699 DO j=1,3
1700 IF(irect(j,ej) <= jrmin)THEN
1701 jrmin=irect(j,ej)
1702 jmin =j
1703 END IF
1704 END DO
1705 irj(1)=irect(jmin,ej)
1706 irj(2)=irect(mod(jmin ,3)+1,ej)
1707 irj(3)=irect(mod(jmin+1,3)+1,ej)
1708 irj(4)=irj(3)
1709 END IF
1710C-----------------------------------------------
1711 iop=0
1712 IF(irj(3)/=irj(4) .AND. iri(3)/=iri(4))THEN
1713 IF(iri(1)==irj(1))THEN
1714 IF(iri(2)==irj(4))THEN
1715 IF(iri(3)==irj(3))THEN
1716 IF(iri(4)==irj(2))THEN
1717 iop=1
1718 END IF
1719 END IF
1720 END IF
1721 END IF
1722 ELSEIF(irj(3)==irj(4) .AND. iri(3)==iri(4))THEN
1723 IF(iri(1)==irj(1))THEN
1724 IF(iri(2)==irj(3))THEN
1725 IF(iri(3)==irj(2))THEN
1726 iop=1
1727 END IF
1728 END IF
1729 END IF
1730 END IF
1731C----6---------------------------------------------------------------7---------8
1732 RETURN