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

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

(define number-of-refinements 5)

(define fine-mg-data
  (MG-build-data (instantiate Obstacle-multigrid-data) coarse-square-triangulation (+ number-of-refinements 1)))

(define mg-data
  (with-access 
   fine-mg-data (Obstacle-multigrid-data triangulations spaces projectors injectors interpolating-projectors)
   (instantiate Obstacle-multigrid-data
		triangulations:           (list->vector (reverse (cdr (reverse (vector->list triangulations)))))
		spaces:                   (list->vector (reverse (cdr (reverse (vector->list spaces)))))
		projectors:               (list->vector (reverse (cdr (reverse (vector->list projectors)))))
		injectors:                (list->vector (reverse (cdr (reverse (vector->list injectors)))))
		interpolating-projectors: (list->vector (reverse (cdr (reverse (vector->list interpolating-projectors))))))))

'(for-each (lambda (data)
	    (with-access
	     data (Obstacle-multigrid-data triangulations spaces projectors injectors interpolating-projectors)
	     (newline)
	     (display triangulations) (newline)
	     (display spaces) (newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-domain-space p)) (cdr (vector->list projectors))))(newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-range-space p)) (cdr (vector->list projectors))))(newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-domain-space p)) (cdr (vector->list injectors))))(newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-range-space p)) (cdr (vector->list injectors))))(newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-domain-space p)) (cdr (vector->list interpolating-projectors))))(newline)
	     (display (map (lambda (p) (Linear-finite-element-operator-range-space p)) (cdr (vector->list interpolating-projectors))))(newline)))
	  (list fine-mg-data mg-data))

(define t (vector-ref (Multigrid-data-triangulations mg-data) number-of-refinements))

(define fes (vector-ref (Multigrid-data-spaces mg-data) number-of-refinements))

(define fine-t (vector-ref (Multigrid-data-triangulations fine-mg-data) (+ number-of-refinements 1)))

(define fine-fes (vector-ref (Multigrid-data-spaces fine-mg-data) (+ number-of-refinements 1)))

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

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

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

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

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

(define coarse->fine-operator (Operator-adjoint (coarse<-fine fes fine-fes)))

;; we've defined the following nonsymmetric tests, and tried to solve them by setting up
;; the normal equations in MG-add-operator-and-rhs-information! (i.e., calculating m-transpose
;; and multiplying both m and rhs by m-transpose), but it didn't work at all, except for the
;; first problem for some reason (and even then it didn't work well).  So support
;; for nonsymmetric problems in multigrid has been disabled.

(define (run-test check solver-thunk)
  
  (define (format-for-display x)
    (if (< x 1.0)
	(let ((x (+ 1.0 x)))
	  (let ((string-x (number->string x)))
	    (substring string-x 1 12)))
	(number->string x)))
  
  (let* ((at-start (##process-statistics))
	 (result (solver-thunk))
	 (at-end (##process-statistics))
	 (user+sys-ns (+ (- (f64vector-ref at-end 0) (f64vector-ref at-start 0))
			 (- (f64vector-ref at-end 1) (f64vector-ref at-start 1))))
	 (gc-ns (+ (- (f64vector-ref at-end 3) (f64vector-ref at-start 3))
		   (- (f64vector-ref at-end 4) (f64vector-ref at-start 4))))
	 (difference (Vector-subtract (Operator-apply coarse->fine-operator result) check)))
    
    (for-each display (list #\tab
			    (format-for-display (sqrt (Vector-inner-product difference L2-inner-product-operator difference)))))
    (for-each display (list #\tab
			    (format-for-display (sqrt (Vector-inner-product difference H1-inner-product-operator difference)))))
    (for-each display (list #\tab
			    (round (* 1000 (- user+sys-ns gc-ns)))))
    (for-each display (list #\tab
			    (round (* 1000 gc-ns))
			    #\newline))))

(define (run-multigrid-tests tests)
  
  (for-each (lambda (problem-data)

	      (define check
		(let ()
		  (MG-add-operator-and-rhs-information! fine-mg-data problem-data)
		  (MG-propagate-level-information! fine-mg-data)
		  (MG-add-solver-information! fine-mg-data richardson-smoother! "  Richardson" 24 10 0 2)
		  (MG-solver fine-mg-data)))
	      
	      (for-each display (list #\newline (Problem-data-description problem-data) #\newline))
  
	      (MG-add-operator-and-rhs-information! mg-data problem-data)
	      (MG-propagate-level-information! mg-data)

	      (for-each (lambda (smoother-maker smoother-name)
			  (for-each (lambda (p)
				      (for-each (lambda (m1)
						  
						  (define (iota m)
						    (let loop ((result '()) (m m))
						      (if (< m 0)
							  result
							  (loop (cons m result) (- m 1)))))
						  
						  (for-each (lambda (m2)
							      (for-each (lambda (r)
									  (MG-add-solver-information! mg-data smoother-maker smoother-name r m1 m2 p)
									  (display (Multigrid-data-solver-description mg-data))
									  (run-test check (lambda ()
												   (MG-solver mg-data))))
									'(21 22 23 24)))
							    (list 0)));(iota m1)))
						'(4 5 6 7 8 9)))
				    '(1 2)))
			(list richardson-smoother!
			      jacobi-smoother!
			      gauss-seidel-smoother!)
			(list "  Richardson"
			      "      Jacobi"
			      "Gauss-Seidel")))
	    tests))
(define symmetric-tests
  (FLOAT
   (list 
    ;; u(x,y) = x^2 - y^2 (is harmonic), Neumann boundary conditions,
    ;; a=c=1, b = 0, f = u, g = 0, 2, or -2 depending on the boundary
    (instantiate Obstacle-problem-data
		 description: "a=1, c=1, u=x^2-y^2, Neumann"
		 a: (lambda (p) 1.0)
		 c: (lambda (p) 1.0)
		 f: (lambda (p) 
		      (let ((x (Point-x p))
			    (y (Point-y p)))
			(- (* x x) (* y y))))
		 g: (lambda (p)
		      (let ((x (Point-x p))
			    (y (Point-y p)))
			(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")))))
		 obstacle: (lambda (p)
			     (let ((x (Point-x p))
				   (y (Point-y p)))
			       (- 0.5
				  (* 16.0 (- x 0.5) (- x 0.5))
				  (* 16.0 (- y 0.5) (- y 0.5)))))
		 u: (lambda (p) 
		      (let ((x (Point-x p))
			    (y (Point-y p)))
			(- (* x x) (* y y))))))))