(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 triangulation (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) ; capacity
(lambda (p) 1.0) ; a
(lambda (p) (make-Point 2.0 1.0)) ; b
(lambda (p) 1.0) ; c
(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
100 ; 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)))