added OpenGL Visualisation for SWE - all classes

src/opengl/camera.cpp

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include "camera.h"
+#include <sstream>
+/**
+    Constructor
+
+	@param view_distance	initial view distance from the origin 
+	@param window_title		title of the current window
+
+*/
+Camera::Camera(float view_distance, const char* window_title) {
+	// Initialize member variables
+	cameraX = -0.3f;
+	cameraY = 1.0f;
+	cameraZ = 2.0f;
+	objectX = 0.0f;
+	objectY = 0.0f;
+	objectZ = 0.0f;
+	angleX = 0.0f;
+	angleY = 0.0f;
+	zoomfactor = 0.85f*view_distance;
+	win_title = window_title;
+	frames = 0;
+	lastTime = 0;
+	oldMouseX = 0;
+	oldMouseY = 0;
+	
+	// Reset framebuffer
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+}
+/**
+    Set the camera via gluLookAt and set the light position afterwards
+
+*/
+void Camera::setCamera() {
+	// Clear our buffers
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+	
+	// Set camera position, center of scene and up vector 
+	gluLookAt(cameraX*zoomfactor, cameraY*zoomfactor, cameraZ*zoomfactor, 
+		0,0,0,   // center of scene
+		0, 1, 0); // up vector
+	
+	// Translate and rotate our object
+	glTranslated(objectX, objectY, objectZ);
+	this->rotateObject();
+	glScalef(1.0f, 1.0f, -1.0f);
+	
+	// Position the light source
+	// Note: Coordinates are multiplied by current modelview-matrix
+	// by OpenGL
+	GLfloat LightPosition[]=	{ 1.0f, 1.0f, 0.0f, 0.0f };
+	glLightfv(GL_LIGHT0, GL_POSITION,LightPosition);	
+}
+
+/**
+    Rotate camera 
+
+*/
+void Camera::rotateObject() {
+	glRotatef(angleX, 1.0, 0.0f, 0.0f);
+	glRotatef(angleY, 0.0f, 1.0f, 0.0f);
+}
+
+/**
+    Increment viewing orientation of the camera
+
+	@param angle_dX		angle relative to the x-axis
+	@param angle_dY		angle relative to the rotated y-axis
+
+*/
+void Camera::orient( float angle_dX, float angle_dY ) 
+{
+	angleX += angle_dX;
+	angleY += angle_dY;
+}
+
+/**
+    Zoom in
+
+	@param scaleFactor		factor which is used for zooming			
+
+*/
+void Camera::zoomIn( float scaleFactor ) 
+{
+	if (zoomfactor > 3.0f) {
+		zoomfactor /= scaleFactor; 
+	}	  
+}
+/**
+    Zoom out
+
+	@param scaleFactor		factor which is used for zooming			
+
+*/
+void Camera::zoomOut( float scaleFactor ) 
+{
+	zoomfactor *= scaleFactor;		  
+}
+
+/**
+    Calculates the current framerate, updates the window title and
+	swaps framebuffers to display the new image
+
+*/
+void Camera::displayImage() {
+    // Set up stringstream for window title
+	std::stringstream ss;
+	ss.setf(std::ios::fixed, std::ios::floatfield);
+	ss.setf(std::ios::showpoint);
+	ss.precision(1);
+	
+	// Calculate framerate
+	float fps;
+	frames++;
+	if (lastTime == 0) {
+		lastTime = SDL_GetTicks();
+	} else if( (SDL_GetTicks() - lastTime) >= 500) {
+		fps = (float)frames/(SDL_GetTicks() - lastTime)*1000.0f;
+		ss << win_title << " (" << fps << " fps)";
+		SDL_WM_SetCaption(ss.str().c_str(), NULL);
+
+		frames = 0;
+        lastTime = SDL_GetTicks();
+    }
+
+	// Swap frame buffer
+    SDL_GL_SwapBuffers();
+}
+
+/**
+    User starts dragging. 
+	Remember the old mouse coordinates.	
+
+*/
+void Camera::startPanning(int xPos, int yPos) {
+	oldMouseX = xPos;
+	oldMouseY = yPos;
+}
+
+/**
+    User drags our object. 
+	Transform screen coordinates into world coordinates
+	and update the objects position
+
+*/
+void Camera::panning(int newX, int newY) {
+	GLint viewport[4];
+	GLdouble modelview[16];
+	GLdouble projection[16];
+	GLfloat winX, winY, winZ;
+	GLdouble posX, posY, posZ;
+	GLdouble pos2X, pos2Y, pos2Z;
+	GLdouble dX, dY, dZ;
+	
+	// Draw invisible fake plane
+	setCamera();
+	glColorMask( GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE );
+	GLfloat dim = 100000.0f;
+	glBegin(GL_QUADS);	
+		glColor3f(1.0f, 1.0f, 0.3f);
+		glVertex3f( dim,0.0f, dim);
+		glVertex3f(-dim,0.0f, dim);		
+		glVertex3f(-dim,0.0f, -dim);	
+		glVertex3f( dim,0.0f, -dim);				
+	glEnd();
+	glColorMask(GL_TRUE,GL_TRUE,GL_TRUE,GL_TRUE); 	
+
+    // Get projection matrices
+	glGetDoublev( GL_MODELVIEW_MATRIX, modelview );
+	glGetDoublev( GL_PROJECTION_MATRIX, projection );
+	glGetIntegerv( GL_VIEWPORT, viewport );
+
+	// Get drag positions in object coordinates
+	winX = (GLfloat)newX;
+	winY = (GLfloat)viewport[3] - (GLfloat)newY;
+	glReadPixels( newX, GLint(winY), 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ );
+	gluUnProject( winX, winY,  winZ, modelview, projection, viewport, &posX, &posY, &posZ);
+	
+	winX = (GLfloat)oldMouseX;
+	winY = (GLfloat)viewport[3] - (GLfloat)oldMouseY;
+	glReadPixels( oldMouseX, GLint(winY), 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ );
+	gluUnProject( winX, winY,  winZ, modelview, projection, viewport, &pos2X, &pos2Y, &pos2Z);
+	
+	glClear(GL_DEPTH_BUFFER_BIT);
+	
+	// Calculate drag distances
+	dX = (posX-pos2X);
+	dY = (posY-pos2Y);
+	dZ = (posZ-pos2Z);
+
+	// Convert drag distances to world coordinates
+	glPushMatrix();
+		glLoadIdentity();
+		rotateObject();
+		glScalef(1.0f, 1.0f, -1.0f);
+		glGetDoublev( GL_MODELVIEW_MATRIX, modelview );
+		GLdouble x = modelview[0]*dX + modelview[4]*dY + modelview[8]*dZ; 
+		GLdouble y = modelview[1]*dX + modelview[5]*dY + modelview[9]*dZ; 
+		GLdouble z = modelview[2]*dX + modelview[6]*dY + modelview[10]*dZ; 
+	glPopMatrix();
+	
+	// Set the new world coordinates of our object
+	objectX += x;
+	objectY += y;
+	objectZ += z;
+
+	// Save our current mouse position
+	oldMouseX = newX;
+	oldMouseY = newY;
+}
\ No newline at end of file

src/opengl/camera.h

+#ifndef CAMERA_H
+#define CAMERA_H
+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include <SDL.h>
+#include <SDL_opengl.h> 
+
+class Camera {
+public:
+	Camera(float view_distance, const char* window_title);
+	// Set modelview matrix
+	void setCamera ();
+
+	// Change viewing properties
+	void orient( float angX, float angY );
+	void zoomIn( float scaleFactor );
+	void zoomOut( float scaleFactor );
+	void startPanning(int xPos, int yPos);
+	void panning(int newX, int newY);
+
+	// Update framebuffer
+	void displayImage();
+
+private:
+	// Position of the camera
+	float cameraX;
+	float cameraY;
+	float cameraZ;
+
+	// Position of the object
+	GLdouble objectX;
+	GLdouble objectY;
+	GLdouble objectZ;
+
+	// Zoom factor
+	float zoomfactor;
+	float angleX, angleY;
+
+	// Helper variables
+	unsigned int frames;
+	unsigned int lastTime;
+	unsigned int oldMouseX, oldMouseY;
+	unsigned int newMouseX, newMouseY;
+	
+	// Window title
+	const char* win_title;
+	void rotateObject();
+};
+
+#endif
\ No newline at end of file

