Reply To: Test the water model only

Home Forums AR Sandbox Forum Test the water model only Reply To: Test the water model only

#102025
Oliver Kreylos
Keymaster

It’s possible, but not trivial. I made a separate application to do just that, but I never bundled it with the SARndbox package. Here is the source:

/***********************************************************************
FlowModel2D - Utility to visualize a 2D simulation of shallow water flow
using a Saint-Venant system.
Copyright (c) 2012-2015 Oliver Kreylos

This file is part of the Augmented Reality Sandbox (SARndbox).

The Augmented Reality Sandbox is free software; you can redistribute it
and/or modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

The Augmented Reality Sandbox 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 the Augmented Reality Sandbox; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
***********************************************************************/

#include <iostream>
#include <fstream>
#include <Misc/SizedTypes.h>
#include <Misc/FileNameExtensions.h>
#include <Misc/CreateNumberedFileName.h>
#include <IO/File.h>
#include <IO/ValueSource.h>
#include <Math/Math.h>
#include <Math/Constants.h>
#include <Math/Noise.h>
#include <GL/gl.h>
#include <GL/GLVertexTemplates.h>
#include <GL/GLNormalTemplates.h>
#include <GL/GLColorTemplates.h>
#include <GL/GLColor.h>
#include <GL/GLColorOperations.h>
#include <GL/GLMaterialTemplates.h>
#include <GL/GLContextData.h>
#include <GL/GLObject.h>
#include <GL/Extensions/GLARBVertexProgram.h>
#include <GLMotif/StyleSheet.h>
#include <GLMotif/PopupWindow.h>
#include <GLMotif/RowColumn.h>
#include <GLMotif/Label.h>
#include <GLMotif/TextFieldSlider.h>
#include <Vrui/Vrui.h>
#include <Vrui/Application.h>
#include <Vrui/DisplayState.h>
#include <Vrui/OpenFile.h>
#include <Kinect/FrameBuffer.h>

#include "Types.h"
#include "DepthImageRenderer.h"
#include "WaterTable2.h"
#include "SurfaceRenderer.h"
#include "WaterRenderer.h"

class FlowModel2D:public Vrui::Application,public GLObject
	{
	/* Embedded classes: */
	private:
	typedef double Scalar;
	typedef GLColor<GLfloat,3> Color;
	
	struct DataItem:public GLObject::DataItem
		{
		/* Elements: */
		public:
		unsigned int frameNumber; // Frame number to allow multi-view rendering
		};
	
	/* Elements: */
	private:
	double gridOrigin[2]; // Origin of the grid in world coordinates
	unsigned int gridSize[2]; // Number of grid cells
	Scalar cellSize[2]; // Width and height of a grid cell
	Kinect::FrameBuffer bathymetry; // Bathymetry grid
	DepthImageRenderer* depthImageRenderer; // Helper object for bathymetry management and rendering
	WaterTable2* waterTable; // Object to simulate water flow
	float inflow; // Amount of water flowing in at the right boundary
	double speed; // Simulation speed relative to real time
	SurfaceRenderer* bathymetryRenderer; // Helper object to render the bathymetry
	WaterRenderer* waterSurfaceRenderer; // Helper object to render the water surface
	unsigned int frameNumber; // Frame counter
	unsigned int saveWaterSurfaceFrame; // Frame number at which to save the current water surface
	
	/* Private methods: */
	void speedSliderCallback(GLMotif::TextFieldSlider::ValueChangedCallbackData* cbData);
	void addWater(GLContextData& contextData) const;
	void saveWaterGrid(const GLfloat* q,const char* bilFileName) const;
	
	/* Constructors and destructors: */
	public:
	FlowModel2D(int& argc,char**& argv);
	virtual ~FlowModel2D(void);
	
	/* Methods from Vrui::Application: */
	virtual void frame(void);
	virtual void display(GLContextData& contextData) const;
	virtual void eventCallback(EventID eventId,Vrui::InputDevice::ButtonCallbackData* cbData);
	
	/* Methods from GLObject: */
	virtual void initContext(GLContextData& contextData) const;
	};

/****************************
Methods of class FlowModel2D:
****************************/

