1 /*
2     Copyright (c) 2005-2019 Intel Corporation
3 
4     Licensed under the Apache License, Version 2.0 (the "License");
5     you may not use this file except in compliance with the License.
6     You may obtain a copy of the License at
7 
8         http://www.apache.org/licenses/LICENSE-2.0
9 
10     Unless required by applicable law or agreed to in writing, software
11     distributed under the License is distributed on an "AS IS" BASIS,
12     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13     See the License for the specific language governing permissions and
14     limitations under the License.
15 */
16 
17 #include "video.h"
18 #include <cmath>
19 
20 using namespace std;
21 
22 #ifdef _MSC_VER
23 // warning C4068: unknown pragma
24 #pragma warning(disable: 4068)
25 // warning C4351: new behavior: elements of array 'array' will be default initialized
26 #pragma warning(disable: 4351)
27 #endif
28 
29 #include "universe.h"
30 
31 const colorcomp_t MaterialColor[4][3] = { // BGR
32     {96,0,0},     // WATER
33     {0,48,48},    // SANDSTONE
34     {32,32,23}    // SHALE
35 };
36 
InitializeUniverse(video const & colorizer)37 void Universe::InitializeUniverse(video const& colorizer) {
38 
39     pulseCounter = pulseTime = 100;
40     pulseX = UniverseWidth/3;
41     pulseY = UniverseHeight/4;
42     // Initialize V, S, and T to slightly non-zero values, in order to avoid denormal waves.
43     for( int i=0; i<UniverseHeight; ++i )
44 #pragma ivdep
45         for( int j=0; j<UniverseWidth; ++j ) {
46             T[i][j] = S[i][j] = V[i][j] = ValueType(1.0E-6);
47         }
48     for( int i=1; i<UniverseHeight-1; ++i ) {
49         for( int j=1; j<UniverseWidth-1; ++j ) {
50             float x = float(j-UniverseWidth/2)/(UniverseWidth/2);
51             ValueType t = (ValueType)i/UniverseHeight;
52             MaterialType m;
53             D[i][j] = 1.0;
54             // Coefficient values are fictitious, and chosen to visually exaggerate
55             // physical effects such as Rayleigh waves.  The fabs/exp line generates
56             // a shale layer with a gentle upwards slope and an anticline.
57             if( t<0.3f ) {
58                 m = WATER;
59                 M[i][j] = 0.125;
60                 L[i][j] = 0.125;
61             } else if( fabs(t-0.7+0.2*exp(-8*x*x)+0.025*x)<=0.1 ) {
62                 m = SHALE;
63                 M[i][j] = 0.5;
64                 L[i][j] = 0.6;
65             } else {
66                 m = SANDSTONE;
67                 M[i][j] = 0.3;
68                 L[i][j] = 0.4;
69             }
70             material[i][j] = m;
71         }
72     }
73     ValueType scale = 2.0f/ColorMapSize;
74     for( int k=0; k<4; ++k ) {
75         for( int i=0; i<ColorMapSize; ++i ) {
76             colorcomp_t c[3];
77             ValueType t = (i-ColorMapSize/2)*scale;
78             ValueType r = t>0 ? t : 0;
79             ValueType b = t<0 ? -t : 0;
80             ValueType g = 0.5f*fabs(t);
81             std::memcpy(c, MaterialColor[k], sizeof(c));
82             c[2] = colorcomp_t(r*(255-c[2])+c[2]);
83             c[1] = colorcomp_t(g*(255-c[1])+c[1]);
84             c[0] = colorcomp_t(b*(255-c[0])+c[0]);
85             ColorMap[k][i] = colorizer.get_color(c[2], c[1], c[0]);
86         }
87     }
88     // Set damping coefficients around border to reduce reflections from boundaries.
89     ValueType d = 1.0;
90     for( int k=DamperSize-1; k>0; --k ) {
91         d *= 1-1.0f/(DamperSize*DamperSize);
92         for( int j=1; j<UniverseWidth-1; ++j ) {
93             D[k][j] *= d;
94             D[UniverseHeight-1-k][j] *= d;
95         }
96         for( int i=1; i<UniverseHeight-1; ++i ) {
97             D[i][k] *= d;
98             D[i][UniverseWidth-1-k] *= d;
99         }
100     }
101     drawingMemory = colorizer.get_drawing_memory();
102 }
UpdatePulse()103 void Universe::UpdatePulse() {
104   if( pulseCounter>0 ) {
105     ValueType t = (pulseCounter-pulseTime/2)*0.05f;
106     V[pulseY][pulseX] += 64*sqrt(M[pulseY][pulseX])*exp(-t*t);
107     --pulseCounter;
108   }
109 }
110 
111 
UpdateStress(Rectangle const & r)112 void Universe::UpdateStress(Rectangle const& r ) {
113     drawing_area  drawing(r.StartX(),r.StartY(),r.Width(),r.Height(),drawingMemory);
114     for( int i=r.StartY(); i<r.EndY() ; ++i ) {
115         drawing.set_pos(1, i-r.StartY());
116 #pragma ivdep
117         for( int j=r.StartX(); j<r.EndX() ; ++j ) {
118             S[i][j] += M[i][j]*(V[i][j+1]-V[i][j]);
119             T[i][j] += M[i][j]*(V[i+1][j]-V[i][j]);
120             int index = (int)(V[i][j]*(ColorMapSize/2)) + ColorMapSize/2;
121             if( index<0 ) index = 0;
122             if( index>=ColorMapSize ) index = ColorMapSize-1;
123             color_t* c = ColorMap[material[i][j]];
124             drawing.put_pixel(c[index]);
125         }
126     }
127 }
128 
SerialUpdateStress()129 void Universe::SerialUpdateStress() {
130     Rectangle area(0, 0, UniverseWidth-1, UniverseHeight-1);
131     UpdateStress(area);
132 }
133 
134 //void Universe::ParallelUpdateStress(tbb::affinity_partitioner &affinity) {
135 //    tbb::parallel_for( tbb::blocked_range<int>( 0, UniverseHeight-1 ), // Index space for loop
136 //                       UpdateStressBody(*this),                             // Body of loop
137 //                       affinity );                                     // Affinity hint
138 //}
139 
UpdateVelocity(Rectangle const & r)140 void Universe::UpdateVelocity(Rectangle const& r) {
141     for( int i=r.StartY(); i<r.EndY(); ++i )
142 #pragma ivdep
143         for( int j=r.StartX(); j<r.EndX(); ++j )
144             V[i][j] = D[i][j]*(V[i][j] + L[i][j]*(S[i][j] - S[i][j-1] + T[i][j] - T[i-1][j]));
145 }
146 
SerialUpdateVelocity()147 void Universe::SerialUpdateVelocity() {
148     UpdateVelocity(Rectangle(1,1,UniverseWidth-1,UniverseHeight-1));
149 }
150 
151 
152 //void Universe::ParallelUpdateVelocity(tbb::affinity_partitioner &affinity) {
153 //    tbb::parallel_for( tbb::blocked_range<int>( 1, UniverseHeight ), // Index space for loop
154 //                       UpdateVelocityBody(*this),                    // Body of loop
155 //                       affinity );                                   // Affinity hint
156 //}
157 
SerialUpdateUniverse()158 void Universe::SerialUpdateUniverse() {
159   UpdatePulse();
160   SerialUpdateStress();
161   SerialUpdateVelocity();
162 }
163 
164 //void Universe::ParallelUpdateUniverse() {
165 //  /** Affinity is an argument to parallel_for to hint that an iteration of a loop
166 //    is best replayed on the same processor for each execution of the loop.
167 //    It is a static object because it must remember where the iterations happened
168 //    in previous executions. */
169 //  static tbb::affinity_partitioner affinity;
170 //  UpdatePulse();
171 //  ParallelUpdateStress(affinity);
172 //  ParallelUpdateVelocity(affinity);
173 //}
174 
TryPutNewPulseSource(int x,int y)175 bool Universe::TryPutNewPulseSource(int x, int y){
176     if(pulseCounter == 0) {
177         pulseCounter = pulseTime;
178         pulseX = x; pulseY = y;
179         return true;
180     }
181     return false;
182 }
183 
SetDrawingMemory(const drawing_memory & dmem)184 void Universe::SetDrawingMemory(const drawing_memory &dmem) {
185     drawingMemory = dmem;
186 }
187