Finish the FEM.R1 solution for Dirichlet boundary conditions.

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index 922469f730dff4af6ed9d8387928c60335acd223..57e2412c83af938f736200bee0ce2b83299cff5f 100644 (file)
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -43,11 +43,19 @@ import Linear.Matrix (
   Row,
   (!!!),
   construct,
+  element_sum2,
   fromList,
   ifoldl2,
+  map2,
   nrows,
-  set_idx )
+  rows2,
+  set_idx,
+  toList,
+  transpose,
+  zip2,
+  zipwith2 )
 import Linear.System ( solve_positive_definite )
+import Piecewise ( Piecewise(..), from_intervals )
 import Polynomials.Orthogonal ( legendre )
 
 -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
@@ -190,7 +198,7 @@ affine_inv (x1,x2) x =
 -- * Load vector
 
 -- | Normalized integrals of orthogonal basis functions over
---   n[-1,1]. The test case below comes from Sage where the
+--   [-1,1]. The test case below comes from Sage where the
 --   orthogonality of the polynomials' derivatives can easily be
 --   tested.
 --
@@ -568,3 +576,59 @@ coefficients pde params =
   where
     matrix = (big_K pde params) + (big_M pde params)
     b      = big_F pde params
+
+
+solution :: forall m n l a.
+            (Arity m, Arity n, Arity l,
+            Algebraic.C a, Eq a, RealField.C a, ToRational.C a, Show a)
+         => PDE a
+         -> Params m n (S l) a
+         -> Piecewise a
+solution pde params =
+  from_intervals $ map head $ toList $ solved_column
+  where
+    global_coeffs :: Col (S l) a
+    global_coeffs = coefficients pde params
+
+    ptr :: Mat m (S n) Int
+    ptr = pointer params
+
+    -- Each mesh element has an associated row in the pointer
+    -- matrix. Stick them together.
+    mesh_with_ptr_rows :: Col m (Interval a, Row (S n) Int)
+    mesh_with_ptr_rows = zip2 (mesh params) (rows2 ptr)
+
+    make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
+    make_local_coeffs (interval, ptr_row) =
+      construct lambda
+        where
+          lambda _ j = if (ptr_row !!! (0,j)) == zero
+                       then zero
+                       else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
+
+    -- Create a column vector for each mesh element containing the global
+    -- coefficients corresponding to that element.
+    local_coeffs :: Col m (Row (S n) a)
+    local_coeffs = map2 make_local_coeffs mesh_with_ptr_rows
+
+    global_basis_functions :: Col (S n) (a -> a)
+    global_basis_functions =
+      construct lambda
+        where lambda i _ = big_N (toInteger i)
+
+    mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
+    mesh_with_coeffs = zip2 (mesh params) local_coeffs
+
+    solved_column :: Col m (Interval a, (a -> a))
+    solved_column = map2 solve_piece $ mesh_with_coeffs
+
+    solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
+    solve_piece (interval, coeffs_row) = (interval, f)
+      where
+        coeffs_col = transpose coeffs_row
+
+        f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
+          where
+            xi = (affine interval) x
+            combine ci ni = ci*(ni xi)
+