;; The following routine builds a Linear-finite-element-space (i.e., a Finite-element-space
;; consisting of linear elements) from a triangulation. Each vertex in the triangulation
;; is associated with a basis element, and the neighbors of the element are the elements
;; associated with the vertices connected to the original vertex by (Vertex-edges v).
;; We use the same numbering scheme for the vertices as for the elements, and rely on that
;; a lot in the following code.
(define (Triangulation->Linear-finite-element-space t)
(let ((element-list
(map (lambda (v)
(instantiate Linear-element
index: (Vertex-index v)
neighbors: (Neighbors<-list
(cons (Vertex-index v) ; v is a neighbor to itelf
(map (lambda (e)
(Vertex-index
(if (eq? (Edge-vertex-1 e) v)
(Edge-vertex-2 e)
(Edge-vertex-1 e))))
(Vertex-edges->list (Vertex-edges v)))))))
(Triangulation-vertices->list (Triangulation-vertices t)))))
(instantiate Linear-finite-element-space
triangulation: t
neighbors: (Finite-element-space-neighbors<-list (map Element-neighbors element-list))
elements: (Finite-element-space-elements<-list element-list))))
;; We anticipate that Finite-element-operators will have matrix implementations
;; of one type or another (Sparse-matrix-operators, most likely).
(define-class Finite-element-operator Operator-with-matrix-implementation
((= offsets immutable:) ; The offset into indices and coefficients for each row
(= indices immutable:) ; a vector of the indices of ...
(= coefficients maybe-uninitialized: immutable:) ; nonzero-coefficients in the matrix representation
(= domain-space immutable:) ; (pointers to) the domain ...
(= range-space immutable:))) ; and range finite-element-spaces
(define-class Linear-finite-element-operator Finite-element-operator ())
(define-handy-method (initialize! (o Linear-finite-element-operator))
;; initialize various fields that depend on the space
;; if coefficients is not defined, we set it to arrays of floating-point
;; zeros associated with the associated indices field (which must be defined).
(if (not (field-defined? o 'coefficients))
(initialize-field-value!
o
(make-f64vector (indices-length indices) 0.0)
'coefficients))
;; these are easy, so we do them here
(if (not (field-defined? o 'domain?))
(initialize-field-value!
o
(lambda (v)
(and (Finite-element-vector? v)
(eq? domain-space (Finite-element-vector-space v))))
'domain?))
(if (not (field-defined? o 'range?))
(initialize-field-value!
o
(lambda (v)
(and (Finite-element-vector? v)
(eq? range-space (Finite-element-vector-space v))))
'range?))
;; Here both indices and coefficients are initialized, so we can rely
;; on them since they're immutable.
(if (not (field-defined? o 'matrix-implementation))
(initialize-field-value!
o
(instantiate Sparse-matrix-operator
domain-dimension: (Finite-element-space-neighbors-length domain-space)
range-dimension: (Finite-element-space-neighbors-length range-space)
indices: indices
coefficients: coefficients
offsets: offsets)
'matrix-implementation))
;; initialize the mapping field to multiplication by the sparse matrix
(initialize-field-value!
o
(lambda (v)
(instantiate Linear-finite-element-vector
data: (Operator-apply matrix-implementation (Finite-element-vector-data v))
space: range-space))
'mapping)
;; call superclass's initialize! method
(call-next-method))
(define-handy-method (->Linear-finite-element-operator (space Linear-finite-element-space))
(let ((neighbors neighbors))
(let ((indices (list->indices (mappend1 (lambda (x)
(vector->list x))
(vector->list neighbors))))
(offsets (let* ((n (fx+ (vector-length neighbors) 1))
(result (make-vector n)))
(vector-set! result 0 0)
(do ((i 1 (fx+ i 1)))
((fx= i n) result)
(vector-set! result i (fx+ (vector-ref result (fx- i 1))
(vector-length (vector-ref neighbors (fx- i 1)))))))))
(instantiate Linear-finite-element-operator
indices: indices
offsets: offsets
domain-space: space
range-space: space))))
(define-handy-method (check-same-class (m Linear-finite-element-operator) n)
(and (call-next-method)
(eq? domain-space
(Linear-finite-element-operator-domain-space n))
(eq? range-space
(Linear-finite-element-operator-range-space n))))
(define (wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space)
(instantiate Linear-finite-element-operator
indices: (Sparse-matrix-operator-indices matrix)
coefficients: (Sparse-matrix-operator-coefficients matrix)
offsets: (Sparse-matrix-operator-offsets matrix)
matrix-implementation: matrix
domain-space: domain-space
range-space: range-space))
(define-handy-method (Vector-add (m Linear-finite-element-operator) n)
(if (check-same-class m n)
(let ((matrix (Vector-add matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
(error "Vector-add: Linear-finite-element-operators are not the same class: " m n)))
(define-handy-method (Vector-subtract (m Linear-finite-element-operator) n)
(if (check-same-class m n)
(let ((matrix (Vector-subtract matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
(error "Vector-subtract: Linear-finite-element-operators are not the same class: " m n)))
(define-handy-method (Vector-scale a (m Linear-finite-element-operator))
(let ((matrix (Vector-scale a matrix-implementation)))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space)))
(define-handy-method (Vector-*axpy a (m Linear-finite-element-operator) n)
(if (check-same-class m n)
(let ((matrix (Vector-*axpy a matrix-implementation (Linear-finite-element-operator-matrix-implementation n))))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
(error "Vector-*axpy: Linear-finite-element-operators are not the same class: " m n)))
(define-handy-method (Vector-Linf-norm (m Operator-with-matrix-implementation))
(Vector-Linf-norm matrix-implementation))
;; we compose two linear-finite-element-operators by composing their matrix implementations
;; and building a new operator with instantiate. The immutability of the fields of the
;; Finite-element-operator gets in the way here, since we may like to define this method
;; on Finite-element-operators and clone one of the arguments to make it work with subclasses
;; too, but we'll let future generations worry about this ;-).
(define-method (Operator-compose (op1 Linear-finite-element-operator) (op2 Linear-finite-element-operator op2))
(if (eq? (Linear-finite-element-operator-range-space op2)
(Linear-finite-element-operator-domain-space op1))
(let ((new-implementation
(Operator-compose (Linear-finite-element-operator-matrix-implementation op1)
(Linear-finite-element-operator-matrix-implementation op2))))
(wrap-matrix-implementation-in-linear-finite-element-operator new-implementation
(Linear-finite-element-operator-domain-space op2)
(Linear-finite-element-operator-range-space op1)))
(error "Operator-compose: Linear-finite-element-operators do not have compatable domains and ranges.")))
(define-handy-method (Operator->identity (m Linear-finite-element-operator))
(if (eq? range-space
domain-space)
(let ((matrix (Operator->identity matrix-implementation)))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix domain-space range-space))
(error "Operator->identity: Linear-finite-element-operator does not have same domain and range." m)))
;; Takes a Linear-finite-element-operator and returns the adjoint operator.
(define-handy-method (Operator-adjoint (m Linear-finite-element-operator))
(let ((matrix (Operator-adjoint matrix-implementation)))
(wrap-matrix-implementation-in-linear-finite-element-operator matrix range-space domain-space)))
;; calculate the gradient of the linear function that is 1.0 at p1 and 0.0 at p2 and p3
(define (grad-phi p1 p2 p3)
(declare (inlining-limit 10000)) ; I want this inlined
(let ((v (Point-subtract p1 p2))
(u (Point-subtract p3 p2)))
(let ((w (Point-subtract v (Point-scale (fl/ (Point-dot-product v u)
(Point-dot-product u u))
u))))
(Point-scale (fl/ (Point-dot-product w w)) w))))
;;; We'll work in barycentric coordinates in triangles.
;;; We'll define an interface for barycentric coordinates
(define make-barycentric f64vector)
(define barycentric-1/2-1/2-0 (make-barycentric 0.5 0.5 0.0))
(define barycentric-1/2-0-1/2 (make-barycentric 0.5 0.0 0.5))
(define barycentric-0-1/2-1/2 (make-barycentric 0.0 0.5 0.5))
(define barycentric-coordinate f64vector-ref)
(define barycentric-dimension f64vector-length)
(define-generic (evaluate-barycentric f coordinates (simplex)))
(define-method (evaluate-barycentric f coordinates (triangle Triangle))
(if (not (fx= (barycentric-dimension coordinates) 3))
(error "evaluate-barycentric: Need three barycentric coordinates for a triangle: " f coordinates triangle)
(evaluate-barycentric-triangle f coordinates triangle)))
(define-generic (evaluate-barycentric-triangle (f) coordinates triangle)
(cond ((number? f)
f)
((procedure? f)
(with-access
triangle (Triangle vertex-1 vertex-2 vertex-3)
(f (Point-*axpy (barycentric-coordinate coordinates 0)
(Vertex-point vertex-1)
(Point-*axpy (barycentric-coordinate coordinates 1)
(Vertex-point vertex-2)
(Point-scale (barycentric-coordinate coordinates 2)
(Vertex-point vertex-3)))))))
(else
(error "evaluate-barycentric-triangle: No method for this object: " f coordinates triangle))))
(define-method (evaluate-barycentric-triangle (f Linear-finite-element-vector) coordinates triangle)
(with-access
f (Linear-finite-element-vector data)
(with-access
triangle (Triangle vertex-1 vertex-2 vertex-3)
(fl+ (fl* (barycentric-coordinate coordinates 0) (f64vector-ref data (Vertex-index vertex-1)))
(fl* (barycentric-coordinate coordinates 1) (f64vector-ref data (Vertex-index vertex-2)))
(fl* (barycentric-coordinate coordinates 2) (f64vector-ref data (Vertex-index vertex-3)))))))
(define-method (evaluate-barycentric f coordinates (edge Edge))
(if (not (fx= (barycentric-dimension coordinates) 2))
(error "evaluate-barycentric: Need three barycentric coordinates for an edge: " f coordinates edge)
(evaluate-barycentric-edge f coordinates edge)))
(define-generic (evaluate-barycentric-edge (f) coordinates edge)
(cond ((number? f)
f)
((procedure? f)
(with-access
edge (Edge vertex-1 vertex-2)
(f (Point-*axpy (barycentric-coordinate coordinates 0)
(Vertex-point vertex-1)
(Point-scale (barycentric-coordinate coordinates 1)
(Vertex-point vertex-2))))))
(else
(error "evaluate-barycentric-edge: No method for this object: " f coordinates edge))))
(define-method (evaluate-barycentric-edge (f Linear-finite-element-vector) coordinates edge)
(with-access
f (Linear-finite-element-vector data)
(with-access
edge (Edge vertex-1 vertex-2)
(fl+ (fl* (barycentric-coordinate coordinates 0) (f64vector-ref data (Vertex-index vertex-1)))
(fl* (barycentric-coordinate coordinates 1) (f64vector-ref data (Vertex-index vertex-2)))))))
;; Evaluate the function a at the midpoint of the Edge e. A method for
;; curved-edges will come later.
(define-generic (Edge-midpoint (e)))
(define-handy-method (Edge-midpoint (e Edge))
(declare (inlining-limit 10000)) ; I want this inlined
(Point-scale 0.5 (Point-add (Vertex-point vertex-1)
(Vertex-point vertex-2))))
;; OK, the following functions should be multi-methods based not only on the class
;; of the triangle, but also on the type of the space. And the f64vector of as
;; should be encapsulated, too. But at least this allows us to have a method for
;; Triangle-with-curved-side later.
;; as is an f64vector of the values of a evaluated at the midpoint of the edges
;; in a triangulation.
(define-generic (Triangle-stiffness (t Triangle) a m)
;; integrate a \grad \phi_j\cdot\grad\phi_i on t, given a's at midpoints of edges
(with-access
t (Triangle vertex-1 vertex-2 vertex-3)
(with-access
vertex-1 (Vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Vertex point index)
(let ((point-2 point)
(index-2 index))
(with-access
vertex-3 (Vertex point index)
(let ((point-3 point)
(index-3 index))
;; here we rely crucially on the elements having the same numbering as the
;; vertices in the Triangulation. So grad-1 is the gradient on the triangle
;; of the element that is 1 at vertex-1 and zero at the other vertices, etc.
(let ((grad-1 (grad-phi point-1 point-2 point-3))
(grad-2 (grad-phi point-2 point-3 point-1))
(grad-3 (grad-phi point-3 point-1 point-2))
;; We're going to use the quadrature rule that the integral of a function
;; f on a triangle is approximately equal to the average of f at the midpoints
;; of the three sides of the triangle times the area of the triangle. This
;; is exact for quadratic functions f.
;; Here, the function we're integrating is the dot-product of the gradients
;; of elements, which is constant on the triangle, times a, so we precompute
;; the term involving a.
(weight (fl* (fl+ (evaluate-barycentric a barycentric-1/2-1/2-0 t)
(evaluate-barycentric a barycentric-1/2-0-1/2 t)
(evaluate-barycentric a barycentric-0-1/2-1/2 t))
0.333333333333333333333333333333333333
(Triangle-area t))))
(Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (Point-dot-product grad-1 grad-1)))
(Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (Point-dot-product grad-2 grad-2)))
(Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight (Point-dot-product grad-3 grad-3)))
(let ((val (fl* weight (Point-dot-product grad-1 grad-2))))
(Sparse-matrix-operator-increment-entry m index-1 index-2 val)
(Sparse-matrix-operator-increment-entry m index-2 index-1 val))
(let ((val (fl* weight (Point-dot-product grad-1 grad-3))))
(Sparse-matrix-operator-increment-entry m index-1 index-3 val)
(Sparse-matrix-operator-increment-entry m index-3 index-1 val))
(let ((val (fl* weight (Point-dot-product grad-2 grad-3))))
(Sparse-matrix-operator-increment-entry m index-2 index-3 val)
(Sparse-matrix-operator-increment-entry m index-3 index-2 val)))))))))))
;; Add the stiffness terms, the integral of a times \grad \phi_i\cdot\grad\phi_j to the Linear-finite-element
;; matrix m.
(define (add-stiffness-terms! m a)
(with-access
m (Linear-finite-element-operator matrix-implementation domain-space)
(with-access
domain-space (Finite-element-space triangulation)
(let ((m matrix-implementation))
;; calculate the contribution on each triangle to the terms involving
;; the elements that overlap that triangle nontrivially.
(Triangulation-for-each-triangle
(lambda (t)
(Triangle-stiffness t a m))
triangulation)))))
;; The same principle as Triangle-stiffness, but here we're working with
;; vectors of Points (bs).
(define-generic (Triangle-transport (t Triangle) b m)
;; integrate b \cdot \grad\phi_j \phi_i on t given b values
(with-access
t (Triangle vertex-1 vertex-2 vertex-3)
(with-access
vertex-1 (Vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Vertex point index)
(let ((point-2 point)
(index-2 index))
(with-access
vertex-3 (Vertex point index)
(let ((point-3 point)
(index-3 index))
(let ((grad-1 (grad-phi point-1 point-2 point-3))
(grad-2 (grad-phi point-2 point-3 point-1))
(grad-3 (grad-phi point-3 point-1 point-2))
(edge-1-b (evaluate-barycentric b barycentric-1/2-1/2-0 t))
(edge-2-b (evaluate-barycentric b barycentric-1/2-0-1/2 t))
(edge-3-b (evaluate-barycentric b barycentric-0-1/2-1/2 t))
(weight (fl* (Triangle-area t)
;; average of 3 points time 1/2, phi_i at midpoint of each edge
0.166666666666666666666666)))
(Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight
(fl+ (Point-dot-product grad-1 edge-3-b)
(Point-dot-product grad-1 edge-1-b))))
(Sparse-matrix-operator-increment-entry m index-2 index-1 (fl* weight
(fl+ (Point-dot-product grad-1 edge-1-b)
(Point-dot-product grad-1 edge-2-b))))
(Sparse-matrix-operator-increment-entry m index-3 index-1 (fl* weight
(fl+ (Point-dot-product grad-1 edge-2-b)
(Point-dot-product grad-1 edge-3-b))))
(Sparse-matrix-operator-increment-entry m index-1 index-2 (fl* weight
(fl+ (Point-dot-product grad-2 edge-3-b)
(Point-dot-product grad-2 edge-1-b))))
(Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight
(fl+ (Point-dot-product grad-2 edge-1-b)
(Point-dot-product grad-2 edge-2-b))))
(Sparse-matrix-operator-increment-entry m index-3 index-2 (fl* weight
(fl+ (Point-dot-product grad-2 edge-2-b)
(Point-dot-product grad-2 edge-3-b))))
(Sparse-matrix-operator-increment-entry m index-1 index-3 (fl* weight
(fl+ (Point-dot-product grad-3 edge-3-b)
(Point-dot-product grad-3 edge-1-b))))
(Sparse-matrix-operator-increment-entry m index-2 index-3 (fl* weight
(fl+ (Point-dot-product grad-3 edge-1-b)
(Point-dot-product grad-3 edge-2-b))))
(Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight
(fl+ (Point-dot-product grad-3 edge-2-b)
(Point-dot-product grad-3 edge-3-b)))))))))))))
(define (add-transport-terms! m b)
(with-access
m (Linear-finite-element-operator matrix-implementation domain-space)
(with-access
domain-space (Finite-element-space triangulation)
(let ((m matrix-implementation))
(Triangulation-for-each-triangle
(lambda (t)
(Triangle-transport t b ))
triangulation)))))
(define-generic (Triangle-mass (t Triangle) c m)
;; integrates c * \phi_i * \phi_j on t given c's on edge midpoints.
(with-access
t (Triangle vertex-1 vertex-2 vertex-3)
(with-access
vertex-1 (Vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Vertex point index)
(let ((point-2 point)
(index-2 index))
(with-access
vertex-3 (Vertex point index)
(let ((point-3 point)
(index-3 index))
(let ((edge-1-c (evaluate-barycentric c barycentric-1/2-1/2-0 t))
(edge-2-c (evaluate-barycentric c barycentric-1/2-0-1/2 t))
(edge-3-c (evaluate-barycentric c barycentric-0-1/2-1/2 t)))
;; the value of \phi_i\phi_j at the midpoint of an edge is 1/4 or zero;
;; here we precalculate the weight 1/4 * 1/3 * area of t
(let ((weight (fl* 0.08333333333333333333333333333333333333
(Triangle-area t))))
(Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (fl+ edge-1-c edge-2-c)))
(Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (fl+ edge-2-c edge-3-c)))
(Sparse-matrix-operator-increment-entry m index-3 index-3 (fl* weight (fl+ edge-3-c edge-2-c)))
(let ((val (fl* weight edge-1-c)))
(Sparse-matrix-operator-increment-entry m index-1 index-2 val)
(Sparse-matrix-operator-increment-entry m index-2 index-1 val))
(let ((val (fl* weight edge-3-c)))
(Sparse-matrix-operator-increment-entry m index-1 index-3 val)
(Sparse-matrix-operator-increment-entry m index-3 index-2 val))
(let ((val (fl* weight edge-2-c)))
(Sparse-matrix-operator-increment-entry m index-2 index-3 val)
(Sparse-matrix-operator-increment-entry m index-3 index-1 val))))))))))))
(define (add-mass-terms! m c)
(with-access
m (Linear-finite-element-operator matrix-implementation domain-space)
(with-access
domain-space (Finite-element-space triangulation)
(let ((m matrix-implementation))
(Triangulation-for-each-triangle
(lambda (t)
(Triangle-mass t c m))
triangulation)))))
;; For integration on the boundary, we use the two-point Guass-rule on [0,1],
;; which has the advantage of not evaluating functions at vertices on the boundary,
;; where, e.g., there may be a discontinuity in the normal derivative of u
;; because the normal to the boundary points in different directions on each side
;; of a vertex.
(define left-gauss-point .21132486540518708) ; (* 0.5 (- 1.0 (/ (sqrt 3.0))))
(define right-gauss-point .7886751345948129) ; (* 0.5 (+ 1.0 (/ (sqrt 3.0))))
(define-generic (Edge-boundary-term (e Edge) gamma m)
;; integrates gamma * phi_i * phi_j on an edge, increment the correct coefficient
(with-access
e (Edge vertex-1 vertex-2 index)
(with-access
vertex-1 (Boundary-vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Boundary-vertex point index)
(let ((point-2 point)
(index-2 index))
;; We're considering point-1 to be the left point and point-2 to be the right (it doesn't matter).
;; So the gauss point near point-1 is right-gauss-point * point-1 + left-gauss-point * point-2.
(let ((val-1 (gamma (Point-add (Point-scale right-gauss-point point-1)
(Point-scale left-gauss-point point-2))))
(val-2 (gamma (Point-add (Point-scale left-gauss-point point-1)
(Point-scale right-gauss-point point-2))))
(weight (fl* (Edge-length e) 0.5)))
;; The basis element that is 1 at point-1 takes the value right-gauss-point at left-gauss-point, and
;; the value left-gauss-point at right-gauss-point. (Yes, I know it's confusing.)
(Sparse-matrix-operator-increment-entry m index-1 index-1 (fl* weight (fl+ (fl* right-gauss-point right-gauss-point val-1)
(fl* left-gauss-point left-gauss-point val-2))))
(Sparse-matrix-operator-increment-entry m index-2 index-2 (fl* weight (fl+ (fl* right-gauss-point right-gauss-point val-2)
(fl* left-gauss-point left-gauss-point val-1))))
(Sparse-matrix-operator-increment-entry m index-1 index-2 (fl* weight left-gauss-point right-gauss-point (fl+ val-1 val-2)))
(Sparse-matrix-operator-increment-entry m index-2 index-1 (fl* weight left-gauss-point right-gauss-point (fl+ val-1 val-2))))))))))
(define (add-boundary-terms! m gamma)
(with-access
m (Linear-finite-element-operator matrix-implementation domain-space)
(let ((m matrix-implementation))
(with-access
domain-space (Finite-element-space triangulation)
(let ((triangulation triangulation)
(n (Triangulation-boundary-vertices-length triangulation)))
;; go through the boundary edges, calculate the length of each edge, and then increment
;; the contribution of that edge to the integrals
(Triangulation-for-each-boundary-edge
(lambda (e)
(Edge-boundary-term e gamma m))
triangulation))))))
;; The functions to calculate the right-hand-side of the matrix problems.
(define-generic (Triangle-interior-integral (t Triangle) f data)
;; integrate f \phi_i on t, given f values at edge midpoints.
(with-access
t (Triangle vertex-1 vertex-2 vertex-3 edge-1 edge-2 edge-3)
(with-access
vertex-1 (Vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Vertex point index)
(let ((point-2 point)
(index-2 index))
(with-access
vertex-3 (Vertex point index)
(let ((point-3 point)
(index-3 index))
(let ((edge-1-f (evaluate-barycentric f barycentric-1/2-1/2-0 t))
(edge-2-f (evaluate-barycentric f barycentric-1/2-0-1/2 t))
(edge-3-f (evaluate-barycentric f barycentric-0-1/2-1/2 t)))
(let ((weight (fl* 0.1666666666666666666666666666666666666
(Triangle-area t))))
(f64vector-set! data index-1 (fl+ (f64vector-ref data index-1)
(fl* weight (fl+ edge-2-f edge-1-f))))
(f64vector-set! data index-2 (fl+ (f64vector-ref data index-2)
(fl* weight (fl+ edge-1-f edge-3-f))))
(f64vector-set! data index-3 (fl+ (f64vector-ref data index-3)
(fl* weight (fl+ edge-2-f edge-3-f))))))))))))))
(define (add-interior-integral-terms! v f)
(with-access
v (Linear-finite-element-vector data
space)
(with-access
space (Finite-element-space triangulation)
(let ((data data)
(triangulation triangulation))
(Triangulation-for-each-triangle
(lambda (t)
(Triangle-interior-integral t f data))
triangulation)))))
;; Again, here we use the two-point Gauss rule.
(define-generic (Edge-boundary-integral (e Edge) g data)
;; integrate g * \phi_i on the edge e.
(with-access
e (Edge vertex-1 vertex-2 index)
(with-access
vertex-1 (Vertex point index)
(let ((point-1 point)
(index-1 index))
(with-access
vertex-2 (Vertex point index)
(let ((point-2 point)
(index-2 index))
(let ((val-1 (g (Point-add (Point-scale right-gauss-point point-1)
(Point-scale left-gauss-point point-2))))
(val-2 (g (Point-add (Point-scale left-gauss-point point-1)
(Point-scale right-gauss-point point-2))))
(weight (fl* (Edge-length e) 0.5)))
(f64vector-set!
data
index-1
(fl+ (f64vector-ref data index-1)
(fl* weight (fl+ (fl* right-gauss-point val-1) ; this element has the value right-gauss-point at the Gauss point near point-1
(fl* left-gauss-point val-2))))) ; and the value left-gauss-point at the Gauss point near point-2
(f64vector-set!
data
index-2
(fl+ (f64vector-ref data index-2)
(fl* weight (fl+ (fl* left-gauss-point val-1) ; this element has the value left-gauss-point at the Gauss point near point-1
(fl* right-gauss-point val-2)))))))))))) ; and the value right-gauss-point at the Gauss point near point-2
(define (add-boundary-integral-terms! v g)
(with-access
v (Linear-finite-element-vector space data)
(with-access
space (Finite-element-space triangulation)
(let ((triangulation triangulation)
(data data))
;; go through the boundary edges, calculate the length of each edge, and then increment
;; the contribution of that edge to the integrals
(Triangulation-for-each-boundary-edge
(lambda (e)
(Edge-boundary-integral e g data))
triangulation)))))