From: Michael Orlitzky Date: Sun, 13 Apr 2014 08:12:52 +0000 (-0400) Subject: Fix bugs in big_K, big_M. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=7ac2f5fdd7db6c998949af5efdb45501dfb051da;p=numerical-analysis.git Fix bugs in big_K, big_M. --- diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs index ae7548a..3e68695 100644 --- 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 =