void FlowModel2D::speedSliderCallback(GLMotif::TextFieldSlider::ValueChangedCallbackData* cbData)
	{
	/* Set the simulation speed: */
	speed=cbData->value;
	}

void FlowModel2D::addWater(GLContextData& contextData) const
	{
	#if 1
	glBegin(GL_LINES);
	glVertexAttrib1fARB(1,inflow);
	glVertex2d((double(gridSize[0])-0.5)*cellSize[0],469.5);
	glVertex2d((double(gridSize[0])-0.5)*cellSize[0],double(gridSize[1])*cellSize[1]);
	glEnd();
	#endif
	
	#if 0
	glBegin(GL_POINTS);
	glVertexAttrib1fARB(1,inflow);
	glVertex2d(257.5,1075.5);
	glEnd();
	#endif
	}

namespace {

/****************
Helper functions:
****************/

inline GLfloat writeWaterHeight(IO::File& bilFile,GLfloat water,GLfloat bathymetry)
	{
	bilFile.write<Misc::Float32>(water>bathymetry+1.0e-8f?water:-9999.0f);
	}

}

void FlowModel2D::saveWaterGrid(const GLfloat* q,const char* bilFileName) const
	{
	/* Construct the header file name: */
	const char* bilFileNameExt=Misc::getExtension(bilFileName);
	std::string hdrFileName(bilFileName,bilFileNameExt);
	hdrFileName.append(".hdr");
	
	/* Write the header file: */
	{
	std::ofstream hdrFile(hdrFileName.c_str());
	hdrFile<<"BYTEORDER I"<<std::endl;
	hdrFile<<"LAYOUT BIL"<<std::endl;
	hdrFile<<"NBANDS 1"<<std::endl;
	hdrFile<<"NBITS 32"<<std::endl;
	hdrFile<<"NCOLS "<<gridSize[0]<<std::endl;
	hdrFile<<"NROWS "<<gridSize[1]<<std::endl;
	hdrFile<<"BANDROWBYTES "<<gridSize[0]*4<<std::endl;
	hdrFile<<"TOTALROWBYTES "<<gridSize[0]*4<<std::endl;
	hdrFile<<"XDIM "<<cellSize[0]<<std::endl;
	hdrFile<<"YDIM "<<cellSize[1]<<std::endl;
	hdrFile<<"ULXMAP "<<gridOrigin[0]<<std::endl;
	hdrFile<<"ULYMAP "<<gridOrigin[1]<<std::endl;
	hdrFile<<"NODATA_VALUE -9999.0"<<std::endl;
	}
	
	/* Open the bil file: */
	IO::FilePtr bilFile(Vrui::openFile(bilFileName,IO::File::WriteOnly));
	bilFile->setEndianness(Misc::LittleEndian);
	
	/* Write the water surface, setting all "dry" vertices to the invalid data value: */
	const float* bathy=bathymetry.getData<float>();
	const GLfloat* qRowPtr=q+(gridSize[1]-1)*gridSize[0]*3;
	const float* bRowPtr=bathy+(gridSize[1]-2)*(gridSize[0]-1);
	
	const GLfloat* qPtr=qRowPtr;
	writeWaterHeight(*bilFile,qPtr[0],bRowPtr[0]);
	qPtr+=3;
	for(unsigned int x=1;x<gridSize[0]-1;++x,qPtr+=3)
		writeWaterHeight(*bilFile,qPtr[0],(bRowPtr[x-1]+bRowPtr[x])*0.5f);
	writeWaterHeight(*bilFile,qPtr[0],bRowPtr[gridSize[0]-2]);
	qPtr+=3;
	
	qRowPtr-=gridSize[0]*3;
	
	for(unsigned int y=1;y<gridSize[1]-1;++y,qRowPtr-=gridSize[0]*3,bRowPtr-=gridSize[0]-1)
		{
		const float* bNextRowPtr=bRowPtr-(gridSize[0]-1);
		qPtr=qRowPtr;
		writeWaterHeight(*bilFile,qPtr[0],(bRowPtr[0]+bNextRowPtr[0])*0.5f);
		qPtr+=3;
		for(unsigned int x=1;x<gridSize[0]-1;++x,qPtr+=3)
			writeWaterHeight(*bilFile,qPtr[0],(bRowPtr[x-1]+bRowPtr[x]+bNextRowPtr[x-1]+bNextRowPtr[x])*0.25f);
		writeWaterHeight(*bilFile,qPtr[0],(bRowPtr[gridSize[0]-2]+bNextRowPtr[gridSize[0]-2])*0.5f);
		qPtr+=3;
		}
	
	qPtr=qRowPtr;
	writeWaterHeight(*bilFile,qPtr[0],bRowPtr[0]);
	qPtr+=3;
	for(unsigned int x=1;x<gridSize[0]-1;++x,qPtr+=3)
		writeWaterHeight(*bilFile,qPtr[0],(bRowPtr[x-1]+bRowPtr[x])*0.5f);
	writeWaterHeight(*bilFile,qPtr[0],bRowPtr[gridSize[0]-2]);
	
	std::cout<<"Saved water surface as "<<bilFileName<<std::endl;
	}

