libsim Versione 7.1.11

◆ contng_simc()

integer function contng_simc ( type (xy), dimension(:), intent(in)  co,
integer, intent(out)  nt,
integer, dimension(:), intent(out)  ipt,
integer, intent(out)  nl,
integer, dimension(:), intent(out)  ipl 
)
private

THIS SUBROUTINE PERFORMS TRIANGULATION.

IT DIVIDES THE X-Y PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS CORRESPONDING TO THE BORDER LINE SEGMENTS. AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE ARE LISTED COUNTER-CLOCKWISE. POINT NUMBERS OF THE END POINTS OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE, LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.

Return 0 if all right return 1 if IDENTICAL INPUT DATA POINTS FOUND return 2 if ALL DATA ARE COLLINEAR DATA POINTS

Parametri
[in]coARRAY OF DIMENSION NDP CONTAINING THE COORDINATES OF THE DATA POINTS
[out]ntNUMBER OF TRIANGLES
[out]nlNUMBER OF BORDER LINE SEGMENTS
[out]iptARRAY OF DIMENSION 6*NDP-15, WHERE THE POINT NUMBERS OF THE VERTEXES OF THE (IT)TH TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND, (3*IT-1)ST, AND (3*IT)TH ELEMENTS, IT=1,2,...,NT
[out]iplARRAY OF DIMENSION 6*NDP, WHERE THE POINT NUMBERS OF THE END POINTS OF THE (IL)TH BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE NUMBER ARE TO BE STORED AS THE (3*IL-2)ND, (3*IL-1)ST, AND (3*IL)TH ELEMENTS, IL=1,2,..., NL.

Definizione alla linea 367 del file space_utilities.F90.

