libsim  Versione 7.1.8

◆ contng_simc()

integer function space_utilities::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.

368 nl0 = 3
369 nlt3 = 9
370 ipl(1) = ip1
371 ipl(2) = ip2
372 ipl(3) = 1
373 ipl(4) = ip2
374 ipl(5) = ip3
375 ipl(6) = 1
376 ipl(7) = ip3
377 ipl(8) = ip1
378 ipl(9) = 1
379 
380 call l4f_log(l4f_debug,"end triangulation first triangle")
381 
382  !
383  ! ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
384  !
385 l400 : 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
417 230 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  !
496 300 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  !
549 330 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)
555 340 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  !
564 350 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
611 end DO l400
612 
613 call 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  !
619 DO 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
626 end DO
627 call l4f_log(l4f_debug,"end triangulation rearranging")
628 
629 nt = nt0
630 nl = nl0
631 
632 call l4f_log(l4f_debug,"end triangulation")
633 
634 RETURN
635 
636 contains
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 
645 elemental double precision function vdsqf(co1,co2)
646 type(xy),intent(in) :: co1,co2
647 
648 vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649 if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
650 end 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 
659 double precision function side(co1,co2,co3)
660 type(xy),intent(in):: co1,co2,co3
661 
662 side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
663 end function side
664 
665 end function contng_simc
666 
667 
673 INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
674 
675 double precision,intent(in) :: X(:), Y(:)
678 integer,intent(in) :: I1,I2,I3,I4
679 
680 double precision :: X0(4), Y0(4)
681 double precision :: C2SQ,C1SQ,A3SQ,B2SQ,B3SQ,A1SQ,A4SQ,B1SQ,B4SQ,A2SQ,C4SQ,C3SQ
682 integer :: IDX
683 double precision :: S1SQ,S2SQ,S3SQ,S4SQ,U1,U2,U3,U4
684 
685 equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686  (b4sq,a2sq),(c4sq,c3sq)
687 
688 x0(1) = x(i1)
689 y0(1) = y(i1)
690 x0(2) = x(i2)
691 y0(2) = y(i2)
692 x0(3) = x(i3)
693 y0(3) = y(i3)
694 x0(4) = x(i4)
695 y0(4) = y(i4)
696 idx = 0
697 u3 = (y0(2)-y0(3))*(x0(1)-x0(3))-(x0(2)-x0(3))*(y0(1)-y0(3))
698 u4 = (y0(1)-y0(4))*(x0(2)-x0(4))-(x0(1)-x0(4))*(y0(2)-y0(4))
699 IF (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
713 end IF
714 conxch_simc = idx
715 RETURN
716 
717 END FUNCTION conxch_simc
718 
719 end module space_utilities
Space utilities, derived from NCAR software.

Generated with Doxygen.