#|
Our strategy for Dirichlet boundary conditions is to ensure that all
Linear-finite-element-vectors we work with have zero coefficients on the boundary.
I.e., we don't define a separate space with a different (lesser) dimension.
A Dirichlet linear finite element vector is just a regular linear finite element vector
for which all the boundary coefficients are zero.
I suppose we could define a new class for this, and define check-same-class to check
for this, but I don't see the point.
We do define a separate Dirichlet-linear-finite-element-operator class, because
we can change the operator mapping to always return a (virtual)
Dirichlet-linear-finite-element-vector and change the operator's domain? and
range? to actually check that the boundary coefficients are zero. This doesn't
take so long (if there are $O(N)$ vertices, there are $O(\sqrt N)$ boundary
vertices) and it catches a lot of errors.
|#
(define-class Dirichlet-linear-finite-element-operator Linear-finite-element-operator ())
(define-handy-method (initialize! (o Dirichlet-linear-finite-element-operator))
;; set up the usual operator with Neumann boundary conditions
(call-next-method)
;; The Dirichlet operator mapping is the same, except we set the boundary
;; element coefficients to zero.
(let ((neumann-mapping mapping))
(set! mapping (lambda (v)
(->Dirichlet-linear-finite-element-vector! (neumann-mapping v)))))
;; To test the domain, we need to check that the argument's in the right space, but
;; also that the boundary elements have zero coefficients.
(let ((neumann-domain? domain?))
(set! domain?
(lambda (v)
(and (neumann-domain? v)
(Dirichlet-linear-finite-element-vector? v)))))
;; To test the range, we need to check that the result is in the right space, but
;; also that the boundary elements have zero coefficients.
(let ((neumann-range? range?))
(set! range?
(lambda (v)
(and (neumann-range? v)
(Dirichlet-linear-finite-element-vector? v)))))
;; Now we return the modified object
o)
#|
The next thing is a bit tricky.
It takes a Linear-finite-element-operator, makes a new object in the subclass
Dirichlet-linear-finite-element-operator, and copies over all the slot values,
(including the mapping), and calls the initialize! method for Dirichlet-linear-finite-element-operator.
The first thing that does is call the initialize! method for the superclass,
Linear-finite-element-operator (which is why we had to change the initialize! method on
Linear-finite-element-operator to not assume the mapping is uninitialized).
The initialize! method for Dirichlet-linear-finite-element-operator then modifies
the mapping, range?, and domain? leaving all the other structure intact (so we can get
an overestimate of the $L_\infty$ norm of the operator for the Richardson smoother).
(Which is why we had to make those slots mutable in the Operator superclass.)
And we have to leave the call to (call-next-method) in initialize! because someone might
try to build a Dirichlet-linear-finite-element-operator with instantiate.
Whew!
|#
(define-handy-method (->Dirichlet-linear-finite-element-operator (o Linear-finite-element-operator))
(duplicate (o Dirichlet-linear-finite-element-operator)))
(define (->Dirichlet-linear-finite-element-vector! v)
(with-access
v (Linear-finite-element-vector data space)
(with-access
space (Linear-finite-element-space triangulation)
(let ((data data)
(triangulation triangulation))
(Triangulation-for-each-boundary-vertex
(lambda (v)
(with-access
v (Vertex index)
(f64vector-set! data index 0.)))
triangulation)
v))))
(define (Dirichlet-linear-finite-element-vector? v)
(and (Linear-finite-element-vector? v)
(let* ((fes (Finite-element-vector-space v))
(data (Finite-element-vector-data v))
(t (Finite-element-space-triangulation fes)))
(call-with-current-continuation
(lambda (k)
(Triangulation-for-each-boundary-vertex
(lambda (v)
(with-access
v (Vertex index)
(or (flzero? (f64vector-ref data index))
(k #f))))
t)
#t)))))
(define-handy-method (gauss-seidel-smoother! (m Dirichlet-linear-finite-element-operator))
(let* ((t (Finite-element-space-triangulation (Linear-finite-element-operator-range-space m)))
(vertices (Triangulation-vertices t))
(boundary-index? (lambda (i) (Boundary-vertex? (vector-ref vertices i)))) ; we're breaking the abstraction here, but let's fix it later.
(matrix-smoother! (Sparse-matrix-operator-gauss-seidel-smoother! matrix-implementation boundary-index?)))
(lambda (rhs #!optional (result (instantiate Linear-finite-element-vector space: (Linear-finite-element-operator-range-space m))))
(matrix-smoother! (Linear-finite-element-vector-data rhs)
(Linear-finite-element-vector-data result))
result)))
(define (CG-add-Dirichlet-operator-and-rhs-information! cg-data problem-data) ; builds the operator and the right-hand-side from problem-data
(with-access
cg-data (Conjugate-gradient-data triangulation space operator right-hand-side)
(with-access
problem-data (Problem-data a c b gamma f g)
(let ((m (->Linear-finite-element-operator space))
(rhs (instantiate Linear-finite-element-vector
space: space)))
(cond (a => (lambda (a) (add-stiffness-terms! m a))))
(cond (c => (lambda (c) (add-mass-terms! m c))))
(cond (b => (lambda (b) (add-transport-terms! m b))))
(cond (gamma => (lambda (gamma) (add-boundary-terms! m gamma))))
(cond (f => (lambda (f) (add-interior-integral-terms! rhs f))))
(cond (g => (lambda (g) (add-boundary-integral-terms! rhs g))))
(set! operator (->Dirichlet-linear-finite-element-operator m))
(set! right-hand-side (->Dirichlet-linear-finite-element-vector! rhs))))))