368nl0 = 3
369nlt3 = 9
370ipl(1) = ip1
371ipl(2) = ip2
372ipl(3) = 1
373ipl(4) = ip2
374ipl(5) = ip3
375ipl(6) = 1
376ipl(7) = ip3
377ipl(8) = ip1
378ipl(9) = 1
379
380call l4f_log(l4f_debug,"end triangulation first triangle")
381
382 !
383 ! ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
384 !
385l400 : DO jp1=4,ndp
386 ip1 = iwp(jp1)
387 x1 = co(ip1)%x
388 y1 = co(ip1)%y
389 !
390 ! - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
391 !
392 ip2 = ipl(1)
393 jpmn = 1
394 dxmn = co(ip2)%x-x1
395 dymn = co(ip2)%y-y1
396 dsqmn = dxmn**2+dymn**2
397 armn = dsqmn*ratio
398 jpmx = 1
399 dxmx = dxmn
400 dymx = dymn
401 dsqmx = dsqmn
402 armx = armn
403 DO jp2=2,nl0
404 ip2 = ipl(3*jp2-2)
405 dx = co(ip2)%x-x1
406 dy = co(ip2)%y-y1
407 ar = dy*dxmn-dx*dymn
408 IF (ar <= armn) then
409 dsqi = dx**2+dy**2
410 IF (ar.GE.(-armn) .AND. dsqi.GE.dsqmn) GO TO 230
411 jpmn = jp2
412 dxmn = dx
413 dymn = dy
414 dsqmn = dsqi
415 armn = dsqmn*ratio
416 end IF
417230 ar = dy*dxmx-dx*dymx
418 IF (ar .LT. (-armx)) cycle
419 dsqi = dx**2+dy**2
420 IF (ar.LE.armx .AND. dsqi.GE.dsqmx) cycle
421 jpmx = jp2
422 dxmx = dx
423 dymx = dy
424 dsqmx = dsqi
425 armx = dsqmx*ratio
426 end DO
427 IF (jpmx .LT. jpmn) jpmx = jpmx+nl0
428 nsh = jpmn-1
429 IF (nsh > 0) then
430 !
431 ! - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
432 ! - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
433 !
434 nsht3 = nsh*3
435 DO jp2t3=3,nsht3,3
436 jp3t3 = jp2t3+nlt3
437 ipl(jp3t3-2) = ipl(jp2t3-2)
438 ipl(jp3t3-1) = ipl(jp2t3-1)
439 ipl(jp3t3) = ipl(jp2t3)
440 end DO
441 DO jp2t3=3,nlt3,3
442 jp3t3 = jp2t3+nsht3
443 ipl(jp2t3-2) = ipl(jp3t3-2)
444 ipl(jp2t3-1) = ipl(jp3t3-1)
445 ipl(jp2t3) = ipl(jp3t3)
446 end DO
447 jpmx = jpmx-nsh
448 !
449 ! - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
450 ! - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
451 ! - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
452 !
453 end IF
454
455 jwl = 0
456 l310 : DO jp2=jpmx,nl0
457 jp2t3 = jp2*3
458 ipl1 = ipl(jp2t3-2)
459 ipl2 = ipl(jp2t3-1)
460 it = ipl(jp2t3)
461 !
462 ! - - ADDS A TRIANGLE TO THE IPT ARRAY.
463 !
464 nt0 = nt0+1
465 ntt3 = ntt3+3
466 ipt(ntt3-2) = ipl2
467 ipt(ntt3-1) = ipl1
468 ipt(ntt3) = ip1
469 !
470 ! - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
471 !
472 IF (jp2 == jpmx) then
473 ipl(jp2t3-1) = ip1
474 ipl(jp2t3) = nt0
475 end IF
476 IF (jp2 == nl0) then
477 nln = jpmx+1
478 nlnt3 = nln*3
479 ipl(nlnt3-2) = ip1
480 ipl(nlnt3-1) = ipl(1)
481 ipl(nlnt3) = nt0
482 end IF
483 !
484 ! - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
485 ! - - LINE SEGMENTS.
486 !
487 itt3 = it*3
488 ipti = ipt(itt3-2)
489 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
490 ipti = ipt(itt3-1)
491 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
492 ipti = ipt(itt3)
493 !
494 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
495 !
496300 IF (conxch_simc(co%X,co%Y,ip1,ipti,ipl1,ipl2) .EQ. 0) cycle l310
497 !
498 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
499 !
500 ipt(itt3-2) = ipti
501 ipt(itt3-1) = ipl1
502 ipt(itt3) = ip1
503 ipt(ntt3-1) = ipti
504 IF (jp2 .EQ. jpmx) ipl(jp2t3) = it
505 IF (jp2.EQ.nl0 .AND. ipl(3).EQ.it) ipl(3) = nt0
506 !
507 ! - - SETS FLAGS IN THE IWL ARRAY.
508 !
509 jwl = jwl+4
510 iwl(jwl-3) = ipl1
511 iwl(jwl-2) = ipti
512 iwl(jwl-1) = ipti
513 iwl(jwl) = ipl2
514 end DO l310
515 nl0 = nln
516 nlt3 = nlnt3
517 nlf = jwl/2
518 IF (nlf .EQ. 0) cycle l400
519 !
520 ! - IMPROVES TRIANGULATION.
521 !
522 ntt3p3 = ntt3+3
523 DO irep=1,nrep
524 l370 : DO ilf=1,nlf
525 ilft2 = ilf*2
526 ipl1 = iwl(ilft2-1)
527 ipl2 = iwl(ilft2)
528 !
529 ! - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
530 ! - - THE FLAGGED LINE SEGMENT.
531 !
532 ntf = 0
533 DO itt3r=3,ntt3,3
534 itt3 = ntt3p3-itt3r
535 ipt1 = ipt(itt3-2)
536 ipt2 = ipt(itt3-1)
537 ipt3 = ipt(itt3)
538 IF (ipl1.NE.ipt1 .AND. ipl1.NE.ipt2 .AND. ipl1.NE.ipt3) cycle
539 IF (ipl2.NE.ipt1 .AND. ipl2.NE.ipt2 .AND. ipl2.NE.ipt3) cycle
540 ntf = ntf+1
541 itf(ntf) = itt3/3
542 IF (ntf .EQ. 2) GO TO 330
543 end DO
544 IF (ntf .LT. 2) cycle l370
545 !
546 ! - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
547 ! - - ON THE LINE SEGMENT.
548 !
549330 it1t3 = itf(1)*3
550 ipti1 = ipt(it1t3-2)
551 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
552 ipti1 = ipt(it1t3-1)
553 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
554 ipti1 = ipt(it1t3)
555340 it2t3 = itf(2)*3
556 ipti2 = ipt(it2t3-2)
557 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
558 ipti2 = ipt(it2t3-1)
559 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
560 ipti2 = ipt(it2t3)
561 !
562 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
563 !
564350 IF (conxch_simc(co%X,co%Y,ipti1,ipti2,ipl1,ipl2) .EQ. 0) cycle l370
565 !
566 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
567 !
568 ipt(it1t3-2) = ipti1
569 ipt(it1t3-1) = ipti2
570 ipt(it1t3) = ipl1
571 ipt(it2t3-2) = ipti2
572 ipt(it2t3-1) = ipti1
573 ipt(it2t3) = ipl2
574 !
575 ! - - SETS NEW FLAGS.
576 !
577 jwl = jwl+8
578 iwl(jwl-7) = ipl1
579 iwl(jwl-6) = ipti1
580 iwl(jwl-5) = ipti1
581 iwl(jwl-4) = ipl2
582 iwl(jwl-3) = ipl2
583 iwl(jwl-2) = ipti2
584 iwl(jwl-1) = ipti2
585 iwl(jwl) = ipl1
586 DO jlt3=3,nlt3,3
587 iplj1 = ipl(jlt3-2)
588 iplj2 = ipl(jlt3-1)
589 IF ((iplj1.EQ.ipl1 .AND. iplj2.EQ.ipti2) .OR. &
590 (iplj2.EQ.ipl1 .AND. iplj1.EQ.ipti2)) ipl(jlt3) = itf(1)
591 IF ((iplj1.EQ.ipl2 .AND. iplj2.EQ.ipti1) .OR. &
592 (iplj2.EQ.ipl2 .AND. iplj1.EQ.ipti1)) ipl(jlt3) = itf(2)
593 end DO
594 end DO l370
595 nlfc = nlf
596 nlf = jwl/2
597 IF (nlf .EQ. nlfc) cycle l400
598 !
599 ! - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
600 !
601 jwl = 0
602 jwl1mn = (nlfc+1)*2
603 nlft2 = nlf*2
604 DO jwl1=jwl1mn,nlft2,2
605 jwl = jwl+2
606 iwl(jwl-1) = iwl(jwl1-1)
607 iwl(jwl) = iwl(jwl1)
608 end DO
609 nlf = jwl/2
610 end DO
611end DO l400
612
613call l4f_log(l4f_debug,"end triangulation appending")
614
615 !
616 ! REARRANGE THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
617 ! ARE LISTED COUNTER-CLOCKWISE.
618 !
619DO itt3=3,ntt3,3
620 ip1 = ipt(itt3-2)
621 ip2 = ipt(itt3-1)
622 ip3 = ipt(itt3)
623 IF (side(co(ip1),co(ip2),co(ip3)) .GE. 10.0) cycle
624 ipt(itt3-2) = ip2
625 ipt(itt3-1) = ip1
626end DO
627call l4f_log(l4f_debug,"end triangulation rearranging")
628
629nt = nt0
630nl = nl0
631
632call l4f_log(l4f_debug,"end triangulation")
633
634RETURN
635
636contains
637
638!!$double precision function DSQF(U1,V1,U2,V2)
639!!$double precision,intent(in) :: U1,V1,U2,V2
640!!$
641!!$DSQF = (U2-U1)**2+(V2-V1)**2
642!!$end function DSQF
643
644
645elemental double precision function vdsqf(co1,co2)
646type(xy),intent(in) :: co1,co2
647
648vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
650end function vdsqf
651
652
653!!$double precision function SIDE(U1,V1,U2,V2,U3,V3)
654!!$double precision,intent(in):: U1,V1,U2,V2,U3,V3
655!!$
656!!$SIDE = (V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
657!!$end function SIDE
658
659double precision function side(co1,co2,co3)
660type(xy),intent(in):: co1,co2,co3
661
662side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
663end function side
664
665end function contng_simc
666
667
673INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
674
675double precision,intent(in) :: X(:), Y(:)
678integer,intent(in) :: I1,I2,I3,I4
679
680double precision :: X0(4), Y0(4)
681double precision :: C2SQ,C1SQ,A3SQ,B2SQ,B3SQ,A1SQ,A4SQ,B1SQ,B4SQ,A2SQ,C4SQ,C3SQ
682integer :: IDX
683double precision :: S1SQ,S2SQ,S3SQ,S4SQ,U1,U2,U3,U4
684
685equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686 (b4sq,a2sq),(c4sq,c3sq)
687
688x0(1) = x(i1)
689y0(1) = y(i1)
690x0(2) = x(i2)
691y0(2) = y(i2)
692x0(3) = x(i3)
693y0(3) = y(i3)
694x0(4) = x(i4)
695y0(4) = y(i4)
696idx = 0
697u3 = (y0(2)-y0(3))*(x0(1)-x0(3))-(x0(2)-x0(3))*(y0(1)-y0(3))
698u4 = (y0(1)-y0(4))*(x0(2)-x0(4))-(x0(1)-x0(4))*(y0(2)-y0(4))
699IF (u3*u4 > 0.0) then
700 u1 = (y0(3)-y0(1))*(x0(4)-x0(1))-(x0(3)-x0(1))*(y0(4)-y0(1))
701 u2 = (y0(4)-y0(2))*(x0(3)-x0(2))-(x0(4)-x0(2))*(y0(3)-y0(2))
702 a1sq = (x0(1)-x0(3))**2+(y0(1)-y0(3))**2
703 b1sq = (x0(4)-x0(1))**2+(y0(4)-y0(1))**2
704 c1sq = (x0(3)-x0(4))**2+(y0(3)-y0(4))**2
705 a2sq = (x0(2)-x0(4))**2+(y0(2)-y0(4))**2
706 b2sq = (x0(3)-x0(2))**2+(y0(3)-y0(2))**2
707 c3sq = (x0(2)-x0(1))**2+(y0(2)-y0(1))**2
708 s1sq = u1*u1/(c1sq*max(a1sq,b1sq))
709 s2sq = u2*u2/(c2sq*max(a2sq,b2sq))
710 s3sq = u3*u3/(c3sq*max(a3sq,b3sq))
711 s4sq = u4*u4/(c4sq*max(a4sq,b4sq))
712 IF (min(s1sq,s2sq) < min(s3sq,s4sq)) idx = 1
713end IF
714conxch_simc = idx
715RETURN
716
717END FUNCTION conxch_simc
718
719end module space_utilities
Space utilities, derived from NCAR software.

Generated with Doxygen.