FlowModel2D::FlowModel2D(int& argc,char**& argv)
	:Vrui::Application(argc,argv),
	 depthImageRenderer(0),waterTable(0),inflow(1.0f/16.0f),speed(0.0),bathymetryRenderer(0),waterSurfaceRenderer(0),
	 frameNumber(0U),
	 saveWaterSurfaceFrame(~0U)
	{
	/* Parse the command line: */
	gridSize[0]=256U;
	gridSize[1]=512U;
	cellSize[0]=cellSize[1]=Scalar(1.0);
	const char* demFileName=0;
	float boundaryHeight=0.0f;
	for(int i=1;i<argc;++i)
		{
		if(argv[i][0]=='-')
			{
			if(strcasecmp(argv[i]+1,"gridSize")==0)
				{
				for(int j=0;j<2;++j)
					{
					++i;
					gridSize[j]=(unsigned int)(atoi(argv[i]));
					}
				}
			else if(strcasecmp(argv[i]+1,"cellSize")==0)
				{
				for(int j=0;j<2;++j)
					{
					++i;
					cellSize[j]=Scalar(atof(argv[i]));
					}
				}
			else if(strcasecmp(argv[i]+1,"inflow")==0)
				{
				++i;
				inflow=float(atof(argv[i]));
				}
			else if(strcasecmp(argv[i]+1,"save")==0)
				{
				++i;
				saveWaterSurfaceFrame=atoi(argv[i]);
				}
			}
		else if(demFileName==0)
			demFileName=argv[i];
		}
	
	if(demFileName!=0)
		{
		/* Load a digital elevation model for the bathymetry: */
		const char* demFileNameExt=Misc::getExtension(demFileName);
		std::string hdrFileName(demFileName,demFileNameExt);
		hdrFileName.append(".hdr");
		IO::ValueSource hdrFile(Vrui::openFile(hdrFileName.c_str()));
		hdrFile.setPunctuation('\n',true);
		hdrFile.skipWs();
		
		Misc::Endianness demEndianness=Misc::HostEndianness;
		Scalar demGridOrigin[2]={Scalar(0),Scalar(0)};
		unsigned int demGridSize[2]={0U,0U};
		Scalar demCellSize[2]={Scalar(0),Scalar(0)};
		unsigned int demValueSize=0;
		while(!hdrFile.eof())
			{
			std::string token=hdrFile.readString();
			if(token!="\n")
				{
				if(token=="BYTEORDER")
					{
					std::string byteOrder=hdrFile.readString();
					if(byteOrder=="I")
						demEndianness=Misc::LittleEndian;
					else if(byteOrder=="M")
						demEndianness=Misc::BigEndian;
					else
						std::cerr<<"Unknown byte order "<<byteOrder<<" in DEM file "<<demFileName<<std::endl;
					}
				else if(token=="LAYOUT")
					{
					std::string layout=hdrFile.readString();
					if(layout!="BIL"&&layout!="BIP"&&layout!="BSQ")
						std::cerr<<"Unknown layout "<<layout<<" in DEM file "<<demFileName<<std::endl;
					}
				else if(token=="NBANDS")
					{
					if(hdrFile.readUnsignedInteger()!=1)
						std::cerr<<"Unsupported number of bands in DEM file "<<demFileName<<std::endl;
					}
				else if(token=="NBITS")
					{
					unsigned int numBits=hdrFile.readUnsignedInteger();
					demValueSize=numBits/8;
					if(numBits%8!=0U||demValueSize<1U||demValueSize==3U||demValueSize>4U)
						{
						std::cerr<<"Unsupported value size "<<numBits<<" in DEM file "<<demFileName<<std::endl;
						demValueSize=0U;
						}
					}
				else if(token=="NCOLS")
					demGridSize[0]=hdrFile.readUnsignedInteger();
				else if(token=="NROWS")
					demGridSize[1]=hdrFile.readUnsignedInteger();
				else if(token=="BANDROWBYTES")
					{
					unsigned int bandRowBytes=hdrFile.readUnsignedInteger();
					if(bandRowBytes!=demGridSize[0]*demValueSize)
						std::cerr<<"Mismatching band row size "<<bandRowBytes<<" in DEM file "<<demFileName<<std::endl;
					}
				else if(token=="TOTALROWBYTES")
					{
					unsigned int totalRowBytes=hdrFile.readUnsignedInteger();
					if(totalRowBytes!=demGridSize[0]*demValueSize)
						std::cerr<<"Mismatching total row size "<<totalRowBytes<<" in DEM file "<<demFileName<<std::endl;
					}
				else if(token=="ULXMAP")
					demGridOrigin[0]=hdrFile.readNumber();
				else if(token=="ULYMAP")
					demGridOrigin[1]=hdrFile.readNumber();
				else if(token=="XDIM")
					demCellSize[0]=hdrFile.readNumber();
				else if(token=="YDIM")
					demCellSize[1]=hdrFile.readNumber();
				else
					std::cerr<<"Unknown token "<<token<<" in DEM file "<<demFileName<<std::endl;
				}
			}
		if(demGridSize[0]!=0U&&demGridSize[1]!=0U&&demCellSize[0]>Scalar(0)&&demCellSize[1]>Scalar(0)&&demValueSize!=0U)
			{
			/* Load the DEM: */
			IO::FilePtr demFile=Vrui::openFile(demFileName);
			demFile->setEndianness(demEndianness);
			
			for(int i=0;i<2;++i)
				{
				gridOrigin[i]=demGridOrigin[i]-Math::div2(demCellSize[i]);
				gridSize[i]=demGridSize[i]+1;
				cellSize[i]=demCellSize[i];
				}
			bathymetry=Kinect::FrameBuffer(demGridSize[0],demGridSize[1],demGridSize[1]*demGridSize[0]*sizeof(float));
			float* b=bathymetry.getData<float>();
			if(demValueSize==1U)
				{
				Misc::UInt8* row=new Misc::UInt8[demGridSize[0]];
				for(unsigned int y=0;y<demGridSize[1];++y)
					{
					demFile->read<Misc::UInt8>(row,demGridSize[0]);
					float* bPtr=b+(demGridSize[1]-1-y)*demGridSize[0];
					for(unsigned int x=0;x<demGridSize[0];++x,++bPtr)
						*bPtr=float(row[x]);
					}
				delete[] row;
				}
			else if(demValueSize==2U)
				{
				Misc::SInt16* row=new Misc::SInt16[demGridSize[0]];
				for(unsigned int y=0;y<demGridSize[1];++y)
					{
					demFile->read<Misc::SInt16>(row,demGridSize[0]);
					float* bPtr=b+(demGridSize[1]-1-y)*demGridSize[0];
					for(unsigned int x=0;x<demGridSize[0];++x,++bPtr)
						*bPtr=float(row[x]);
					}
				delete[] row;
				}
			else
				{
				Misc::SInt32* row=new Misc::SInt32[demGridSize[0]];
				for(unsigned int y=0;y<demGridSize[1];++y)
					{
					demFile->read<Misc::SInt32>(row,demGridSize[0]);
					float* bPtr=b+(demGridSize[1]-1-y)*demGridSize[0];
					for(unsigned int x=0;x<demGridSize[0];++x,++bPtr)
							*bPtr=float(row[x])*0.001f;
					}
				delete[] row;
				}
			}
		}
	else
		{
		/* Create custom bathymetry: */
		bathymetry=Kinect::FrameBuffer(gridSize[0]-1,gridSize[1]-1,(gridSize[1]-1)*(gridSize[0]-1)*sizeof(float));
		float* b=bathymetry.getData<float>();
		
		GLfloat* bPtr=b;
		for(int y=0;y<gridSize[1]-1;++y)
			for(int x=0;x<gridSize[0]-1;++x,++bPtr)
				{
				#if 0
				
				/* Flat bathymetry: */
				*bPtr=0.0f;
				
				#elif 1
				
				/* Swimming pool: */
				if(x>0&&x<gridSize[0]-2&&y>0&&y<gridSize[1]-2)
					{
					#if 0
					/* Gaussian blob island: */
					GLfloat cx=GLfloat(gridSize[0])*0.5f;
					GLfloat cy=GLfloat(gridSize[1])*0.5f;
					GLfloat arg=Math::exp(-(Math::sqr(GLfloat(x)-cx)+Math::sqr(GLfloat(y)-cy))/Math::sqr(20.0f))*50.0f;
					*bPtr=arg;
					#else
					*bPtr=0.0f;
					#endif
					}
				else
					*bPtr=100.0f;
				
				#elif 0
				
				/* Gaussian blob island: */
				GLfloat cx=GLfloat(gridSize[0])*0.5f;
				GLfloat cy=GLfloat(gridSize[1])*0.5f;
				GLfloat arg=Math::exp(-(Math::sqr(GLfloat(x)-cx)+Math::sqr(GLfloat(y)-cy))/Math::sqr(20.0f))*25.0f;
				*bPtr=arg;
				
				#else
				
				/* Reservoir with outflow channel: */
				if(x==0||x==gridSize[0]-2||y==0||y==gridSize[1]-2)
					*bPtr=50.0f;
				else if(x>=5&&x<=gridSize[0]-7&&y>=5&&y<gridSize[1]/4)
					*bPtr=0.0f;
				else if(x>=gridSize[0]/2-15&&x<gridSize[0]/2+35&&y>=gridSize[1]/4+5)
					*bPtr=0.0f;
				else if(y>=gridSize[1]/4+5)
					*bPtr=5.0f;
				else if(x>=gridSize[0]/2-10&&x<gridSize[0]/2+30&&y>=5)
					*bPtr=0.0f;
				else
					*bPtr=50.0f;
				
				#endif
				}
		}
	
	/*********************************************************************
	Visualize GPU-based water flow simulation:
	*********************************************************************/
	
	/* Create a depth image renderer: */
	unsigned int bathySize[2];
	for(int i=0;i<2;++i)
		bathySize[i]=gridSize[i]-1;
	depthImageRenderer=new DepthImageRenderer(bathySize);
	
	/* Create a transformation from depth image space into water table space: */
	PTransform depthProjection=PTransform::identity;
	PTransform::Matrix& dpm=depthProjection.getMatrix();
	dpm(0,0)=cellSize[0];
	dpm(0,3)=Math::div2(cellSize[0]);
	dpm(1,1)=cellSize[1];
	dpm(1,3)=Math::div2(cellSize[1]);
	depthImageRenderer->setDepthProjection(depthProjection);
	
	/* Create a base plane equation for elevation rendering: */
	depthImageRenderer->setBasePlane(Plane(Plane::Vector(0,0,1),0));
	if(bathymetry.isValid())
		depthImageRenderer->setDepthImage(bathymetry);
	
	/* Create the water flow simulator: */
	GLfloat wtCellSize[2];
	for(int i=0;i<2;++i)
		wtCellSize[i]=GLfloat(cellSize[i]);
	waterTable=new WaterTable2(gridSize[0],gridSize[1],wtCellSize);
	waterTable->setDryBoundary(false);
	
	/* Ensure that the water table's OpenGL state is initialized before this object's: */
	dependsOn(waterTable);
	
	/* Add a rendering function to add water at the right edge: */
	waterTable->addRenderFunction(Misc::createFunctionCall(this,&FlowModel2D::addWater));
	
	/* Create a renderer for the bathymetry: */
	bathymetryRenderer=new SurfaceRenderer(depthImageRenderer);
	bathymetryRenderer->setDrawContourLines(false);
	bathymetryRenderer->setIlluminate(true);
	// bathymetryRenderer->setWaterTable(waterTable);
	// bathymetryRenderer->setWaterOpacity(5.0f);
	
	/* Create a renderer for the water surface: */
	waterSurfaceRenderer=new WaterRenderer(waterTable);
	
	/* Create a GUI: */
	GLMotif::PopupWindow* controlPopup=new GLMotif::PopupWindow("ControlPopup",Vrui::getWidgetManager(),"Simulation Control");
	GLMotif::RowColumn* control=new GLMotif::RowColumn("Control",controlPopup,false);
	control->setOrientation(GLMotif::RowColumn::HORIZONTAL);
	control->setPacking(GLMotif::RowColumn::PACK_TIGHT);
	
	new GLMotif::Label("SpeedLabel",control,"Speed");
	
	GLMotif::TextFieldSlider* speedSlider=new GLMotif::TextFieldSlider("SpeedSlider",control,5,10.0f*Vrui::getUiStyleSheet()->fontHeight);
	speedSlider->setSliderMapping(GLMotif::TextFieldSlider::LINEAR);
	speedSlider->setValueType(GLMotif::TextFieldSlider::FLOAT);
	speedSlider->setValueRange(0.0,2.0,0.125);
	speedSlider->getSlider()->addNotch(1.0);
	speedSlider->setValue(speed);
	speedSlider->getValueChangedCallbacks().add(this,&FlowModel2D::speedSliderCallback);
	
	control->manageChild();
	
	Vrui::popupPrimaryWidget(controlPopup);
	
	/* Create a tool class to save the current water surface as a BIL file: */
	addEventTool("Save Water Surface",0,0);
	
	/* Initialize the navigation transformation: */
	Vrui::Scalar s0=Vrui::Scalar(gridSize[0]/2)*Vrui::Scalar(cellSize[0]);
	Vrui::Scalar s1=Vrui::Scalar(gridSize[1]/2)*Vrui::Scalar(cellSize[1]);
	Vrui::setNavigationTransformation(Vrui::Point(s0,s1,Vrui::Scalar(10)),Math::sqrt(Math::sqr(s0)+Math::sqr(s1)),Vrui::Vector(0,1,0));
	}

