1 {- 2 convection term 3 4 XZ, 24/10/91 5 -} 6 7 {- 8 Modified to employ S_array. 9 The way in which the c_mat is constructed is also changed. 10 The new implementation is more efficient. 11 12 XZ, 7/2/92 13 -} 14 15 {- 16 Removed a bug 17 18 XZ, 25/2/92 19 -} 20 21 module C_matrix ( c_mat ) where 22 23 import Defs 24 import S_Array -- not needed w/ proper module handling 25 import Norm -- ditto 26 import Asb_routs 27 28 ----------------------------------------------------------- 29 -- Element convection matrix. -- 30 -- Used in "c_mat". -- 31 ----------------------------------------------------------- 32 33 cc_list :: 34 [ 35 ( 36 [Frac_type], 37 [Frac_type], 38 [Frac_type], 39 [(Frac_type,Frac_type)], 40 [(Frac_type,Frac_type)], 41 [(Frac_type,Frac_type)] 42 ) 43 ] 44 cc_list = 45 [ 46 select [a11,a12,a13,a14,a15,a16], 47 select [a12,a22,a23,a24,a25,a26], 48 select [a13,a23,a33,a34,a35,a36], 49 select [a14,a24,a34,a44,a45,a46], 50 select [a15,a25,a35,a45,a55,a56], 51 select [a16,a26,a36,a46,a56,a66] 52 ] 53 where 54 select = \tup_list -> 55 ( 56 map select1 tup_list, 57 map select2 tup_list, 58 map select3 tup_list, 59 map select4 tup_list, 60 map select5 tup_list, 61 map select6 tup_list 62 ) 63 select1 = \(x,_,_,_,_,_) -> x 64 select2 = \(_,x,_,_,_,_) -> x 65 select3 = \(_,_,x,_,_,_) -> x 66 select4 = \(_,_,_,x,_,_) -> x 67 select5 = \(_,_,_,_,x,_) -> x 68 select6 = \(_,_,_,_,_,x) -> x 69 70 c4 = 4 :: Frac_type 71 c9 = (-9) :: Frac_type 72 c11 = 11 :: Frac_type 73 c12 = 12 :: Frac_type 74 c16 = (-16) :: Frac_type 75 c18 = (-18) :: Frac_type 76 c20 = (-20) :: Frac_type 77 c24 = 24 :: Frac_type 78 c32 = (-32) :: Frac_type 79 c48 = (-48) :: Frac_type 80 c78 = 78 :: Frac_type 81 c80 = 80 :: Frac_type 82 c96 = (-96) :: Frac_type 83 c120 = 120 :: Frac_type 84 c128 = 128 :: Frac_type 85 c160 = 160 :: Frac_type 86 c192 = 192 :: Frac_type 87 c384 = 384 :: Frac_type 88 89 a11 = (c78,c18,c18,(c24,c24),(c24,c120),(c24,c120)) 90 91 a12 = (c9,c9,c11,(c4,c16),(c4,c16),(c16,c16)) 92 93 a13 = (c9,c11,c9,(c16,c4),(c16,c16),(c4,c16)) 94 95 a14 = (c12,c20,c20,(c48,c48),(c48,c16),(c48,c16)) 96 97 a15 = ((-c48),c16,c32,(c32,c16),(c32,(-c48)),(c16,(-c48))) 98 99 a16 = ((-c48),c32,c16,(c16,c32),(c16,(-c48)),(c32,(-c48))) 100 101 a22 = (c18,c78,c18,(c24,c120),(c24,c24),(c120,c24)) 102 103 a23 = (c11,c9,c9,(c16,c16),(c16,c4),(c16,c4)) 104 105 a24 = (c16,(-c48),c32,(c32,(-c48)),(c32,c16),((-c48),c16)) 106 107 a25 = (c20,c12,c20,(c48,c16),(c48,c48),(c16,c48)) 108 109 a26 = (c32,(-c48),c16,(c16,(-c48)),(c16,c32),((-c48),c32)) 110 111 a33 = (c18,c18,c78,(c120,c24),(c120,c24),(c24,c24)) 112 113 a34 = (c16,c32,(-c48),((-c48),c32),((-c48),c16),(c32,c16)) 114 115 a35 = (c32,c16,(-c48),((-c48),c16),((-c48),c32),(c16,c32)) 116 117 a36 = (c20,c20,c12,(c16,c48),(c16,c48),(c48,c48)) 118 119 a44 = (c96,c160,c160,(c384,c384),(c384,c128),(c384,c128)) 120 121 a45 = ((-c16),(-c16),c80,(c192,c128),(c192,c128),(c128,c128)) 122 123 a46 = ((-c16),c80,(-c16),(c128,c192),(c128,c128),(c192,c128)) 124 125 a55 = (c160,c96,c160,(c384,c128),(c384,c384),(c128,c384)) 126 127 a56 = (c80,(-c16),(-c16),(c128,c128),(c128,c192),(c128,c192)) 128 129 a66 = (c160,c160,c96,(c128,c384),(c128,c384),(c384,c384)) 130 131 ----------------------------------------------------------- 132 -- Converction term. -- 133 -- Used in "get_rh1". -- 134 -- Calls "cc_mat". -- 135 ----------------------------------------------------------- 136 137 c_mat 138 :: My_Array Int (((Frac_type,Frac_type,Frac_type), 139 (Frac_type,Frac_type,Frac_type)) 140 -> ([Frac_type],[Frac_type]) -> [Frac_type]) 141 142 c_mat = -- parameter: element_factors velocities 143 s_listArray bnds 144 ( map 145 ( 146 \ (col1,col2,col3,col4,col5,col6) -> 147 148 ( \ (gd1,gd2,gd3) u -> 149 map (\ col->list_match_prod col u) 150 [ 151 map ((*) gd1) col1, 152 map ((*) gd2) col2, 153 map ((*) gd3) col3, 154 map (f' gd2 gd3) col4, 155 map (f' gd1 gd3) col5, 156 map (f' gd1 gd2) col6 157 ] ) `bindTo` ( \ f -> 158 159 (\ (y1,y2) (u1,u2) -> zipWith (+) (f y1 u1) (f y2 u2)) ) 160 ) cc_list 161 ) 162 where 163 f' = \g1 g2 x -> (fst x) * g1 + (snd x) * g2 164 bnds = (1,v_nodel) 165 bindTo x k = k x -- essentially Haskell 1.0 "let"