| 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 |
|---|