src/opengl/controller.cpp

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include "controller.h"
+
+#include "../scenarios/SWE_simple_scenarios_vis.h"
+
+/**
+    Constructor
+
+	@param sim				instance of simulation class
+	@param vis				instance of visualization class
+
+*/	
+Controller::Controller(Simulation* sim, Visualization* vis) {
+	simulation = sim;
+	visualization = vis;
+	isActive = true;
+	done = false;
+	paused = false;
+	allowStep = false;
+}
+
+/**
+    Process all user events in a loop
+	Returns true, when user wants to quit
+
+*/	
+bool Controller::handleEvents() {
+	// Handle events
+	SDL_Event event;
+	allowStep = false;
+
+	while ( SDL_PollEvent(&event) ) 
+	{
+		switch( event.type )
+		{
+		case SDL_MOUSEMOTION:
+			if (event.motion.state & SDL_BUTTON(SDL_BUTTON_LEFT)) {
+				// Left mouse button dragging
+				visualization->camera->orient(0.5f*event.motion.yrel, 0.5f*event.motion.xrel);
+			}
+			else if (event.motion.state & SDL_BUTTON(SDL_BUTTON_RIGHT)) {
+				// Panning happens here
+				visualization->camera->panning(event.motion.x, event.motion.y);			
+			}
+			break;
+		case SDL_MOUSEBUTTONDOWN:
+			// Scroll wheel down
+			if ( event.button.button == SDL_BUTTON_WHEELDOWN ) {
+				// Zoom out
+				visualization->camera->zoomOut(1.2f);
+			} else if (event.button.button == SDL_BUTTON_RIGHT) {
+				visualization->camera->startPanning(event.button.x, event.button.y);
+			}
+			break;
+		case SDL_MOUSEBUTTONUP:
+			// Scroll wheel up
+			if ( event.button.button == SDL_BUTTON_WHEELUP ) {
+				// Zoom in
+				visualization->camera->zoomIn(1.15f);
+			}
+			break;
+		case SDL_ACTIVEEVENT:
+			// Don't draw anything if we're getting minimized
+			if ( event.active.state & SDL_APPACTIVE ) {
+				isActive = (event.active.gain != 0);	
+			}
+			break;			    
+		case SDL_VIDEORESIZE:
+			// Our window was resized
+			visualization->resizeWindow(event.resize.w, event.resize.h);
+			break;
+		case SDL_KEYDOWN:
+			// User has pressed a key
+			done = handleKeyPress( &event.key.keysym);
+			break;
+		case SDL_QUIT:
+			// Window close request
+			done = true;
+			break;
+		default:
+			break;
+		}
+	}
+	return done;
+}
+
+/**
+	Returns true, when window has focus 
+
+*/	
+bool Controller::hasFocus() {
+	return isActive && !done;
+}
+
+/**
+    Return whether program is currently paused
+
+*/	
+bool Controller::isPaused() {
+	return paused && !allowStep;
+}
+
+/**
+    Process single keyboard event
+	@param *keysym			pointer to the sdl keyevent structure
+
+*/	
+bool Controller::handleKeyPress( SDL_keysym *keysym) {
+
+	switch ( keysym->sym ) {
+		case SDLK_ESCAPE:
+			// ESC key was pressed -> exit
+			return true;
+			break;
+		case SDLK_r:
+			// Restart simulation
+			allowStep = paused;
+			simulation->restart();
+			break;
+		case SDLK_w:
+			// Switch rendering mode
+			visualization->toggleRenderingMode();
+			break;
+		case SDLK_s:
+			// Save simulation data to file
+			simulation->saveToFile();
+			break;
+		case SDLK_RIGHT:
+			// Advance single timestep when paused
+			allowStep = paused;
+			break;
+		case SDLK_SPACE:
+			// Pause/Resume
+			paused = !paused;
+			break;
+		case SDLK_1:
+			// Load scenario 1
+			{
+			  allowStep = paused;
+                          SWE_RadialDamBreakScenarioVisInfo* 
+                            newScene = new SWE_RadialDamBreakScenarioVisInfo();
+			  SWE_VisInfo* visInfo = newScene;
+
+			  // define grid size and initial time step
+			  float dx = (newScene->getBoundaryPos(BND_RIGHT) - newScene->getBoundaryPos(BND_LEFT) )/SWE_Block::getNx();
+			  float dy = (newScene->getBoundaryPos(BND_TOP) - newScene->getBoundaryPos(BND_BOTTOM) )/SWE_Block::getNy();
+			  SWE_Block::initGridData(SWE_Block::getNx(),SWE_Block::getNy(),dx,dy);
+
+			  simulation->loadNewScenario(newScene, visInfo);
+			  visualization->updateBathymetryVBO(simulation);
+			}
+			break;
+		case SDLK_2:
+			// Load scenario 2
+			{
+			  allowStep = paused;
+			  SWE_Scenario* newScene = new SWE_BathymetryDamBreakScenario;
+			  // define grid size and initial time step
+			  float dx = (newScene->getBoundaryPos(BND_RIGHT) - newScene->getBoundaryPos(BND_LEFT) )/SWE_Block::getNx();
+			  float dy = (newScene->getBoundaryPos(BND_TOP) - newScene->getBoundaryPos(BND_BOTTOM) )/SWE_Block::getNy();
+			  SWE_Block::initGridData(SWE_Block::getNx(),SWE_Block::getNy(),dx,dy);
+
+			  simulation->loadNewScenario(newScene, NULL);
+			  visualization->updateBathymetryVBO(simulation);
+			}
+			break;
+		case SDLK_3:
+			// Load scenario 3
+			{
+			  allowStep = paused;
+			  SWE_SplashingPoolScenarioVisInfo* 
+                            newScene = new SWE_SplashingPoolScenarioVisInfo;
+
+			  // define grid size and initial time step
+			  float dx = (newScene->getBoundaryPos(BND_RIGHT) - newScene->getBoundaryPos(BND_LEFT) )/SWE_Block::getNx();
+			  float dy = (newScene->getBoundaryPos(BND_TOP) - newScene->getBoundaryPos(BND_BOTTOM) )/SWE_Block::getNy();
+			  SWE_Block::initGridData(SWE_Block::getNx(),SWE_Block::getNy(),dx,dy);
+
+			  simulation->loadNewScenario(newScene, newScene);
+			  visualization->updateBathymetryVBO(simulation);
+			}
+			break;
+		default:
+			break;
+	}
+	return false;
+}

src/opengl/controller.h

+#ifndef CONTROLLER_H
+#define CONTROLLER_H
+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include <SDL.h>
+#include "simulation.h"
+#include "visualization.h"
+
+class Controller {
+public:
+	Controller(Simulation* sim, Visualization* vis);
+	
+	// Process new events
+	bool handleEvents();
+
+	// Return if window is enabled
+	bool hasFocus();
+	// Return if program is paused
+	bool isPaused();
+private:
+	// Internal functions
+	bool isActive;
+	bool done;
+	bool paused;
+	bool allowStep;
+
+	// References to other classes
+	Simulation* simulation;
+	Visualization* visualization;
+	
+	// Handle keyboard events
+	bool handleKeyPress( SDL_keysym *keysym);
+};
+
+#endif
\ No newline at end of file

src/opengl/fragment.glsl

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+
+varying vec3 V;
+varying vec3 N;
+varying vec4 ambient;
+
+/* Shading of water surfaces with fresnel effect */
+
+void main()
+{
+	// Helper variable for color
+	vec4 color = ambient;
+	
+	// Get normal and view vector
+	vec3 view = normalize(V);
+	vec3 normal = normalize(N); ;
+	if (!gl_FrontFacing) {
+		// Handle backface rendering
+		normal = -normal;
+	}
+	
+	// Compute fresnel factor
+	float fresnelFactor = 0.1+0.1*pow(1+dot(normal,view), 3.0);
+		
+	// Combine color with fresnel factor
+	color += gl_LightSource[0].specular * fresnelFactor;
+	color = mix(color, vec4(0.2, 0.4, 0.9, 1.0), 0.2);
+    
+    // Set constant transparency 
+	gl_FragColor = vec4(color.xyz, 0.7);
+}

