(declare (standard-bindings) (extended-bindings) (not safe) (block) (not interrupts-enabled))

(declare (not core))    ; include the file that contains the definitions of points
(include "points.scm")  ; so those definitions can be inlined here.
(declare (core))

(define triangulation (refine (refine (refine (refine (refine (refine coarse-square-triangulation)))))))

(define space (Triangulation->Linear-finite-element-space triangulation))

(define (square x)
  (FLOAT (* x x)))

(define u (solve-transport-problem (lambda (p) 1.0) ; c
				   (lambda (p) 1.0) ; a
				   (lambda (p) (make-Point 2.0 1.0)) ; b
				   (lambda (p) 0.0) ; u0 
				   (FLOAT
				    (let ((pi (* 4.0 (atan 1.0))))
				      (lambda (p) 
					(let ((x (Point-x p))
					      (y (Point-y p)))
					  (- (* pi (+ (* 2.0 (sin (* 2.0 pi x)) (square (sin (* pi y))))
						      (* (sin (* 2.0 pi y)) (square (sin (* pi x))))))
					     (* 2.0 pi pi (+ (* (cos (* 2.0 pi x)) (square (sin (* pi y))))
							     (* (cos (* 2.0 pi y)) (square (sin (* pi x)))))))))))
				   0.0025 ; delta-t
				   space
				   10000 ; N
				   ))

(define L2-inner-product-operator (->Linear-finite-element-operator space))

(add-mass-terms! L2-inner-product-operator (lambda (p) 1.0))

(define H1-inner-product-operator (->Linear-finite-element-operator space))

(add-stiffness-terms! H1-inner-product-operator (lambda (p) 1.0))

(add-mass-terms! H1-inner-product-operator (lambda (p) 1.0))

(define solution (instantiate Linear-finite-element-vector
			      space: space
			      data: (FLOAT
				     (let ((pi (* 4.0 (atan 1.0))))
				       (list->f64vector 
					(map (lambda (v)
					       (with-access
						v (Vertex point)
						(* (square (sin (* pi (Point-x point))))
						   (square (sin (* pi (Point-y point)))))))
					     (Triangulation-vertices->list (Triangulation-vertices triangulation))))))))

(define difference (Vector-subtract u solution))

(define L2-error (sqrt (Vector-inner-product difference L2-inner-product-operator difference)))

(define H1-error (sqrt (Vector-inner-product difference H1-inner-product-operator difference)))