FlowModel2D::~FlowModel2D(void)
	{
	delete waterSurfaceRenderer;
	delete bathymetryRenderer;
	delete waterTable;
	delete depthImageRenderer;
	}

void FlowModel2D::frame(void)
	{
	/* Request another simulation step: */
	++frameNumber;
	
	/* Request another frame: */
	Vrui::scheduleUpdate(Vrui::getNextAnimationTime());
	}

void FlowModel2D::display(GLContextData& contextData) const
	{
	/* Get the data item: */
	DataItem* dataItem=contextData.retrieveDataItem<DataItem>(this);
	
	if(dataItem->frameNumber!=frameNumber)
		{
		/* Run a simulation step: */
		#if 1
		
		GLfloat totalTimeStep=GLfloat(Vrui::getFrameTime()*speed);
		unsigned int numSteps=0;
		while(numSteps<50&&totalTimeStep>1.0e-8f)
			{
			waterTable->setMaxStepSize(totalTimeStep);
			totalTimeStep-=waterTable->runSimulationStep(false,contextData);
			++numSteps;
			}
		if(totalTimeStep>1.0e-8f)
			std::cout<<"Ran out of time by "<<totalTimeStep<<std::endl;
		
		#else
		
		for(int i=0;i<100;++i)
			waterTable->runSimulationStep(false,contextData);
		
		#endif
		
		if(saveWaterSurfaceFrame==frameNumber)
			{
			/* Bind the quantity texture: */
			glActiveTextureARB(GL_TEXTURE0_ARB);
			waterTable->bindQuantityTexture(contextData);
			
			/* Read the current quantity texture: */
			GLfloat* q=new GLfloat[gridSize[1]*gridSize[0]*3];
			glGetTexImage(GL_TEXTURE_RECTANGLE_ARB,0,GL_RGB,GL_FLOAT,q);
			
			/* Save the water surface as a DEM file: */
			char bilFileName[256];
			saveWaterGrid(q,Misc::createNumberedFileName("WaterSurface.bil",4,bilFileName));
			
			/* Clean up: */
			delete[] q;
			}
		
		dataItem->frameNumber=frameNumber;
		}
	
	/* Get the display state object: */
	const Vrui::DisplayState& ds=Vrui::getDisplayState(contextData);
	
	/* Draw the bathymetry surface: */
	glMaterialAmbientAndDiffuse(GLMaterialEnums::FRONT,GLColor<GLfloat,4>(0.7f,0.5f,0.3f));
	glMaterialSpecular(GLMaterialEnums::FRONT,GLColor<GLfloat,4>(0.0f,0.0f,0.0f));
	glMaterialShininess(GLMaterialEnums::FRONT,0.0f);
	bathymetryRenderer->renderSinglePass(ds.viewport,ds.projection,ds.modelviewNavigational,contextData);
	
	/* Draw the water surface: */
	glMaterialAmbientAndDiffuse(GLMaterialEnums::FRONT,GLColor<GLfloat,4>(0.4f,0.5f,0.8f));
	glMaterialSpecular(GLMaterialEnums::FRONT,GLColor<GLfloat,4>(1.0f,1.0f,1.0f));
	glMaterialShininess(GLMaterialEnums::FRONT,64.0f);
	waterSurfaceRenderer->render(ds.projection,ds.modelviewNavigational,contextData);
	}

