00001 #ifndef OPENTISSUE_DYNAMICS_SWE_SWE_SHALLOW_WATER_EQUATIONS_H
00002 #define OPENTISSUE_DYNAMICS_SWE_SWE_SHALLOW_WATER_EQUATIONS_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/utility/gl/gl_util.h>
00013 #include <OpenTissue/core/math/math_precision.h>
00014 #include <OpenTissue/core/math/big/big_types.h>
00015 #include <OpenTissue/core/math/big/big_conjugate_gradient.h>
00016 #include <OpenTissue/utility/utility_timer.h>
00017
00018 #include <boost/multi_array.hpp>
00019
00020 #include <cmath>
00021 #include <cassert>
00022 #include <iostream>
00023
00024 namespace OpenTissue
00025 {
00026 namespace swe
00027 {
00028
00040 template<typename real_type_ = double>
00041 class ShallowWaterEquations
00042 {
00043 public:
00044
00045 typedef real_type_ real_type;
00046 typedef boost::multi_array<real_type, 2> array_type;
00047 typedef unsigned int index_type;
00048 typedef ublas::compressed_matrix<real_type> matrix_type;
00049 typedef ublas::vector<real_type> vector_type;
00050
00051 protected:
00052
00053 real_type g;
00054 real_type timestep;
00055 unsigned int X;
00056 unsigned int Y;
00057 real_type dx;
00058 real_type dy;
00059 array_type water;
00060 array_type b;
00061 array_type h;
00062 array_type u;
00063 array_type v;
00064 array_type d;
00065 array_type h_tilde;
00066 array_type u_tilde;
00067 array_type v_tilde;
00068 array_type x_tilde;
00069 array_type y_tilde;
00070 matrix_type A;
00071 vector_type hnew;
00072 vector_type rhs;
00073
00074 public:
00075
00076 ShallowWaterEquations()
00077 : g(9.81)
00078 , timestep(0)
00079 , X(0)
00080 , Y(0)
00081 , dx(0)
00082 , dy(0)
00083 , water(boost::extents[0][0])
00084 , b(boost::extents[0][0])
00085 , h(boost::extents[0][0])
00086 , u(boost::extents[0][0])
00087 , v(boost::extents[0][0])
00088 , d(boost::extents[0][0])
00089 , h_tilde(boost::extents[0][0])
00090 , u_tilde(boost::extents[0][0])
00091 , v_tilde(boost::extents[0][0])
00092 , x_tilde(boost::extents[0][0])
00093 , y_tilde(boost::extents[0][0])
00094 {}
00095
00096 ~ShallowWaterEquations() { clear(); }
00097
00098 public:
00099
00107 void init( const unsigned int M, const unsigned int N, const real_type width, const real_type height )
00108 {
00109 clear();
00110
00111 g = 9.81;
00112
00113 this->X = M;
00114 this->Y = N;
00115 this->dx = width / X;
00116 this->dy = height / Y;
00117
00118 water.resize(boost::extents[X][Y]);
00119 b.resize(boost::extents[X][Y]);
00120 h.resize(boost::extents[X][Y]);
00121 u.resize(boost::extents[X][Y]);
00122 v.resize(boost::extents[X][Y]);
00123 d.resize(boost::extents[X][Y]);
00124 h_tilde.resize(boost::extents[X][Y]);
00125 u_tilde.resize(boost::extents[X][Y]);
00126 v_tilde.resize(boost::extents[X][Y]);
00127 x_tilde.resize(boost::extents[X][Y]);
00128 y_tilde.resize(boost::extents[X][Y]);
00129
00130 for ( index_type x = 0; x < X; ++x )
00131 for ( index_type y = 0; y < Y; ++y )
00132 {
00133 water[ x ][ y ] = 1;
00134 b[ x ][ y ] = 0;
00135 h[ x ][ y ] = 2;
00136 u[ x ][ y ] = 0;
00137 v[ x ][ y ] = 0;
00138 d[ x ][ y ] = 0;
00139 h_tilde[ x ][ y ] = 0;
00140 u_tilde[ x ][ y ] = 0;
00141 v_tilde[ x ][ y ] = 0;
00142 x_tilde[ x ][ y ] = 0;
00143 y_tilde[ x ][ y ] = 0;
00144 }
00145
00146 int nvars = X * Y;
00147
00148
00149 A.resize( nvars, nvars, false );
00150 A.clear();
00151
00152 hnew.resize( nvars );
00153 hnew.clear();
00154 rhs.resize( nvars );
00155 rhs.clear();
00156 }
00157
00158 void clear()
00159 {
00160 X = 0;
00161 Y = 0;
00162 dx = 0;
00163 dy = 0;
00164 }
00165
00166
00173 void draw( const bool debugMode ) const
00174 {
00175
00176 if(X<=1 || Y<=1)
00177 return;
00178
00179 if ( debugMode )
00180 {
00181 for ( index_type i = 0; i < X - 1; ++i )
00182 {
00183 for ( index_type j = 0; j < Y - 1; ++j )
00184 {
00185 real_type x = i * dx;
00186 real_type y = j * dx;
00187 real_type z00 = h[ i ][ j ];
00188 real_type z10 = h[ i + 1 ][ j ];
00189 real_type z01 = h[ i ][ j + 1 ];
00190 real_type z11 = h[ i + 1 ][ j + 1 ];
00191 gl::ColorPicker( 0.0, 0.0, 0.5 );
00192 glBegin( GL_LINE_LOOP );
00193 glVertex3f( x, y, z00 );
00194 glVertex3f( x + dx, y, z10 );
00195 glVertex3f( x + dx, y + dy, z11 );
00196 glVertex3f( x, y + dy, z01 );
00197 glEnd();
00198 real_type z00_tilde = h_tilde[ i ][ j ];
00199 real_type z10_tilde = h_tilde[ i + 1 ][ j ];
00200 real_type z01_tilde = h_tilde[ i ][ j + 1 ];
00201 real_type z11_tilde = h_tilde[ i + 1 ][ j + 1 ];
00202 gl::ColorPicker( 0., 1., 0. );
00203 glBegin( GL_LINE_LOOP );
00204 glVertex3f( x, y, z00_tilde );
00205 glVertex3f( x + dx, y, z10_tilde );
00206 glVertex3f( x + dx, y + dy, z11_tilde );
00207 glVertex3f( x, y + dy, z01_tilde );
00208 glEnd();
00209 real_type b00 = b[ i ][ j ];
00210 real_type b10 = b[ i + 1 ][ j ];
00211 real_type b01 = b[ i ][ j + 1 ];
00212 real_type b11 = b[ i + 1 ][ j + 1 ];
00213 gl::ColorPicker( 0.5, 0.1, 0 );
00214 glBegin( GL_LINE_LOOP );
00215 glVertex3f( x, y, b00 );
00216 glVertex3f( x + dx, y, b10 );
00217 glVertex3f( x + dx, y + dy, b11 );
00218 glVertex3f( x, y + dy, b01 );
00219 glEnd();
00220 }
00221 }
00222 for ( index_type i = 0; i < X; ++i )
00223 {
00224 for ( index_type j = 0; j < Y; ++j )
00225 {
00226 real_type x = i * dx;
00227 real_type y = j * dx;
00228 real_type z00 = h[ i ][ j ];
00229 real_type z00_tilde = h_tilde[ i ][ j ];
00230 gl::ColorPicker( 0, 0.5, 0.5 );
00231 glBegin( GL_LINES );
00232 glVertex3f( x, y, z00 );
00233 glVertex3f( x, y, z00_tilde );
00234 glEnd();
00235 real_type vx = u[ i ][ j ];
00236 real_type vy = v[ i ][ j ];
00237 gl::ColorPicker( 0.5, 0.5, 0 );
00238 glBegin( GL_LINES );
00239 glVertex3f( x, y, z00 );
00240 glVertex3f( x + vx, y + vy, z00 );
00241 glEnd();
00242 real_type px = x_tilde[ i ][ j ];
00243 real_type py = y_tilde[ i ][ j ];
00244 gl::ColorPicker( 0, 0.5, 0 );
00245 glBegin( GL_LINES );
00246 glVertex3f( x, y, z00 );
00247 glVertex3f( px, py, z00 );
00248 glEnd();
00249 real_type wx = u_tilde[ i ][ j ];
00250 real_type wy = v_tilde[ i ][ j ];
00251 gl::ColorPicker( 0.5, 0, 0.5 );
00252 glBegin( GL_LINES );
00253 glVertex3f( x, y, z00 );
00254 glVertex3f( x + wx, y + wy, z00 );
00255 glEnd();
00256 }
00257 }
00258 }
00259 else
00260 {
00261 glBegin( GL_QUADS );
00262 gl::ColorPicker( 0.3, 0.2, 0.1, 1.0 );
00263 for ( index_type i = 0; i < X - 1; ++i )
00264 {
00265 for ( index_type j = 0; j < Y - 1; ++j )
00266 {
00267 real_type x = i * dx;
00268 real_type y = j * dx;
00269 real_type b00 = b[ i ][ j ];
00270 real_type b10 = b[ i + 1 ][ j ];
00271 real_type b01 = b[ i ][ j + 1 ];
00272 real_type b11 = b[ i + 1 ][ j + 1 ];
00273
00274
00275
00276
00277 float nx = -dy * ( b[ i + 1 ][ j ] - b[ i ][ j ] );
00278 float ny = -( dx * ( b[ i ][ j + 1 ] - b[ i ][ j ] ) );
00279 float nz = dx * dy;
00280
00281 glNormal3f( nx, ny, nz );
00282 glVertex3f( x, y, b00 );
00283 glNormal3f( nx, ny, nz );
00284 glVertex3f( x + dx, y, b10 );
00285 glNormal3f( nx, ny, nz );
00286 glVertex3f( x + dx, y + dy, b11 );
00287 glNormal3f( nx, ny, nz );
00288 glVertex3f( x, y + dy, b01 );
00289 }
00290 }
00291 gl::ColorPicker( 0.2, 0.2, 0.2, 1.0 );
00292 for ( index_type i = 0; i < X - 1; ++i )
00293 {
00294 for ( index_type j = 0; j < Y - 1; ++j )
00295 {
00296 real_type x = i * dx;
00297 real_type y = j * dx;
00298
00299 bool s00 = ( water[ i ][ j ] == 0 );
00300 bool s10 = ( water[ i + 1 ][ j ] == 0 );
00301 bool s01 = ( water[ i ][ j + 1 ] == 0 );
00302 bool s11 = ( water[ i + 1 ][ j + 1 ] == 0 );
00303 if ( s00 && s10 && s01 && s11 )
00304 {
00305 glNormal3f( 0, 0, 1 );
00306 glVertex3f( x, y, 4 );
00307 glNormal3f( 0, 0, 1 );
00308 glVertex3f( x + dx, y, 4 );
00309 glNormal3f( 0, 0, 1 );
00310 glVertex3f( x + dx, y + dy, 4 );
00311 glNormal3f( 0, 0, 1 );
00312 glVertex3f( x, y + dy, 4 );
00313 }
00314 else
00315 {
00316 if ( s00 && s01 )
00317 {
00318 glNormal3f( 1, 0, 0 );
00319 glVertex3f( x, y, 4 );
00320 glNormal3f( 1, 0, 0 );
00321 glVertex3f( x, y, 0 );
00322 glNormal3f( 1, 0, 0 );
00323 glVertex3f( x, y + dy, 0 );
00324 glNormal3f( 1, 0, 0 );
00325 glVertex3f( x, y + dy, 4 );
00326 }
00327 if ( s10 && s11 )
00328 {
00329 glNormal3f( -1, 0, 0 );
00330 glVertex3f( x + dx, y, 0 );
00331 glNormal3f( -1, 0, 0 );
00332 glVertex3f( x + dx, y, 4 );
00333 glNormal3f( -1, 0, 0 );
00334 glVertex3f( x + dx, y + dy, 4 );
00335 glNormal3f( -1, 0, 0 );
00336 glVertex3f( x + dx, y + dy, 0 );
00337 }
00338 if ( s00 && s10 )
00339 {
00340 glNormal3f( 0, 1, 0 );
00341 glVertex3f( x, y, 0 );
00342 glNormal3f( 0, 1, 0 );
00343 glVertex3f( x, y, 4 );
00344 glNormal3f( 0, 1, 0 );
00345 glVertex3f( x + dx, y, 4 );
00346 glNormal3f( 0, 1, 0 );
00347 glVertex3f( x + dx, y, 0 );
00348 }
00349 if ( s01 && s11 )
00350 {
00351 glNormal3f( 0, -1, 0 );
00352 glVertex3f( x, y + dy, 0 );
00353 glNormal3f( 0, -1, 0 );
00354 glVertex3f( x + dx, y + dy, 0 );
00355 glNormal3f( 0, -1, 0 );
00356 glVertex3f( x + dx, y + dy, 4 );
00357 glNormal3f( 0, -1, 0 );
00358 glVertex3f( x, y + dy, 4 );
00359 }
00360 }
00361 }
00362 }
00363 real_type x = 0;
00364 for ( index_type j = 0; j < Y - 1; ++j )
00365 {
00366 real_type y = j * dx;
00367
00368 bool s00 = ( water[ 0 ][ j ] == 0 );
00369 bool s01 = ( water[ 0 ][ j + 1 ] == 0 );
00370 if ( s00 && s01 )
00371 {
00372 glNormal3f( 1, 0, 0 );
00373 glVertex3f( x, y, 4 );
00374 glNormal3f( 1, 0, 0 );
00375 glVertex3f( x, y, 0 );
00376 glNormal3f( 1, 0, 0 );
00377 glVertex3f( x, y + dy, 0 );
00378 glNormal3f( 1, 0, 0 );
00379 glVertex3f( x, y + dy, 4 );
00380 }
00381 }
00382 x = ( X - 2 ) * dx;
00383 for ( index_type j = 0; j < Y - 1; ++j )
00384 {
00385 real_type y = j * dx;
00386
00387 bool s10 = ( water[ X - 1 ][ j ] == 0 );
00388 bool s11 = ( water[ X - 1 ][ j + 1 ] == 0 );
00389 if ( s10 && s11 )
00390 {
00391 glNormal3f( -1, 0, 0 );
00392 glVertex3f( x + dx, y, 0 );
00393 glNormal3f( -1, 0, 0 );
00394 glVertex3f( x + dx, y, 4 );
00395 glNormal3f( -1, 0, 0 );
00396 glVertex3f( x + dx, y + dy, 4 );
00397 glNormal3f( -1, 0, 0 );
00398 glVertex3f( x + dx, y + dy, 0 );
00399 }
00400 }
00401 real_type y = 0;
00402 for ( index_type i = 0; i < X - 1; ++i )
00403 {
00404 real_type x = i * dx;
00405
00406 bool s00 = ( water[ i ][ 0 ] == 0 );
00407 bool s10 = ( water[ i + 1 ][ 0 ] == 0 );
00408 if ( s00 && s10 )
00409 {
00410 glNormal3f( 0, 1, 0 );
00411 glVertex3f( x, y, 0 );
00412 glNormal3f( 0, 1, 0 );
00413 glVertex3f( x, y, 4 );
00414 glNormal3f( 0, 1, 0 );
00415 glVertex3f( x + dx, y, 4 );
00416 glNormal3f( 0, 1, 0 );
00417 glVertex3f( x + dx, y, 0 );
00418 }
00419 }
00420 y = ( Y - 2 ) * dx;
00421 for ( index_type i = 0; i < X - 1; ++i )
00422 {
00423 real_type x = i * dx;
00424
00425 bool s01 = ( water[ i ][ Y - 1 ] == 0 );
00426 bool s11 = ( water[ i + 1 ][ Y - 1 ] == 0 );
00427 if ( s01 && s11 )
00428 {
00429 glBegin( GL_POLYGON );
00430 glNormal3f( 0, -1, 0 );
00431 glVertex3f( x, y + dy, 0 );
00432 glNormal3f( 0, -1, 0 );
00433 glVertex3f( x + dx, y + dy, 0 );
00434 glNormal3f( 0, -1, 0 );
00435 glVertex3f( x + dx, y + dy, 4 );
00436 glNormal3f( 0, -1, 0 );
00437 glVertex3f( x, y + dy, 4 );
00438 glEnd();
00439 }
00440 }
00441
00442 bool disableBlend = false;
00443 if ( !glIsEnabled( GL_BLEND ) )
00444 {
00445 glEnable( GL_BLEND );
00446 glBlendFunc( GL_ONE, GL_ONE_MINUS_SRC_ALPHA );
00447 disableBlend = true;
00448 }
00449
00450 gl::ColorPicker( 0.0, 0.15, 0.5, 0.35 );
00451
00452 for ( index_type i = 0; i < X - 1; ++i )
00453 {
00454 for ( index_type j = 0; j < Y - 1; ++j )
00455 {
00456 real_type x = i * dx;
00457 real_type y = j * dx;
00458 real_type z00 = h[ i ][ j ];
00459 real_type z10 = h[ i + 1 ][ j ];
00460 real_type z01 = h[ i ][ j + 1 ];
00461 real_type z11 = h[ i + 1 ][ j + 1 ];
00462
00463
00464
00465 float nx = -dy * ( h[ i + 1 ][ j ] - h[ i ][ j ] );
00466 float ny = -( dx * ( h[ i ][ j + 1 ] - h[ i ][ j ] ) );
00467 float nz = dx * dy;
00468 glNormal3f( nx, ny, nz );
00469 glVertex3f( x, y, z00 );
00470 glNormal3f( nx, ny, nz );
00471 glVertex3f( x + dx, y, z10 );
00472 glNormal3f( nx, ny, nz );
00473 glVertex3f( x + dx, y + dy, z11 );
00474 glNormal3f( nx, ny, nz );
00475 glVertex3f( x, y + dy, z01 );
00476 }
00477 }
00478 glEnd();
00479 if ( disableBlend )
00480 {
00481 glDisable( GL_BLEND );
00482 }
00483
00484 }
00485 }
00486
00487
00492 void run( const real_type time_step, bool statistics = false )
00493 {
00494 assert( time_step > 0 && "timestep must be positive");
00495
00496 OpenTissue::utility::Timer<double> timer;
00497 double timeSetup = 0;
00498 double timeAssemble = 0;
00499 double timeCG = 0;
00500 double timeUpdate = 0;
00501 timer.start();
00502
00503 this->timestep = time_step;
00504
00505 compute_departure_points();
00506 compute_value_at_departure_points( h, h_tilde );
00507 compute_value_at_departure_points( u, u_tilde );
00508 compute_value_at_departure_points( v, v_tilde );
00509 compute_depth();
00510
00511 timer.stop();
00512 timeSetup = timer();
00513 timer.start();
00514
00515 assemble( A, rhs );
00516
00517 timer.stop();
00518 timeAssemble = timer();
00519 timer.start();
00520
00521 math::big::conjugate_gradient(A, hnew, rhs);
00522
00523 timer.stop();
00524 timeCG = timer();
00525 timer.start();
00526
00527 set_height( hnew );
00528 update_u();
00529 update_v();
00530
00531 apply_velocity_constraints();
00532 nullify();
00533
00534 timer.stop();
00535 timeUpdate = timer();
00536
00537 if (statistics)
00538 {
00539 std::cout << "CFD_ShallowWaterEquations::run(" << timestep << ") took (in milisecs.):" << std::endl;
00540 std::cout << " Setup : " << timeSetup << std::endl;
00541 std::cout << " Assemble : " << timeAssemble << std::endl;
00542 std::cout << " ConjGrad : " << timeCG << std::endl;
00543 std::cout << " Update : " << timeUpdate << std::endl;
00544 }
00545 }
00546
00556 void setSeaVelocity( const unsigned int i, const unsigned int j, const real_type vx, const real_type vy )
00557 {
00558 assert( i < X && "index out of bounds");
00559 assert( j < Y && "index out of bounds");
00560 u[ i ][ j ] = vx;
00561 v[ i ][ j ] = vy;
00562 }
00563
00572 void setSeaHeight( const unsigned int i, const unsigned int j, const real_type height )
00573 {
00574 assert( i < X && "index out of bounds" );
00575 assert( j < Y && "index out of bounds" );
00576 h[ i ][ j ] = height;
00577 if ( h[ i ][ j ] < b[ i ][ j ] )
00578 b[ i ][ j ] = h[ i ][ j ];
00579 }
00580
00589 void setSeaBottom( const unsigned int i, const unsigned int j, const real_type bottom )
00590 {
00591 assert( i < X && "index out of bounds" );
00592 assert( j < Y && "index out of bounds" );
00593 b[ i ][ j ] = bottom;
00594 if ( b[ i ][ j ] > h[ i ][ j ] )
00595 h[ i ][ j ] = b[ i ][ j ];
00596 }
00597
00606 void setShore( const unsigned int i, const unsigned int j, const bool shore )
00607 {
00608 assert( i < X && "index out of bounds" );
00609 assert( j < Y && "index out of bounds" );
00610 water[ i ][ j ] = !shore;
00611 }
00612
00613 protected:
00614
00615 void set_height( vector_type & h_new )
00616 {
00617 typename vector_type::iterator hval = h_new.begin();
00618 for ( index_type x = 0; x < X; ++x )
00619 {
00620 for ( index_type y = 0; y < Y; ++y )
00621 {
00622 h[ x ][ y ] = *hval++;
00623
00624
00625
00626
00627
00628 if ( h[ x ][ y ] < b[ x ][ y ] )
00629 h[ x ][ y ] = b[ x ][ y ];
00630 }
00631 }
00632 }
00633
00634
00635 void assemble( matrix_type & A, vector_type & rhs ) const
00636 {
00637 A.clear();
00638
00639
00640
00641
00642
00643
00644 real_type dxdx = dx * dx;
00645 real_type dydy = dy * dy;
00646 real_type dt = timestep;
00647 real_type dtdt = dt * dt;
00648 real_type dt_2dx = dt / ( 2 * dx );
00649 real_type dt_2dy = dt / ( 2 * dy );
00650 real_type dtdtgdxdx = ( dtdt * g ) / ( dxdx );
00651 real_type dtdtgdydy = ( dtdt * g ) / ( dydy );
00652 real_type dtdtgdxdx_4 = 0.25 * dtdtgdxdx;
00653 real_type dtdtgdydy_4 = 0.25 * dtdtgdydy;
00654
00655 for ( index_type i = 0;i < X;++i )
00656 {
00657 for ( index_type j = 0;j < Y;++j )
00658 {
00659 unsigned int row = i * Y + j;
00660
00661 unsigned int idx_ip1_j = ( i + 1 ) * Y + j;
00662 unsigned int idx_i_j = i * Y + j;
00663 unsigned int idx_i_jp1 = i * Y + ( j + 1 );
00664 unsigned int idx_im1_j = ( i - 1 ) * Y + j;
00665
00666
00667
00668 unsigned int idx_i_jm1 = i * Y + ( j - 1 );
00669
00670 A( row, idx_i_j ) += 1;
00671
00672 real_type dBi = 0;
00673 if ( ( i + 1 ) < X )
00674 dBi = b[ i + 1 ][ j ];
00675 else
00676 dBi = b[ i ][ j ];
00677 if ( i >= 1 )
00678 dBi -= b[ i - 1 ][ j ];
00679 else
00680 dBi -= b[ i ][ j ];
00681
00682
00683 if ( dBi )
00684 {
00685
00686 real_type k = dtdtgdxdx_4 * dBi;
00687 if ( ( i + 1 ) < X )
00688 A( row, idx_ip1_j ) += k;
00689 else
00690 A( row, idx_i_j ) += 2 * k;
00691
00692 if ( i >= 1 )
00693 A( row, idx_im1_j ) -= k;
00694 else
00695 A( row, idx_i_j ) -= 2 * k;
00696 }
00697
00698 real_type dBj = 0;
00699 if ( ( j + 1 ) < Y )
00700 dBj = b[ i ][ j + 1 ];
00701 else
00702 dBj = b[ i ][ j ];
00703 if ( j >= 1 )
00704 dBj -= b[ i ][ j - 1 ];
00705 else
00706 dBj -= b[ i ][ j ];
00707
00708
00709 if ( dBj )
00710 {
00711
00712 real_type h = dtdtgdydy_4 * dBj;
00713 if ( ( j + 1 ) < Y )
00714 A( row, idx_i_jp1 ) += h;
00715 else
00716 A( row, idx_i_j ) += 2 * h;
00717
00718 if ( j >= 1 )
00719 A( row, idx_i_jm1 ) -= h;
00720 else
00721 A( row, idx_i_j ) -= 2 * h;
00722 }
00723
00724
00725 if ( d[ i ][ j ] )
00726 {
00727
00728 real_type w = dtdtgdxdx * d[ i ][ j ];
00729 if ( i >= 1 )
00730 A( row, idx_im1_j ) -= w;
00731 else
00732 A( row, idx_i_j ) -= w;
00733 A( row, idx_i_j ) += 2 * w;
00734 if ( ( i + 1 ) < X )
00735 A( row, idx_ip1_j ) -= w;
00736 else
00737 A( row, idx_i_j ) -= w;
00738
00739
00740 real_type q = dtdtgdydy * d[ i ][ j ];
00741 if ( j >= 1 )
00742 A( row, idx_i_jm1 ) -= q;
00743 else
00744 A( row, idx_i_j ) -= q;
00745 A( row, idx_i_j ) += 2 * q;
00746 if ( ( j + 1 ) < Y )
00747 A( row, idx_i_jp1 ) -= q;
00748 else
00749 A( row, idx_i_j ) -= q;
00750 }
00751
00752
00753 rhs( row ) = h_tilde[ i ][ j ] + u_tilde[ i ][ j ] * dBi * dt_2dx + v_tilde[ i ][ j ] * dBj * dt_2dy;
00754
00755 real_type dU = 0;
00756 if ( ( i + 1 ) < X )
00757 dU = u_tilde[ i + 1 ][ j ];
00758 if ( i >= 1 )
00759 dU -= u_tilde[ i - 1 ][ j ];
00760
00761 real_type dV = 0;
00762 if ( ( j + 1 ) < Y )
00763 dV = v_tilde[ i ][ j + 1 ];
00764 if ( j >= 1 )
00765 dV -= v_tilde[ i ][ j - 1 ];
00766
00767
00768 rhs( row ) -= d[ i ][ j ] * ( dU * dt_2dx + dV * dt_2dy );
00769 }
00770 }
00771 }
00772
00773
00777 void compute_depth()
00778 {
00779 static real_type epsilon = math::working_precision<real_type>();
00780
00781 for ( unsigned int i = 0; i < X; ++i )
00782 {
00783 for ( unsigned int j = 0; j < Y; ++j )
00784 {
00785 d[ i ][ j ] = h[ i ][ j ] - b[ i ][ j ];
00786
00787
00788
00789
00790
00791 if ( d[ i ][ j ] < epsilon )
00792 d[ i ][ j ] = 0;
00793 }
00794 }
00795 }
00796
00800 void compute_departure_points()
00801 {
00802 real_type x = 0;
00803 for ( unsigned int i = 0;i < X;++i )
00804 {
00805 real_type y = 0;
00806 for ( unsigned int j = 0;j < Y;++j )
00807 {
00808 x_tilde[ i ][ j ] = x - timestep * u[ i ][ j ];
00809 y_tilde[ i ][ j ] = y - timestep * v[ i ][ j ];
00810 y += dy;
00811 }
00812 x += dx;
00813 }
00814 }
00815
00816
00821 void compute_value_at_departure_points( array_type & value, array_type & value_tilde )
00822 {
00823 using std::floor;
00824
00825 for ( unsigned int i = 0; i < X; ++i )
00826 {
00827 for ( unsigned int j = 0; j < Y; ++j )
00828 {
00829 real_type tmp_x = x_tilde[ i ][ j ] / dx;
00830 real_type tmp_y = y_tilde[ i ][ j ] / dy;
00831 int i00 = static_cast<int>( floor( tmp_x ) );
00832 int j00 = static_cast<int>( floor( tmp_y ) );
00833 real_type x_low = i00 * dx;
00834 real_type y_low = j00 * dy;
00835 real_type x_high = x_low + dx;
00836 real_type y_high = y_low + dy;
00837
00838 real_type val00 = 0;
00839 real_type val10 = 0;
00840 real_type val01 = 0;
00841 real_type val11 = 0;
00842
00843 if ( i00 >= 0 && j00 >= 0 && static_cast<unsigned int>(i00) < X && static_cast<unsigned int>(j00) < Y )
00844 val00 = value[ i00 ][ j00 ];
00845 if ( ( i00 + 1 ) >= 0 && j00 >= 0 && static_cast<unsigned int>( i00 + 1 ) < X && static_cast<unsigned int>(j00) < Y )
00846 val10 = value[ i00 + 1 ][ j00 ];
00847 if ( i00 >= 0 && ( j00 + 1 ) >= 0 && static_cast<unsigned int>(i00) < X && static_cast<unsigned int>( j00 + 1 ) < Y )
00848 val01 = value[ i00 ][ j00 + 1 ];
00849 if ( ( i00 + 1 ) >= 0 && ( j00 + 1 ) >= 0 && static_cast<unsigned int>( i00 + 1 ) < X && static_cast<unsigned int>( j00 + 1 ) < Y )
00850 val11 = value[ i00 + 1 ][ j00 + 1 ];
00851
00852 real_type s1 = x_high - x_tilde[ i ][ j ];
00853 real_type s0 = x_tilde[ i ][ j ] - x_low;
00854
00855 real_type tmp0 = ( s0 * val10 + s1 * val00 ) / dx;
00856 real_type tmp1 = ( s0 * val11 + s1 * val01 ) / dx;
00857
00858 real_type t1 = y_high - y_tilde[ i ][ j ];
00859 real_type t0 = y_tilde[ i ][ j ] - y_low;
00860
00861 value_tilde[ i ][ j ] = ( t0 * tmp1 + t1 * tmp0 ) / dy;
00862 }
00863 }
00864 }
00865
00869 void update_u()
00870 {
00871 for ( unsigned int i = 0; i < X; ++i )
00872 {
00873 for ( unsigned int j = 0; j < Y; ++j )
00874 {
00875 real_type Dh = 0;
00876 if ( ( i + 1 ) < X && water[ i + 1 ][ j ] )
00877 Dh = h[ i + 1 ][ j ];
00878 else
00879 Dh = h[ i ][ j ];
00880 if ( i >= 1 && water[ i - 1 ][ j ] )
00881 Dh -= h[ i - 1 ][ j ];
00882 else
00883 Dh -= h[ i ][ j ];
00884 Dh /= ( 2 * dx );
00885 u[ i ][ j ] = u_tilde[ i ][ j ] - timestep * g * Dh;
00886 }
00887 }
00888 }
00889
00893 void update_v()
00894 {
00895 for ( unsigned int i = 0; i < X; ++i )
00896 {
00897 for ( unsigned int j = 0; j < Y; ++j )
00898 {
00899 real_type Dh = 0;
00900 if ( ( j + 1 ) < Y && water[ i ][ j + 1 ] )
00901 Dh = h[ i ][ j + 1 ];
00902 else
00903 Dh = h[ i ][ j ];
00904 if ( j >= 1 && water[ i ][ j - 1 ] )
00905 Dh -= h[ i ][ j - 1 ];
00906 else
00907 Dh -= h[ i ][ j ];
00908 Dh /= ( 2 * dy );
00909 v[ i ][ j ] = v_tilde[ i ][ j ] - timestep * g * Dh;
00910 }
00911 }
00912 }
00913
00918 void apply_velocity_constraints()
00919 {
00920 for ( unsigned int i = 0;i < X;++i )
00921 {
00922 for ( unsigned int j = 0;j < Y;++j )
00923 {
00924 if ( water[ i ][ j ] )
00925 {
00926
00927 bool wall_left = false;
00928 bool wall_right = false;
00929 bool wall_up = false;
00930 bool wall_down = false;
00931 bool wall_left_up = false;
00932 bool wall_right_up = false;
00933 bool wall_left_down = false;
00934 bool wall_right_down = false;
00935
00936 if ( i >= 1 )
00937 wall_left = ( water[ i - 1 ][ j ] == 0 );
00938 if ( ( i + 1 ) < X )
00939 wall_right = ( water[ i + 1 ][ j ] == 0 );
00940 if ( ( j + 1 ) < Y )
00941 wall_up = ( water[ i ][ j + 1 ] == 0 );
00942 if ( j >= 1 )
00943 wall_down = ( water[ i ][ j - 1 ] == 0 );
00944 if ( i >= 1 && ( j + 1 ) < Y )
00945 wall_left_up = ( water[ i - 1 ][ j + 1 ] == 0 );
00946 if ( ( i + 1 ) < X && ( j + 1 ) < Y )
00947 wall_right_up = ( water[ i + 1 ][ j + 1 ] == 0 );
00948 if ( i >= 1 && j >= 1 )
00949 wall_left_down = ( water[ i - 1 ][ j - 1 ] == 0 );
00950 if ( ( i + 1 ) < X && j >= 1 )
00951 wall_right_down = ( water[ i + 1 ][ j - 1 ] == 0 );
00952
00953
00954 if ( wall_left || wall_right )
00955 u[ i ][ j ] = 0;
00956
00957 if ( wall_up || wall_down )
00958 v[ i ][ j ] = 0;
00959
00960
00961 if ( wall_left_down && !wall_left && !wall_down )
00962 {
00963 real_type nx = 1;
00964 real_type ny = 1;
00965 real_type ux = u[ i ][ j ];
00966 real_type uy = v[ i ][ j ];
00967 ux = ux - ( ux * nx + uy * ny ) * nx;
00968 uy = uy - ( ux * nx + uy * ny ) * ny;
00969 u[ i ][ j ] = ux;
00970 v[ i ][ j ] = uy;
00971
00972 }
00973 if ( wall_right_down && !wall_right && !wall_down )
00974 {
00975 real_type nx = -1;
00976 real_type ny = 1;
00977 real_type ux = u[ i ][ j ];
00978 real_type uy = v[ i ][ j ];
00979 ux = ux - ( ux * nx + uy * ny ) * nx;
00980 uy = uy - ( ux * nx + uy * ny ) * ny;
00981 u[ i ][ j ] = ux;
00982 v[ i ][ j ] = uy;
00983 }
00984 if ( wall_left_up && !wall_left && !wall_up )
00985 {
00986 real_type nx = 1;
00987 real_type ny = -1;
00988 real_type ux = u[ i ][ j ];
00989 real_type uy = v[ i ][ j ];
00990 ux = ux - ( ux * nx + uy * ny ) * nx;
00991 uy = uy - ( ux * nx + uy * ny ) * ny;
00992 u[ i ][ j ] = ux;
00993 v[ i ][ j ] = uy;
00994 }
00995 if ( wall_right_up && !wall_right && !wall_up )
00996 {
00997 real_type nx = -1;
00998 real_type ny = -1;
00999 real_type ux = u[ i ][ j ];
01000 real_type uy = v[ i ][ j ];
01001 ux = ux - ( ux * nx + uy * ny ) * nx;
01002 uy = uy - ( ux * nx + uy * ny ) * ny;
01003 u[ i ][ j ] = ux;
01004 v[ i ][ j ] = uy;
01005 }
01006 }
01007 else
01008 {
01009 u[ i ][ j ] = 0;
01010 v[ i ][ j ] = 0;
01011 }
01012
01013 }
01014 }
01015 }
01016
01017
01018 void nullify()
01019 {
01020 real_type threshold = 10e-5;
01021
01022 for ( unsigned int i = 0;i < X;++i )
01023 {
01024 for ( unsigned int j = 0;j < Y;++j )
01025 {
01026
01027 if ( d[ i ][ j ] < threshold )
01028 {
01029 bool left = true;
01030 bool right = true;
01031 bool up = true;
01032 bool down = true;
01033 bool left_up = true;
01034 bool right_up = true;
01035 bool left_down = true;
01036 bool right_down = true;
01037
01038 if ( i >= 1 )
01039 left = ( d[ ( i - 1 ) ][ j ] < threshold );
01040 if ( ( i + 1 ) < X )
01041 right = ( d[ ( i + 1 ) ][ j ] < threshold );
01042 if ( ( j + 1 ) < Y )
01043 up = ( d[ i ][ ( j + 1 ) ] < threshold );
01044 if ( j >= 1 )
01045 down = ( d[ i ][ ( j - 1 ) ] < threshold );
01046 if ( i >= 1 && ( j + 1 ) < Y )
01047 left_up = ( d[ ( i - 1 ) ][ ( j + 1 ) ] < threshold );
01048 if ( ( i + 1 ) < X && ( j + 1 ) < Y )
01049 right_up = ( d[ ( i + 1 ) ][ ( j + 1 ) ] < threshold );
01050 if ( i >= 1 && j >= 1 )
01051 left_down = ( d[ ( i - 1 ) ][ ( j - 1 ) ] < threshold );
01052 if ( ( i + 1 ) < X && j >= 1 )
01053 right_down = ( d[ ( i + 1 ) ][ ( j - 1 ) ] < threshold );
01054
01055 if ( left && right && up && down && left_down && right_down && left_up && right_up )
01056 {
01057 u[ i ][ j ] = 0;
01058 v[ i ][ j ] = 0;
01059 }
01060 }
01061 }
01062 }
01063 }
01064
01065 };
01066
01067 }
01068 }
01069
01070
01071 #endif
