;; The following functions use the Sparse-vector and Sparse-matrix{-operator}
;; machinery in sparse-linear-algebra, 

;; This function defines the Operator that projects an element of the
;; dual of a Linear-finite-element-vector from fine-fes to a
;; Linear-finite-element-vector of coarse-fes.
;; This function uses  the fact that
;; (Vertex-fine vertex) of the coarse triangulation contains
;; the underlying vertex of the fine triangulation, and (Edge-fine e)
;; of the coarse triangulation contains the fine vertex at the
;; midpoint of e in the coarse triangulation.

;; This should be a generic function on any Finite-element-space.

(define (coarse<-fine coarse-space fine-space)
  (with-access
   coarse-space (Finite-element-space triangulation)
   (let ((coarse-triangulation triangulation))
     (let ((pseudo-diagonal-matrix
	    (->Sparse-matrix-operator
	     (instantiate Sparse-matrix                                    ; has Sparse-vector rows
			  rows: (map                                       ; for each vertex
				 (lambda (v)                               ; make a Sparse-vector
				   (with-access
				    v (Vertex fine)
				    (list (make-Entry (Vertex-index fine) ; underlying vertex has coefficient 1.0
						      1.0))))
				 (Triangulation-vertices->list             ; elements have same indices as vertices
				  (Triangulation-vertices coarse-triangulation)))
			  domain-dimension: (Finite-element-space-neighbors-length fine-space)
			  range-dimension: (Finite-element-space-neighbors-length coarse-space))))
	   (off-diagonal-matrix
	    (->Sparse-matrix-operator                                      ; has indices and coefficients
	     (instantiate Sparse-matrix                                    ; has Sparse-vector rows
			  rows: (map                                       ; for each vertex
				 (lambda (v)                               ; make a Sparse-vector
				   (with-access
				    v (Vertex fine edges)
				    (sort!                                 ; sort by index
				     (map (lambda (e)                ; other vertices have coefficient 0.5
					     (make-Entry (Vertex-index (Edge-fine e))
							 0.5))
					   (Vertex-edges->list edges))
				     (lambda (x y)
				       (fx< (Entry-index x) (Entry-index y))))))
				 (Triangulation-vertices->list             ; elements have same indices as vertices
				  (Triangulation-vertices coarse-triangulation)))
			  domain-dimension: (Finite-element-space-neighbors-length fine-space)
			  range-dimension: (Finite-element-space-neighbors-length coarse-space)))))
       (let ((matrix (Vector-add pseudo-diagonal-matrix off-diagonal-matrix)))
	 (with-access
	  off-diagonal-matrix (Sparse-matrix-operator indices offsets range-dimension)
	  (let ((pseudo-diagonal-indices (Sparse-matrix-operator-indices pseudo-diagonal-matrix))
		(indices indices)
		(offsets offsets)
		(range-dimension range-dimension))
	    (Sparse-matrix-operator-mapping-set!
	     matrix
	     (lambda (v)
	       (let ((pseudo-diagonal-indices pseudo-diagonal-indices)
		     (indices    indices)  ; off-diagonal indices
		     (offsets    offsets)  ; off-diagonal-offsets
		     (data v))
		 (let* ((n range-dimension)
			(result (make-f64vector n)))
		   (let outer ((i (fx- n 1))
			       (next-offset (vector-ref offsets n)))
		     (if (fx< i 0)
			 result
			 (let ((offset (vector-ref offsets i)))
			   (cond ((fx= (fx- next-offset offset) 6) ; interior vertices in a refined triangulation
				  (f64vector-set! 
				   result 
				   i 
				   (fl* 0.5 
					(fl+ (fl+ (fl+ (f64vector-ref data (indices-ref indices offset))
						       (f64vector-ref data (indices-ref indices (fx+ offset 1))))
						  (fl+ (f64vector-ref data (indices-ref indices (fx+ offset 2)))
						       (f64vector-ref data (indices-ref indices (fx+ offset 3)))))
					     (fl+ (fl+ (f64vector-ref data (indices-ref indices (fx+ offset 4)))
						       (f64vector-ref data (indices-ref indices (fx+ offset 5))))
						  (fl* 2.0 (f64vector-ref data (indices-ref pseudo-diagonal-indices i)))))))
				  (outer (fx- i 1) offset))
				 ((fx= (fx- next-offset offset) 4) ; boundary vertices in a refined triangulation
				  (f64vector-set! 
				   result 
				   i 
				   (fl+ (f64vector-ref data (indices-ref pseudo-diagonal-indices i))
					(fl* 0.5 
					     (fl+ (fl+ (f64vector-ref data (indices-ref indices offset))
						       (f64vector-ref data (indices-ref indices (fx+ offset 1))))
						  (fl+ (f64vector-ref data (indices-ref indices (fx+ offset 2)))
						       (f64vector-ref data (indices-ref indices (fx+ offset 3))))))))
				  (outer (fx- i 1) offset))
				 (else
				  (do ((j (fx- next-offset 1) (fx- j 1))
				       (inner-result 0.0 (fl+ inner-result 
							      (f64vector-ref data (indices-ref indices j)))))
				      ((fx< j offset) 
				       (f64vector-set! result i (fl+ (f64vector-ref data (indices-ref pseudo-diagonal-indices i))
								     (fl* 0.5 inner-result)))
				       (outer (fx- i 1) offset))))))))))))
	    (wrap-matrix-implementation-in-linear-finite-element-operator matrix fine-space coarse-space))))))))

(define (projector->injector m)
  (with-access
   m (Linear-finite-element-operator domain-space range-space matrix-implementation)
   (let ((matrix (Operator-adjoint matrix-implementation)))
     (with-access
      matrix (Sparse-matrix-operator indices offsets)
      (let ((indices indices)
	    (offsets offsets))
	(Sparse-matrix-operator-mapping-set! 
	 matrix
	 (lambda (data)
	   (let ((indices indices)
		 (offsets offsets))
	     (let* ((n (vector-length offsets))
		    (result (make-f64vector (fx- n 1))))
	       (let loop ((i (fx- n 2)) 
			  (previous-offset (vector-ref offsets (fx- n 1))))
		 (if (fx< i 0)
		     result
		     (let ((offset (vector-ref offsets i)))
		       (cond ((fx= (fx- previous-offset offset) 2)
			      (f64vector-set! result i (fl* 0.5 (fl+ (f64vector-ref data (indices-ref indices offset))
								     (f64vector-ref data (indices-ref indices (fx+ offset 1))))))
			      (loop (fx- i 1) offset))
			     (else
			      (f64vector-set! result i (f64vector-ref data (indices-ref indices offset)))
			      (loop (fx- i 1) offset))))))))))))
     (wrap-matrix-implementation-in-linear-finite-element-operator matrix range-space domain-space))))