void FlowModel2D::eventCallback(EventID eventId,Vrui::InputDevice::ButtonCallbackData* cbData)
	{
	/* Save the water surface on this frame: */
	saveWaterSurfaceFrame=frameNumber+1; // Frame number will be incremented right after this
	}

void FlowModel2D::initContext(GLContextData& contextData) const
	{
	/* Initialize required OpenGL extensions: */
	GLARBVertexProgram::initExtension();
	
	/* Create a new data item: */
	DataItem* dataItem=new DataItem;
	dataItem->frameNumber=0U;
	
	/* Store the data item: */
	contextData.addDataItem(this,dataItem);
	
	if(bathymetry.isValid())
		{
		/* Upload the bathymetry grid into the water table: */
		waterTable->updateBathymetry(bathymetry.getData<float>(),contextData);
		}
	
	/* Initialize the water table's water level surface: */
	GLfloat cx=GLfloat(gridSize[0])*0.25f;
	GLfloat cy=GLfloat(gridSize[1])*0.333f;
	GLfloat* q=new GLfloat[gridSize[1]*gridSize[0]];
	GLfloat* qPtr=q;
	for(int y=0;y<gridSize[1];++y)
		for(int x=0;x<gridSize[0];++x,++qPtr)
			{
			#if 0
			
			/* Dam failure: */
			if(y<gridSize[1]/4)
				qPtr[0]=40.0f;
			else
				qPtr[0]=0.0f;
			
			#elif 0
			
			/* Gaussian water blob: */
			GLfloat arg=Math::exp(-(Math::sqr(GLfloat(x)-cx)+Math::sqr(GLfloat(y)-cy))/Math::sqr(16.0f))*40.0f+10.0f;
			qPtr[0]=arg;
			
			#elif 0
			
			/* Rectangular water blob: */
			if(x>=gridSize[0]/4-10&&x<gridSize[0]/4+10&&y>=gridSize[1]/3-10&&y<gridSize[1]/3+10)
				qPtr[0]=40.0f;
			else
				qPtr[0]=0.0f;
			
			#else
			
			/* Flat surface: */
			qPtr[0]=-1000.0f;
			
			#endif
			}
	
	waterTable->setWaterLevel(q,contextData);
	delete[] q;
	}

VRUI_APPLICATION_RUN(FlowModel2D)

It’s a Vrui application like SARndbox; building it is left as an exercise for the reader.

Comments are closed.