[project @ 2000-03-30 14:06:46 by simonpj]
[nofib.git] / real / pic / Utils.hs
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