Fix bugs in big_K, big_M.

diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
index ae7548a86d75d42311691c8bdcfa780384b3938c..3e6869557d66f5601c769ff3863448a666729d60 100644 (file)
--- a/src/FEM/R1.hs
+++ b/src/FEM/R1.hs
@@ -370,7 +370,7 @@ big_K_elem :: forall m n l a b.
       -> b
       -> Mat l l a
 big_K_elem pde params _ k cur_K _ =
-  ifoldl2 accum cur_K (big_N's_matrix :: Mat m (S n) (a -> a))
+  ifoldl2 accum cur_K (big_N's_matrix :: Mat (S n) (S n) (a -> a))
   where
     accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
     accum i j prev_K these_N's =
@@ -380,8 +380,11 @@ big_K_elem pde params _ k cur_K _ =
         (x1,x2) = (mesh params) !!! (k,0)
         q = affine_inv (x1,x2)
         integrand x = ((big_A pde) (q x)) * (these_N's x)
-        -- The pointer matrix numbers from 1 so subtract one here to
-        -- get the right index.
+        -- The pointer matrix numbers from 1 so subtract one below to
+        -- get the right index. The indices i,j have upper bounds
+        -- dependent on the element k. Since we statically create the
+        -- matrix of basis function derivatives, we have to check here
+        -- whether or not i,j exceed the max index.
         row_idx = ((pointer params) !!! (k,i)) - 1
         col_idx = ((pointer params) !!! (k,j)) - 1
         integral = (two/(x2 - x1))* (gaussian integrand)
@@ -436,7 +439,7 @@ big_M_elem :: forall m n l a b.
       -> b
       -> Mat l l a
 big_M_elem pde params _ k cur_M _ =
-  ifoldl2 accum cur_M (big_Ns_matrix :: Mat m (S n) (a -> a))
+  ifoldl2 accum cur_M (big_Ns_matrix :: Mat (S n) (S n) (a -> a))
   where
     accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
     accum i j prev_M these_Ns =