src/opengl/shader.cpp

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+
+#include "shader.h"
+
+/**
+    Constructor.
+	Check whether shaders are supported. If yes, load vertex and fragment 
+	shader from textfile into memory and compile
+
+	@param  vertexShaderFile	name of the text file containing the vertex 
+								shader code
+	@param	fragmentShaderFile  name of the text file containing the fragment
+								shader code
+	
+*/	
+Shader::Shader(char const * vertexShaderFile,  char const * fragmentShaderFile) {
+	shdrSupport = isExtensionSupported("GL_ARB_vertex_shader") 
+					&& isExtensionSupported("GL_ARB_shader_objects")
+					&& isExtensionSupported("GL_ARB_fragment_shader");
+
+	shdrLoaded = false;
+
+	if (shdrSupport) {
+		// Load extensions
+		glCreateShader = (PFNGLCREATESHADERPROC) SDL_GL_GetProcAddress("glCreateShader");
+		glCreateProgram = (PFNGLCREATEPROGRAMPROC) SDL_GL_GetProcAddress("glCreateProgram");
+		glAttachShader = (PFNGLATTACHSHADERPROC)SDL_GL_GetProcAddress("glAttachShader");
+		glCompileShader = (PFNGLCOMPILESHADERPROC)SDL_GL_GetProcAddress("glCompileShader");
+		glUseProgram = (PFNGLUSEPROGRAMPROC)SDL_GL_GetProcAddress("glUseProgram");
+		glDetachShader = (PFNGLDETACHSHADERPROC)SDL_GL_GetProcAddress("glDetachShader");
+		glDeleteShader = (PFNGLDELETESHADERPROC) SDL_GL_GetProcAddress("glDeleteShader");
+		glLinkProgram = (PFNGLLINKPROGRAMPROC)SDL_GL_GetProcAddress("glLinkProgram");
+		glShaderSource = (PFNGLSHADERSOURCEPROC)SDL_GL_GetProcAddress("glShaderSource");
+		glDeleteProgram = (PFNGLDELETEPROGRAMPROC)SDL_GL_GetProcAddress("glDeleteProgram");
+		glGetObjectParameterivARB = (PFNGLGETOBJECTPARAMETERIVARBPROC) SDL_GL_GetProcAddress("glGetObjectParameterivARB");
+		glGetShaderiv = (PFNGLGETSHADERIVPROC) SDL_GL_GetProcAddress("glGetShaderiv");
+		glGetShaderInfoLog = (PFNGLGETSHADERINFOLOGPROC) SDL_GL_GetProcAddress("glGetShaderInfoLog");
+		glGetProgramInfoLog = (PFNGLGETPROGRAMINFOLOGPROC) SDL_GL_GetProcAddress("glGetProgramInfoLog");
+
+		// Read shader files
+		bool readSuccess = false;
+		vertexShaderSource = NULL;
+		fragmentShaderSource = NULL;
+		readSuccess = readShaderFile(vertexShaderFile, vertexShaderSource, vertexShaderLength);
+		readSuccess &= readShaderFile(fragmentShaderFile, fragmentShaderSource, fragmentShaderLength);
+
+		// Check if reading was completed
+		if (readSuccess) {
+			// Create Vertex-Shader 
+			vertexShader = glCreateShader(GL_VERTEX_SHADER);
+			glShaderSource(vertexShader, 1, const_cast<GLchar const **>(&vertexShaderSource), &vertexShaderLength);
+			glCompileShader(vertexShader);
+
+			// Create Fragment-Shader 
+			fragmentShader = glCreateShader(GL_FRAGMENT_SHADER);
+			glShaderSource(fragmentShader, 1, const_cast<GLchar const **>(&fragmentShaderSource), &fragmentShaderLength);
+			glCompileShader(fragmentShader);
+
+			if (isShaderCompiled(vertexShader, "Vertex-Shader")
+				&& isShaderCompiled(fragmentShader, "Fragment-Shader"))
+			{
+				// Create Shader program
+				program = glCreateProgram();
+				glAttachShader(program, vertexShader);
+				glAttachShader(program, fragmentShader);
+				glLinkProgram(program);
+				shdrLoaded = true;
+			} else {
+				// Errors while compiling shaders
+				shdrSupport = false;
+				glDeleteShader(vertexShader);
+				delete[] vertexShaderSource;
+				glDeleteShader(fragmentShader);
+				delete[] fragmentShaderSource;
+			}
+		} else {
+			if (vertexShaderSource != NULL) delete[] vertexShaderSource;
+			if (fragmentShaderSource != NULL) delete[] fragmentShaderSource;
+		}
+	}
+}
+/**
+    Destructor.
+	Unload shaders and free resources.
+	
+*/	
+Shader::~Shader() {
+	// Unload Shaders
+	if (shdrLoaded) {
+		glUseProgram(0);
+		glDeleteProgram(program);
+
+		glDetachShader(program, vertexShader);
+		glDeleteShader(vertexShader);
+		delete[] vertexShaderSource;
+
+		glDetachShader(program, fragmentShader);
+		glDeleteShader(fragmentShader);
+		delete[] fragmentShaderSource;
+	}
+}
+/**
+    Replaces OpenGL shaders by our custom shaders
+
+*/	
+void Shader::enableShader() {
+	if (shdrLoaded) {
+		glUseProgram(program);
+	}
+}
+/**
+    Restores OpenGL default shaders 
+
+*/
+void Shader::disableShader() {
+	if (shdrLoaded) {
+		glUseProgram(0);
+	}
+}
+
+/**
+	Returns, whether shaders are supported by graphics hardware
+*/
+bool Shader::shadersSupported() {
+	return shdrSupport;
+}
+
+/**
+    Returns, whether shaders could by loaded successfully 
+
+*/
+bool Shader::shadersLoaded() {
+	return shdrLoaded;
+}
+
+/**
+    Returns, whether a special extension is supported by the current 
+	graphics card
+
+    @param	szTargetExtention	string describing the extension to look for	
+
+*/	
+bool Shader::isExtensionSupported(const char* szTargetExtension )
+{
+	const unsigned char *pszExtensions = NULL;
+	const unsigned char *pszStart;
+	unsigned char *pszWhere, *pszTerminator;
+	
+	// Extension names should not have spaces
+	pszWhere = (unsigned char *) strchr( szTargetExtension, ' ' );
+	if( pszWhere || *szTargetExtension == '\0' )
+		return false;
+
+	// Get Extensions String
+	pszExtensions = glGetString( GL_EXTENSIONS );
+
+	// Search The Extensions String For An Exact Copy
+	pszStart = pszExtensions;
+	for(;;)
+	{
+		pszWhere = (unsigned char *) strstr( (const char *) pszStart, szTargetExtension );
+		if( !pszWhere )
+			break;
+		pszTerminator = pszWhere + strlen( szTargetExtension );
+		if( pszWhere == pszStart || *( pszWhere - 1 ) == ' ' )
+			if( *pszTerminator == ' ' || *pszTerminator == '\0' )
+				return true;
+		pszStart = pszTerminator;
+	}
+	return false;
+}
+
+/**
+    Returns, whether a special extension is supported by the current 
+	graphics card
+
+    @param	filename		shader file
+	@return shaderSource 	file content
+	@return length			length of the file
+	@returns bool	true, if file has been read successfully
+
+*/	
+bool Shader::readShaderFile(char const * filename, GLchar * & shaderSource, GLint & length)
+{
+	using namespace std;
+	string fullFileName = "";
+
+	// On unix: determine directory of exec 
+#ifdef __unix__
+	// Code taken from: 
+	// http://www.gamedev.net/community/forums/topic.asp?topic_id=459511
+	std::string path = "";
+	pid_t pid = getpid();
+	char buf[20] = {0};
+	sprintf(buf,"%d",pid);
+	std::string _link = "/proc/";
+	_link.append( buf );
+	_link.append( "/exe");
+	char proc[512];
+	int ch = readlink(_link.c_str(),proc,512);
+	if (ch != -1) {
+		proc[ch] = 0;
+		path = proc;
+		std::string::size_type t = path.find_last_of("/");
+		path = path.substr(0,t);
+	}
+
+	fullFileName = path + string("/");
+#endif
+
+	fullFileName.append(filename);
+	ifstream file(fullFileName.c_str(), ios::in);
+
+	length = 0;
+
+	if (!file.good())
+	{
+		return false;
+	}
+
+	file.seekg(0, ios::end);
+	length = file.tellg();
+	file.seekg(ios::beg);
+
+	if (length == 0 )
+	{
+		return false;
+	}
+
+	shaderSource = new GLchar[length];
+
+	if (shaderSource == 0)
+	{
+		return false;
+	}
+
+	memset(shaderSource, 0, length);
+	file.read(shaderSource, length);
+
+	file.close();
+
+	return true;
+}
+
+/**
+    Returns, whether a shader has been compiled with success
+
+    @param	shader		shader id
+	@param prefix 		custom error message
+
+*/	
+bool Shader::isShaderCompiled(GLuint shader, char const * prefix)
+{
+	using namespace std;
+	GLint compiled(0);
+
+	glGetObjectParameterivARB(shader, GL_COMPILE_STATUS, &compiled);
+
+	if (!compiled)
+	{
+		GLint maxLength(0);
+
+		glGetShaderiv(shader, GL_INFO_LOG_LENGTH , &maxLength);
+
+		cerr << prefix
+			<< " -- Compiler Log: "
+			<< endl
+			;
+
+		if (maxLength > 1)
+		{
+			GLchar * log = new GLchar[maxLength];
+			glGetShaderInfoLog(shader, maxLength, 0, log);
+			cerr << log
+				;
+			delete[] log;
+		}
+	}
+
+	return (compiled != 0);
+}
+/**
+    Returns, whether a shader has been linked with success
+
+    @param	shader		shader id
+	@param prefix 		custom error message
+
+*/	
+bool Shader::isProgramLinked(GLuint program, char const * prefix)
+{
+	using namespace std;
+	GLint linked(0);
+
+	glGetObjectParameterivARB(program, GL_LINK_STATUS, &linked);
+
+	if (!linked)
+	{
+		GLint maxLength(0);
+
+		glGetShaderiv(program, GL_INFO_LOG_LENGTH , &maxLength);
+
+		cerr << prefix
+			<< " -- Linker Log: "
+			<< endl
+			;
+
+		if (maxLength > 1)
+		{
+			GLchar * log = new GLchar[maxLength];
+			glGetProgramInfoLog(program, maxLength, 0, log);
+			cerr << log
+				;
+			delete[] log;
+		}
+	}
+
+	return (linked != 0);
+}

src/opengl/shader.h

+#ifndef SHADER_H
+#define SHADER_H
+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include <SDL.h>
+#include <SDL_opengl.h>
+#include <iostream>
+#include <fstream>
+
+class Shader {
+public:
+	// Constructor and Destructor
+	// Load vertex shader file and fragment shader file 
+	Shader(char const * vertexShaderFile,  char const * fragmentShaderFile) ;
+	~Shader();
+
+	// Check if shaders are supported
+	bool shadersSupported();
+	bool shadersLoaded();
+
+	// Shader control
+	void enableShader();
+	void disableShader();
+
+private:	
+	// State flags
+	bool shdrSupport;
+	bool shdrLoaded;
+
+	// Helper functions
+	bool readShaderFile(char const * filename, GLchar * & shaderSource, GLint & length);
+	bool isExtensionSupported(const char* szTargetExtension );
+	bool isShaderCompiled(GLuint shader, char const * prefix);
+	bool isProgramLinked(GLuint program, char const * prefix);
+
+	// Vertex-Shader
+	GLuint   vertexShader;
+	GLchar * vertexShaderSource;
+	GLint    vertexShaderLength;
+	
+	// Fragment-Shader
+	GLuint   fragmentShader;
+	GLchar * fragmentShaderSource;
+	GLint    fragmentShaderLength;
+
+	// Shaders id
+	GLuint program;
+
+	// Shader extension function pointers
+	PFNGLCREATESHADERPROC glCreateShader;
+	PFNGLCREATEPROGRAMPROC glCreateProgram;
+	PFNGLATTACHSHADERPROC glAttachShader;
+	PFNGLCOMPILESHADERPROC glCompileShader;
+	PFNGLUSEPROGRAMPROC glUseProgram;
+	PFNGLDETACHSHADERPROC glDetachShader;
+	PFNGLDELETESHADERPROC glDeleteShader;
+	PFNGLLINKPROGRAMPROC glLinkProgram;
+	PFNGLSHADERSOURCEPROC glShaderSource;
+	PFNGLDELETEPROGRAMPROC glDeleteProgram;
+	
+	// Shader objects extension pointers
+	PFNGLGETOBJECTPARAMETERIVARBPROC glGetObjectParameterivARB;
+	PFNGLGETSHADERIVPROC glGetShaderiv;
+	PFNGLGETSHADERINFOLOGPROC glGetShaderInfoLog;
+	PFNGLGETPROGRAMINFOLOGPROC glGetProgramInfoLog;
+
+};
+#endif

