A..7 6_Turing_2D_ColorMap.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 }



桂田 祐史