;; The following routine builds a Linear-finite-element-space (i.e., a Finite-element-space
;; consisting of linear elements) from a triangulation.  Each vertex in the triangulation
;; is associated with a basis element, and the neighbors of the element are the elements
;; associated with the vertices connected to the original vertex by (Vertex-edges v).
;; We use the same numbering scheme for the vertices as for the elements, and rely on that
;; a lot in the following code.

(define (Triangulation->Linear-finite-element-space t)
  (let ((element-list 
	 (map (lambda (v)
		(instantiate Linear-element
			     index: (Vertex-index v)
			     neighbors: (Neighbors<-list
					 (cons (Vertex-index v)  ; v is a neighbor to itelf
					       (map (lambda (e)
						      (Vertex-index
						       (if (eq? (Edge-vertex-1 e) v)
							   (Edge-vertex-2 e)
							   (Edge-vertex-1 e))))
						    (Vertex-edges->list (Vertex-edges v)))))))
	      (Triangulation-vertices->list (Triangulation-vertices t)))))
    (instantiate Linear-finite-element-space
		 triangulation: t
		 neighbors: (Finite-element-space-neighbors<-list (map Element-neighbors element-list))
		 elements: (Finite-element-space-elements<-list element-list))))

;; We anticipate that Finite-element-operators will have matrix implementations
;; of one type or another (Sparse-matrix-operators, most likely).

(define-class Finite-element-operator Operator-with-matrix-implementation
  ((= offsets immutable:)                           ; The offset into indices and coefficients for each row
   (= indices immutable:)                           ; a vector of the indices of ...
   (= coefficients maybe-uninitialized: immutable:) ; nonzero-coefficients in the matrix representation
   (= domain-space immutable:)                      ; (pointers to) the domain ... 
   (= range-space immutable:)))                     ; and range finite-element-spaces

(define-class Linear-finite-element-operator Finite-element-operator ())


