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

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


(define t coarse-square-triangulation)

'(set! t (time (refine (refine (refine (refine (refine (refine (refine (refine (refine t)))))))))))

(set! t (time (refine (refine (refine (refine (refine (refine (refine t)))))))))

'(set! t (time (refine (refine (refine (refine (refine (refine t))))))))

'(set! t (refine (refine t)))


(time (Triangulation-consistent? t))

(define fes (time (Triangulation->Linear-finite-element-space t)))


(define m (time (->Linear-finite-element-operator fes)))


(time (add-stiffness-terms! m (lambda (p) 1.0)))


(time (add-mass-terms! m (lambda (p) 1.0)))


(time (add-transport-terms! m (lambda (p) (make-Point 0.0 0.0))))


(time (add-boundary-terms! m (lambda (p) 0.0)))


(define vec1 (time (instantiate Linear-finite-element-vector
			  space: fes)))


(time (add-interior-integral-terms! vec1 
				   (lambda (p)
				     (let ((x (Point-x p))
					   (y (Point-y p)))
				       (declare (flonum))
				       (- (* x x) (* y y))))))


(time (add-boundary-integral-terms! vec1 
				   (lambda (p)
				     (let ((x (Point-x p))
					   (y (Point-y p)))
				       (declare (flonum))
				       (cond ((= x 0.0)
					      0.0)
					     ((= x 1.0)
					      2.0)
					     ((= y 0.0)
					      0.0)
					     ((= y 1.0)
					      -2.0)
					     (else
					      (error "shouldn't happen")))))))


(define result (time (conjugate-gradient
		      m
		      vec1
		      (lambda (A preconditioner z p r x c m)
			(or (< c 1.0e-20)
			    (= m 200)))
		      (gauss-seidel-preconditioner m))))

(define (iter i)
  (do ((j 0 (+ j 1))
       (vec vec1 (Operator-apply m vec1)))
      ((= j i) vec)))

(time (iter 100))

(define check (make-f64vector (Vector-dimension result)))


(vector-for-each-with-index-1 (lambda (i v)
				(with-access
				 v (Vertex point)
				 (f64vector-set! check i (- (* (Point-x point)
							       (Point-x point))
							    (* (Point-y point)
							       (Point-y point))))))
			      (Triangulation-vertices t))

(pp (time (Vector-L2-norm (Vector-subtract check (Finite-element-vector-data result)))))