;;; An Element is a basis element of a finite-element space.
;;; All the elements in a space are ordered and indexed;
;;; it's useful if the element knows its index (I think).
;;; An element's neighbors are the indices of the other elements
;;; whose support has nontrivial intersection with this element.
(define-class Element Object
((= index)
(= neighbors))) ; including itself
;; we ensure that the index is the *last* element of the neighbors list.
(define-handy-method (initialize! (e Element))
(let loop ((nabors (Neighbors->list neighbors)) (result (list index)))
(if (null? nabors)
(set! neighbors (Neighbors<-list result))
(loop (cdr nabors)
(if (eq? (car nabors) index)
result
(cons (car nabors) result)))))
(call-next-method))
;;; Neighbors will be used in various parts of the the system, so
;;; we turn it into an abstract data type with
;;; constructors, deconstructors, accessors, and modifiers.
(define (Neighbors->list n)
(vector->list n))
(define (Neighbors<-list l)
(list->vector l))
(define (Neighbors-ref n i)
(vector-ref n i))
(define (Neighbors-set! n i val)
(vector-set! n i val))
(define (Neighbors-length n)
(vector-length n))
;;; I anticipate there being Cubic-elements, mixed elements, etc.
;;; Now we will use just Linear-element (piecewise linear, continuous
;;; on triangles).
(define-class Linear-element Element ())
;;; A finite element space of any kind is associated with a triangulation,
;;; and has elements, each of which is an Element.
;;; Note that we're going to layer this as much as possible; the
;;; Finite-element-space is an object that doesn't contain the
;;; numerical coordinates of any vector in the space, that comes in
;;; the Finite-element-vector.
(define-class Finite-element-space Object
((= triangulation)
(= elements) ; the elements of the finite element space
(= neighbors))) ; the elements' neighbors
;;; Neighbors and elements are stored in a vector
(define (Finite-element-space-elements->list e)
(vector->list e))
(define (Finite-element-space-elements<-list l)
(list->vector l))
(define (Finite-element-space-neighbors->list e)
(vector->list e))
(define (Finite-element-space-neighbors<-list l)
(list->vector l))
(define (Finite-element-space-neighbors-length space)
(vector-length (Finite-element-space-neighbors space)))
;;; A linear finite element space, i.e., a space containing Linear-elements.
(define-class Linear-finite-element-space Finite-element-space ())
;;; A Finite-element-vector contains a (pointer to) a Finite-element-space
;;; and an f64vector containing the numerical data of the vector.
;;; Note: this class is treated as a parameterized class, with the parameter
;;; being the space.
(define-class Finite-element-vector Object
((= space)
(= data maybe-uninitialized: immutable:)))
(define-handy-method (initialize! (v Finite-element-vector))
;; set it to zero if not otherwise initialized
(if (not (field-defined? v 'data))
(initialize-field-value!
v
(make-f64vector (Finite-element-space-neighbors-length space) 0.0)
'data))
(call-next-method))
;;; Since it's a subclass of Vector, we need to define all the
;;; methods for that space. We use the generic methods
;;; on the data fields.
(define-method (Vector-dimension (x Finite-element-vector))
(Vector-dimension (Finite-element-vector-data x)))
(define-handy-method (check-same-class (x Finite-element-vector) y)
(and (call-next-method)
(eq? space (Finite-element-vector-space y))))
(define-method (Vector-add (x Finite-element-vector) y)
(if (check-same-class x y)
(instantiate-from (x Finite-element-vector)
data: (Vector-add (Finite-element-vector-data x)
(Finite-element-vector-data y)))
(error "Vector-add---underlying finite element spaces are not the same")))
(define-method (Vector-subtract (x Finite-element-vector) y)
(if (check-same-class x y)
(instantiate-from (x Finite-element-vector)
data: (Vector-subtract (Finite-element-vector-data x)
(Finite-element-vector-data y)))
(error "Vector-subtract---underlying finite element spaces are not the same")))
(define-handy-method (Vector-scale x (y Finite-element-vector))
(instantiate-from (y Finite-element-vector)
data: (Vector-scale x data)))
(define-method (Vector-*axpy a (x Finite-element-vector) y)
(if (check-same-class x y)
(instantiate-from (x Finite-element-vector)
data: (Vector-*axpy a
(Finite-element-vector-data x)
(Finite-element-vector-data y)))
(error "Vector-*axpy---underlying finite element spaces are not the same")))
(define-method (Vector-dot-product (x Finite-element-vector) y)
(if (check-same-class x y)
(Vector-dot-product (Finite-element-vector-data x)
(Finite-element-vector-data y))
(error "Vector-dot-product---underlying finite element spaces are not the same")))
(define-handy-method (Vector-zero (x Finite-element-vector))
(instantiate-from (x Finite-element-vector)
data: (Vector-zero data)))
(define-handy-method (Vector-copy (x Finite-element-vector))
;; new copy of data, don't copy space
(instantiate-from (x Finite-element-vector)
data: (Vector-copy data)))
;;; I thought it would be useful to define a subclass for Linear-finite-element-vector
(define-class Linear-finite-element-vector Finite-element-vector ())
;;; It did prove to be useful---I'll define max and min for
;;; Linear-finite-element-vectors since it makes sense there (somewhat)
;;; with the Lagrange basis elements; I won't define it for general
;;; Finite-element-vectors.
(define-method (Lattice-min (x Linear-finite-element-vector) y)
(if (check-same-class x y)
(instantiate-from (x Finite-element-vector)
data: (Lattice-min (Finite-element-vector-data x)
(Finite-element-vector-data y)))
(error "Lattice-min---underlying finite element spaces are not the same")))
(define-method (Lattice-max (x Linear-finite-element-vector) y)
(if (check-same-class x y)
(instantiate-from (x Finite-element-vector)
data: (Lattice-max (Finite-element-vector-data x)
(Finite-element-vector-data y)))
(error "Lattice-max---underlying finite element spaces are not the same")))
;;; it's useful for testing purposes to be able to interpolate a given function to
;;; a linear finite element space
(define (interpolate-procedure->Linear-finite-element-vector u fes)
;; u is a function of Point p
;; fes is a Linear-finite-element-space
(let ((t (Finite-element-space-triangulation fes)))
(instantiate Linear-finite-element-vector
data: (let ((result (make-f64vector (Triangulation-vertices-length t))))
;; we've constructed the f64vector, now we fill it with the values
;; of u at (Vertex-point v) of each vertex in t
(Triangulation-for-each-vertex
(lambda (v)
(with-access
v (Vertex point index)
(f64vector-set! result index (u point))))
t)
;; return it so it can be placed in Linear-finite-element-vector-data
result)
space: fes)))