source: trunk/SRC/ToBeReviewed/TRIANGULATION/definetri.pro

;+
;
; @file_comments
; Define a triangulation array like <proidl>TRIANGULATE</proidl>.
;
; But in a VERY SIMPLE CASE:
; the points are regularly-gridded on nx*ny array.
; Find a Delaunay triangulation for this set of points is easy:
; Points define (nx-1)*(ny-1) rectangles which we can cut in 2 triangles.
;
; cf. figure above
;
; <fixe>
;      ny-1*---*---*. . . . . .*---*---*
;          |  +|  +|           |  +|  +|
;          | + | + |           | + | + |
;          |+  |+  |           |+  |+  |
;      ny-2*---*---*. . . . . .*---*---*
;          .       .           .       .
;          .       .           .       .
;          .       .           .       .
;         1*---*---*. . . . . .*---*---*
;          |  +|  +|           |  +|  +|
;          | + | + |           | + | + |
;          |+  |+  |           |+  |+  |
;         0*---*---*. . . . . .*---*---*
;           0   1   2        nx-3  nx-2 nx-1
; </fixe>
;
;  You have 2 ways to cut a rectangle:
;      1) the upward diagonal       2) the downward diagonal
;
; <fixe>
;          *---*                        *---*
;          |  +|                        |+  |
;          | + |                        | + |
;          |+  |                        |  +|
;          *---*                        *---*
; </fixe>
;
; @categories
; Utilities
;
; @param nx {in}{required}
; The x dimension array
;
; @param ny {in}{required}
; The y dimension array
;
; @param downward {in}{optional}
; When downward is undefined all rectangles are cut in using the upward
; diagonal.
; downward is a vector which contains the rectangles numbers which are cut in
; using the downward diagonal.
; The rectangle number is defined by the index (in a nx*ny vector) of the
; lower-left corner of the rectangle.
;
; @returns
; triangles is a 2d array and its dimensions are 3 and 2*(nx-1)*(ny-1).
; triangles is defined like in the TRIANGULATE procedure.
;
; @examples
;
;   IDL> triangles=definetri(3,3,[1,3])
;
; triangles will be this kind of triangulation:
;
; <fixe>
;          *---*---*
;          |+  |  +|
;          | + | + |
;          |  +|+  |
;          *---*---*
;          |  +|+  |
;          | + | + |
;          |+  |  +|
;          *---*---*
;
; <fixe>
;
; @history
; sebastien Masson (smlod\@ipsl.jussieu.fr)
;                       4/3/1999
;
; @version
; $Id$
;-
FUNCTION definetri, nx, ny, downward
;
  compile_opt idl2, strictarrsubs
;
   nx = long(nx)
   ny = long(ny)
   if n_elements(downward) NE 0 THEN BEGIN
      if n_elements(downward) GT (nx-1)*(ny-1) then begin
         print, 'downward a trop d''elements par rapport a nx et ny!'
         return,  -1
      endif
      downward = long(downward)
   ENDIF
; we define triangles
   triangles = lonarr(3, 2*(nx-1)*(ny-1))
;----------------------------------------------------------------------------------
; we cut the rectangles with the upward diagonal
;----------------------------------------------------------------------------------
   if n_elements(downward) NE (nx-1)*(ny-1) then BEGIN ; there is some rectangle to cut.
; we define upward: upward is a vector which contains the rectangles
; numbers which are cut in using the upward diagonal.
; The rectangle number is defined by the index (in a nx*ny vector) of
; the lower-left corner of the rectangle.
      upward = bytarr(nx, ny)+1
      upward[*, ny-1] = 0
      upward[nx-1, *] = 0
      if n_elements(downward) NE 0 then upward[downward] = 0
      upward = where(upward EQ 1)
      n1 = n_elements(upward)
;
; 4 corners indexes of a rectangle number i are
;
;       i+nx  i+nx+1
;          *---*
;          |  +|
;          | + |
;          |+  |
;          *---*
;          i   i+1
;
      trinumber = 2*(upward-upward/nx)
; we define the right triangles
      triangles[0, trinumber] = upward
      triangles[1, trinumber] = upward+1
      triangles[2, trinumber] = upward+1+nx
; we define the left triangles
      triangles[0, trinumber+1] = upward+1+nx
      triangles[1, trinumber+1] = upward+nx
      triangles[2, trinumber+1] = upward
   ENDIF ELSE n1 = 0
;----------------------------------------------------------------------------------
; we cut the rectangles with the downward diagonal
;----------------------------------------------------------------------------------
   if n_elements(downward) NE 0 then BEGIN
      n2 = n_elements(downward)
      trinumber = 2*(downward-downward/nx)
; we define the right triangles
      triangles[0, trinumber] = downward+1
      triangles[1, trinumber] = downward+nx+1
      triangles[2, trinumber] = downward+nx
; we define the left triangles
      triangles[0, trinumber+1] = downward+nx
      triangles[1, trinumber+1] = downward
      triangles[2, trinumber+1] = downward+1
   endif

   return, triangles
end