1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
#include <blitz/array.h>
#include <fstream>
using namespace blitz;
void makeLogo();
int main()
{
makeLogo();
return 0;
}
void setInitialConditions(Array<float,2>& c, Array<float,2>& P1,
Array<float,2>& P2, Array<float,2>& P3, int N, int M);
void snapshot(const Array<float,2>& P, const Array<float,2>& c);
void makeLogo()
{
const int N = 300, M = 900;
int niters = 3000;
Array<float,2> P1, P2, P3, c;
allocateArrays(shape(N,M), P1, P2, P3, c);
Range I(1,N-2), J(1,M-2);
setInitialConditions(c, P1, P2, P3, N, M);
for (int iter=0; iter < niters; ++iter)
{
P3(I,J) = (2-4*c(I,J)) * P2(I,J)
+ c(I,J)*(P2(I-1,J) + P2(I+1,J) + P2(I,J-1) + P2(I,J+1))
- P1(I,J);
cycleArrays(P1,P2,P3);
snapshot(P2, c);
}
}
void setInitialConditions(Array<float,2>& c, Array<float,2>& P1,
Array<float,2>& P2, Array<float,2>& P3, int N, int M)
{
// Set the velocity field
c = 0.3;
ifstream ifs("blitz3.pgm");
char tmpBuf[128];
int pixel;
ifs.getline(tmpBuf, 128);
ifs.getline(tmpBuf, 128);
ifs.getline(tmpBuf, 128);
for (int pi=0; pi < 199; ++pi)
{
for (int pj=0; pj < 798; ++pj)
{
ifs >> pixel;
if (pixel)
c(pi+50,pj+56) = 0.02;
}
}
// Initial pressure distribution: gaussian pulse
using namespace blitz::tensor;
int cr = N/6-1;
// int cc = 7.0*M/8.0-1;
float s2 = 64.0 * 9.0 / pow2(N/2.0);
P1 = 0.0;
// P2 = exp(-(pow2(i-cr)+pow2(j-cc)) * s2);
P2 = exp(-(pow2(i-cr)) * s2);
P3 = 0.0;
}
void snapshot(const Array<float,2>& P, const Array<float,2>& c)
{
static int count = 0, snapshotNum = 0;
if (++count < 50)
return;
count = 0;
++snapshotNum;
char filename[128];
sprintf(filename, "logo%03d.m", snapshotNum);
ofstream ofs(filename);
int N = P.length(firstDim);
int M = P.length(secondDim);
float Pmin = -0.6;
float PScale = 1.0/1.2;
float VScale = 1.0;
ofs << "P" << snapshotNum << " = [ ";
for (int i=0; i < N; ++i)
{
for (int j=0; j < M; ++j)
{
float value1 = (P(i,j)-Pmin)*PScale;
float value2 = c(i,j)*VScale;
int r1 = value1 * 4096;
int r2 = value2 * 4096;
ofs << r1 << " " << r2 << " ";
}
if (i < N-1)
ofs << ";" << endl;
}
ofs << "];" << endl;
}
|
