(define coarse-square-triangulation
  
  ;; This defines a triangulation of the unit square with two triangles,
  ;; five edges, and four vertices.
  
  (let ((pll (make-Point 0.0 0.0)) ; lower-left
	(plr (make-Point 1.0 0.0)) ; lower-right
	(pul (make-Point 0.0 1.0)) ; upper-left
	(pur (make-Point 1.0 1.0))); upper-right
    (with-co-instantiation
     ((vll Boundary-vertex                  ; lower-left
	   point: pll
	   edges: (Vertex-edges<-list (list bottom diagonal left)) 
	   boundary-edges: (Boundary-vertex-boundary-edges<-list (list bottom left)) 
	   )
      (vlr Boundary-vertex                  ; lower-right
	   point: plr 
	   edges: (Vertex-edges<-list (list right bottom))         
	   boundary-edges: (Boundary-vertex-boundary-edges<-list (list right bottom))
	   )
      (vul Boundary-vertex                  ; upper-left
	   point: pul 
	   edges: (Vertex-edges<-list (list left top))             
	   boundary-edges: (Boundary-vertex-boundary-edges<-list (list left top))
	   )
      (vur Boundary-vertex                  ; upper-right
	   point: pur 
	   edges: (Vertex-edges<-list (list top diagonal right))   
	   boundary-edges: (Boundary-vertex-boundary-edges<-list (list top right))   
	   )
      (bottom   Boundary-edge 
		vertex-1: vll vertex-2: vlr 
		)
      (top      Boundary-edge 
		vertex-1: vul vertex-2: vur
		)
      (left     Boundary-edge 
		vertex-1: vll vertex-2: vul
		)
      (right    Boundary-edge 
		vertex-1: vlr vertex-2: vur 
		)
      (diagonal Edge 
		vertex-1: vll vertex-2: vur
		)
      (tul Triangle              ; upper-left
	   vertex-1: vll vertex-2: vur vertex-3: vul 
	   edge-1: diagonal edge-2: top edge-3: left)
      (tlr Triangle              ; lower-right
	   vertex-1: vll vertex-2: vlr vertex-3: vur 
	   edge-1: bottom edge-2: right edge-3: diagonal))
     (instantiate Triangulation
		  vertices:  (Triangulation-vertices<-list 
			      (list vll vlr vul vur))
		  edges:     (Triangulation-edges<-list     
			      (list bottom top left right diagonal))
		  triangles: (Triangulation-triangles<-list 
			      (list tul tlr))
		  has-polygonal-boundary?: #t))))

(define-generic (refine (o)))

;; We have temporarily replaced with-co-instantiation with instantiate
;; and explicitly initialized all fields so that fill-other-fields!
;; is not called.  This doubles (!) the speed of (refine triangulation)
;; even with the space-filling-curve ordering of all the vertices
;; edges and triangles.  I hope Christian can change the macro expansion
;; to be smarter so that we can use with-co-instantiation again.

;; Each old vertex introduces a new vertex at the same location.
;; We'll store this new vertex in Vertex-fine of each old vertex.

(define-handy-method (refine (v Vertex))
  (set! fine (instantiate Vertex
			  index: #f
			  fine: '()
			  edges: '()
			  point: point)))

(define-handy-method (refine (v Boundary-vertex))
  (set! fine (instantiate Boundary-vertex
			  index: #f
			  fine: '()
			  edges: '()
			  boundary-index: #f
			  boundary-edges: '()
		  point: point)))  ; in the refine methods only

;; Each old edge has a new vertex in the middle, and 
;; introduces two new edges.  We'll store a list of the new vertex
;; and the two new edges in Edge-fine

(define-handy-method (refine (e Edge))
  (let ((vertex-1 (Vertex-fine vertex-1))
	(vertex-2 (Vertex-fine vertex-2)))
    (let ((point-1 (Vertex-point vertex-1))
	  (point-2 (Vertex-point vertex-2)))
      (let ((middle (instantiate Vertex 
				 index: #f
				 fine: '()
				 edges: '()
				 point: (Point-scale 0.5 
						     (Point-add point-1 
								point-2)))))
	(let ((left (instantiate Edge
				 index: #f
				 fine: '()
				 vertex-1: vertex-1      ; existing vertex is vertex-1
				 vertex-2: middle))
	      (right (instantiate Edge
				  index: #f
				  fine: '()
				  vertex-1: vertex-2     ; existing vertex is vertex-1
				  vertex-2: middle)))
	  (Edge-fine-set! e (list middle left right))
	  
	  ;; add edge information to new vertices
	  
	  (for-each 
	   (lambda (e)
	     (Edge-for-each-vertex 
	      (lambda (v)
		(with-access
		 v (Vertex edges)
		 (set! edges (cons e edges))))
	      e))
	   (cdr (Edge-fine e))))))))

(define-handy-method (refine (e Boundary-edge))
  (let ((vertex-1 (Vertex-fine vertex-1))
	(vertex-2 (Vertex-fine vertex-2)))
    (let ((point-1 (Vertex-point vertex-1))
	  (point-2 (Vertex-point vertex-2)))
      (let ((middle (instantiate Boundary-vertex 
				 index: #f
				 fine: '()
				 edges: '()
				 boundary-index: #f
				 boundary-edges: '()
				 point: (Point-scale 0.5 
						     (Point-add point-1 
								point-2)))))
	(let ((left (instantiate Boundary-edge
				 index: #f
				 fine: '()
				 boundary-index: #f
				 vertex-1: vertex-1      ; existing vertex is vertex-1
				 vertex-2: middle))
	      (right (instantiate Boundary-edge
				  index: #f
				  fine: '()
				  boundary-index: #f
				  vertex-1: vertex-2     ; existing vertex is vertex-1
				  vertex-2: middle)))
	  (Edge-fine-set! e (list middle left right))
	  
	  ;; add edge information to new vertices
	  
	  (for-each 
	   (lambda (e)
	     (Edge-for-each-vertex 
	      (lambda (v)
		(with-access
		 v (Boundary-vertex edges boundary-edges)
		 (set! edges (cons e edges))
		 (set! boundary-edges (cons e boundary-edges))))
	      e))
	   (cdr (Edge-fine e))))))))

;; each old triangle is decomposed into four new triangles and 
;; generates three new edges

(define-handy-method (refine (t Triangle))
  (let ((old-vertex-1 vertex-1)
	(old-vertex-2 vertex-2)
	(old-vertex-3 vertex-3)
	(old-edge-1 edge-1)
	(old-edge-2 edge-2)
	(old-edge-3 edge-3))
    (let ((vertex-1 (Vertex-fine old-vertex-1))
	  (vertex-2 (Vertex-fine old-vertex-2))
	  (vertex-3 (Vertex-fine old-vertex-3))
	  (mid-vertex-1 (car (Edge-fine old-edge-1)))
	  (mid-vertex-2 (car (Edge-fine old-edge-2)))
	  (mid-vertex-3 (car (Edge-fine old-edge-3))))
      (let ((mid1 (instantiate Edge
			       index: #f
			       fine: '()
			       vertex-1: mid-vertex-1
			       vertex-2: mid-vertex-2))
	    (mid2 (instantiate Edge
			       index: #f
			       fine: '()
			       vertex-1: mid-vertex-2
			       vertex-2: mid-vertex-3))
	    (mid3 (instantiate Edge
			       index: #f
			       fine: '()
			       vertex-1: mid-vertex-3
			       vertex-2: mid-vertex-1)))
	(let ((t1 (instantiate Triangle 
			       index: #f
			       fine: '()
			       vertex-1: vertex-1 
			       vertex-2: mid-vertex-1 
			       vertex-3: mid-vertex-3
			       edge-1: (if (eq? (Edge-vertex-1 old-edge-1) old-vertex-1)
					   (cadr (Edge-fine old-edge-1))
					   (caddr (Edge-fine old-edge-1)))
			       edge-2: mid3
			       edge-3: (if (eq? (Edge-vertex-1 old-edge-3) old-vertex-1)
					   (cadr (Edge-fine old-edge-3))
					   (caddr (Edge-fine old-edge-3)))))
	      (t2 (instantiate Triangle
			       index: #f
			       fine: '()
			       vertex-1: vertex-2
			       vertex-2: mid-vertex-2
			       vertex-3: mid-vertex-1
			       edge-1: (if (eq? (Edge-vertex-1 old-edge-2) old-vertex-2)
					   (cadr (Edge-fine old-edge-2))
					   (caddr (Edge-fine old-edge-2)))
			       edge-2: mid1
			       edge-3: (if (eq? (Edge-vertex-1 old-edge-1) old-vertex-2)
					   (cadr (Edge-fine old-edge-1))
					   (caddr (Edge-fine old-edge-1)))))
	      (t3 (instantiate Triangle
			       index: #f
			       fine: '()
			       vertex-1: vertex-3
			       vertex-2: mid-vertex-3
			       vertex-3: mid-vertex-2
			       edge-1: (if (eq? (Edge-vertex-1 old-edge-3) old-vertex-3)
					   (cadr (Edge-fine old-edge-3))
					   (caddr (Edge-fine old-edge-3)))
			       edge-2: mid2
			       edge-3: (if (eq? (Edge-vertex-1 old-edge-2) old-vertex-3)
					   (cadr (Edge-fine old-edge-2))
					   (caddr (Edge-fine old-edge-2)))
			       ))
	      (t4 (instantiate Triangle
			       index: #f
			       fine: '()
			       vertex-1: mid-vertex-1
			       vertex-2: mid-vertex-2
			       vertex-3: mid-vertex-3
			       edge-1: mid1
			       edge-2: mid2
			       edge-3: mid3)))
	  
	  (let ((new-triangles (list t1 t2 t3 t4))
		(new-edges (list mid1 mid2 mid3)))
	    (Triangle-fine-set! t (list new-triangles new-edges))
	    
	    ;; add new edge information to points; none of the new 
	    ;;edges are on the boundary
	    
	    (for-each (lambda (e)
			(Edge-for-each-vertex 
			 (lambda (v)
			   (with-access
			    v (Vertex edges)
			    ;; accumulate in a list
			    (set! edges (cons e edges))))
			 e))
		      new-edges)))))))

(define-method (refine (t Triangulation))
  (let ((old-vertices (Triangulation-vertices->list 
		       (Triangulation-vertices t)))
	(old-edges (Triangulation-edges->list 
		    (Triangulation-edges t)))
	(old-triangles (Triangulation-triangles->list 
			(Triangulation-triangles t))))
    
    (for-each refine old-vertices)
    
    (for-each refine old-edges)
    
    (for-each refine old-triangles)
    
    (let ((new-Triangulation
	   (instantiate Triangulation
			triangles: (Triangulation-triangles<-list 
				    (mappend1 (lambda (t) 
					       (car (Triangle-fine t))) 
					     old-triangles))
			edges: (Triangulation-edges<-list 
				(append (mappend1 (lambda (e) 
						   (cdr (Edge-fine e))) 
						 old-edges)
					(mappend1 (lambda (t) 
						   (cadr (Triangle-fine t))) 
						 old-triangles)))
			vertices: (Triangulation-vertices<-list 
				   (append (map Vertex-fine 
						old-vertices)
					   (map (lambda (e) 
						  (car (Edge-fine e))) 
						old-edges)))
			has-polygonal-boundary?: (Triangulation-has-polygonal-boundary? t)
			;; if the triangulation has a polygonal boundary,
			;; then the refined triangulation has the same limits
			;; as the original triangulation
			limits: (and (Triangulation-has-polygonal-boundary? t)
				     (Triangulation-limits t)))))
      
      ;; convert temporary lists to internal data structures
      
      (Triangulation-for-each-vertex 
       (lambda (v)
	 (with-access
	  v (Vertex edges)
	  (set! edges (Vertex-edges<-list edges))))
       new-Triangulation)
      (Triangulation-for-each-boundary-vertex
       (lambda (v)
	 (with-access
	  v (Boundary-vertex boundary-edges)
	  (set! boundary-edges (Boundary-vertex-boundary-edges<-list boundary-edges))))
       new-Triangulation)
      
      ;; convert fine fields on coarse triangulation
      ;; to what we need outside refine.

      ;; Vertex-fine is OK.

      (Triangulation-for-each-edge
       (lambda (e)
	 (with-access
	  e (Edge fine)
	  (set! fine (car fine))))   ; keep only the new vertex
       t)
      
      (Triangulation-for-each-triangle
       (lambda (t)
	 (Triangle-fine-set! t '()))   ; I don't believe we need anything here
       t)
	  
      new-Triangulation)))

(define (Triangulation-consistent? t)
  (Triangulation-for-each-triangle 
   (lambda (t)
     (with-access
      t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
      (for-each 
       (lambda (edge vertex-1 vertex-2)
	 (if (not (and (Edge-contains-vertex? edge vertex-1)
		       (Edge-contains-vertex? edge vertex-2)))
	     (for-each 
	      display 
	      (list "edge-k is not between vertex-k and vertex-k+1: ("
		    (Point-x (Vertex-point (Edge-vertex-1 edge))) ","
		    (Point-y (Vertex-point (Edge-vertex-1 edge))) ") ("
		    (Point-x (Vertex-point (Edge-vertex-2 edge))) ","
		    (Point-y (Vertex-point (Edge-vertex-2 edge))) ")\n"))))
       (list edge-1 edge-2 edge-3)
       (list vertex-1 vertex-2 vertex-3) 
       (list vertex-2 vertex-3 vertex-1))))
   t)
  (Triangulation-for-each-edge 
   (lambda (e)
     (Edge-for-each-vertex 
      (lambda (v)
	(if (not (Vertex-contains-edge? v e))
	    (for-each display (list "edge is not in vertex-edges: ("
				    (Point-x (Vertex-point v)) ","
				    (Point-y (Vertex-point v)) ")\n"))))
      e))
   t)
  (Triangulation-for-each-boundary-edge
   (lambda (e)
     (Edge-for-each-vertex 
      (lambda (v)
	(if (not (Boundary-vertex? v))
	    (for-each display (list "Boundary edge has a non-boundary vertex: "
				    (Point-x (Vertex-point v)) ","
				    (Point-y (Vertex-point v)) ")\n"))
	    (if (not (Boundary-vertex-contains-boundary-edge? v e))
		(for-each display (list "Boundary edge is not in Vertex-boundary-edges: "
					(Point-x (Vertex-point v)) ","
					(Point-y (Vertex-point v)) ")\n")))))
      e))
   t)
  (Triangulation-for-each-vertex 
   (lambda (v)
     (Vertex-for-each-edge 
      (lambda (e)
	(if (not (Edge-contains-vertex? e v))
	    (for-each 
	     display 
	     (list 
	      "edge is in vertex-edges, but does not contain this vertex: ("
	      (Point-x (Vertex-point v)) ","
	      (Point-y (Vertex-point v)) ")\n"))))
      v))
   t)
  (Triangulation-for-each-boundary-vertex
   (lambda (v)
     (Boundary-vertex-for-each-boundary-edge
      (lambda (e)
	(if (not (Boundary-edge? e))
	    (for-each display (list "Boundary vertex has a non-boundary edge in boundary-edges: "
				    (Point-x (Vertex-point v)) ","
				    (Point-y (Vertex-point v)) ")\n"))))
      v))
   t))

(define (lists->Triangulation coordinates elements)

  ;; we assume that coordinates is a
  ;; list of
  ;;   lists of two inexact reals
  ;; and the i'th vertex will have a point with coordinates equal to the i'th
  ;; element in this list (starting at 0)

  ;; we assume that elements is a
  ;; list of
  ;;   lists of three indices of vertices in the triangle making up the element

  (define (coordinates->vertices coordinates)

    ;; builds a vector of vertices from the list of coordinates
    ;; doesn't try to figure out which vertices are on the boundary,
    ;; elements->edges will change vertices on the boundary to boundary-vertices
    ;; (now I am beginning to doubt the utility of the concept of boundary vertices ...)

    (let ((result
	   (list->vector
	    (map (lambda (point-coordinates)
		   (instantiate Vertex
				point: (apply make-Point point-coordinates)))
		 coordinates))))
      (vector-for-each-with-index (lambda (i v) (Vertex-index-set! v i)) result)
      result))

  (define (elements->edges! elements vertices)   ;; convertes some of the vertices to boundary-vertices
    
    (define (Vertex->boundary-vertex! i)
      (let ((v (vector-ref vertices i)))
	(if (not (Boundary-vertex? v))
	    (let ((new-v (instantiate Boundary-vertex
				      point: (Vertex-point v)
				      index: (Vertex-index v)
				      )))
	      (vector-set! vertices i new-v)))))
    
    (let* ((redundant-edge-vertex-indices
	    (sort! (let loop ((elements elements))
		     
		     ;; calculate the indices of the vertices of all the edges in the triangle
		     ;; we assume that the indices in each element are sorted in
		     ;; increasing order, so the indices in each edge are also
		     ;; sorted in increasing order
		     
		     (if (null? elements)
			 '()
			 (let ((element (car elements)))
			   `(,(list (car element) (cadr element))
			     ,(list (car element) (caddr element))
			     ,(list (cadr element) (caddr element))
			     ,@(loop (cdr elements))))))
		   
		   ;; we sort all the edges in lexicographical order of the indices
		   ;; of their vertices
		   
		   (lambda (l1 l2)
		     (or (fx< (car l1) (car l2))
			 (and (fx= (car l1) (car l2))
			      (fx< (cadr l1) (cadr l2)))))))
	   (edge-vertex-indices+counts
	    (unique-with-count equal? redundant-edge-vertex-indices))
	   (ignore
	    (for-each (lambda (count+vertex-indices)
			(let ((count (car count+vertex-indices))
			      (vertex-indices (cdr count+vertex-indices)))
			  (if (fx= count 1)
			      ;; this edge is part of only one triangle, so it is a boundary edge
			      ;; so we convert to boundary vertices the vertices of the edge
			      (for-each Vertex->boundary-vertex! vertex-indices))))
		      edge-vertex-indices+counts))
	   
	   (edges
	    (list->vector
	     (map (lambda (count+vertex-indices)
		    (let ((count (car count+vertex-indices))
			  (vertex-indices (cdr count+vertex-indices)))
		      (case count
			((1)
			 ;; this edge is part of only one triangle, so it
			 ;; is a boundary edge
			 (apply (lambda (index-1 index-2)
				  (instantiate Boundary-edge
					       vertex-1: (vector-ref vertices index-1)
					       vertex-2: (vector-ref vertices index-2)))
				vertex-indices))
			((2)
			 (apply (lambda (index-1 index-2)
				  (instantiate Edge
					       vertex-1: (vector-ref vertices index-1)
					       vertex-2: (vector-ref vertices index-2)))
				vertex-indices))
			(else
			 (error "should not happen")))))
		  edge-vertex-indices+counts)))
	   )
      
      ;; add edges to vertices
      
      (vector-for-each (lambda (e)
			 (Edge-for-each-vertex (lambda (v)
						 (Vertex-edges-set! v (cons e (Vertex-edges v))))
					       e))
		       edges)

      ;; turn them into proper Vertex-edges

      (vector-for-each (lambda (v)
			 (Vertex-edges-set! v (Vertex-edges<-list (Vertex-edges v))))
		       vertices)
      
      ;; add boundary edges to boundary vertices
      
      (vector-for-each (lambda (v)
			 (if (Boundary-vertex? v)
			     (begin
			       (Vertex-for-each-edge (lambda (e)
						       (if (Boundary-edge? e)
							   (Boundary-vertex-boundary-edges-set! v (cons e (Boundary-vertex-boundary-edges v)))))
						     v)
			       (Boundary-vertex-boundary-edges-set! v (Boundary-vertex-boundary-edges<-list (Boundary-vertex-boundary-edges v))))))
		       vertices)
      
      edges))

  (define (elements->triangles elements vertices edges)

    (define (find-edge v1 v2)
      (let loop ((edges (Vertex-edges->list (Vertex-edges v1))))
	(if (null? edges)
	    (error "elements->triangles!: ran out of edges: " v1 v2)
	    (let ((edge (car edges)))
	      (if (Edge-contains-vertex? edge v2)
		  edge
		  (loop (cdr edges)))))))
    
    (let ((result
	   (list->vector
	    (map
	     (lambda (element)
	       (let ((element-vertices (map (lambda (i) (vector-ref vertices i)) element)))
		 (let ((element-vertices
			(if (apply Point-counter-clockwise? (map Vertex-point element-vertices))
			    element-vertices
			    (reverse element-vertices))))
		   (let ((v1 (car element-vertices))
			 (v2 (cadr element-vertices))
			 (v3 (caddr element-vertices)))
		     (let ((e1 (find-edge v1 v2))
			   (e2 (find-edge v2 v3))
			   (e3 (find-edge v3 v1)))
		       (instantiate Triangle
				    vertex-1: v1
				    vertex-2: v2
				    vertex-3: v3
				    edge-1: e1
				    edge-2: e2
				    edge-3: e3))))))
	     elements))))

      ;; I don't see that any more processing is needed
      
      result))
  
  (let* ((elements  (map (lambda (element) (sort element fx<)) elements))
	 (vertices  (coordinates->vertices coordinates))
	 (edges     (elements->edges! elements vertices))
	 (triangles (elements->triangles elements vertices edges)))
    (instantiate Triangulation
		 vertices: vertices
		 edges: edges
		 triangles: triangles
		 has-polygonal-boundary?: #t)))