src/opengl/simulation.cu

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Michael Bader, Kaveh Rahnema, Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include "simulation.h"
+#include <stdlib.h>
+#include "../SWE_BlockCUDA.hh"
+
+/**
+    Constructor. 
+	Initializes SWE_BlockCUDA and creates a new instance of it.
+
+    @param nx			number of grid cells in x-direction
+	@param ny			number of grid cells in y-direction
+	@param dx			size of one cell in x-direction
+	@param dy			size of one cell in y-direction
+	@param scene		instance of SWE_Scenario describing the boundaries
+
+*/
+Simulation::Simulation (int nx, int ny, float dx, float dy, 
+			SWE_Scenario* scene, SWE_VisInfo* visInfo, 
+                        SWE_BlockCUDA* _splash) 
+: myScenario(scene),
+  maxDim( (nx > ny) ? nx : ny ),
+  maxCellSize( (dx > dy) ? dx : dy ),
+  splash(_splash),
+  curTime(0.0f),
+  isFirstStep(1),
+  fileNumber(0),
+  useFileInput(false)
+{
+
+        Simulation::initBoundaries(myScenario);
+
+	if (visInfo != NULL) {
+		wAverage = visInfo->waterHeightAtRest();
+		wScale = visInfo->waterVerticalScaling()*(maxDim/20.0f);
+		wOffset = visInfo->waterDistanceFromGround()*(maxDim/20.0f);
+		bAverage = visInfo->bathyVerticalCenter();
+		bScale = visInfo->bathyVerticalScaling()*(maxDim/20.0f);
+		bOffset = visInfo->bathyDistanceFromGround()*(maxDim/20.0f);
+	} else {
+		getScalingApproximation(splash->getWaterHeight(), splash->getBathymetry());
+	}
+}
+/** 
+	Destructor. 
+	Delete current SWE_BlockCUDA and SWE_Scenario instance.
+*/
+Simulation::~Simulation () {
+}
+
+void Simulation::loadNewScenario(SWE_Scenario* scene, SWE_VisInfo* visInfo) {
+	delete myScenario;
+	myScenario = scene;
+	curTime = 0.0f;
+	isFirstStep = 1;
+	useFileInput = false;
+	Simulation::initBoundaries(myScenario);
+	if (visInfo != NULL) {
+		wAverage = visInfo->waterHeightAtRest();
+		wScale = visInfo->waterVerticalScaling()*(maxDim/20.0f);
+		wOffset = visInfo->waterDistanceFromGround()*(maxDim/20.0f);
+		bAverage = visInfo->bathyVerticalCenter();
+		bScale = visInfo->bathyVerticalScaling()*(maxDim/20.0f);
+		bOffset = visInfo->bathyDistanceFromGround()*(maxDim/20.0f);
+	} else {
+		getScalingApproximation(splash->getWaterHeight(), splash->getBathymetry());
+	}
+}
+
+/**
+    This is the main simulation procedure. Simulates one timestep and computes
+	normals afterwards.
+	
+    @param vbo_resource			cuda resource holding the vertex positions
+	@param vbo_normals			cuda resource holding the normals
+*/
+void Simulation::runCuda(struct cudaGraphicsResource **vbo_resource, struct cudaGraphicsResource **vbo_normals)
+{
+    // map OpenGL buffer object for writing from CUDA
+    float3 *dptr, *dptr2;
+	size_t num_bytes, num_bytes2; 
+    cudaGraphicsMapResources(1, vbo_resource, 0);
+    
+    cudaGraphicsResourceGetMappedPointer((void **)&dptr, &num_bytes,  
+						       *vbo_resource);
+
+	// run simulation and fill our vbo
+	calculateWaterSurface(dptr);
+
+    // unmap buffer object
+    cudaGraphicsUnmapResources(1, vbo_resource, 0);
+	checkCUDAError("Fehler bei Simulations-VBO\n");
+
+	//Calculate normals
+    cudaGraphicsMapResources(1, vbo_resource, 0);
+    cudaGraphicsResourceGetMappedPointer((void **)&dptr, &num_bytes,  
+						       *vbo_resource);
+	cudaGraphicsMapResources(1, vbo_normals, 0);
+    cudaGraphicsResourceGetMappedPointer((void **)&dptr2, &num_bytes2,  
+						       *vbo_normals);
+
+	// calculate normals of our new water surface
+	calculateNormals(dptr, dptr2);
+
+    // unmap buffer objects
+	cudaGraphicsUnmapResources(1, vbo_resource, 0);
+	checkCUDAError("Fehler bei Simulations-VBO\n");
+    cudaGraphicsUnmapResources(1, vbo_normals, 0);
+	checkCUDAError("Fehler bei Normalen-VBO\n");
+}
+
+/**
+    Init new boundaries defined by the scene
+
+    @param scene			instance of SWE_Scenario
+*/
+void Simulation::initBoundaries(SWE_Scenario* scene) {
+std::cout << "Init Scenario\n" << flush;
+	splash->initScenario(*scene);
+
+	splash->setGhostLayer();
+}
+
+
+/**
+    Restarts the simulation. Restores the initial bondaries.
+
+*/
+void Simulation::restart() {
+	curTime = 0.0f;
+	isFirstStep = 1;
+	initBoundaries(myScenario);
+}
+
+
+/**
+    Sets the bathymetry buffer. Buffer contains vertex position and 
+    vertex normal in sequence.
+    @param bath				float array in which computed values will be stored
+*/
+void Simulation::setBathBuffer(float* bath) {
+	// splash->setBathBuffer(output, bAverage, bScale, bOffset);
+	
+        const Float2D& b = splash->getBathymetry();
+        int nx = b.getRows()-2;
+        int ny = b.getCols()-2;
+        
+        // Set vertex coordinates
+	for (int j=0; j<ny+1;j++) {
+		for (int i=0;i<nx+1;i++) {
+				bath[(j*(ny+1) + i)*6] = (float) i;
+				bath[(j*(ny+1) + i)*6 + 1]= scaleFunction(0.25f*(b[i][j]+b[i+1][j]+b[i][j+1]+b[i+1][j+1]), 
+                                                                          bAverage, bScale, bOffset);
+				bath[(j*(ny+1) + i)*6 + 2] = (float) j;
+				bath[(j*(ny+1) + i)*6 + 3] = 0.0f;
+				bath[(j*(ny+1) + i)*6 + 4] = 0.0f;
+				bath[(j*(ny+1) + i)*6 + 5] = 0.0f;
+		}
+	}
+	// Calculate normals
+	for(int j=0; j < ny; j++)
+	{
+		for(int i=0; i < nx; i++)
+		{
+			// Calculate normal vectors for each triangle
+			float normal1[3];
+			float normal2[3];
+
+			calculateNormal(&bath[(j*(ny+1) + i)*6], 
+				&bath[((j+1)*(ny+1) + i + 1)*6],
+				&bath[((j+1)*(ny+1) + i)*6],
+				normal1);
+
+			calculateNormal(&bath[(j*(ny+1) + i)*6],
+				&bath[(j*(ny+1) + i + 1)*6],
+				&bath[((j+1)*(ny+1) + i + 1)*6],
+				normal2);
+			// Copy normals to array
+			for (int k=0; k < 3; k++) {
+				bath[(j*(ny+1) + i)*6 + 3 + k] = (normal1[k]+normal2[k])*0.5f;
+			}
+		}
+	}
+
+	// Fill boundary regions
+	for(int x=0; x < ny; x++) {
+		for (int i=0; i < 3; i++) {
+			bath[(ny*(ny+1) + x)*6 + 3 + i] = bath[((ny-1)*(ny+1) + x)*6 + 3 + i];
+		}
+	}
+	for(int y=0; y < nx; y++) {
+		for (int i=0; i < 3; i++) {
+			bath[(y*(ny+1) + nx)*6 + 3 + i] = bath[(y*(ny+1) + nx - 1)*6 + 3 + i];
+		}
+	}
+	for (int i=0; i < 3; i++) {
+			bath[(ny*(ny+1) + nx)*6 + 3 + i] = bath[((ny-1)*(ny+1) + nx - 1)*6 + 3 + i];
+	}
+}
+
+/**
+    Used for debugging purposes only. Write the current simulation and 
+	visualization buffer to a file.
+*/
+void Simulation::writeDebugOutput(float3* destBuffer) {
+	splash->writeVTKFile3D(generateFileName("debugoutput",0));
+	if (destBuffer != NULL) {
+		debugVisBuffer(destBuffer);
+	}
+}
+/**
+    Save the current simulation state to a vtk file
+*/
+void Simulation::saveToFile() {
+	splash->writeVTKFile3D(generateFileName("simdata", fileNumber));
+	std::cout << "Saved data to: " << generateFileName("simdata", fileNumber) << std::endl;
+	fileNumber++;
+}
+
+/**
+    Advance single step in simulation on graphics card 
+	and copy results to visualization buffer
+
+    @param destBuffer			vertex buffer object holding all water surface data
+*/
+void Simulation::calculateWaterSurface(float3* destBuffer) {
+	if (isFirstStep == 1) {
+		isFirstStep = 0;
+	} else {
+		//splash->simulateConstTimestep();
+		curTime = splash->simulate(curTime, curTime + 0.001f);
+	}
+	// splash->updateVisBuffer(destBuffer, wAverage, wScale, wOffset);
+	updateVisBuffer(destBuffer);
+}
+
+/**
+    Computes a first approximation of the scaling values needed
+	for visualization.
+	Gets called before simulation starts and determines the average,
+	mininimum and maximum values of the bathymetry and water surface data.
+	Uses latter values to estimate the scaling factors.
+	@param h	Fload2D array that holds water height
+	@param b	Fload2D array that holds bathymetry data
+*/	
+void Simulation::getScalingApproximation(const Float2D& h, const Float2D& b) {
+	// Averages of B and (B + H)
+	float avgB, avgH; 
+	// Minimum values
+	float minB, minH;
+	// Maximum values
+	float maxB, maxH;
+        
+        int nx = h.getRows()-2;
+        int ny = h.getCols()-2;
+	int maxDim = (nx > ny) ? nx : ny;
+
+	avgH = 0.0f;
+	avgB = 0.0f;
+	minB = b[1][1];
+	minH = h[1][1];
+	maxB = b[1][1];
+	maxH = h[1][1];
+
+	for(int i=1; i<=nx; i++) {
+		for(int j=1; j<=ny; j++) {
+			// Update averages
+			avgH += (h[i][j] + b[i][j])/(nx*ny);
+			avgB += b[i][j]/(nx*ny);
+			// Update minima
+			if ((h[i][j] + b[i][j]) < minH)
+				minH = (h[i][j] + b[i][j]);
+			if (b[i][j] < minB)
+				minB = b[i][j];
+			// Update maxima
+			if ((h[i][j] + b[i][j]) > maxH)
+				maxH = (h[i][j] + b[i][j]);
+			if (b[i][j] > maxB)
+				maxB = b[i][j];			
+		}
+	}
+	// cout << "max and min bathymetry: " << maxB << " ; " << minB << endl;
+	// bAverage = avgB;
+	bAverage = (maxB+minB)/2;
+	cout << "Average bathymetry: " << bAverage << endl;
+	if (abs(maxB - minB) > 0.0001f) {
+		bScale = (maxDim/20.0f)/(maxB - minB);
+	} else {
+		bScale = (maxDim/20.0f);
+	}
+	cout << "Scaling of bathymetry: " << bScale << endl;
+	bOffset = bScale*(avgB - minB) + maxDim/15.0f;
+	// bOffset = minB+100;
+	cout << "bathymetry offset: " << bOffset << endl;
+
+	wAverage = avgH;
+	cout << "Average water height: " << avgH << endl;
+	if ((maxH - minH) < 0.0001f) {
+		wScale = 1.0f/(maxH- minH);
+	} else {
+		wScale = 1.0f;
+	}
+	wScale = (maxDim/40.0)*wScale;
+	cout << "Scaling of water level: " << wScale << endl;
+	wOffset = bOffset + (maxDim/5.0f);
+	cout << "water level offset: " << wOffset << endl;
+
+}
+
+/**
+    Compute normal of a triangle
+	@param fVert1-fVert3	vertices of the triangle
+	@param fNormal			resulting normal
+*/	
+void Simulation::calculateNormal(float fVert1[], float fVert2[],
+                                 float fVert3[], float fNormal[]) {
+   float Qx, Qy, Qz, Px, Py, Pz;
+
+   Qx = fVert2[0]-fVert1[0];
+   Qy = fVert2[1]-fVert1[1];
+   Qz = fVert2[2]-fVert1[2];
+   Px = fVert3[0]-fVert1[0];
+   Py = fVert3[1]-fVert1[1];
+   Pz = fVert3[2]-fVert1[2];
+
+   fNormal[0] = Py*Qz - Pz*Qy;
+   fNormal[1] = Pz*Qx - Px*Qz;
+   fNormal[2] = Px*Qy - Py*Qx;
+}
+
+//==================================================================
+// member functions for visualisation and output
+//==================================================================
+/**
+    Scale function used for visualization 
+*/	
+__device__ __host__ 
+float scaleFunction(float val, float average, float scale, float offset) {
+	return (val - average)*scale + offset;
+}
+
+// /**
+//     Return coordinates of visBuffer[x][y]
+// 	@param x,y		coordinates
+// 	@param ny		length in y-direction
+// */	
+// __device__
+// int getVisBufCoord(int x, int y, int ny) {
+// 	return y*(ny+1) + x;
+// }
+
+/**
+	Transform cell-centered values into node-centered values
+	for visualization purposes
+	
+	@param visBuffer		pointer to visualization buffer
+	@param hd				h-data on device
+	@param bd				b-data on device
+	@param nx, ny			x- and y-dimensions
+	@param dx, dy			cellsize in x- and y-direction
+	@param average			center (average) of h-data (water surface) 
+							for visualization
+	@param scale			scaling factor which is applied for visualization
+	@param offset			offset factor for visualization
+ */
+__global__
+void kernelCalcVisBuffer(float3* visBuffer, const float* hd, const float* bd, 
+                         int nx, int ny, float average, float scale, float offset)
+{
+   int i = TILE_SIZE*blockIdx.x + threadIdx.x;
+   int j = TILE_SIZE*blockIdx.y + threadIdx.y;
+   if ((i <= nx) && (j <= ny)) {
+	   int index = i*(nx+2) + j;
+	   int index2 = (i+1)*(nx+2) + j;	
+	   visBuffer[j*(ny+1) + i] = make_float3(
+		   i, 
+		   scaleFunction(0.25*(
+				hd[index]+hd[index + 1]+ hd[index2] + hd[index2 + 1] +
+				bd[index]+bd[index + 1]+ bd[index2] + bd[index2 + 1])
+				, average, scale, offset),
+		   j);
+   } 
+   //Corresponding C-Code:
+   //	for (int j=0; j<ny+1;j++)
+   //		for (int i=0;i<nx+1;i++)
+   //			Vtk_file << i*dx<<" "<<j*dy<<" "
+   //			         << 0.25*(h[i][j]+h[i+1][j]+h[i][j+1]+h[i+1][j+1]
+   //				         +b[i][j]+b[i+1][j]+b[i][j+1]+b[i+1][j+1]) 
+}
+
+/**
+	copy h-values into visualization buffer
+	@param _visBuffer		visualization buffer
+
+*/
+void Simulation::updateVisBuffer(float3* _visBuffer) {
+	
+        const float* hd = splash->getCUDA_waterHeight();
+        const float* bd = splash->getCUDA_bathymetry();
+        int nx = splash->getNx();
+        int ny = splash->getNy();
+//         // Fill ghost layer corner cells
+// 	kernelHdBufferEdges<<<1,1>>>(hd, nx, ny);
+	// Interpolate cell centered h-values
+	dim3 dimBlock(TILE_SIZE,TILE_SIZE);
+	dim3 dimGrid((nx+TILE_SIZE)/TILE_SIZE,(ny+TILE_SIZE)/TILE_SIZE);
+	kernelCalcVisBuffer<<<dimGrid,dimBlock>>>(_visBuffer, hd, bd, nx, ny, 
+                                                  wAverage, wScale, wOffset);
+}
+/**
+	Function used for debugging. Outputs the current visBuffer 
+	@param visBuffer		current visualization buffer
+	
+*/
+void Simulation::debugVisBuffer(float3* _visBuffer) {
+        int nx = splash->getNx();
+        int ny = splash->getNy();
+	float* dbgBuf = new float[(nx+1)*(ny+1)*3];
+	int size = (nx+1)*(ny+1)*sizeof(float3);
+	cudaMemcpy(dbgBuf,_visBuffer, size, cudaMemcpyDeviceToHost);
+    checkCUDAError("memory of visualization buffer not transferred");
+	for (int i = 0; i < 3*(nx+1)*(ny+1); i = i + 3) {
+		std::cout << "(" << dbgBuf[i] << ", " << dbgBuf[i+1] << ", " << dbgBuf[i+2] << ") \n";
+	}
+	delete[] dbgBuf;
+}
+
+/**
+    Compute normal of a triangle
+	@param fVert1-fVert3	vertices of the triangle
+	@return	resulting normal
+*/	
+inline __device__
+float3 calculateNormal(float3 fVert1, float3 fVert2,
+                             float3 fVert3)
+   {
+   float Qx, Qy, Qz, Px, Py, Pz;
+
+   Qx = fVert2.x-fVert1.x;
+   Qy = fVert2.y-fVert1.y;
+   Qz = fVert2.z-fVert1.z;
+   Px = fVert3.x-fVert1.x;
+   Py = fVert3.y-fVert1.y;
+   Pz = fVert3.z-fVert1.z;
+
+   return make_float3(Py*Qz - Pz*Qy, Pz*Qx - Px*Qz, Px*Qy - Py*Qx);
+}
+
+/**
+    Compute new normals resulting from updated water surface
+	@param visBuffer			vertex buffer holding water surface vertices
+	@param normalBuffer			buffer which will contain new normals
+	@param nx, ny				
+*/	
+__global__
+void kernelCalcNormals(float3* visBuffer, float3* normalBuffer, int nx, int ny)
+{
+   int i = TILE_SIZE*blockIdx.x + threadIdx.x;
+   int j = TILE_SIZE*blockIdx.y + threadIdx.y;
+   if ((i < nx + 1) && (j < ny + 1)) {
+	   int _i = i;
+	   int _j = j;
+	   // Handle boundaries
+	   if (i >= nx) _i = i - 1;
+	   if (j >= ny) _j = j - 1;
+	   normalBuffer[j*(ny+1) + i] = calculateNormal(visBuffer[(_j*(ny+1) + _i)], 
+				visBuffer[(_j+1)*(ny+1) + _i + 1],
+				visBuffer[(_j+1)*(ny+1) + _i]);
+   } 
+}
+
+/**
+    Compute normals of the new water surface on graphics card
+
+    @param vertexBuffer			cuda resource holding the vertex positions
+	@param destBuffer			cuda resource holding the new normals
+*/
+void Simulation::calculateNormals(float3* vertexBuffer, float3* destBuffer) {
+// 	splash->calculateNormals(vertexBuffer, destBuffer);
+	int nx = splash->getNx();
+	int ny = splash->getNy();
+	dim3 dimBlock(TILE_SIZE,TILE_SIZE);
+	dim3 dimGrid((nx+TILE_SIZE)/TILE_SIZE,(ny+TILE_SIZE)/TILE_SIZE);
+	kernelCalcNormals<<<dimGrid,dimBlock>>>(vertexBuffer, destBuffer, nx, ny);
+}
+
+/**
+	Helper function to check for errors in CUDA calls
+	param msg			custom error message
+	source: NVIDIA
+
+ */

