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

(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)
  (fl* 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 
				   (let ((pi (fl* 4.0 (atan 1.0))))
				     (lambda (p) 
				       (let ((x (Point-x p))
					     (y (Point-y p)))
					 (fl- (fl* pi (fl+ (fl* 2.0 (flsin (fl* 2.0 pi x)) (square (flsin (fl* pi y))))
							   (fl* (flsin (fl* 2.0 pi y)) (square (flsin (fl* pi x))))))
					      (fl* 2.0 pi pi (fl+ (fl* (flcos (fl* 2.0 pi x)) (square (flsin (fl* pi y))))
								  (fl* (flcos (fl* 2.0 pi y)) (square (flsin (fl* pi x))))))))))
				   0.0025 ; delta-t
				   space
				   1000 ; 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: (let ((pi (fl* 4.0 (flatan 1.0))))
				      (list->f64vector 
				       (map (lambda (v)
					      (with-access
					       v (Vertex point)
					       (fl* (square (flsin (fl* pi (Point-x point))))
						    (square (flsin (fl* 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)))