1 -- 2 -- Patricia Fasel 3 -- Los Alamos National Laboratory 4 -- 1990 August 5 -- 6 module Utils (applyOpToMesh, coarseMesh, fineMesh, genRand, log2) where 7 8 import PicType 9 import Array -- 1.3 10 infix 1 =: 11 (=:) a b = (a,b) 12 13 -- apply the given operator to a mesh of given size 14 -- operator is applied to the position and to the 8 surrounding positions 15 -- so value(i,j) is decided by (i-1,j-1), (i,j-1), (i+1,j-1), etc. 16 -- corners and edges are handled as if the mesh was a torus 17 18 applyOpToMesh :: (Mesh -> Range -> Value) -> Mesh -> Indx -> Mesh 19 applyOpToMesh operator mesh n' = 20 array ((0,0), (n,n)) 21 22 (-- corners (nw, ne, sw, se) 23 [(0,0) =: operator mesh [n, 0, 1, n, 0, 1]] ++ 24 [(0,n) =: operator mesh [n, 0, 1, n1, n, 0]] ++ 25 [(n,0) =: operator mesh [n1, n, 0, n, 0, 1]] ++ 26 [(n,n) =: operator mesh [n1, n, 0, n1, n, 0]] ++ 27 28 -- edges (top row, bottom row, left column, right column) 29 [(0,j) =: operator mesh [n, 0, n1, (j-1), j, (j+1)] | j<-[1..n1]] ++ 30 [(n,j) =: operator mesh [n1, n, 0, (j-1), j, (j+1)] | j<-[1..n1]] ++ 31 [(i,0) =: operator mesh [(i-1), i, (i+1), n, 0, 1] | i<-[1..n1]] ++ 32 [(i,n) =: operator mesh [(i-1), i, (i+1), n1, n, 0] | i<-[1..n1]] ++ 33 34 -- internal 35 [(i,j) =: operator mesh [(i-1), i, (i+1), (j-1), j, (j+1)] 36 | i <- [1..n1], j <- [1..n1]]) 37 where 38 n = n'-1 39 n1 = n'-2 40 41 42 -- project a mesh onto a mesh of half the rank 43 -- a b c d e f a c e 44 -- g h i j k l => m o q 45 -- m n o p q r y 0 2 46 -- s t u v w x 47 -- y z 0 1 2 3 48 -- 4 5 6 7 8 9 49 50 coarseMesh :: Mesh -> Indx -> Mesh 51 coarseMesh mesh n = 52 array ((0,0), (nHalf,nHalf)) 53 [(i,j) =: mesh!(i*2, j*2) | i <- [0..nHalf], j <- [0..nHalf]] 54 where 55 nHalf = n `div` 2 - 1 56 57 58 -- interpolate a mesh of half rank onto a full mesh 59 -- values aren't just copied but are a function of the letter shown 60 -- a b c d e f A B C 61 -- g h i j k l <= D E F 62 -- m n o p q r G H I 63 -- s t u v w x 64 -- y z 0 1 2 3 65 -- 4 5 6 7 8 9 66 -- 67 -- a = A, c = B, e = C, m = D, o = E, q = F, y = G, 0 = H, 2 = I 68 -- g = .5(A+D), i = .5(B+E), k = .5(C+F), b = .5(A+B), d = .5(B+C), f = .5(C+A) 69 -- s = .5(D+G), u = .5(E+H), w = .5(F+I), n = .5(D+E), p = .5(E+F), r = .5(F+D) 70 -- 4 = .5(G+A), 6 = .5(H+B), 8 = .5(I+C), z = .5(G+H), 1 = .5(H+I), 3 = .5(I+G) 71 -- h = .25(A+B+D+E), j = .25(B+C+E+F), l = .25(C+A+F+D) ... 72 73 fineMesh :: Mesh -> Indx -> Mesh 74 fineMesh mesh nHalf' = 75 array ((0,0), (n,n)) 76 77 (-- corners (nw, ne, sw, se) 78 [(0,0) =: 3] ++ 79 [(0,n) =: 3] ++ 80 [(n,0) =: 3] ++ 81 [(n,n) =: 3] ++ 82 83 -- edges (north, south) 84 [(0,2*j) =: 4| j<-[1..nHalf]] ++ 85 [(0,2*j-1) =: 4| j<-[1..nHalf]] ++ 86 [(n,2*j) =: 4| j<-[1..nHalf]] ++ 87 [(n,2*j-1) =: 4| j<-[1..nHalf]] ++ 88 89 -- edges (west, east) 90 [(2*i,0) =: 5| i<-[1..nHalf]] ++ 91 [(2*i-1,0) =: 5| i<-[1..nHalf]] ++ 92 [(2*i,n) =: 5| i<-[1..nHalf]] ++ 93 [(2*i-1,n) =: 5| i<-[1..nHalf]] ++ 94 95 -- interior 96 [(2*i,2*j) =: 6| i<-[1..nHalf], j<-[1..nHalf]] ++ 97 [(2*i,2*j-1) =: 6| i<-[1..nHalf], j<-[1..nHalf]] ++ 98 [(2*i-1,2*j) =: 6| i<-[1..nHalf], j<-[1..nHalf]] ++ 99 [(2*i-1,2*j-1) =: 6| i<-[1..nHalf], j<-[1..nHalf]]) 100 {- 101 array ((0,0), (n,n)) 102 103 (-- corners (nw, ne, sw, se) 104 [(0,0) =: mesh!(0,0)] ++ 105 [(0,n) =: 0.5*(mesh!(0,0) + mesh!(0,nHalf))] ++ 106 [(n,0) =: 0.5*(mesh!(0,0) + mesh!(nHalf,0))] ++ 107 [(n,n) =: 0.25*(mesh!(0,0) + mesh!(0,nHalf) + mesh!(nHalf,0) + 108 mesh!(nHalf,nHalf))] ++ 109 110 -- edges (north, south) 111 [(0,2*j) =: mesh!(0,j)| j<-[1..nHalf]] ++ 112 [(0,2*j-1) =: 0.5*(mesh!(0,j) + mesh!(0,j-1)) | j<-[1..nHalf]] ++ 113 [(n,2*j) =: 0.5*(mesh!(0,j) + mesh!(nHalf,j)) | j<-[1..nHalf]] ++ 114 [(n,2*j-1) =: 0.25*(mesh!(0,j) + mesh!(0,j-1) + mesh!(nHalf,j) + 115 mesh!(nHalf,j-1)) | j<-[1..nHalf]] ++ 116 117 -- edges (west, east) 118 [(2*i,0) =: mesh!(i,0) | i<-[1..nHalf]] ++ 119 [(2*i-1,0) =: 0.5*(mesh!(i,0) + mesh!(i,nHalf)) | i<-[1..nHalf]] ++ 120 [(2*i,n) =: 0.5*(mesh!(i,0) + mesh!(i,nHalf)) | i<-[1..nHalf]] ++ 121 [(2*i-1,n) =: 0.25*(mesh!(i,0) + mesh!(i,nHalf) + mesh!(i-1,0) + 122 mesh!(i-1,nHalf)) | i<-[1..nHalf]] ++ 123 124 -- interior 125 [(2*i,2*j) =: mesh!(i,j) | i<-[1..nHalf], j<-[1..nHalf]] ++ 126 [(2*i,2*j-1) =: 0.5*(mesh!(i,j) + mesh!(i,j-1)) 127 | i<-[1..nHalf], j<-[1..nHalf]] ++ 128 [(2*i-1,2*j) =: 0.5*(mesh!(i,j) + mesh!(i-1,j)) 129 | i<-[1..nHalf], j<-[1..nHalf]] ++ 130 [(2*i-1,2*j-1) =: 0.25*(mesh!(i,j) + mesh!(i,j-1) + mesh!(i-1,j) + 131 mesh!(i-1,j-1)) | i<-[1..nHalf], j<-[1..nHalf]]) 132 -} 133 where 134 nHalf = nHalf'-1 135 n = 2 * nHalf' - 1 136 137 138 -- random number generator 139 140 genRand :: Value -> Value 141 genRand seed = 142 r1 / 655357 143 where 144 r1 = (31453257*seed + 271829) `fiRem` 655357 145 x `fiRem` m = x - fromIntegral ((truncate x `div` m) * m) 146 147 148 log2 :: Int -> Int 149 log2 n = log2' n 0 150 where 151 log2' n accum 152 | n > 1 = log2' (n `div` 2) (accum+1) 153 | otherwise = accum