source: trunk/ToBeReviewed/TRIANGULATION/definetri_e.pro @ 29

function numtri, index, nx, ny

   y=index/nx
   x=index-y*nx
   numtri = (y NE 0)*(nx-1)*(2*(y-1)+1) + (2-(y EQ (ny-1) OR y EQ (ny-1)))*x

   return, numtri
end


;+
; NAME:definetri
;
; PURPOSE:Define a triangulation array like TRIANGULATE but for a
; E-grid type
;
; CATEGORY: make contours with E-grid type
;
; CALLING SEQUENCE:triangles=definetri(nx, ny [,vertical])
; 
; INPUTS: nx and ny are the array dimensions 
;
; OPTIONAL INPUTS:
;
;        vertical: When vertical is undefine all rectangles are cut
;        in using the horizontal diagonal. Vertical is a vector which
;        contains the rectangles numbers which are cut in using the
;        vertical diagonal.
;        The rectangle number is define by the index (in a nx*ny
;        vector) of the lower-left corner of the rectangle.
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;        triangles is a 2d array and is dimensions are 3 and
;        2*(nx-1)*(ny-1)
;        triangles is define like in the TRIANGULATE procedure.
;        
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; EXAMPLE:
;
;
; MODIFICATION HISTORY: sebastien Masson (smlod@ipsl.jussieu.fr)
;                       June 2001
;
;-
FUNCTION definetri_e, nx, ny, singular, SHIFTED = shifted
   nx = long(nx)
   ny = long(ny)
   triangles = lonarr(3, 2*(nx-1)*(ny-1))
;
; build the base triangulation with the diamond cut in two triangles
; by the vertical diagonal
;
; first line
   index = lindgen(nx-1)
   trinumber = index
   triangles[0, trinumber] = index
   triangles[1, trinumber] = index+1
   triangles[2, trinumber] = index+(nx+1-shifted)
; last line
   index = (ny-1)*nx+lindgen(nx-1)
   trinumber = numtri(index, nx, ny)
   triangles[0, trinumber] = index
   triangles[1, trinumber] = index+(-nx+((index/nx+1-shifted) MOD 2))
   triangles[2, trinumber] = index+1
; other lines
   if ny GT 2 then begin
      index = lindgen(nx, ny)
      index = index[0:nx-2, 1:ny-2]
      index = index[*]
      oddeven = (index/nx+1-shifted) MOD 2
      trinumber = numtri(index, nx, ny)
      triangles[0, trinumber] = index
      triangles[1, trinumber] = index-nx+oddeven
      triangles[2, trinumber] = index+nx+oddeven
      triangles[0, trinumber+1] = index+nx+oddeven
      triangles[1, trinumber+1] = index-nx+oddeven
      triangles[2, trinumber+1] = index+1
   endif
;
; cut the diamond specified by singular in two triangles
; by the horizontal diagonal
;
   IF keyword_set(singular) then BEGIN
      yindex = singular/nx
      otherline = where(yindex NE 0 AND yindex NE (ny-1))
      if otherline[0] NE -1 then begin
         index = singular[otherline]
         oddeven = (index/nx+1-shifted) MOD 2
         trinumber = numtri(index, nx, ny)
         triangles[0, trinumber] = index
         triangles[1, trinumber] = index-nx+oddeven
         triangles[2, trinumber] = index+1
         triangles[0, trinumber+1] = index
         triangles[1, trinumber+1] = index+1
         triangles[2, trinumber+1] = index+nx+oddeven
      endif

   endif
   return, triangles
end