;;; An Element is a basis element of a finite-element space.
;;; All the elements in a space are ordered and indexed;
;;; it's useful if the element knows its index (I think).
;;; An element's neighbors are the indices of the other elements
;;; whose support has nontrivial intersection with this element.

(define-class Element Object
  ((= index)
   (= neighbors))) ; including itself

;; we  ensure that the index is the *last* element of the neighbors list.

(define-handy-method (initialize! (e Element))

  (let loop ((nabors (Neighbors->list neighbors)) (result (list index)))
    (if (null? nabors)
	(set! neighbors (Neighbors<-list result))
	(loop (cdr nabors) 
	      (if (eq? (car nabors) index)
		  result
		  (cons (car nabors) result)))))
  
  (call-next-method))


;;; Neighbors will be used in various parts of the the system, so
;;; we turn it into an abstract data type with
;;; constructors, deconstructors, accessors, and modifiers.

(define (Neighbors->list n)
  (vector->list n))

(define (Neighbors<-list l)
  (list->vector l))

(define (Neighbors-ref n i)
  (vector-ref n i))

(define (Neighbors-set! n i val)
  (vector-set! n i val))

(define (Neighbors-length n)
  (vector-length n))

;;; I anticipate there being Cubic-elements, mixed elements, etc.
;;; Now we will use just Linear-element (piecewise linear, continuous
;;; on triangles).

(define-class Linear-element Element ())

;;; A finite element space of any kind is associated with a triangulation,
;;; and has elements, each of which is an Element.
;;; Note that we're going to layer this as much as possible; the
;;; Finite-element-space is an object that doesn't contain the
;;; numerical coordinates of any vector in the space, that comes in
;;; the Finite-element-vector.

(define-class Finite-element-space Object
  ((= triangulation)
   (= elements)        ; the elements of the finite element space
   (= neighbors)))     ; the elements' neighbors

;;; Neighbors and elements are stored in a vector

(define (Finite-element-space-elements->list e)
  (vector->list e))

(define (Finite-element-space-elements<-list l)
  (list->vector l))

(define (Finite-element-space-neighbors->list e)
  (vector->list e))

(define (Finite-element-space-neighbors<-list l)
  (list->vector l))

(define (Finite-element-space-neighbors-length space)
  (vector-length (Finite-element-space-neighbors space)))

;;; A linear finite element space, i.e., a space containing Linear-elements.

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

;;; A Finite-element-vector contains a (pointer to) a Finite-element-space
;;; and an f64vector containing the numerical data of the vector.
;;; Note: this class is treated as a parameterized class, with the parameter
;;; being the space.

(define-class Finite-element-vector Object
  ((= space)
   (= data maybe-uninitialized: immutable:)))

(define-handy-method (initialize! (v Finite-element-vector))
  ;; set it to zero if not otherwise initialized
  (if (not (field-defined? v 'data))
      (initialize-field-value!
       v
       (make-f64vector (Finite-element-space-neighbors-length space) 0.0)
       'data))
  (call-next-method))

;;; Since it's a subclass of Vector, we need to define all the
;;; methods for that space.  We mainly use the generic methods
;;; on the data fields.

(define-method (Vector-dimension (x Finite-element-vector))
  (f64vector-length (Finite-element-vector-data x)))

(define-handy-method (check-same-class (x Finite-element-vector) y)
  (and (call-next-method)
       (eq? space (Finite-element-vector-space y))))

(define-method (Vector-add (x Finite-element-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Finite-element-vector)
			data: (Vector-add (Finite-element-vector-data x)
					  (Finite-element-vector-data y)))
      (error "Vector-add---underlying finite element spaces are not the same")))

(define-method (Vector-subtract (x Finite-element-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Finite-element-vector)
			data: (Vector-subtract (Finite-element-vector-data x)
					       (Finite-element-vector-data y)))
      (error "Vector-subtract---underlying finite element spaces are not the same")))

(define-handy-method (Vector-scale x (y Finite-element-vector))
  (instantiate-from (y Finite-element-vector)
		    data: (Vector-scale x data)))

(define-method (Vector-*axpy a (x Finite-element-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Finite-element-vector)
			data: (Vector-*axpy a
					    (Finite-element-vector-data x)
					    (Finite-element-vector-data y)))
      (error "Vector-*axpy---underlying finite element spaces are not the same")))

(define-method (Vector-dot-product (x Finite-element-vector) y)
  (if (check-same-class x y)
      (Vector-dot-product (Finite-element-vector-data x)
			  (Finite-element-vector-data y))
      (error "Vector-dot-product---underlying finite element spaces are not the same")))

(define-handy-method (Vector-zero (x Finite-element-vector))
  (instantiate-from (x Finite-element-vector)
		    data: (Vector-zero data)))

(define-handy-method (Vector-copy (x Finite-element-vector))
  ;; new copy of data, don't copy space
  (instantiate-from (x Finite-element-vector)
		    data: (Vector-copy data)))

;;; I thought it would be useful to define a subclass for Linear-finite-element-vector

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

;;; It did prove to be useful---I'll define max and min for
;;; Linear-finite-element-vectors since it makes sense there (somewhat)
;;; with the Lagrange basis elements; I won't define it for general
;;; Finite-element-vectors.

(define-method (Lattice-min (x Linear-finite-element-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Finite-element-vector)
			data: (Lattice-min (Finite-element-vector-data x)
					   (Finite-element-vector-data y)))
      (error "Lattice-min---underlying finite element spaces are not the same")))

(define-method (Lattice-max (x Linear-finite-element-vector) y)
  (if (check-same-class x y)
      (instantiate-from (x Finite-element-vector)
			data: (Lattice-max (Finite-element-vector-data x)
					   (Finite-element-vector-data y)))
      (error "Lattice-max---underlying finite element spaces are not the same")))

;;; it's useful for testing purposes to be able to interpolate a given function to
;;; a linear finite element space

(define (interpolate-procedure->Linear-finite-element-vector u fes)

  ;; u is a function of Point p
  ;; fes is a Linear-finite-element-space

  (let ((t (Finite-element-space-triangulation fes)))

    (instantiate Linear-finite-element-vector
		 data: (let ((result (make-f64vector (Triangulation-vertices-length t))))
			 ;; we've constructed the f64vector, now we fill it with the values
			 ;; of u at (Vertex-point v) of each vertex in t
			 (Triangulation-for-each-vertex
			  (lambda (v)
			    (with-access
			     v (Vertex point index)
			     (f64vector-set! result index (u point))))
			  t)
			 ;; return it so it can be placed in Linear-finite-element-vector-data
			 result)
		 space: fes)))