src/opengl/simulation.h

+#ifndef SIMULATION_H
+#define SIMULATION_H
+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Michael Bader, Kaveh Rahnema, Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include <math.h>
+#include <cuda_runtime.h>
+
+#include "../SWE_BlockCUDA.hh"
+#include "../scenarios/SWE_Scenario.h"
+#include "../scenarios/SWE_VisInfo.h"
+
+void checkCUDAError(const char *msg);
+
+class Simulation {
+
+  public:
+    // Constructor + Destructor 
+    Simulation (int nx, int ny, float dx, float dy, 
+                SWE_Scenario* scene, SWE_VisInfo* visInfo, SWE_BlockCUDA* init_splash);
+    ~Simulation();
+
+    // Restart simulation
+    void restart();
+    // Load new scenario after initialization
+    void loadNewScenario(SWE_Scenario* scene, SWE_VisInfo* visInfo);
+    // Save simulation state to file
+    void saveToFile();
+    // Return the bathymetry data
+    void setBathBuffer(float* output);
+    // Simulate single timestep on graphics card
+    void runCuda(struct cudaGraphicsResource **vbo_resource, struct cudaGraphicsResource **vbo_normals);
+
+    // Debugging
+    void writeDebugOutput(float3* destBuffer = NULL);
+
+  protected:
+
+    // Instance of SWE_BlockCUDA 
+    SWE_BlockCUDA* splash;
+    // Current scenario
+    SWE_Scenario* myScenario;
+
+    // Current simulation time
+    float curTime;
+    // Is this our first simulation step?
+    int isFirstStep;
+    // Maximum of cell sizes
+    float maxCellSize;
+    // Scaling factors used by visualization
+    float bOffset, bScale, wOffset, wScale, bAverage, wAverage;
+    // Initalize boundaries defined by the scene
+    void initBoundaries(SWE_Scenario* scene);
+    // Initalize boundaries defined by an input file 
+
+    int fileNumber;
+    // Use file input as boundary data
+    bool useFileInput;
+    // Holds maximum dimension
+    int maxDim;
+    // Compute new water surface
+    void calculateWaterSurface(float3* destBuffer);
+    // Compute normals of the water surface for shading
+    void calculateNormals(float3* vertexBuffer, float3* destBuffer);
+
+    void getScalingApproximation(const Float2D& h, const Float2D& b);
+
+    void updateVisBuffer(float3* _visBuffer);
+    void debugVisBuffer(float3* _visBuffer);
+
+    /**
+        Scale function used for visualization 
+    */	
+    __device__ __host__ 
+    static float scaleFunction(float val, float average, float scale, float offset) {
+       return (val - average)*scale + offset;
+    }
+
+    static void calculateNormal(float fVert1[], float fVert2[],
+				float fVert3[], float fNormal[]);
+
+};
+
+#endif

