(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))) 
	      (FLOAT
	       (+ (+ (* t (* t (* (* x (- 1.0 x)) (* y (- 1.0 y)))))    
		     (* (* 2.0 (* t (* x (- 1.0 x)) (* (* y (- 1.0 y)))))))
		  (+ (* 2.0 (* t (* t (* y (- 1.0 y)))))
		     (* 2.0 (* t (* t (* x (- 1.0 x)))))))))))
  
(define g (lambda (t p)
	    (let ((x (Point-x p))
		  (y (Point-y p)))
                  (FLOAT
		   (+  (* t (* t (* y (- y 1.0)))) 
		       (* t (* t (* x (- x 1.0)))))))))
      
(define gamma (lambda (p) 0.0))

(define u0 (lambda (p) 0.0))

(define (u t p)
  (FLOAT
   (let ((x (Point-x p))
	 (y (Point-y p)))
     (* t t x (- 1.0 x) y (- 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 (FLOAT (< c 1.0e-20))
				 (FIX (= m 1000)))))

(define (test-solution)
  (let ((solution (time (parabolic-equation-solution a c f g gamma u0 fes T N K))))
    (do ((i 0 (+ i 1)))
	((> i N))
      (let ((rhs-vector (instantiate Linear-finite-element-vector
				     space: fes)))
	(add-interior-integral-terms! rhs-vector (lambda (p) ( u (FLOAT (* delta-t (##flonum.<-fixnum 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))))))