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

(declare (not core))
(include "points.scm")
(declare (core))

(define a (lambda (p) 1.0))

(define c (lambda (p) 1.0))

(define f (lambda (t p) 
	    (let ((x (Point-x p))
		  (y (Point-y p))) 
	      (fl+ (fl+ (fl* t (fl* t (fl* (fl* x (fl- 1.0 x)) (fl* y (fl- 1.0 y)))))    
			(fl* (fl* 2.0 (fl* t (fl* x (fl- 1.0 x)) (fl* (fl* y (fl- 1.0 y)))))))
		   (fl+ (fl* 2.0 (fl* t (fl* t (fl* y (fl- 1.0 y)))))
			(fl* 2.0 (fl* t (fl* t (fl* x (fl- 1.0 x))))))))))
  
(define g (lambda (t p)
	    (let ((x (Point-x p))
		  (y (Point-y p)))
	      (fl+  (fl* t (fl* t (fl* y (fl- y 1.0)))) 
		    (fl* t (fl* t (fl* x (fl- x 1.0))))))))
      
(define gamma (lambda (p) 0.0))

(define u0 (lambda (p) 0.0))

(define (u t p)
  (let ((x (Point-x p))
	(y (Point-y p)))
    (fl* t t x (fl- 1.0 x) y (fl- 1.0 y))))

(define N 128)

(define T 2.0)

(define K 5)

(define delta-t (/ T N))

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

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

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

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

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

(define stopping-criterion (lambda (A preconditioner z p r x c m)
			     (or  (fl< c 1.0e-20)
				  (fx= m 1000))))

(define (test-solution)
  (let ((solution (time (parabolic-equation-solution a c f g gamma u0 fes T N K))))
    (do ((i 0 (fx+ i 1)))
	((fx> i N))
      (let ((rhs-vector (instantiate Linear-finite-element-vector
				     space: fes)))
	(add-interior-integral-terms! rhs-vector (lambda (p) ( u (fl* delta-t (fixnum->flonum i)) p)))
	(let* ((u-at-time-t (conjugate-gradient L2-inner-product-operator
						rhs-vector 
						stopping-criterion))
	       (difference (Vector-subtract (vector-ref solution i) u-at-time-t)))
	  (display (sqrt (Vector-inner-product difference L2-inner-product-operator difference)))
	  (newline)
	  (newline))))))