Commit whatever I was working on last and have forgotten about.

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index 57e2412c83af938f736200bee0ce2b83299cff5f..217abd3009acaac271e2a17fe99eccd6e5df9050 100644 (file)
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -632,3 +632,61 @@ solution pde params =
             xi = (affine interval) x
             combine ci ni = ci*(ni xi)
 
+
+
+
+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)
+         => Col (S l) a
+         -> PDE a
+         -> Params m n (S l) a
+         -> Piecewise a
+solution' global_coeffs 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)
+