(define-handy-method (initialize! (o Linear-finite-element-operator))
  
  ;; initialize various fields that depend on the space
  
  ;; if coefficients is not defined, we set it to arrays of floating-point
  ;; zeros associated with the associated indices field (which must be defined).
  
  (if (not (field-defined? o 'coefficients))
      (initialize-field-value! 
       o
       (make-f64vector (indices-length indices) 0.0)
       'coefficients))

  ;; these are easy, so we do them here
  
  (if (not (field-defined? o 'domain?))
      (initialize-field-value!
       o
       (lambda (v)
	 (and (Finite-element-vector? v)
	      (eq? domain-space (Finite-element-vector-space v))))
       'domain?))
  
  (if (not (field-defined? o 'range?))
      (initialize-field-value!
       o
       (lambda (v)
	 (and (Finite-element-vector? v)
	       (eq? range-space (Finite-element-vector-space v))))
       'range?))
  
  ;; Here both indices and coefficients are initialized, so we can rely
  ;; on them since they're immutable.
  
  (if (not (field-defined? o 'matrix-implementation))
      (initialize-field-value!
       o
       (instantiate Sparse-matrix-operator
		    domain-dimension: (Finite-element-space-neighbors-length domain-space)
		    range-dimension: (Finite-element-space-neighbors-length range-space)
		    indices: indices
		    coefficients: coefficients
		    offsets: offsets)
       'matrix-implementation))
  
  ;; initialize the mapping field to multiplication by the sparse matrix
  
  (initialize-field-value!
   o
   (lambda (v)
     (instantiate Linear-finite-element-vector
		  data: (Operator-apply matrix-implementation (Finite-element-vector-data v))
		  space: range-space))
   'mapping)
  
  ;; call superclass's initialize! method
  
  (call-next-method))

(define-handy-method (->Linear-finite-element-operator (space Linear-finite-element-space))
  (let ((neighbors neighbors))
    (let ((indices (list->indices (mappend1 (lambda (x)
					      (vector->list x))
					    (vector->list neighbors))))
	  (offsets (let* ((n (fx+ (vector-length neighbors) 1))
			  (result (make-vector n)))
		     (vector-set! result 0 0)
		     (do ((i 1 (fx+ i 1)))
			 ((fx= i n) result)
		       (vector-set! result i (fx+ (vector-ref result (fx- i 1))
						  (vector-length (vector-ref neighbors (fx- i 1)))))))))
      (instantiate Linear-finite-element-operator
		   indices: indices
		   offsets: offsets
		   domain-space: space
		   range-space: space))))

(define-handy-method (check-same-class (m Linear-finite-element-operator) n)
  (and (call-next-method)
       (eq? domain-space
	    (Linear-finite-element-operator-domain-space n))
       (eq? range-space
	    (Linear-finite-element-operator-range-space n))))

(define (wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space)
  (instantiate Linear-finite-element-operator
	       indices: (Sparse-matrix-operator-indices matrix)
	       coefficients: (Sparse-matrix-operator-coefficients matrix)
	       offsets: (Sparse-matrix-operator-offsets matrix)
	       matrix-implementation: matrix
	       domain-space: domain-space
	       range-space: range-space))

(define-handy-method (Vector-add (m Linear-finite-element-operator) n)
  (if (check-same-class m n)
      (let ((matrix (Vector-add matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
	(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
      (error "Vector-add: Linear-finite-element-operators are not the same class: " m n)))

(define-handy-method (Vector-subtract (m Linear-finite-element-operator) n)
  (if (check-same-class m n)
      (let ((matrix (Vector-subtract matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
	(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
      (error "Vector-subtract: Linear-finite-element-operators are not the same class: " m n)))

(define-handy-method (Vector-scale a (m Linear-finite-element-operator))
  (let ((matrix (Vector-scale a matrix-implementation)))
    (wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space)))

(define-handy-method (Vector-*axpy a (m Linear-finite-element-operator) n)
  (if (check-same-class m n)
      (let ((matrix (Vector-*axpy a matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
	(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
      (error "Vector-*axpy: Linear-finite-element-operators are not the same class: " m n)))

(define-handy-method (Vector-Linf-norm (m Operator-with-matrix-implementation))
  (Vector-Linf-norm matrix-implementation))

;; we compose two linear-finite-element-operators by composing their matrix implementations
;; and building a new operator with instantiate.  The immutability of the fields of the
;; Finite-element-operator gets in the way here, since we may like to define this method
;; on Finite-element-operators and clone one of the arguments to make it work with subclasses
;; too, but we'll let future generations worry about this ;-).

(define-method (Operator-compose (op1 Linear-finite-element-operator) (op2 Linear-finite-element-operator op2))
  (if (eq? (Linear-finite-element-operator-range-space op2)
	   (Linear-finite-element-operator-domain-space op1))
      (let ((new-implementation
	     (Operator-compose (Linear-finite-element-operator-matrix-implementation op1)
			       (Linear-finite-element-operator-matrix-implementation op2))))
	(wrap-matrix-implementation-in-linear-finite-element-operator new-implementation 
								      (Linear-finite-element-operator-domain-space op2)
								      (Linear-finite-element-operator-range-space op1)))
      (error "Operator-compose: Linear-finite-element-operators do not have compatable domains and ranges.")))

(define-handy-method (Operator->identity (m Linear-finite-element-operator))
  (if (eq? range-space
	   domain-space)
      (let ((matrix (Operator->identity matrix-implementation)))
	(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
      (error "Operator->identity: Linear-finite-element-operator does not have same domain and range." m)))

;; Takes a Linear-finite-element-operator and returns the adjoint operator.

(define-handy-method (Operator-adjoint (m Linear-finite-element-operator))
  (let ((matrix (Operator-adjoint matrix-implementation)))
    (wrap-matrix-implementation-in-linear-finite-element-operator matrix range-space domain-space)))

;; calculate the gradient of the linear function that is 1.0 at p1 and 0.0 at p2 and p3

(define (grad-phi p1 p2 p3)
  (declare (inlining-limit 10000))  ; I want this inlined
  (let ((v (Point-subtract p1 p2))
	(u (Point-subtract p3 p2)))
    (let ((w (Point-subtract v (Point-scale (fl/ (Point-dot-product v u)
						 (Point-dot-product u u))
					    u))))
      (Point-scale (fl/ (Point-dot-product w w)) w))))

;; Evaluate the function a at the midpoint of the Edge e.  A method for
;; curved-edges will come later.

(define-generic (Edge-midpoint (e)))

(define-handy-method (Edge-midpoint (e Edge))
  (declare (inlining-limit 10000))  ; I want this inlined
  (Point-scale 0.5 (Point-add (Vertex-point vertex-1) 
			      (Vertex-point vertex-2))))

;; OK, the following functions should be multi-methods based not only on the class
;; of the triangle, but also on the type of the space.  And the f64vector of as
;; should be encapsulated, too.  But at least this allows us to have a method for
;; Triangle-with-curved-side later.

;; as is an f64vector of the values of a evaluated at the midpoint of the edges
;; in a triangulation.

(define-generic (Triangle-stiffness (t Triangle) as m)

  ;; integrate a \grad \phi_j\cdot\grad\phi_i on t, given a's at midpoints of edges

  (with-access
   t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
   (with-access
    vertex-1 (Vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 (with-access
	  vertex-3 (Vertex point index)
	  (let ((point-3 point)
		(index-3 index))
	    
	    ;; here we rely crucially on the elements having the same numbering as the
	    ;; vertices in the Triangulation.  So grad-1 is the gradient on the triangle
	    ;; of the element that is 1 at vertex-1 and zero at the other vertices, etc.
	    
	    (let ((grad-1 (grad-phi point-1 point-2 point-3))
		  (grad-2 (grad-phi point-2 point-3 point-1))
		  (grad-3 (grad-phi point-3 point-1 point-2))
		  
		  ;; We're going to use the quadrature rule that the integral of a function
		  ;; f on a triangle is approximately equal to the average of f at the midpoints
		  ;; of the three sides of the triangle times the area of the triangle.  This
		  ;; is exact for quadratic functions f.

		  ;; Here, the function we're integrating is the dot-product of the gradients
		  ;; of elements, which is constant on the triangle, times a, so we precompute
		  ;; the term involving a.

		  (weight (fl* (fl+ (f64vector-ref as (Edge-index edge-1))
				    (f64vector-ref as (Edge-index edge-2))
				    (f64vector-ref as (Edge-index edge-3)))
			       0.333333333333333333333333333333333333
			       (Triangle-area t))))
	      (Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (Point-dot-product grad-1 grad-1)))
	      (Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (Point-dot-product grad-2 grad-2)))
	      (Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight (Point-dot-product grad-3 grad-3)))
	      (let ((val (fl* weight (Point-dot-product grad-1 grad-2))))
		(Sparse-matrix-operator-increment-entry m index-1 index-2 val)
		(Sparse-matrix-operator-increment-entry m index-2 index-1 val))
	      (let ((val (fl* weight (Point-dot-product grad-1 grad-3))))
		(Sparse-matrix-operator-increment-entry m index-1 index-3 val)
		(Sparse-matrix-operator-increment-entry m index-3 index-1 val))
	      (let ((val (fl* weight (Point-dot-product grad-2 grad-3))))
		(Sparse-matrix-operator-increment-entry m index-2 index-3 val)
		(Sparse-matrix-operator-increment-entry m index-3 index-2 val)))))))))))

;; Add the stiffness terms, the integral of a times \grad \phi_i\cdot\grad\phi_j to the Linear-finite-element
;; matrix m.

(define (add-stiffness-terms! m a)
  (with-access
   m (Linear-finite-element-operator matrix-implementation domain-space)
   (with-access
    domain-space (Finite-element-space triangulation)
    (let ((m matrix-implementation)
	  (triangulation triangulation)
	  (as 	  ;; make vector to hold the values of a evaluated at all the edges
	   (make-f64vector (Triangulation-edges-length triangulation))))
      
      ;; evaluate a at all the edge mid-points
      
      (Triangulation-for-each-edge
       (lambda (e)
	 (f64vector-set! as (Edge-index e) (a (Edge-midpoint e))))
       triangulation)
      
      ;; calculate the contribution on each triangle to the terms involving
      ;; the elements that overlap that triangle nontrivially.
      
      (Triangulation-for-each-triangle
       (lambda (t)
	 (Triangle-stiffness t as m))
       triangulation)))))

;; The same principle as Triangle-stiffness, but here we're working with
;; vectors of Points (bs).

(define-generic (Triangle-transport (t Triangle) bs m)

  ;; integrate b \cdot \grad\phi_j \phi_i on t given b values

  (with-access
   t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
   (with-access
    vertex-1 (Vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 (with-access
	  vertex-3 (Vertex point index)
	  (let ((point-3 point)
		(index-3 index))
	    (let ((grad-1 (grad-phi point-1 point-2 point-3))
		  (grad-2 (grad-phi point-2 point-3 point-1))
		  (grad-3 (grad-phi point-3 point-1 point-2))
		  (edge-1-b (vector-ref bs (Edge-index edge-1)))
		  (edge-2-b (vector-ref bs (Edge-index edge-2)))
		  (edge-3-b (vector-ref bs (Edge-index edge-3)))
		  (weight (fl* (Triangle-area t)               
			       ;; average of 3 points time 1/2, phi_i at midpoint of each edge
			       0.166666666666666666666666)))
	      (Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight 
									     (fl+ (Point-dot-product grad-1 edge-3-b)
										  (Point-dot-product grad-1 edge-1-b))))
	      (Sparse-matrix-operator-increment-entry m index-2 index-1 (fl* weight 
									     (fl+ (Point-dot-product grad-1 edge-1-b)
										  (Point-dot-product grad-1 edge-2-b))))
	      (Sparse-matrix-operator-increment-entry m index-3 index-1 (fl* weight 
									     (fl+ (Point-dot-product grad-1 edge-2-b)
										  (Point-dot-product grad-1 edge-3-b))))
	      (Sparse-matrix-operator-increment-entry m index-1 index-2 (fl* weight 
									     (fl+ (Point-dot-product grad-2 edge-3-b)
										  (Point-dot-product grad-2 edge-1-b))))
	      (Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight 
									     (fl+ (Point-dot-product grad-2 edge-1-b)
										  (Point-dot-product grad-2 edge-2-b))))
	      (Sparse-matrix-operator-increment-entry m index-3 index-2 (fl* weight 
									     (fl+ (Point-dot-product grad-2 edge-2-b)
										  (Point-dot-product grad-2 edge-3-b))))
	      (Sparse-matrix-operator-increment-entry m index-1 index-3 (fl* weight 
									     (fl+ (Point-dot-product grad-3 edge-3-b)
										  (Point-dot-product grad-3 edge-1-b))))
	      (Sparse-matrix-operator-increment-entry m index-2 index-3 (fl* weight 
									     (fl+ (Point-dot-product grad-3 edge-1-b)
										  (Point-dot-product grad-3 edge-2-b))))
	      (Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight 
									     (fl+ (Point-dot-product grad-3 edge-2-b)
										  (Point-dot-product grad-3 edge-3-b)))))))))))))


(define (add-transport-terms! m b)
  (with-access
   m (Linear-finite-element-operator matrix-implementation domain-space)
   (with-access
    domain-space (Finite-element-space triangulation)
    (let ((m matrix-implementation)
	  (triangulation triangulation)
	  (bs 	  ;; make vector to hold the values of b evaluated at all the edges
	   (make-vector (Triangulation-edges-length triangulation))))
      
      
      (Triangulation-for-each-edge
       (lambda (e)
	 (vector-set! bs (Edge-index e) (b (Edge-midpoint e))))
       triangulation)
      
      (Triangulation-for-each-triangle
       (lambda (t)
	 (Triangle-transport t bs m))
       triangulation)))))

(define-generic (Triangle-mass (t Triangle) cs m)
  
  ;; integrates c * \phi_i * \phi_j on t given c's on edge midpoints.

  (with-access
   t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
   (with-access
    vertex-1 (Vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 (with-access
	  vertex-3 (Vertex point index)
	  (let ((point-3 point)
		(index-3 index))
	    (let ((edge-1-c (f64vector-ref cs (Edge-index edge-1)))
		  (edge-2-c (f64vector-ref cs (Edge-index edge-2)))
		  (edge-3-c (f64vector-ref cs (Edge-index edge-3))))
	      
	      ;; the value of \phi_i\phi_j at the midpoint of an edge is 1/4 or zero;
	      ;; here we precalculate the weight 1/4 * 1/3 * area of t
	      
	      (let ((weight (fl* 0.08333333333333333333333333333333333333
				 (Triangle-area t))))
		(Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (fl+ edge-1-c edge-3-c)))
		(Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (fl+ edge-2-c edge-1-c)))
		(Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight (fl+ edge-3-c edge-2-c)))
		(let ((val (fl* weight edge-1-c)))
		  (Sparse-matrix-operator-increment-entry m index-1 index-2 val)
		  (Sparse-matrix-operator-increment-entry m index-2 index-1 val))
		(let ((val (fl* weight edge-3-c)))
		  (Sparse-matrix-operator-increment-entry m index-1 index-3 val)
		  (Sparse-matrix-operator-increment-entry m index-3 index-1 val))
		(let ((val (fl* weight edge-2-c)))
		  (Sparse-matrix-operator-increment-entry m index-2 index-3 val)
		  (Sparse-matrix-operator-increment-entry m index-3 index-2 val))))))))))))



(define (add-mass-terms! m c)
  (with-access
   m (Linear-finite-element-operator matrix-implementation  domain-space)
   (with-access
    domain-space (Finite-element-space triangulation)
    (let ((m matrix-implementation)
	  (triangulation triangulation)
	  (cs 	  ;; make vector to hold the value of c evaluated at all the edges
	   (make-f64vector (Triangulation-edges-length triangulation))))
      
      ;; evaluate c at all the edge mid-points
      
      (Triangulation-for-each-edge
       (lambda (e)
	 (f64vector-set! cs (Edge-index e) (c (Edge-midpoint e))))
       triangulation)
            
      (Triangulation-for-each-triangle
       (lambda (t)
	 (Triangle-mass t cs m))
       triangulation)))))

;; For integration on the boundary, we use the two-point Guass-rule on [0,1],
;; which has the advantage of not evaluating functions at vertices on the boundary,
;; where, e.g., there may be a discontinuity in the normal derivative of u
;; because the normal to the boundary points in different directions on each side
;; of a vertex.

(define left-gauss-point  .21132486540518708) ; (* 0.5 (- 1.0 (/ (sqrt 3.0))))

(define right-gauss-point .7886751345948129)  ; (* 0.5 (+ 1.0 (/ (sqrt 3.0))))

(define-generic (Edge-boundary-term (e Edge) gamma m)

  ;; integrates gamma * phi_i * phi_j on an edge, increment the correct coefficient

  (with-access
   e (Edge vertex-1 vertex-2 index)
   (with-access
    vertex-1 (Boundary-vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Boundary-vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 ;; We're considering point-1 to be the left point and point-2 to be the right (it doesn't matter).
	 ;; So the gauss point near point-1 is right-gauss-point * point-1 + left-gauss-point * point-2.
	 
	 (let ((val-1 (gamma (Point-add (Point-scale right-gauss-point point-1)
					(Point-scale left-gauss-point point-2))))
	       (val-2 (gamma (Point-add (Point-scale left-gauss-point point-1)
					(Point-scale right-gauss-point point-2))))
	       (weight (fl* (Edge-length e) 0.5)))
	   
	   ;; The basis element that is 1 at point-1 takes the value right-gauss-point at left-gauss-point, and
	   ;; the value left-gauss-point at right-gauss-point.  (Yes, I know it's confusing.)
	   
	   (Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (fl+ (fl* right-gauss-point right-gauss-point val-1)
										      (fl* left-gauss-point left-gauss-point val-2))))
	   (Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (fl+ (fl* left-gauss-point left-gauss-point val-2)
										      (fl* right-gauss-point right-gauss-point val-2))))
	   (Sparse-matrix-operator-increment-entry m index-1 index-2 (fl* weight left-gauss-point right-gauss-point (fl+ val-1 val-2)))
	   (Sparse-matrix-operator-increment-entry m index-2 index-1 (fl* weight left-gauss-point right-gauss-point (fl+ val-1 val-2))))))))))

(define (add-boundary-terms! m gamma)
  (with-access
   m (Linear-finite-element-operator matrix-implementation domain-space)
   (let ((m matrix-implementation))
     (with-access
      domain-space (Finite-element-space triangulation)
      (let ((triangulation triangulation)
	    (n (Triangulation-boundary-vertices-length triangulation)))
	(let ((boundary-vertex-gammas (make-f64vector n)))
	  
	  ;; go through the boundary edges, calculate the length of each edge, and then increment
	  ;; the contribution of that edge to the integrals
	  
	  (Triangulation-for-each-boundary-edge
	   (lambda (e)
	     (Edge-boundary-term e gamma m))
	   triangulation)))))))

;; The functions to calculate the right-hand-side of the matrix problems.

(define-generic (Triangle-interior-integral (t Triangle) fs data)

  ;; integrate f \phi_i on t, given f values at edge midpoints.

  (with-access
   t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
   (with-access
    vertex-1 (Vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 (with-access
	  vertex-3 (Vertex point index)
	  (let ((point-3 point)
		(index-3 index))
	    (let ((edge-1-f (f64vector-ref fs (Edge-index edge-1)))
		  (edge-2-f (f64vector-ref fs (Edge-index edge-2)))
		  (edge-3-f (f64vector-ref fs (Edge-index edge-3))))
	      (let ((weight (fl* 0.1666666666666666666666666666666666666
				 (Triangle-area t))))
		(f64vector-set! data index-1 (fl+ (f64vector-ref data index-1)
						  (fl* weight (fl+ edge-3-f edge-1-f))))
		(f64vector-set! data index-2 (fl+ (f64vector-ref data index-2)
						  (fl* weight (fl+ edge-1-f edge-2-f))))
		(f64vector-set! data index-3 (fl+ (f64vector-ref data index-3)
						  (fl* weight (fl+ edge-2-f edge-3-f))))))))))))))


(define (add-interior-integral-terms! v f)
  (with-access
   v (Linear-finite-element-vector data
				   space)
   (with-access
    space (Finite-element-space triangulation)
    (let ((fs 	  ;; make vector to hold the value of f evaluated at all the edges
	   (make-f64vector (Triangulation-edges-length triangulation)))
	  (data data)
	  (triangulation triangulation))
      
      ;; evaluate f at all the edge mid-points
      
      (Triangulation-for-each-edge
       (lambda (e)
	 (f64vector-set! fs (Edge-index e) (f (Edge-midpoint e))))
       triangulation)
      
      (Triangulation-for-each-triangle
       (lambda (t)
	 (Triangle-interior-integral t fs data))
       triangulation)))))

;; Again, here we use the two-point Gauss rule.

(define-generic (Edge-boundary-integral (e Edge) g data)

  ;; integrate g * \phi_i on the edge e.

  (with-access
   e (Edge vertex-1 vertex-2 index)
   (with-access
    vertex-1 (Vertex point index)
    (let ((point-1 point)
	  (index-1 index))
      (with-access
       vertex-2 (Vertex point index)
       (let ((point-2 point)
	     (index-2 index))
	 (let ((val-1 (g (Point-add (Point-scale right-gauss-point point-1)
				    (Point-scale left-gauss-point point-2))))
	       (val-2 (g (Point-add (Point-scale left-gauss-point point-1)
				    (Point-scale right-gauss-point point-2))))
	       (weight (fl* (Edge-length e) 0.5)))
	   (f64vector-set!
	    data 
	    index-1 
	    (fl+ (f64vector-ref data index-1)
		 (fl* weight (fl+ (fl* right-gauss-point val-1)            ; this element has the value right-gauss-point at the Gauss point near point-1
				  (fl* left-gauss-point val-2)))))         ; and the value left-gauss-point at the Gauss point near point-2
	   (f64vector-set!
	    data 
	    index-2
	    (fl+ (f64vector-ref data index-2)
		 (fl* weight (fl+ (fl* left-gauss-point val-1)             ; this element has the value left-gauss-point at the Gauss point near point-1
				  (fl* right-gauss-point val-2)))))))))))) ; and the value right-gauss-point at the Gauss point near point-2




(define (add-boundary-integral-terms! v g)
  (with-access
   v (Linear-finite-element-vector space data)
   (with-access
    space (Finite-element-space triangulation)
    (let ((triangulation triangulation)
	  (data data))
      
      ;; go through the boundary edges, calculate the length of each edge, and then increment
      ;; the contribution of that edge to the integrals
      
      (Triangulation-for-each-boundary-edge
       (lambda (e)
	 (Edge-boundary-integral e g data))
       triangulation)))))