(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)))