src/opengl/visualization.cpp

+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+
+#include "visualization.h"
+
+/**
+	Constructor. All dimensions are node-based, this means a grid consisting of 
+	2x2 cells would have 3x3 nodes.
+
+	@param: window_title	title of the window created
+	@param: _grid_x_size	number of nodes of the grid (in x-direction)
+	@param: _grid_y_size	number of nodes of the grid (in y-direction)
+
+*/
+Visualization::Visualization(int windowWidth, int windowHeight, const char* window_title, 
+							 int _grid_xsize, int _grid_ysize) {
+	// Initialize member variables
+	grid_xsize = _grid_xsize;
+	grid_ysize = _grid_ysize;
+	renderMode = SHADED;
+
+	// Initialize rendering
+	initSDL(windowWidth, windowHeight);
+	initGLWindow(windowWidth, windowHeight);
+	initGLDefaults();
+	initCUDA();
+	
+	// Load camera and shaders
+	camera = new Camera(_grid_xsize*1.5f, window_title);
+	shaders = new Shader("vertex.glsl", "fragment.glsl");
+	if (shaders->shadersSupported()) {
+		printf("Shaders supported!\n");
+	} else {
+		printf("Shaders are NOT supported! Normal rendering mode\n");
+	}
+	if (shaders->shadersLoaded()) {
+		printf("Shaders successfully loaded!\n");
+		renderMode = WATERSHADER;
+	} else {
+		printf("Shaders error while loading shaders\n");
+	}
+}
+
+/** 
+	Destructor (see note below)
+*/
+Visualization::~Visualization() {
+	delete camera;
+	delete shaders;
+	SDL_Quit();
+}
+
+/**
+	Frees all memory we used for geometry data 
+	Needs to be called before destructor gets called 
+	in order to work correctly
+
+*/
+void Visualization::cleanUp() {
+	deleteVBO(&vboWaterSurface, cuda_vbo_watersurface);
+	deleteVBO(&vboNormals, cuda_vbo_normals);
+	deleteVBO(&verticesIndex);
+	deleteVBO(&vboBathymetry);
+}
+
+
+
+/**
+    Main rendering function. Draws the scene and updates screen
+
+*/	
+void Visualization::renderDisplay() {
+	// Set camera
+	camera->setCamera();
+	
+	// Draw Scene
+	DrawBottom();
+	DrawBathymetry(vboBathymetry, verticesIndex);
+	
+	// Shaded pass
+	DrawWaterSurface(vboWaterSurface, vboNormals, verticesIndex);
+
+	// Update framebuffer
+	camera->displayImage();  
+
+	// Check for errors
+	if (glGetError() != GL_NO_ERROR) {
+		printf("OpenGL error occured..\n");
+	}
+}
+
+/**
+    Draws the water surface geometry (triangles)
+
+    @param vboID			id of the vertex buffer object storing 
+							the vertex positions
+	@param vboID			id of the array buffer object storing
+							the corresponding normals for each vertex
+	@param verticesIndex	id of the array buffer object storing the 
+							vertex indices in rendering sequence 
+*/
+void Visualization::DrawWaterSurface(GLuint vboID, GLuint vboNormals, GLuint verticesIndex) {
+	if (renderMode == WATERSHADER) {
+		// Depth pass first
+		glColorMask( GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE );
+		renderMode = SHADED;
+		DrawWaterSurface(vboID, vboNormals, verticesIndex);
+		renderMode = WATERSHADER;
+		glColorMask(GL_TRUE,GL_TRUE,GL_TRUE,GL_TRUE);
+
+		// Now render all visible geometry
+		shaders->enableShader();
+		glEnable(GL_BLEND);
+		glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+	}
+	else if (renderMode == SHADED) {
+		// Enable lighting
+		glEnable(GL_COLOR_MATERIAL);
+		glEnable(GL_LIGHTING);
+	} else {
+		glDisable(GL_LIGHTING);
+	}
+	
+	// Enable array rendering
+	glEnableClientState(GL_NORMAL_ARRAY);
+	glEnableClientState(GL_VERTEX_ARRAY);
+
+	// Set rendering to VBO mode
+	glBindBuffer(GL_ARRAY_BUFFER,vboNormals);
+	glNormalPointer(GL_FLOAT, 0, 0);
+	glBindBuffer(GL_ARRAY_BUFFER, vboID);
+    glVertexPointer(3, GL_FLOAT, 0, 0);
+	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, verticesIndex);
+	
+	// Enable VBO access and render triangles
+	glPushMatrix();
+		glTranslatef(-(grid_xsize-1)/2.0f,0.0f,-(grid_ysize-1)/2.0f);
+		glColor3f(0.3f*1.2f, 0.45f*1.2f, 0.9f*1.2f);
+		if (renderMode == WIREFRAME)
+			glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
+		glDrawElements(GL_TRIANGLES, 6*(grid_xsize - 1)*(grid_ysize - 1), GL_UNSIGNED_INT, NULL);
+		glPolygonMode( GL_FRONT_AND_BACK, GL_FILL );
+	glPopMatrix();
+	
+	// Disable array rendering
+	glDisableClientState(GL_VERTEX_ARRAY);
+	glDisableClientState(GL_NORMAL_ARRAY);
+	glDisable(GL_LIGHTING);
+	glDisable(GL_BLEND);
+	shaders->disableShader();
+
+}
+
+/**
+    Draws the bathymetry geometry
+
+    @param vboID			id of the vertex buffer object storing 
+							the vertex positions
+	@param verticesIndex	id of the array buffer object storing the 
+							vertex indices in rendering sequence 
+*/
+
+void Visualization::DrawBathymetry(GLuint vboID, GLuint verticesIndex) {
+	GLsizei stride = 6*sizeof(float);
+	
+	// Enable lighting
+	glEnable(GL_COLOR_MATERIAL);
+	glEnable(GL_LIGHTING);
+	
+	// Enable array rendering
+	glEnableClientState(GL_NORMAL_ARRAY);
+	glEnableClientState(GL_VERTEX_ARRAY);
+	
+	// Bind buffers
+	glBindBuffer(GL_ARRAY_BUFFER, vboID);
+	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, verticesIndex);
+	
+	// Set VBO pointers
+	const GLvoid* normalOffset = (GLvoid*) 12;
+	glNormalPointer(GL_FLOAT, stride, normalOffset);
+	glVertexPointer(3, GL_FLOAT, stride, 0);
+
+	// Render triangles
+	glPushMatrix();
+		glTranslatef(-(grid_xsize-1)/2.0f,0.0f,-(grid_ysize-1)/2.0f);
+		glColor3f(0.4f*1.1f, 0.36f*1.1f, 0.3f*1.1f);
+		glDrawElements(GL_TRIANGLES, 6*(grid_xsize - 1)*(grid_ysize - 1), GL_UNSIGNED_INT, NULL);
+	glPopMatrix();
+
+	// Disable array rendering
+	glDisableClientState(GL_VERTEX_ARRAY);
+	glDisableClientState(GL_NORMAL_ARRAY);
+	glDisable(GL_LIGHTING);
+}
+
+/**
+    Draws the bottom plane in our scene
+
+*/
+void Visualization::DrawBottom()
+{
+	glBegin(GL_QUADS);	
+		// Draw A Quad
+		glNormal3f(0.0f,1.0f,0.0f);
+		glColor3f(0.3f,0.3f,0.3f);
+		glVertex3f( (grid_xsize - 1)/2.0f,0.0f, (grid_ysize - 1)/2.0f);					// Top Right Of The Quad (Bottom)
+		glVertex3f(-(grid_xsize - 1)/2.0f,0.0f, (grid_ysize - 1)/2.0f);					// Top Left Of The Quad (Bottom)
+		glVertex3f(-(grid_xsize - 1)/2.0f,0.0f,-(grid_ysize - 1)/2.0f);					// Bottom Left Of The Quad (Bottom)
+		glVertex3f( (grid_xsize - 1)/2.0f,0.0f,-(grid_ysize - 1)/2.0f);					// Bottom Right Of The Quad (Bottom)
+	glEnd();				
+}
+
+
+
+/**
+	Allocates memory for vertices and other geometry data.
+	
+	@param sim				instance of the simulation class
+
+*/
+void Visualization::init(Simulation* sim) {
+	// Create VBOs
+	createBathymetryVBO(&vboBathymetry, grid_xsize * grid_ysize * 6 * sizeof(float), sim);
+	createIndicesVBO(&verticesIndex, grid_xsize, grid_ysize);
+	createVertexVBO(&vboWaterSurface, grid_xsize * grid_ysize * 3 * sizeof(float), 
+					&cuda_vbo_watersurface, cudaGraphicsMapFlagsNone);	
+	createVertexVBO(&vboNormals ,grid_xsize * grid_ysize * 3 * sizeof(float), 
+					&cuda_vbo_normals, cudaGraphicsMapFlagsWriteDiscard);
+}
+
+/**
+	Returns a pointer to the cuda memory object holding the
+	vertex normals
+*/
+cudaGraphicsResource** Visualization::getCudaNormalsPtr() {
+	return &cuda_vbo_normals;
+}
+/**
+	Returns a pointer to the cuda memory object holding the
+	vertex positions
+*/
+cudaGraphicsResource** Visualization::getCudaWaterSurfacePtr() {
+	return &cuda_vbo_watersurface;
+}
+
+/**
+    Returns, whether a special extension is supported by the current 
+	graphics card
+
+    @param	szTargetExtention	string describing the extension to look for	
+
+*/	
+bool Visualization::IsExtensionSupported(const char* szTargetExtension )
+{
+	const unsigned char *pszExtensions = NULL;
+	const unsigned char *pszStart;
+	unsigned char *pszWhere, *pszTerminator;
+
+	// Extension names should not have spaces
+	pszWhere = (unsigned char *) strchr( szTargetExtension, ' ' );
+	if( pszWhere || *szTargetExtension == '\0' )
+		return false;
+
+	// Get Extensions String
+	pszExtensions = glGetString( GL_EXTENSIONS );
+
+	// Search The Extensions String For An Exact Copy
+	pszStart = pszExtensions;
+	for(;;)
+	{
+		pszWhere = (unsigned char *) strstr( (const char *) pszStart, szTargetExtension );
+		if( !pszWhere )
+			break;
+		pszTerminator = pszWhere + strlen( szTargetExtension );
+		if( pszWhere == pszStart || *( pszWhere - 1 ) == ' ' )
+			if( *pszTerminator == ' ' || *pszTerminator == '\0' )
+				return true;
+		pszStart = pszTerminator;
+	}
+	return false;
+}
+
+/**
+    Initialize the simple directmedia layer (SDL) which is a wrapper library 
+	for OpenGL
+
+    @param	windowWidth		width in pixels of the window to create
+	@param  windowHeight	height in pixels
+
+*/	
+void Visualization::initSDL(int windowWidth, int windowHeight) {
+	// Initialize SDL system
+	if ( SDL_Init(SDL_INIT_VIDEO) < 0 ) {
+		fprintf(stderr, "Unable to initialize SDL: %s\n", SDL_GetError());
+		exit(1);
+	}
+	
+	// Enable double buffering and center window
+	SDL_GL_SetAttribute( SDL_GL_DOUBLEBUFFER, 1 );
+	SDL_putenv(const_cast<char *>("SDL_VIDEO_CENTERED=center"));
+	SDL_WM_SetCaption("Loading...", NULL);
+
+	// Initialize window with OpenGL support
+	if ( SDL_SetVideoMode(windowWidth, windowHeight, 0, SDL_OPENGL | SDL_RESIZABLE | SDL_ANYFORMAT) == NULL ) {
+		fprintf(stderr, "Unable to create OpenGL screen: %s\n", SDL_GetError());
+		SDL_Quit();
+		exit(2);
+	}
+
+	// Check OpenGL extension(s)
+	if (!IsExtensionSupported( "GL_ARB_vertex_buffer_object" )) {
+		printf("Vertex Buffer Objects Extension not supported! Exit..\n");
+		SDL_Quit();
+		exit(1);
+	}
+	// Load Vertex Buffer Extension
+	glGenBuffers = (PFNGLGENBUFFERSARBPROC) SDL_GL_GetProcAddress("glGenBuffersARB");
+	glBindBuffer = (PFNGLBINDBUFFERARBPROC) SDL_GL_GetProcAddress("glBindBufferARB");
+	glBufferData = (PFNGLBUFFERDATAARBPROC) SDL_GL_GetProcAddress("glBufferDataARB");
+	glDeleteBuffers = (PFNGLDELETEBUFFERSARBPROC) SDL_GL_GetProcAddress("glDeleteBuffersARB");
+	
+}
+
+/**
+    Initializes OpenGL projection and viewport settings
+
+    @param	width		width in pixels of the window to create
+	@param  height		height in pixels
+
+*/	
+void Visualization::initGLWindow(int width, int height) {	
+	GLfloat ratio;
+
+	// Protect against a divide by zero
+	if ( height == 0 )
+		height = 1;
+
+	ratio = ( GLfloat )width / ( GLfloat )height;
+	// Setup our viewport.
+	glViewport( 0, 0, ( GLint )width, ( GLint )height );
+
+    // change to the projection matrix and set our viewing volume.
+    glMatrixMode( GL_PROJECTION );
+    glLoadIdentity( );
+	float Zoom = 0.00005f;
+	glFrustum(-Zoom * width, Zoom * width, -Zoom * height, Zoom * height, 0.1f, 10.0f*grid_ysize);
+
+    // Switch back to the modelview
+    glMatrixMode( GL_MODELVIEW );
+}
+
+/**
+    Gets called when window gets resized
+
+    @param	newWidth		new window width in pixels 
+	@param  newHeight		height in pixels
+
+*/
+int Visualization::resizeWindow(int newWidth, int newHeight) {
+	if ( SDL_SetVideoMode( newWidth, newHeight, 0, 
+		SDL_OPENGL | SDL_RESIZABLE | SDL_ANYFORMAT ) == NULL ) 
+	{
+		fprintf( stderr, "Could not get a surface after resize: %s\n", SDL_GetError( ) );
+		return 1;
+	} else {
+		initGLWindow ( newWidth, newHeight );
+		return 0;
+	}
+	initGLWindow(newWidth, newHeight);
+}
+
+/**
+    Initializes various OpenGL settings
+*/
+void Visualization::initGLDefaults() {
+	// Set background color to white
+	glClearColor(1.0f, 1.0f, 1.0f, 0.0f);
+
+    // Smooth (interpolated) shading
+    glShadeModel( GL_SMOOTH );
+
+    // Reset the depth buffer
+    glClearDepth( 1.0f );
+
+    // Enable depth testing
+    glEnable( GL_DEPTH_TEST );
+
+    // Comparisons in depth buffer via <=
+    glDepthFunc( GL_LEQUAL );
+
+    // Some perspective corrections
+    glHint( GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST );
+
+	// Ensure that all of the following transformations will only
+	// apply to our modelview matrix (objects + camera)
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+
+	// Lighting
+	GLfloat LightAmbient[]=		{ 0.0f, 0.0f, 0.0f, 1.0f };
+	GLfloat LightDiffuse[]=		{ 1.0f, 1.0f, 1.0f, 1.0f };
+	GLfloat LightPosition[]=	{ 1.0f, 1.0f, 0.0f, 0.0f };
+	glLightfv(GL_LIGHT0, GL_POSITION,LightPosition);	// Position The Light
+	glLightfv(GL_LIGHT0, GL_AMBIENT, LightAmbient);		// Setup The Ambient Light
+	glLightfv(GL_LIGHT0, GL_DIFFUSE, LightDiffuse);		// Setup The Diffuse Light
+	glEnable(GL_LIGHT0);	
+	glEnable(GL_NORMALIZE);
+	glDisable(GL_LIGHTING);
+}
+
+/**
+    Calculate 1-D array offset for 2-D coordinates
+
+    @param	x			x-position
+	@param	y			y-position
+	@param	width		size of a line
+
+*/	
+int Visualization::coord(int x, int y, int width) {
+	if (width == -1) width = grid_xsize;
+	return (y*width + x);
+}
+
+/**
+    Creates a vertex buffer object in OpenGL 
+
+	@param size				size in bytes to allocate
+	@return	vboID			pointer to an integer depicting the id of the 
+							new vertex buffer object 
+
+*/	
+void Visualization::createVertexVBO(GLuint* vboID, int size)
+{
+	// Create 1 buffer object
+	glGenBuffers(1, vboID);
+	glBindBuffer(GL_ARRAY_BUFFER, *vboID);
+	
+	// Initialize buffer object
+	glBufferData(GL_ARRAY_BUFFER, size, NULL, GL_DYNAMIC_DRAW);
+	
+	// Switch to default buffer
+	glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+/**
+    Creates a vertex buffer object in OpenGL and loads it with the bathymetry
+	data from simulation
+
+
+	@param size				size in bytes to allocate
+	@param sim				pointer to an instance of the simulation class
+	@return	vboID			pointer to an integer depicting the id of the 
+							new vertex buffer object 
+
+*/	
+void Visualization::createBathymetryVBO(GLuint* vboID, int size, Simulation* sim) {
+	GLfloat* vBathy = new GLfloat[size/sizeof(GLfloat)];
+
+	// Create buffer object for vertex indices
+	glGenBuffers(1, vboID);
+	glBindBuffer(GL_ARRAY_BUFFER, *vboID);
+	
+	// Initialize buffer object
+	sim->setBathBuffer(vBathy);
+	glBufferData(GL_ARRAY_BUFFER, size, vBathy, GL_STATIC_DRAW);
+	
+	// Switch to default buffer
+	glBindBuffer(GL_ARRAY_BUFFER, 0);
+	delete[] vBathy;
+}
+
+/**
+    Updates vertex buffer object with new bathymetry
+	data from simulation
+
+	@param sim				pointer to an instance of the simulation class
+
+*/	
+void Visualization::updateBathymetryVBO(Simulation* sim) {
+	int size = grid_xsize * grid_ysize * 6 * sizeof(float);
+	GLfloat* vBathy = new GLfloat[size/sizeof(GLfloat)];
+	
+	// Select buffer
+	glBindBuffer(GL_ARRAY_BUFFER, vboBathymetry);
+	
+	// Update buffer object
+	sim->setBathBuffer(vBathy);
+	glBufferData(GL_ARRAY_BUFFER, size, vBathy, GL_STATIC_DRAW);
+	
+	// Switch to default buffer
+	glBindBuffer(GL_ARRAY_BUFFER, 0);
+	delete[] vBathy;
+}
+/**
+    Creates a vertex buffer object in OpenGL and an associated CUDA resource
+
+	@param size				size in bytes to allocate 
+	@param	vbo_res_flags   cuda flags for memory access
+	@return	vboID			pointer to an integer depicting the id of the 
+							new vertex buffer object 
+	@return	vbo_res			cuda structure created by this function
+	
+*/	
+void Visualization::createVertexVBO(GLuint* vboID, int size, struct cudaGraphicsResource **vbo_res, 
+	       unsigned int vbo_res_flags)
+{
+	createVertexVBO(vboID, size);
+	cudaGraphicsGLRegisterBuffer(vbo_res, *vboID, vbo_res_flags);
+	checkCUDAError("Couldn't register GL buffer");
+}
+
+/**
+    Create an array buffer object which holds a list of vertex indices.
+	Groups of 3 consecutive vertices corresponding to the indices will be  
+	rendered to a triangle.
+	This array buffer object therefore describes how single grid points (nodes)
+	get transformed into a triangle mesh (tesselation).
+
+    @param	vbo				pointer to an integer depicting the id of a 
+							vertex buffer object 
+	@param	xsize			number of grid nodes (in x-direction)
+	@params ysize			number of grid nodes (in y-direction)
+    @return void			
+*/	
+void Visualization::createIndicesVBO(GLuint* vboID, int xsize, int ysize)
+{
+	// Create an array describing the vertex indices to be drawn
+	xsize = xsize;
+	ysize = ysize;
+
+	int noVertices = (xsize-1)*(ysize-1)*6;
+	if ((xsize < 1) || (ysize < 1)) {
+		noVertices = 0;
+	}
+	GLuint* vIndices = new GLuint[noVertices];
+
+	for(int y=0; y < ysize - 1; y++)
+	{
+		for(int x=0; x < xsize - 1; x++)
+		{
+			// Create tessellation of the grid
+			vIndices[coord(x,y, xsize - 1)*6] = coord(x,y);
+			vIndices[coord(x,y, xsize - 1)*6 + 1] = coord(x+1,y+1);
+			vIndices[coord(x,y, xsize - 1)*6 + 2] = coord(x,y+1);
+			vIndices[coord(x,y, xsize - 1)*6 + 3] = coord(x,y);
+			vIndices[coord(x,y, xsize - 1)*6 + 4] = coord(x+1,y);
+			vIndices[coord(x,y, xsize - 1)*6 + 5] = coord(x+1,y+1);
+			
+		}
+	}
+
+	// Create buffer object for vertex indices
+	glGenBuffers(1, vboID);
+	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, *vboID);
+	
+	// Initialize buffer object
+	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint)*noVertices, vIndices, GL_STATIC_DRAW);
+	
+	// Switch to default buffer
+	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+	delete[] vIndices;
+}
+
+/**
+    Frees memory used by a vertex buffer object
+
+    @param	vbo				pointer to an integer depicting the id of a 
+							vertex buffer object 
+*/	
+void Visualization::deleteVBO(GLuint* vbo)
+{
+    if (vbo) {
+		glBindBuffer(1, *vbo);
+		glDeleteBuffers(1, vbo);
+		*vbo = 0;
+    }
+}
+/**
+    Frees memory used by a vertex buffer object and a CUDA resource
+
+    @param	vbo				pointer to an integer depicting the id of a 
+							vertex buffer object 
+	@param	vbo_res			pointer to a CUDA resource structure		
+*/	
+void Visualization::deleteVBO(GLuint* vbo, struct cudaGraphicsResource *vbo_res)
+{
+    if (vbo) {
+		cudaGraphicsUnregisterResource(vbo_res);
+		deleteVBO(vbo);
+    }
+}
+
+/**
+    Initialize CUDA device
+*/	
+void Visualization::initCUDA() {
+	int driver, runtime;
+	cudaDeviceProp  prop;
+    int dev;
+
+	// Display CUDA Versions
+	cudaDriverGetVersion(&driver);
+	cudaRuntimeGetVersion(&runtime);
+	printf("CUDA Driver Version: %d, Runtime Version: %d\n", driver, runtime);
+	
+	// Select GPU for CUDA and OpenGL
+    memset( &prop, 0, sizeof( cudaDeviceProp ) );
+    prop.major = 1;
+    prop.minor = 0;
+    cudaChooseDevice( &dev, &prop );
+    cudaGLSetGLDevice( dev ) ;
+}
+
+
+/**
+    Sets current rendering mode
+
+    @param mode			rendering mode			
+*/	
+void Visualization::setRenderingMode(RenderMode mode) {
+	renderMode = mode;
+}
+
+/**
+    Switches between 3 different rendering modes:
+	- Shaded: Use OpenGL shading
+	- Wireframe: Only render edges of each triangle
+	- Watershader: Use custom GLSL shader for water surface
+*/		
+void Visualization::toggleRenderingMode() {
+	switch( renderMode ) {
+		case SHADED:
+			renderMode = WIREFRAME;
+			break;
+		case WIREFRAME:
+			// Skip watershader if shaders not loaded
+			if (shaders->shadersLoaded()) {
+				renderMode = WATERSHADER;
+			} else {
+				renderMode = SHADED;
+			}
+			break;
+		case WATERSHADER:
+			renderMode = SHADED;
+			break;
+	}
+}
+

