1 | MODULE u_polyg |
---|
2 | !--------------------------------------------------------------------- |
---|
3 | !- utilitaires relatifs aux polygones |
---|
4 | !--------------------------------------------------------------------- |
---|
5 | USE poly_types |
---|
6 | USE m_g2d |
---|
7 | USE mt_polyg |
---|
8 | !- |
---|
9 | IMPLICIT NONE |
---|
10 | !- |
---|
11 | !$ real,parameter :: v_eps = 1.e-4 |
---|
12 | REAL (kind=rp) :: v_eps = 100.0_rp *EPSILON(v_eps) |
---|
13 | !--------------------------------------------------------------------- |
---|
14 | CONTAINS |
---|
15 | != |
---|
16 | SUBROUTINE p_dec (pol,npd,pd) |
---|
17 | !--------------------------------------------------------------------- |
---|
18 | !- decomposition d'un polygone quelconque en elements convexes |
---|
19 | !--------------------------------------------------------------------- |
---|
20 | TYPE(polyg) :: pol |
---|
21 | INTEGER :: npd |
---|
22 | TYPE(polyg),DIMENSION(:) :: pd |
---|
23 | !- |
---|
24 | TYPE(polyg) :: pv |
---|
25 | INTEGER :: nppol |
---|
26 | INTEGER,DIMENSION(pol%n) :: kpt |
---|
27 | INTEGER,DIMENSION(pol%n) :: tabpt |
---|
28 | REAL (kind=rp), DIMENSION(pol%n) :: tabcr |
---|
29 | !- elements relatifs aux croisements |
---|
30 | TYPE :: croisement |
---|
31 | REAL (kind=rp) :: c1,c2 |
---|
32 | INTEGER :: k1,k2,kp |
---|
33 | END TYPE croisement |
---|
34 | INTEGER :: nvcr |
---|
35 | TYPE(croisement),DIMENSION(3*pol%n) :: tvcr |
---|
36 | !- elements relatifs aux segments |
---|
37 | TYPE :: segment |
---|
38 | INTEGER :: ideb,ifin |
---|
39 | END TYPE segment |
---|
40 | INTEGER :: nsgm |
---|
41 | TYPE(segment),DIMENSION(7*pol%n) :: tsgm |
---|
42 | !- |
---|
43 | INTEGER :: npt |
---|
44 | TYPE(c2d),DIMENSION(5*pol%n) :: tpt |
---|
45 | !- |
---|
46 | INTEGER :: n1,n2,n3,ii,ki,kp |
---|
47 | INTEGER :: ipd,ipc,ib,npw |
---|
48 | REAL (kind=rp) :: c1,c2,cb,alfa |
---|
49 | TYPE(c2d) :: d1,d2,ptn |
---|
50 | INTRINSIC huge |
---|
51 | !--------------------------------------------------------------------- |
---|
52 | npd = 0; pv = pol; |
---|
53 | !- reduction du polygone --------------------------------------------- |
---|
54 | CALL p_red (pv) |
---|
55 | !- obtention des elements constitutifs du polygone ------------------- |
---|
56 | nppol = 0; npt = 0; |
---|
57 | IF (pv%n > 2) THEN |
---|
58 | DO n1=1,pv%n |
---|
59 | nppol = nppol+1; ptn = pv%pt(n1); |
---|
60 | CALL p_apt (ptn,npt,tpt,kp) |
---|
61 | kpt(nppol) = kp; |
---|
62 | ENDDO |
---|
63 | ENDIF |
---|
64 | !- recherche des croisements ----------------------------------------- |
---|
65 | nvcr = 0; |
---|
66 | DO n1=1,nppol-2 |
---|
67 | d1 = tpt(kpt(n1+1))-tpt(kpt(n1)); |
---|
68 | DO n2=n1+2,nppol+MIN(n1-2,0) |
---|
69 | d2 = tpt(kpt(MOD(n2,nppol)+1))-tpt(kpt(n2)); |
---|
70 | CALL p_i2d (tpt(kpt(n1)),d1,tpt(kpt(n2)),d2,c1,c2,ii) |
---|
71 | IF (ii > 0) THEN |
---|
72 | IF ( ( (c1 > v_eps).AND.(c1 < (1.-v_eps)) & |
---|
73 | .AND.(c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) & |
---|
74 | .OR.( (c1 > -v_eps).AND.(c1 < (1.+v_eps)) & |
---|
75 | .AND.(c2 > v_eps).AND.(c2 < (1.-v_eps)) ) ) THEN |
---|
76 | nvcr = nvcr+1; |
---|
77 | ptn = tpt(kpt(n2))+c2*d2 |
---|
78 | CALL p_apt (ptn,npt,tpt,kp) |
---|
79 | tvcr(nvcr) = croisement(c1,c2,n1,n2,kp) |
---|
80 | ENDIF |
---|
81 | ENDIF |
---|
82 | ENDDO |
---|
83 | ENDDO |
---|
84 | !- creation des segments avec inclusion des croisements -------------- |
---|
85 | nsgm = 0; |
---|
86 | DO n1=1,nppol |
---|
87 | !--- creation des points |
---|
88 | ki = 1; tabcr(ki) = 0.; tabpt(ki) = kpt(n1); |
---|
89 | DO n2=1,nvcr |
---|
90 | IF (n1 == tvcr(n2)%k1) THEN |
---|
91 | ki = ki+1; tabcr(ki) = tvcr(n2)%c1; tabpt(ki) = tvcr(n2)%kp; |
---|
92 | ELSE IF (n1 == tvcr(n2)%k2) THEN |
---|
93 | ki = ki+1; tabcr(ki) = tvcr(n2)%c2; tabpt(ki) = tvcr(n2)%kp; |
---|
94 | ENDIF |
---|
95 | ENDDO |
---|
96 | ki = ki+1; tabcr(ki) = 1.; tabpt(ki) = kpt(MOD(n1,nppol)+1); |
---|
97 | !--- ordonnancement des points |
---|
98 | DO n3=3,ki-1 |
---|
99 | DO n2=2,n3-1 |
---|
100 | IF (tabcr(n3) < tabcr(n2)) THEN |
---|
101 | tabcr(n2:n3) = (/ tabcr(n3), tabcr(n2:n3-1) /) |
---|
102 | tabpt(n2:n3) = (/ tabpt(n3), tabpt(n2:n3-1) /) |
---|
103 | EXIT |
---|
104 | ENDIF |
---|
105 | ENDDO |
---|
106 | ENDDO |
---|
107 | !--- creation des segments |
---|
108 | DO n2=1,ki-1 |
---|
109 | IF (ABS(tabcr(n2+1)-tabcr(n2)) > v_eps) THEN |
---|
110 | nsgm = nsgm+1; |
---|
111 | alfa = 0.5_rp*(tabcr(n2)+tabcr(n2+1)) |
---|
112 | IF (p_qos(nppol,kpt,tpt,n1,alfa)) THEN |
---|
113 | tsgm(nsgm) = segment(tabpt(n2),tabpt(n2+1)); |
---|
114 | ELSE |
---|
115 | tsgm(nsgm) = segment(tabpt(n2+1),tabpt(n2)); |
---|
116 | ENDIF |
---|
117 | ENDIF |
---|
118 | ENDDO |
---|
119 | ENDDO |
---|
120 | !- creation des polygones sans croisement par chainage des segments |
---|
121 | l_p: & |
---|
122 | DO n1=1,nsgm |
---|
123 | !--- recherche du prochain segment disponible |
---|
124 | IF (tsgm(n1)%ideb <= 0) CYCLE l_p |
---|
125 | !--- creation d'un nouveau polygone |
---|
126 | ipd = tsgm(n1)%ideb; ipc = tsgm(n1)%ifin; |
---|
127 | tsgm(n1)%ideb = -ipd; |
---|
128 | pv%n = 2; pv%pt(1:2) = (/ tpt(ipd),tpt(ipc) /); |
---|
129 | DO |
---|
130 | !----- recherche du segment suivant |
---|
131 | ib = 0; cb = HUGE(cb); |
---|
132 | DO n2=n1+1,nsgm |
---|
133 | IF (tsgm(n2)%ideb == ipc) THEN |
---|
134 | c1 = (tpt(tsgm(n2)%ifin)-pv%pt(pv%n)) & |
---|
135 | .a.(pv%pt(pv%n-1)-pv%pt(pv%n)) |
---|
136 | IF (ABS(c1) > v_eps) THEN |
---|
137 | IF ( (ib == 0) & |
---|
138 | .OR.((cb < 0).AND.(c1 > cb)) & |
---|
139 | .OR.((cb > 0).AND.(c1 > 0).AND.(c1 < cb)) ) THEN |
---|
140 | ib = n2; cb = c1; |
---|
141 | ENDIF |
---|
142 | ENDIF |
---|
143 | ENDIF |
---|
144 | ENDDO |
---|
145 | IF (ib == 0) THEN |
---|
146 | WRITE (unit=*,fmt='(" chainage interrompu !!!")'); |
---|
147 | CYCLE l_p; |
---|
148 | ENDIF |
---|
149 | !----- |
---|
150 | tsgm(ib)%ideb = -ipc; ipc = tsgm(ib)%ifin; |
---|
151 | IF (ipc /= ipd) THEN |
---|
152 | !------- extension du polygone |
---|
153 | pv%n = pv%n+1; pv%pt(pv%n) = tpt(ipc); |
---|
154 | ELSE |
---|
155 | !------- traitement des concavites |
---|
156 | CALL p_dcv (pv,npw,pd(npd+1:)) |
---|
157 | npd = npd+npw |
---|
158 | EXIT |
---|
159 | ENDIF |
---|
160 | ENDDO |
---|
161 | ENDDO l_p |
---|
162 | !---------------------- |
---|
163 | END SUBROUTINE p_dec |
---|
164 | != |
---|
165 | SUBROUTINE p_dcv (pg,npd,pd) |
---|
166 | !--------------------------------------------------------------------- |
---|
167 | !- decomposition d'un polygone sans croisement en polygones convexes |
---|
168 | !--------------------------------------------------------------------- |
---|
169 | TYPE(polyg) :: pg |
---|
170 | INTEGER :: npd |
---|
171 | TYPE(polyg),DIMENSION(:) :: pd |
---|
172 | !- |
---|
173 | TYPE(polyg) :: pw |
---|
174 | TYPE(polyg),DIMENSION((pg%n+1)/2) :: pp |
---|
175 | INTEGER :: jc,jp,js,np |
---|
176 | INTEGER :: ks,k0,k1,k2 |
---|
177 | INTEGER :: jw,ii |
---|
178 | REAL (kind=rp) :: c1,c2,cb |
---|
179 | TYPE(c2d) :: d1,d2,ptn |
---|
180 | INTRINSIC huge |
---|
181 | !--------------------------------------------------------------------- |
---|
182 | npd = 0; |
---|
183 | !- memorisation et reduction du polygone |
---|
184 | np = 1; pp(1)%pt(1:pg%n) = pg%pt(1:pg%n); pp(1)%n = pg%n; |
---|
185 | CALL p_red (pp(1)) |
---|
186 | !- recherche et des traitements des concavites |
---|
187 | DO |
---|
188 | IF (np <= 0) EXIT |
---|
189 | jc = 0; |
---|
190 | DO |
---|
191 | pw = pp(np) |
---|
192 | jc = jc+1; |
---|
193 | IF ( (pw%n <= 3).OR.(jc > pw%n) ) EXIT |
---|
194 | jp = MOD(jc-2+pw%n,pw%n)+1; js = MOD(jc,pw%n)+1; |
---|
195 | IF ( ((pw%pt(js)-pw%pt(jc)).a.(pw%pt(jp)-pw%pt(jc))) & |
---|
196 | < -v_eps) THEN |
---|
197 | !------- traitement d'un rebroussement |
---|
198 | !------- partage du polygone pw par le segment (jc,jc+1) |
---|
199 | d1 = pw%pt(MOD(jc,pw%n)+1)-pw%pt(jc); |
---|
200 | !------- recherche de l'intersection la plus proche |
---|
201 | ks = 0; cb = HUGE(cb); |
---|
202 | DO jw=2,pw%n-2 |
---|
203 | k0 = MOD(jc+jw-1,pw%n)+1 |
---|
204 | d2 = pw%pt(MOD(k0,pw%n)+1)-pw%pt(k0) |
---|
205 | CALL p_i2d (pw%pt(jc),d1,pw%pt(k0),d2,c1,c2,ii) |
---|
206 | IF (ii > 0) THEN |
---|
207 | IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
208 | IF (c1 < v_eps) THEN |
---|
209 | IF (ABS(c1) < ABS(cb)) THEN |
---|
210 | ks = k0; ptn = pw%pt(k0)+c2*d2; cb = c1; |
---|
211 | ENDIF |
---|
212 | ENDIF |
---|
213 | ENDIF |
---|
214 | ENDIF |
---|
215 | ENDDO |
---|
216 | IF (ks /= 0) THEN |
---|
217 | !--------- polygones resultants |
---|
218 | k1 = MIN(jc,ks); k2 = MAX(jc,ks); |
---|
219 | pp(np)%n = k2-k1+1; |
---|
220 | pp(np)%pt(1:pp(np)%n) = (/ ptn, pw%pt(k1+1:k2) /); |
---|
221 | CALL p_red (pp(np)); |
---|
222 | np = np+1; |
---|
223 | pp(np)%n = k1+1+pw%n-k2; |
---|
224 | pp(np)%pt(1:pp(np)%n) = & |
---|
225 | (/ pw%pt(1:k1), ptn, pw%pt(k2+1:pw%n) /); |
---|
226 | CALL p_red (pp(np)); |
---|
227 | jc = 0; |
---|
228 | ELSE |
---|
229 | WRITE (unit=*,fmt='(" echec !!!")'); STOP |
---|
230 | ENDIF |
---|
231 | ENDIF |
---|
232 | ENDDO |
---|
233 | npd = npd+1; pd(npd) = pw; np = np-1; |
---|
234 | ENDDO |
---|
235 | !---------------------- |
---|
236 | END SUBROUTINE p_dcv |
---|
237 | != |
---|
238 | SUBROUTINE p_apt (point,npt,tpt,kp) |
---|
239 | !--------------------------------------------------------------------- |
---|
240 | !- recherche de l'index d'un point dans un tableau |
---|
241 | !- creation de l'element si il n'existe pas |
---|
242 | !--------------------------------------------------------------------- |
---|
243 | TYPE(c2d) :: point |
---|
244 | INTEGER :: npt |
---|
245 | TYPE(c2d),DIMENSION(*) :: tpt |
---|
246 | INTEGER :: kp |
---|
247 | !- |
---|
248 | INTEGER :: i |
---|
249 | TYPE(c2d) :: ddp |
---|
250 | !--------------------------------------------------------------------- |
---|
251 | kp = 0 |
---|
252 | DO i=1,npt |
---|
253 | ddp = point-tpt(i) |
---|
254 | IF ((ddp.s.ddp) < v_eps) THEN |
---|
255 | kp = i |
---|
256 | EXIT |
---|
257 | ENDIF |
---|
258 | ENDDO |
---|
259 | IF (kp == 0) THEN |
---|
260 | npt = npt+1; kp = npt; tpt(kp) = point; |
---|
261 | ENDIF |
---|
262 | !---------------------- |
---|
263 | END SUBROUTINE p_apt |
---|
264 | != |
---|
265 | SUBROUTINE p_i2d (p1,d1,p2,d2,c1,c2,ii) |
---|
266 | !--------------------------------------------------------------------- |
---|
267 | !- recherche de l'intersection de deux droites |
---|
268 | !- |
---|
269 | !- 1 droites secantes |
---|
270 | !- ii = 0 droites paralleles distinctes |
---|
271 | !- -1 droites confondues |
---|
272 | !--------------------------------------------------------------------- |
---|
273 | TYPE(c2d) :: p1,d1,p2,d2 |
---|
274 | REAL (kind=rp) :: c1,c2 |
---|
275 | INTEGER :: ii |
---|
276 | !- |
---|
277 | REAL (kind=rp) :: det |
---|
278 | TYPE(c2d) :: vp,vn |
---|
279 | !--------------------------------------------------------------------- |
---|
280 | det = d1.v.d2; vp = p2-p1; |
---|
281 | IF (ABS(det) > v_eps) THEN |
---|
282 | ii = 1; |
---|
283 | c1 = (vp.v.d2)/det; c2 = (vp.v.d1)/det; |
---|
284 | ELSE |
---|
285 | ii = 0; |
---|
286 | vn = c2d(-d1%y,d1%x); |
---|
287 | det = vn.v.d2; |
---|
288 | IF (ABS(det) > v_eps) THEN |
---|
289 | c1 = (vp.v.d2)/det; c2 = (vp.v.vn)/det; |
---|
290 | IF (ABS(c1) < v_eps) THEN |
---|
291 | ii = -1; |
---|
292 | ENDIF |
---|
293 | ENDIF |
---|
294 | ENDIF |
---|
295 | !---------------------- |
---|
296 | END SUBROUTINE p_i2d |
---|
297 | != |
---|
298 | SUBROUTINE p_orp (pw) |
---|
299 | !--------------------------------------------------------------------- |
---|
300 | !- orientation d'un polygone (sans croisement) |
---|
301 | !--------------------------------------------------------------------- |
---|
302 | !- comptage des intersections de la demi-mediatrice orientee |
---|
303 | !- du segment (pw%n,1) avec le polygone. |
---|
304 | !- si ce nombre est pair, la demi-mediatrice est exterieure, |
---|
305 | !- et on inverse le sens du polygone. |
---|
306 | !--------------------------------------------------------------------- |
---|
307 | TYPE(polyg) :: pw |
---|
308 | !- |
---|
309 | INTEGER :: i,ii,nic,issp,issc |
---|
310 | REAL (kind=rp) :: c1,c2 |
---|
311 | TYPE(c2d) :: pm,vt,vn,vw |
---|
312 | |
---|
313 | !$$$ integer :: ntoto = 0 |
---|
314 | !--------------------------------------------------------------------- |
---|
315 | pm = 0.5_rp*(pw%pt(pw%n)+pw%pt(1)); |
---|
316 | vt = pw%pt(1)-pw%pt(pw%n); |
---|
317 | vn = c2d(-vt%y,vt%x); |
---|
318 | !- |
---|
319 | nic = 0; issp = 0; |
---|
320 | DO i=1,pw%n-1 |
---|
321 | vw = pw%pt(i+1)-pw%pt(i); |
---|
322 | CALL p_i2d (pm,vn,pw%pt(i),vw,c1,c2,ii) |
---|
323 | IF (ii > 0) THEN |
---|
324 | !$$$$$ if (c1 > 0.) then |
---|
325 | IF (c1 > -v_eps) THEN |
---|
326 | IF ( (c2 > v_eps).AND.(c2 < (1.-v_eps)) ) THEN |
---|
327 | !--------- la normale passe par l'interieur du cote (i,i+1) |
---|
328 | nic = nic+1 |
---|
329 | issp = 0; |
---|
330 | ELSE IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
331 | !--------- la normale passe par le sommet i |
---|
332 | issc = NINT(SIGN(1.,vw.s.vt)); |
---|
333 | IF ((issc*issp) > 0) THEN |
---|
334 | nic = nic+1 |
---|
335 | ENDIF |
---|
336 | issp = issc; |
---|
337 | ENDIF |
---|
338 | ENDIF |
---|
339 | ENDIF |
---|
340 | ENDDO |
---|
341 | !- |
---|
342 | IF (MOD(nic,2) == 0) THEN |
---|
343 | !$$$ write (unit=*,fmt='(" inversion du sens du polygone")'); |
---|
344 | !$$$ ntoto = ntoto + 1 |
---|
345 | !$$$ write (unit=*,fmt=*) pw%pt(1:pw%n:1) |
---|
346 | !$$$ if ( ntoto > 4 ) stop |
---|
347 | pw%pt(1:pw%n:1) = pw%pt(pw%n:1:-1); |
---|
348 | !$$$ else |
---|
349 | !$$$ write (unit=*,fmt='(" sens correct")'); |
---|
350 | ENDIF |
---|
351 | !---------------------- |
---|
352 | END SUBROUTINE p_orp |
---|
353 | != |
---|
354 | LOGICAL FUNCTION p_qos (nppol,kpt,tpt,isc,alfa) |
---|
355 | !--------------------------------------------------------------------- |
---|
356 | !- test d'orientation d'un segment de polygone |
---|
357 | !--------------------------------------------------------------------- |
---|
358 | !- comptage des intersections de la demi-mediatrice orientee |
---|
359 | !- du segment (isc,isc+alpha*longueur(isc)) avec le polygone. |
---|
360 | !- si ce nombre est pair, la demi-mediatrice est exterieure, |
---|
361 | !- et le polygone est mal oriente. |
---|
362 | !--------------------------------------------------------------------- |
---|
363 | INTEGER :: nppol |
---|
364 | INTEGER,DIMENSION(*) :: kpt |
---|
365 | TYPE(c2d),DIMENSION(*) :: tpt |
---|
366 | INTEGER :: isc |
---|
367 | REAL (kind=rp) :: alfa |
---|
368 | !- |
---|
369 | INTEGER :: i,ii,nic,issp,issc,isd |
---|
370 | REAL (kind=rp) :: c1,c2 |
---|
371 | TYPE(c2d) :: ptm,pts,normale,segment |
---|
372 | !--------------------------------------------------------------------- |
---|
373 | ptm = tpt(kpt(isc))+alfa*(tpt(kpt(MOD(isc,nppol)+1))-tpt(kpt(isc))); |
---|
374 | segment = ptm-tpt(kpt(isc)); |
---|
375 | normale%x = -segment%y; normale%y = segment%x; |
---|
376 | !- |
---|
377 | nic = 0; issp = 0; |
---|
378 | DO i=0,nppol-2 |
---|
379 | isd = MOD(isc+i,nppol)+1 |
---|
380 | pts = tpt(kpt(MOD(isd,nppol)+1))-tpt(kpt(isd)); |
---|
381 | CALL p_i2d (ptm,normale,tpt(kpt(isd)),pts,c1,c2,ii) |
---|
382 | IF (ii > 0) THEN |
---|
383 | IF (c1 > 0.) THEN |
---|
384 | IF ( (c2 > v_eps).AND.(c2 < (1.-v_eps)) ) THEN |
---|
385 | !--------- la normale passe par l'interieur du cote (isd,isd+1) |
---|
386 | nic = nic+1 |
---|
387 | issp = 0; |
---|
388 | ELSE IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
389 | !--------- la normale passe par le sommet isd |
---|
390 | issc = NINT(SIGN(1.,pts.s.segment)); |
---|
391 | IF ((issc*issp) > 0) THEN |
---|
392 | nic = nic+1 |
---|
393 | ENDIF |
---|
394 | issp = issc; |
---|
395 | ENDIF |
---|
396 | ENDIF |
---|
397 | ENDIF |
---|
398 | ENDDO |
---|
399 | !- |
---|
400 | p_qos = (MOD(nic,2) /= 0) |
---|
401 | !-------------------- |
---|
402 | END FUNCTION p_qos |
---|
403 | != |
---|
404 | SUBROUTINE p_pcp (p1,p2,p3) |
---|
405 | !--------------------------------------------------------------------- |
---|
406 | !- polygone commun a 2 polygones convexes orientes |
---|
407 | !--------------------------------------------------------------------- |
---|
408 | TYPE(polyg) :: p1,p2,p3 |
---|
409 | !- |
---|
410 | TYPE(polyg),DIMENSION(2) :: pw |
---|
411 | INTEGER :: j1d,j1f,j2d,j2f |
---|
412 | INTEGER :: kc,kr |
---|
413 | INTEGER :: ipd,ipf |
---|
414 | TYPE(c2d) :: v1,v2 |
---|
415 | REAL (kind=rp) :: ww |
---|
416 | !--------------------------------------------------------------------- |
---|
417 | !- copie du polygone p1 dans le polygone candidat |
---|
418 | kc = 1; pw(kc)%n = p1%n; pw(kc)%pt(1:pw(kc)%n) = p1%pt(1:pw(kc)%n); |
---|
419 | !- decoupage du polygone candidat par chaque cote de p2 |
---|
420 | j2d = p2%n; |
---|
421 | l_r: & |
---|
422 | DO j2f=1,p2%n |
---|
423 | !--- definition du polygone resultat |
---|
424 | kr = 3-kc; pw(kr)%n = 0; |
---|
425 | !--- definition du vecteur decoupant |
---|
426 | v2 = p2%pt(j2f)-p2%pt(j2d); |
---|
427 | !--- traitement des cotes successifs du polygone candidat |
---|
428 | j1d = pw(kc)%n; |
---|
429 | !--- position relative (origine segment candidat)/(segment decoupant) |
---|
430 | ipd = n_sgn(v2.v.(pw(kc)%pt(j1d)-p2%pt(j2d))) |
---|
431 | DO j1f=1,pw(kc)%n |
---|
432 | !----- definition du vecteur candidat |
---|
433 | v1 = pw(kc)%pt(j1f)-pw(kc)%pt(j1d); |
---|
434 | !----- position relative (fin segment candidat)/(segment decoupant) |
---|
435 | ipf = n_sgn(v2.v.(pw(kc)%pt(j1f)-p2%pt(j2d))) |
---|
436 | !----- si l'origine appartient au cote decoupant ou est du bon cote, |
---|
437 | !----- on l'insere dans le nouveau polygone candidat |
---|
438 | IF (ipd >= 0) THEN |
---|
439 | pw(kr)%n = pw(kr)%n+1 |
---|
440 | pw(kr)%pt(pw(kr)%n) = pw(kc)%pt(j1d) |
---|
441 | ENDIF |
---|
442 | !----- si le segment candidat est coupe par le segment decoupant, |
---|
443 | !----- on introduit l'intersection dans le polygone resultat |
---|
444 | IF ((ipd*ipf) < 0) THEN |
---|
445 | pw(kr)%n = pw(kr)%n+1 |
---|
446 | ww = (v2.v.(pw(kc)%pt(j1d)-p2%pt(j2d)))/(v1.v.v2) |
---|
447 | pw(kr)%pt(pw(kr)%n) = pw(kc)%pt(j1d)+ww*v1 |
---|
448 | ENDIF |
---|
449 | !----- l'extremite du segment sera l'origine du segment suivant |
---|
450 | ipd = ipf; j1d = j1f; |
---|
451 | ENDDO |
---|
452 | !--- si le polygone resultat est degenere, il n'y a pas de recouvrement |
---|
453 | IF (pw(kr)%n < 3) EXIT l_r |
---|
454 | !--- le polygone resultat devient polygone candidat |
---|
455 | kc = kr; |
---|
456 | !--- passage au cote suivant de p2 |
---|
457 | j2d = j2f; |
---|
458 | ENDDO l_r |
---|
459 | !- fin de decoupage : rangement du resultat dans p3 |
---|
460 | IF (pw(kr)%n < 3) THEN |
---|
461 | p3%n = 0 |
---|
462 | ELSE |
---|
463 | p3%n = pw(kr)%n |
---|
464 | p3%pt(1:p3%n) = pw(kr)%pt(1:pw(kr)%n) |
---|
465 | ENDIF |
---|
466 | !---------------------- |
---|
467 | END SUBROUTINE p_pcp |
---|
468 | != |
---|
469 | SUBROUTINE p_red (p0) |
---|
470 | !--------------------------------------------------------------------- |
---|
471 | !- reduction d'un polygone |
---|
472 | !- suppression des doublons et alignements |
---|
473 | !--------------------------------------------------------------------- |
---|
474 | TYPE(polyg) :: p0 |
---|
475 | !- |
---|
476 | INTEGER :: ip,ic,is |
---|
477 | TYPE(c2d) :: dp,ds |
---|
478 | !--------------------------------------------------------------------- |
---|
479 | ic = 1; |
---|
480 | DO |
---|
481 | IF ( (p0%n < 3).OR.(ic > p0%n) ) EXIT |
---|
482 | ip = MOD(ic-2+p0%n,p0%n)+1; is = MOD(ic,p0%n)+1; |
---|
483 | dp = p0%pt(ic)-p0%pt(ip); ds = p0%pt(is)-p0%pt(ic); |
---|
484 | !$$$ if ((dp.s.dp) < v_eps) then |
---|
485 | !$$$ p0%pt(ip:p0%n-1) = p0%pt(ip+1:p0%n); p0%n = p0%n-1; |
---|
486 | !$$$ else if ((ds.s.ds) < v_eps) then |
---|
487 | !$$$ p0%pt(is:p0%n-1) = p0%pt(is+1:p0%n); p0%n = p0%n-1; |
---|
488 | !$$$ else if (abs(dp.v.ds) < v_eps) then |
---|
489 | IF (ABS(dp.v.ds) < v_eps) THEN |
---|
490 | p0%pt(ic:p0%n-1) = p0%pt(ic+1:p0%n); p0%n = p0%n-1; |
---|
491 | ic = MAX(ic-1,1); |
---|
492 | ELSE |
---|
493 | ic = ic+1; |
---|
494 | ENDIF |
---|
495 | ENDDO |
---|
496 | !---------------------- |
---|
497 | END SUBROUTINE p_red |
---|
498 | != |
---|
499 | INTEGER FUNCTION n_sgn (v) |
---|
500 | !--------------------------------------------------------------------- |
---|
501 | REAL (kind=rp):: v |
---|
502 | !--------------------------------------------------------------------- |
---|
503 | IF (v > v_eps) THEN |
---|
504 | n_sgn = 1 |
---|
505 | ELSE IF (v < -v_eps) THEN |
---|
506 | n_sgn = -1 |
---|
507 | ELSE |
---|
508 | n_sgn = 0 |
---|
509 | ENDIF |
---|
510 | !-------------------- |
---|
511 | END FUNCTION n_sgn |
---|
512 | != |
---|
513 | FUNCTION p_srf (pw) |
---|
514 | !--------------------------------------------------------------------- |
---|
515 | !- calcul de la surface d'un polygone (sans croisement) |
---|
516 | !--------------------------------------------------------------------- |
---|
517 | !- sommation des surfaces des triangles orientes composant le polygone |
---|
518 | !--------------------------------------------------------------------- |
---|
519 | REAL (kind=rp) :: p_srf |
---|
520 | TYPE(polyg) :: pw |
---|
521 | !- |
---|
522 | INTEGER :: j1,j2,j3,jp |
---|
523 | !--------------------------------------------------------------------- |
---|
524 | p_srf = 0.0_rp ; j1 = 1; |
---|
525 | DO jp=1,pw%n-2 |
---|
526 | j2 = jp+1; j3 = jp+2; |
---|
527 | p_srf = p_srf+((pw%pt(j2)-pw%pt(j1)).v.(pw%pt(j3)-pw%pt(j1))) |
---|
528 | ENDDO |
---|
529 | p_srf = 0.5_rp * ABS(p_srf) |
---|
530 | !-------------------- |
---|
531 | END FUNCTION p_srf |
---|
532 | != |
---|
533 | INTEGER FUNCTION p_tap (pq,pw) |
---|
534 | !--------------------------------------------------------------------- |
---|
535 | !- test d'appartenance d'un point a un polygone |
---|
536 | !--------------------------------------------------------------------- |
---|
537 | !- calcul du nombre d'intersections de la demi-droite perpendiculaire |
---|
538 | !- au premier cote avec les autres cotes. |
---|
539 | !- si ce nombre est pair, le point appartient au polygone. |
---|
540 | !--------------------------------------------------------------------- |
---|
541 | TYPE(c2d) :: pq |
---|
542 | TYPE(polyg) :: pw |
---|
543 | !- |
---|
544 | INTEGER :: jp,ii,nic,issp,issc |
---|
545 | REAL (kind=rp) :: c1,c2 |
---|
546 | TYPE(c2d) :: vq,vs |
---|
547 | !--------------------------------------------------------------------- |
---|
548 | nic = 0; issp = 0; |
---|
549 | vq = c2d(pw%pt(1)%y-pw%pt(2)%y,pw%pt(2)%x-pw%pt(1)%x) |
---|
550 | DO jp=1,pw%n |
---|
551 | vs = pw%pt(MOD(jp,pw%n)+1)-pw%pt(jp) |
---|
552 | CALL p_i2d (pq,vq,pw%pt(jp),vs,c1,c2,ii) |
---|
553 | IF (ii > 0) THEN |
---|
554 | !----- la demi-droite coupe un cote |
---|
555 | IF (c1 > v_eps) THEN |
---|
556 | !------- le point est externe a la droite portant le cote |
---|
557 | IF ( (c2 > v_eps).AND.(c2 < (1.-v_eps)) ) THEN |
---|
558 | !--------- la demi-droite passe par l'interieur du cote |
---|
559 | nic = nic+1 |
---|
560 | issp = 0; |
---|
561 | ELSE IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
562 | !--------- la demi-droite passe par un sommet |
---|
563 | issc = NINT(SIGN(1.,vq.s.vs)); |
---|
564 | IF ((issc*issp) < 0) THEN |
---|
565 | nic = nic+1 |
---|
566 | issp = 0; |
---|
567 | ENDIF |
---|
568 | issp = issc; |
---|
569 | ENDIF |
---|
570 | ELSE IF (ABS(c1) < v_eps) THEN |
---|
571 | !------- le point est sur la droite portant le cote |
---|
572 | IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
573 | !--------- le point appartient au cote |
---|
574 | nic = -1 |
---|
575 | EXIT |
---|
576 | ENDIF |
---|
577 | ENDIF |
---|
578 | ELSE IF (ii < 0) THEN |
---|
579 | !----- la demi-droite est colineaire a un cote |
---|
580 | IF ( (c2 > -v_eps).AND.(c2 < (1.+v_eps)) ) THEN |
---|
581 | nic = -1 |
---|
582 | EXIT |
---|
583 | ENDIF |
---|
584 | ENDIF |
---|
585 | ENDDO |
---|
586 | !- |
---|
587 | p_tap = MOD(nic,2) |
---|
588 | !-------------------- |
---|
589 | END FUNCTION p_tap |
---|
590 | != |
---|
591 | !----------------- |
---|
592 | END MODULE u_polyg |
---|