Next: A..8 7_Wave_2D_ColorMap.c Up: A. 講習会のサンプル・プログラムを読む Previous: A..6 5_Turing_2D_Contln.c
| 6_Turing_2D_ColorMap.c |
1 /******************************************************************************************
2
3 This program solves reacrtion-diffusion equations using the explicit method.
4
5 Reacrtion-diffusion equations are ut = du uxx + f(u,v) and vt = dv vxx + g(u,v).
6
7 Time loop in this program consists of
8 "Calculation part" and "Draw Part".
9
10 Calculation part is controlled by the interval parameter INTV.
11 This numerical result is visualized by using "g_color_map" made in thie program.
12
13 *******************************************************************************************/
14
15 #include<stdio.h>
16 #include<stdlib.h>
17 #include<glsc3d_3.h>
18
19 #define N (100)
20 #define M (100)
21
22 //////////// Color map: g_color_map ////////////
23 void g_color_map(double Array[][M+2],
24 double dx_width, double dy_width,
25 double L_bottom_x, double L_bottom_y,
26 double max, double min
27 )
28 {
29 int Fine_Grid = 1;
30 for (int j=0; j < M+2; j++) {
31 for (int i=0; i < N+2; i++) {
32 if (Fine_Grid == 0) {
33 g_area_color((Array[i][j]-min)/(max-min), 0, 1.0-(Array[i][j]-min)/(max-min), 0.5);
34 g_box_center_2D(dx_width * (i) + L_bottom_x, dy_width * (j) + L_bottom_y, dx_width, dy_width, G_NO, G_YES);
35 }
36
37 if (Fine_Grid == 1) {
38 double tmpColor;
39 if (j-1<0) tmpColor = Array[i][j];
40 else tmpColor = Array[i][j]*0.5 + Array[i][j-1]*0.5;
41 g_area_color((tmpColor-min)/(max-min), 0, 1.0-(tmpColor-min)/(max-min), 1);
42 g_triangle_2D(dx_width * (i) + L_bottom_x, dy_width * (j) + L_bottom_y,
43 dx_width * (i) + L_bottom_x - dx_width*0.5, dy_width * (j) + L_bottom_y - dy_width*0.5,
44 dx_width * (i) + L_bottom_x + dx_width*0.5, dy_width * (j) + L_bottom_y - dy_width*0.5,
45 G_NO, G_YES);
46
47 if (i+1>=N+2) tmpColor = Array[i][j];
48 else tmpColor = Array[i][j]*0.5 + Array[i+1][j]*0.5;
49 g_area_color((tmpColor-min)/(max-min), 0, 1.0-(tmpColor-min)/(max-min), 1);
50 g_triangle_2D(dx_width * (i) + L_bottom_x, dy_width * (j) + L_bottom_y,
51 dx_width * (i) + L_bottom_x + dx_width*0.5, dy_width * (j) + L_bottom_y - dy_width*0.5,
52 dx_width * (i) + L_bottom_x + dx_width*0.5, dy_width * (j) + L_bottom_y + dy_width*0.5,
53 G_NO, G_YES);
54
55 if (j+1>=M+2) tmpColor = Array[i][j];
56 else tmpColor = Array[i][j]*0.5 + Array[i][j+1]*0.5;
57 g_area_color((tmpColor-min)/(max-min), 0, 1.0-(tmpColor-min)/(max-min), 1);
58 g_triangle_2D(dx_width * (i) + L_bottom_x, dy_width * (j) + L_bottom_y,
59 dx_width * (i) + L_bottom_x + dx_width*0.5, dy_width * (j) + L_bottom_y + dy_width*0.5,
60 dx_width * (i) + L_bottom_x - dx_width*0.5, dy_width * (j) + L_bottom_y + dy_width*0.5,
61 G_NO, G_YES);
62
63 if (i-1<0) tmpColor = Array[i][j];
64 else tmpColor = Array[i][j]*0.5 + Array[i-1][j]*0.5;
65 g_area_color((tmpColor-min)/(max-min), 0, 1.0-(tmpColor-min)/(max-min), 1);
66 g_triangle_2D(dx_width * (i) + L_bottom_x, dy_width * (j) + L_bottom_y,
67 dx_width * (i) + L_bottom_x - dx_width*0.5, dy_width * (j) + L_bottom_y + dy_width*0.5,
68 dx_width * (i) + L_bottom_x - dx_width*0.5, dy_width * (j) + L_bottom_y - dy_width*0.5,
69 G_NO, G_YES);
70 }
71 }
72 }
73 }
74
75 double f(double u, double v)
76 {
77 return u * (1 - u * u) - v;
78 }
79 double g(double u, double v)
80 {
81 return 3.0 * u - 2.0 * v;
82 // return 3.0 * (u + 0.05) - 2.0 * v;
83 // return 3.0 * (u - 0.05) - 2.0 * v;
84 }
85
86 int main()
87 {
88 double Lx = 100.0;
89 double Ly = 100.0;
90 double dx = Lx / N;
91 double dy = Ly / M;
92 double dt = 0.001;
93 double du = 1.0;
94 double dv = 10.0;
95
96 double lambda_u = du * dt / (dx * dx);
97 double lambda_v = dv * dt / (dy * dy);
98
99 double u0[N + 2][M + 2];
100 double u1[N + 2][M + 2];
101 double v0[N + 2][M + 2];
102 double v1[N + 2][M + 2];
103
104 int INTV = 200;
105
106 g_init("Window", 600, 600); //Pixel Size
107
108 g_def_scale_2D(0, //ID
109 -Lx*0.5, Lx*0.5, //xmin,xmax
110 -Ly*0.5, Ly*0.5, //ymin,ymax
111 20.0, 20.0, //Window (Left, Top) Position
112 560, 560); //Window Size (x,y)
113 //Set initial calculation
114 {
115 //Step1
116 double u_Init_Sum = 0.0;
117 double v_Init_Sum = 0.0;
118
119 unsigned int seed = 2; //Seed for random number generated by rand()
120 srand(seed);
121
122 for(int j = 1;j < M + 1; j ++)
123 {
124 for(int i = 1;i < N + 1; i ++)
125 {
126 double Max = 1.0; //Maximum of random number
127 double Min = 0.0; //Minimum of random number
128 double RandNumber = Min + (int)(rand()*(Max-Min+1.0)/(1.0+RAND_MAX));
129 u0[i][j] = 2.0 * (RandNumber - 0.5);
130 v0[i][j] = 2.0 * (RandNumber - 0.5);
131
132 u_Init_Sum += u0[i][j] * dx * dy;
133 v_Init_Sum += v0[i][j] * dx * dy;
134 }
135 }
136
137 //Step2
138 for(int j = 1;j < M + 1; j ++)
139 {
140 u0[0][j] = u0[N][j];
141 u0[N + 1][j] = u0[1][j];
142 v0[0][j] = v0[N][j];
143 v0[N + 1][j] = v0[1][j];
144 }
145 for(int i = 0;i < N + 2; i ++)
146 {
147 u0[i][0] = u0[i][M];
148 u0[i][M + 1] = u0[i][1];
149 v0[i][0] = v0[i][M];
150 v0[i][M + 1] = v0[i][1];
151 }
152 }
153
154 //////////// Start time loop ////////////
155 for (int i_time = 0; ; i_time++) {
156
157 //////////// Draw Part ////////////
158 if (i_time%INTV == 0) {
159 g_cls(); //Clear window
160 g_sel_scale(0); //Select Virtual scale
161 g_line_color(0.0, 0.0, 0.0, 1.0);
162 g_boundary(); //Draw Boundary
163
164 g_color_map(u0, dx, dy, -Lx*0.5 + dx * 0.5, -Ly*0.5 + dy * 0.5, 0.5, -0.5);
165
166 g_finish(); //flush Draw buffer
167 g_sleep(0.01); //Sleep 0.01 sec
168 }
169
170 //////////// Calculation Part ////////////
171 {
172 //Step3
173 for(int j = 1;j < M + 1; j ++)
174 {
175 for(int i = 1;i < N + 1; i ++)
176 {
177 u1[i][j] = u0[i][j]
178 + lambda_u * (u0[i - 1][j] - 2 * u0[i][j] + u0[i + 1][j])
179 + lambda_u * (u0[i][j - 1] - 2 * u0[i][j] + u0[i][j + 1])
180 + dt * f(u0[i][j], v0[i][j]);
181
182 v1[i][j] = v0[i][j]
183 + lambda_v * (v0[i - 1][j] - 2 * v0[i][j] + v0[i + 1][j])
184 + lambda_v * (v0[i][j - 1] - 2 * v0[i][j] + v0[i][j + 1])
185 + dt * g(u0[i][j], v0[i][j]);
186 }
187 }
188 //Step4
189 for(int j = 1;j < M + 1; j ++)
190 {
191 for(int i = 1;i < N + 1; i ++)
192 {
193 u0[i][j] = u1[i][j];
194 v0[i][j] = v1[i][j];
195 }
196 }
197 //Step5
198 for(int j = 1;j < M + 1; j ++)
199 {
200 u0[0][j] = u0[N][j];
201 u0[N + 1][j] = u0[1][j];
202 v0[0][j] = v0[N][j];
203 v0[N + 1][j] = v0[1][j];
204 }
205 for(int i = 0;i < N + 2; i ++)
206 {
207 u0[i][0] = u0[i][M];
208 u0[i][M + 1] = u0[i][1];
209 v0[i][0] = v0[i][M];
210 v0[i][M + 1] = v0[i][1];
211 }
212 }
213
214 //////////// Check part of maximum and minimum of u and v in order to use g_color_map ////////////
215 // double u_Max = u0[0][0];
216 // double u_Min = u0[0][0];
217 // double v_Max = v0[0][0];
218 // double v_Min = v0[0][0];
219 // for(int j = 1;j < M + 1; j ++)
220 // {
221 // for(int i = 1;i < N + 1; i ++)
222 // {
223 // if (u_Max < u0[i][j]) u_Max = u0[i][j];
224 // if (u_Min > u0[i][j]) u_Min = u0[i][j];
225 // if (v_Max < v0[i][j]) v_Max = v0[i][j];
226 // if (v_Min > v0[i][j]) v_Min = v0[i][j];
227 // }
228 // }
229 // printf("%d th: u_Max = % 2.15f\n", i_time, u_Max);
230 // printf("%d th: u_Min = % 2.15f\n", i_time, u_Min);
231 // printf("%d th: v_Max = % 2.15f\n", i_time, v_Max);
232 // printf("%d th: v_Min = % 2.15f\n", i_time, v_Min);
233 }
234 return 0;
235 }
|