src/opengl/visualization.h

+#ifndef VISUALIZATION_H
+#define VISUALIZATION_H
+// =====================================================================
+// This file is part of SWE_CUDA (see file SWE_Block.cu for details).
+// 
+// Copyright (C) 2010,2011 Tobias Schnabel
+// 
+// SWE_CUDA is free software: you can redristribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// SWE_CUDA is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with SWE_CUDA.  If not, see <http://www.gnu.org/licenses/>.
+// =====================================================================
+#include <SDL.h>
+#include <SDL_opengl.h>
+#include <cuda_runtime.h>
+#include <cuda_gl_interop.h>
+#include "camera.h"
+#include "simulation.h"
+#include "shader.h"
+
+void checkCUDAError(const char *msg);
+typedef enum RenderMode {
+   SHADED, WIREFRAME, WATERSHADER
+} RenderMode;
+
+class Visualization {
+public:
+	// Constructor and Destructor
+	Visualization(int windowWidth, int windowHeight, const char* window_title, int _grid_xsize, int _grid_ysize);
+	~Visualization();
+	void init(Simulation* sim);
+	void cleanUp();
+	Camera* camera;
+
+	// Access to CUDA VBO pointers
+	cudaGraphicsResource** getCudaNormalsPtr();
+	cudaGraphicsResource** getCudaWaterSurfacePtr();
+
+	// Main rendering function
+	void renderDisplay();
+	
+	// Helper functions
+	void setRenderingMode(RenderMode mode);
+	void updateBathymetryVBO(Simulation* sim);
+	void toggleRenderingMode();
+	int resizeWindow(int newWidth, int newHeight);
+private:	
+	// Init helper functions
+	void initSDL(int windowWidth, int windowHeight);
+	void initGLWindow(int width, int height);
+	void initGLDefaults();
+	void initCUDA();
+	bool IsExtensionSupported(const char* szTargetExtension );
+
+	// Drawing functions
+	void DrawWaterSurface(GLuint vboID, GLuint vboNormals, GLuint verticesIndex);
+	void DrawBathymetry(GLuint vboID, GLuint verticesIndex);
+	void DrawBottom();
+	int grid_xsize;
+	int grid_ysize;
+
+	// Vertex Buffer objects
+	GLuint vboBathymetry;
+	GLuint verticesIndex;
+	GLuint vboWaterSurface;
+	GLuint vboNormals;
+
+	struct cudaGraphicsResource* cuda_vbo_watersurface;
+	struct cudaGraphicsResource* cuda_vbo_normals;
+
+	// VBO management functions
+	void createIndicesVBO(GLuint* vboID, int xsize, int ysize);
+	void createVertexVBO(GLuint* vboID, int size);
+	void createBathymetryVBO(GLuint* vboID, int size, Simulation* sim);
+	void createVertexVBO(GLuint* vboID, int size, struct cudaGraphicsResource **vbo_res, 
+		unsigned int vbo_res_flags);
+	void deleteVBO(GLuint* vbo);
+	void deleteVBO(GLuint* vbo, struct cudaGraphicsResource *vbo_res);
+
+	// Rendering mode
+	RenderMode renderMode;
+
+	// Shader helper class
+	Shader* shaders;
+
+	// Helper function
+	int coord(int x, int y, int width = -1);
+
+	// VBO Extension Function Pointers
+	PFNGLGENBUFFERSARBPROC glGenBuffers;					// VBO Name Generation Procedure
+	PFNGLBINDBUFFERARBPROC glBindBuffer;					// VBO Bind Procedure
+	PFNGLBUFFERDATAARBPROC glBufferData;					// VBO Data Loading Procedure
+	PFNGLDELETEBUFFERSARBPROC glDeleteBuffers;			// VBO Deletion Procedure
+};
+#endif