Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Created Fluid3D repository

commit 121d23460a721fb4a6f648a83060cd12a1e17e59 0 parents
@christopherbatty authored
Showing with 9,007 additions and 0 deletions.
  1. +26 −0 Fluid3D.sln
  2. BIN  Fluid3D.suo
  3. +252 −0 Fluid3D.vcproj
  4. +793 −0 array1.h
  5. +280 −0 array2.h
  6. +63 −0 array2_utils.h
  7. +272 −0 array3.h
  8. +72 −0 array3_utils.h
  9. +610 −0 fluidsim.cpp
  10. +72 −0 fluidsim.h
  11. +103 −0 levelset_util.cpp
  12. +7 −0 levelset_util.h
  13. +93 −0 main.cpp
  14. +54 −0 pcgsolver/blas_wrapper.h
  15. +312 −0 pcgsolver/pcg_solver.h
  16. +290 −0 pcgsolver/sparse_matrix.h
  17. +150 −0 sparse/.#sparseilu.h.1.10
  18. +10 −0 sparse/CVS/Entries
  19. +9 −0 sparse/CVS/Entries.Extra
  20. +9 −0 sparse/CVS/Entries.Extra.Old
  21. +10 −0 sparse/CVS/Entries.Old
  22. +1 −0  sparse/CVS/Repository
  23. +1 −0  sparse/CVS/Root
  24. +194 −0 sparse/cgsolver.h
  25. +154 −0 sparse/crsolver.h
  26. +69 −0 sparse/ilu.m
  27. +336 −0 sparse/ilu_mex.c
  28. +152 −0 sparse/incomplete_qr.cpp
  29. +15 −0 sparse/incomplete_qr.h
  30. +182 −0 sparse/sparseilu.h
  31. +429 −0 sparse/sparsematrix.h
  32. +208 −0 sparse/sparsemilu.h
  33. +470 −0 util.h
  34. +472 −0 vec.h
  35. +211 −0 viewpls3D/ViewFluid3D.vcproj
  36. +1,114 −0 viewpls3D/gluvi.cpp
  37. +200 −0 viewpls3D/gluvi.h
  38. +370 −0 viewpls3D/main.cpp
  39. +470 −0 viewpls3D/util.h
  40. +472 −0 viewpls3D/vec.h
26 Fluid3D.sln
@@ -0,0 +1,26 @@
+
+Microsoft Visual Studio Solution File, Format Version 10.00
+# Visual Studio 2008
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Fluid3D", "Fluid3D.vcproj", "{CA7F57AF-7216-4435-A1FB-7679ED1868A2}"
+EndProject
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "ViewFluid3D", "viewpls3D\ViewFluid3D.vcproj", "{1FF2B89E-7C6D-4C00-9A74-80FFBD8D10B2}"
+EndProject
+Global
+ GlobalSection(SolutionConfigurationPlatforms) = preSolution
+ Debug|Win32 = Debug|Win32
+ Release|Win32 = Release|Win32
+ EndGlobalSection
+ GlobalSection(ProjectConfigurationPlatforms) = postSolution
+ {CA7F57AF-7216-4435-A1FB-7679ED1868A2}.Debug|Win32.ActiveCfg = Debug|Win32
+ {CA7F57AF-7216-4435-A1FB-7679ED1868A2}.Debug|Win32.Build.0 = Debug|Win32
+ {CA7F57AF-7216-4435-A1FB-7679ED1868A2}.Release|Win32.ActiveCfg = Release|Win32
+ {CA7F57AF-7216-4435-A1FB-7679ED1868A2}.Release|Win32.Build.0 = Release|Win32
+ {1FF2B89E-7C6D-4C00-9A74-80FFBD8D10B2}.Debug|Win32.ActiveCfg = Debug|Win32
+ {1FF2B89E-7C6D-4C00-9A74-80FFBD8D10B2}.Debug|Win32.Build.0 = Debug|Win32
+ {1FF2B89E-7C6D-4C00-9A74-80FFBD8D10B2}.Release|Win32.ActiveCfg = Release|Win32
+ {1FF2B89E-7C6D-4C00-9A74-80FFBD8D10B2}.Release|Win32.Build.0 = Release|Win32
+ EndGlobalSection
+ GlobalSection(SolutionProperties) = preSolution
+ HideSolutionNode = FALSE
+ EndGlobalSection
+EndGlobal
BIN  Fluid3D.suo
Binary file not shown
252 Fluid3D.vcproj
@@ -0,0 +1,252 @@
+<?xml version="1.0" encoding="Windows-1252"?>
+<VisualStudioProject
+ ProjectType="Visual C++"
+ Version="9.00"
+ Name="Fluid3D"
+ ProjectGUID="{CA7F57AF-7216-4435-A1FB-7679ED1868A2}"
+ RootNamespace="Variational"
+ TargetFrameworkVersion="131072"
+ >
+ <Platforms>
+ <Platform
+ Name="Win32"
+ />
+ </Platforms>
+ <ToolFiles>
+ </ToolFiles>
+ <Configurations>
+ <Configuration
+ Name="Debug|Win32"
+ OutputDirectory="$(SolutionDir)$(ConfigurationName)"
+ IntermediateDirectory="$(ConfigurationName)"
+ ConfigurationType="1"
+ >
+ <Tool
+ Name="VCPreBuildEventTool"
+ />
+ <Tool
+ Name="VCCustomBuildTool"
+ />
+ <Tool
+ Name="VCXMLDataGeneratorTool"
+ />
+ <Tool
+ Name="VCWebServiceProxyGeneratorTool"
+ />
+ <Tool
+ Name="VCMIDLTool"
+ />
+ <Tool
+ Name="VCCLCompilerTool"
+ Optimization="0"
+ AdditionalIncludeDirectories=""
+ PreprocessorDefinitions="USE_AMD_BLAS;WIN32;_CRT_SECURE_NO_DEPRECATE"
+ DebugInformationFormat="4"
+ />
+ <Tool
+ Name="VCManagedResourceCompilerTool"
+ />
+ <Tool
+ Name="VCResourceCompilerTool"
+ />
+ <Tool
+ Name="VCPreLinkEventTool"
+ />
+ <Tool
+ Name="VCLinkerTool"
+ AdditionalLibraryDirectories=""
+ GenerateDebugInformation="true"
+ GenerateMapFile="true"
+ RandomizedBaseAddress="1"
+ DataExecutionPrevention="0"
+ />
+ <Tool
+ Name="VCALinkTool"
+ />
+ <Tool
+ Name="VCManifestTool"
+ />
+ <Tool
+ Name="VCXDCMakeTool"
+ />
+ <Tool
+ Name="VCBscMakeTool"
+ />
+ <Tool
+ Name="VCFxCopTool"
+ />
+ <Tool
+ Name="VCAppVerifierTool"
+ />
+ <Tool
+ Name="VCPostBuildEventTool"
+ />
+ </Configuration>
+ <Configuration
+ Name="Release|Win32"
+ OutputDirectory="$(SolutionDir)$(ConfigurationName)"
+ IntermediateDirectory="$(ConfigurationName)"
+ ConfigurationType="1"
+ WholeProgramOptimization="1"
+ >
+ <Tool
+ Name="VCPreBuildEventTool"
+ />
+ <Tool
+ Name="VCCustomBuildTool"
+ />
+ <Tool
+ Name="VCXMLDataGeneratorTool"
+ />
+ <Tool
+ Name="VCWebServiceProxyGeneratorTool"
+ />
+ <Tool
+ Name="VCMIDLTool"
+ />
+ <Tool
+ Name="VCCLCompilerTool"
+ Optimization="2"
+ FavorSizeOrSpeed="1"
+ EnableFiberSafeOptimizations="true"
+ WholeProgramOptimization="true"
+ AdditionalIncludeDirectories=""
+ PreprocessorDefinitions="WIN32;_CRT_SECURE_NO_DEPRECATE;USE_MKL_CBLAS;USE_MKL_LAPACK;_SCL_SECURE_NO_WARNINGS"
+ WarningLevel="3"
+ Detect64BitPortabilityProblems="false"
+ />
+ <Tool
+ Name="VCManagedResourceCompilerTool"
+ />
+ <Tool
+ Name="VCResourceCompilerTool"
+ />
+ <Tool
+ Name="VCPreLinkEventTool"
+ />
+ <Tool
+ Name="VCLinkerTool"
+ AdditionalLibraryDirectories=""
+ RandomizedBaseAddress="1"
+ DataExecutionPrevention="0"
+ />
+ <Tool
+ Name="VCALinkTool"
+ />
+ <Tool
+ Name="VCManifestTool"
+ />
+ <Tool
+ Name="VCXDCMakeTool"
+ />
+ <Tool
+ Name="VCBscMakeTool"
+ />
+ <Tool
+ Name="VCFxCopTool"
+ />
+ <Tool
+ Name="VCAppVerifierTool"
+ />
+ <Tool
+ Name="VCPostBuildEventTool"
+ />
+ </Configuration>
+ </Configurations>
+ <References>
+ </References>
+ <Files>
+ <Filter
+ Name="Source Files"
+ Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
+ UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"
+ >
+ <File
+ RelativePath=".\fluidsim.cpp"
+ >
+ </File>
+ <File
+ RelativePath=".\levelset_util.cpp"
+ >
+ </File>
+ <File
+ RelativePath=".\main.cpp"
+ >
+ </File>
+ </Filter>
+ <Filter
+ Name="Header Files"
+ Filter="h;hpp;hxx;hm;inl;inc;xsd"
+ UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"
+ >
+ <File
+ RelativePath=".\array1.h"
+ >
+ </File>
+ <File
+ RelativePath=".\array2.h"
+ >
+ </File>
+ <File
+ RelativePath=".\array2_utils.h"
+ >
+ </File>
+ <File
+ RelativePath=".\array3_utils.h"
+ >
+ </File>
+ <File
+ RelativePath=".\fluidsim.h"
+ >
+ </File>
+ <File
+ RelativePath=".\gluvi.h"
+ >
+ </File>
+ <File
+ RelativePath=".\hashtable.h"
+ >
+ </File>
+ <File
+ RelativePath=".\levelset_util.h"
+ >
+ </File>
+ <File
+ RelativePath=".\openglutils.h"
+ >
+ </File>
+ <File
+ RelativePath=".\util.h"
+ >
+ </File>
+ <File
+ RelativePath=".\vec.h"
+ >
+ </File>
+ <Filter
+ Name="pcgsolver"
+ >
+ <File
+ RelativePath=".\pcgsolver\blas_wrapper.h"
+ >
+ </File>
+ <File
+ RelativePath=".\pcgsolver\pcg_solver.h"
+ >
+ </File>
+ <File
+ RelativePath=".\pcgsolver\sparse_matrix.h"
+ >
+ </File>
+ </Filter>
+ </Filter>
+ <Filter
+ Name="Resource Files"
+ Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav"
+ UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"
+ >
+ </Filter>
+ </Files>
+ <Globals>
+ </Globals>
+</VisualStudioProject>
793 array1.h
@@ -0,0 +1,793 @@
+#ifndef ARRAY1_H
+#define ARRAY1_H
+
+#include <algorithm>
+#include <cassert>
+#include <climits>
+#include <cstdlib>
+#include <iostream>
+#include <stdexcept>
+#include <vector>
+
+// In this file:
+// Array1<T>: a dynamic 1D array for plain-old-data (not objects)
+// WrapArray1<T>: a 1D array wrapper around an existing array (perhaps objects, perhaps data)
+// For the most part std::vector operations are supported, though for the Wrap version
+// note that memory is never allocated/deleted and constructor/destructors are never called
+// from within the class, thus only shallow copies can be made and some operations such as
+// resize() and push_back() are limited.
+// Note: for the most part assertions are done with assert(), not exceptions...
+
+// gross template hacking to determine if a type is integral or not
+struct Array1True {};
+struct Array1False {};
+template<typename T> struct Array1IsIntegral{ typedef Array1False type; }; // default: no (specializations to yes follow)
+template<> struct Array1IsIntegral<bool>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<signed char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<short>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned short>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<int>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned int>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<long long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned long long>{ typedef Array1True type; };
+
+//============================================================================
+template<typename T>
+struct Array1
+{
+ // STL-friendly typedefs
+
+ typedef T* iterator;
+ typedef const T* const_iterator;
+ typedef unsigned long size_type;
+ typedef long difference_type;
+ typedef T& reference;
+ typedef const T& const_reference;
+ typedef T value_type;
+ typedef T* pointer;
+ typedef const T* const_pointer;
+ typedef std::reverse_iterator<iterator> reverse_iterator;
+ typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
+
+ // the actual representation
+
+ unsigned long n;
+ unsigned long max_n;
+ T* data;
+
+ // STL vector's interface, with additions, but only valid when used with plain-old-data
+
+ Array1(void)
+ : n(0), max_n(0), data(0)
+ {}
+
+ // note: default initial values are zero
+ Array1(unsigned long n_)
+ : n(0), max_n(0), data(0)
+ {
+ if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ data=(T*)std::calloc(n_, sizeof(T));
+ if(!data) throw std::bad_alloc();
+ n=n_;
+ max_n=n_;
+ }
+
+ Array1(unsigned long n_, const T& value)
+ : n(0), max_n(0), data(0)
+ {
+ if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ data=(T*)std::calloc(n_, sizeof(T));
+ if(!data) throw std::bad_alloc();
+ n=n_;
+ max_n=n_;
+ for(unsigned long i=0; i<n; ++i) data[i]=value;
+ }
+
+ Array1(unsigned long n_, const T& value, unsigned long max_n_)
+ : n(0), max_n(0), data(0)
+ {
+ assert(n_<=max_n_);
+ if(max_n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ data=(T*)std::calloc(max_n_, sizeof(T));
+ if(!data) throw std::bad_alloc();
+ n=n_;
+ max_n=max_n_;
+ for(unsigned long i=0; i<n; ++i) data[i]=value;
+ }
+
+ Array1(unsigned long n_, const T* data_)
+ : n(0), max_n(0), data(0)
+ {
+ if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ data=(T*)std::calloc(n_, sizeof(T));
+ if(!data) throw std::bad_alloc();
+ n=n_;
+ max_n=n_;
+ assert(data_);
+ std::memcpy(data, data_, n*sizeof(T));
+ }
+
+ Array1(unsigned long n_, const T* data_, unsigned long max_n_)
+ : n(0), max_n(0), data(0)
+ {
+ assert(n_<=max_n_);
+ if(max_n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ data=(T*)std::calloc(max_n_, sizeof(T));
+ if(!data) throw std::bad_alloc();
+ max_n=max_n_;
+ n=n_;
+ assert(data_);
+ std::memcpy(data, data_, n*sizeof(T));
+ }
+
+ Array1(const Array1<T> &x)
+ : n(0), max_n(0), data(0)
+ {
+ data=(T*)std::malloc(x.n*sizeof(T));
+ if(!data) throw std::bad_alloc();
+ n=x.n;
+ max_n=x.n;
+ std::memcpy(data, x.data, n*sizeof(T));
+ }
+
+ ~Array1(void)
+ {
+ std::free(data);
+#ifndef NDEBUG
+ data=0;
+ n=max_n=0;
+#endif
+ }
+
+ const T& operator[](unsigned long i) const
+ { return data[i]; }
+
+ T& operator[](unsigned long i)
+ { return data[i]; }
+
+ // these are range-checked (in debug mode) versions of operator[], like at()
+ const T& operator()(unsigned long i) const
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ T& operator()(unsigned long i)
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ Array1<T>& operator=(const Array1<T>& x)
+ {
+ if(max_n<x.n){
+ T* new_data=(T*)std::malloc(x.n*sizeof(T));
+ if(!new_data) throw std::bad_alloc();
+ std::free(data);
+ data=new_data;
+ max_n=x.n;
+ }
+ n=x.n;
+ std::memcpy(data, x.data, n*sizeof(T));
+ return *this;
+ }
+
+ bool operator==(const Array1<T>& x) const
+ {
+ if(n!=x.n) return false;
+ for(unsigned long i=0; i<n; ++i) if(!(data[i]==x.data[i])) return false;
+ return true;
+ }
+
+ bool operator!=(const Array1<T>& x) const
+ {
+ if(n!=x.n) return true;
+ for(unsigned long i=0; i<n; ++i) if(data[i]!=x.data[i]) return true;
+ return false;
+ }
+
+ bool operator<(const Array1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]<x[i]) return true;
+ else if(x[i]<data[i]) return false;
+ }
+ return n<x.n;
+ }
+
+ bool operator>(const Array1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]>x[i]) return true;
+ else if(x[i]>data[i]) return false;
+ }
+ return n>x.n;
+ }
+
+ bool operator<=(const Array1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]<x[i]) return true;
+ else if(x[i]<data[i]) return false;
+ }
+ return n<=x.n;
+ }
+
+ bool operator>=(const Array1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]>x[i]) return true;
+ else if(x[i]>data[i]) return false;
+ }
+ return n>=x.n;
+ }
+
+ void add_unique(const T& value)
+ {
+ for(unsigned long i=0; i<n; ++i) if(data[i]==value) return;
+ if(n==max_n) grow();
+ data[n++]=value;
+ }
+
+ void assign(const T& value)
+ { for(unsigned long i=0; i<n; ++i) data[i]=value; }
+
+ void assign(unsigned long num, const T& value)
+ { fill(num, value); }
+
+ // note: copydata may not alias this array's data, and this should not be
+ // used when T is a full object (which defines its own copying operation)
+ void assign(unsigned long num, const T* copydata)
+ {
+ assert(num==0 || copydata);
+ if(num>max_n){
+ if(num>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ std::free(data);
+ data=(T*)std::malloc(num*sizeof(T));
+ if(!data) throw std::bad_alloc();
+ max_n=num;
+ }
+ n=num;
+ std::memcpy(data, copydata, n*sizeof(T));
+ }
+
+ template<typename InputIterator>
+ void assign(InputIterator first, InputIterator last)
+ { assign_(first, last, typename Array1IsIntegral<InputIterator>::type()); }
+
+ template<typename InputIterator>
+ void assign_(InputIterator first, InputIterator last, Array1True check)
+ { fill(first, last); }
+
+ template<typename InputIterator>
+ void assign_(InputIterator first, InputIterator last, Array1False check)
+ {
+ unsigned long i=0;
+ InputIterator p=first;
+ for(; p!=last; ++p, ++i){
+ if(i==max_n) grow();
+ data[i]=*p;
+ }
+ n=i;
+ }
+
+ const T& at(unsigned long i) const
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ T& at(unsigned long i)
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ const T& back(void) const
+ {
+ assert(data && n>0);
+ return data[n-1];
+ }
+
+ T& back(void)
+ {
+ assert(data && n>0);
+ return data[n-1];
+ }
+
+ const T* begin(void) const
+ { return data; }
+
+ T* begin(void)
+ { return data; }
+
+ unsigned long capacity(void) const
+ { return max_n; }
+
+ void clear(void)
+ {
+ std::free(data);
+ data=0;
+ max_n=0;
+ n=0;
+ }
+
+ bool empty(void) const
+ { return n==0; }
+
+ const T* end(void) const
+ { return data+n; }
+
+ T* end(void)
+ { return data+n; }
+
+ void erase(unsigned long index)
+ {
+ assert(index<n);
+ for(unsigned long i=index; i<n-1; ++i)
+ data[i]=data[i-1];
+ pop_back();
+ }
+
+ void fill(unsigned long num, const T& value)
+ {
+ if(num>max_n){
+ if(num>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ std::free(data);
+ data=(T*)std::malloc(num*sizeof(T));
+ if(!data) throw std::bad_alloc();
+ max_n=num;
+ }
+ n=num;
+ for(unsigned long i=0; i<n; ++i) data[i]=value;
+ }
+
+ const T& front(void) const
+ {
+ assert(n>0);
+ return *data;
+ }
+
+ T& front(void)
+ {
+ assert(n>0);
+ return *data;
+ }
+
+ void grow(void)
+ {
+ unsigned long new_size=(max_n*sizeof(T)<ULONG_MAX/2 ? 2*max_n+1 : ULONG_MAX/sizeof(T));
+ T *new_data=(T*)std::realloc(data, new_size*sizeof(T));
+ if(!new_data) throw std::bad_alloc();
+ data=new_data;
+ max_n=new_size;
+ }
+
+ void insert(unsigned long index, const T& entry)
+ {
+ assert(index<=n);
+ push_back(back());
+ for(unsigned long i=n-1; i>index; --i)
+ data[i]=data[i-1];
+ data[index]=entry;
+ }
+
+ unsigned long max_size(void) const
+ { return ULONG_MAX/sizeof(T); }
+
+ void pop_back(void)
+ {
+ assert(n>0);
+ --n;
+ }
+
+ void push_back(const T& value)
+ {
+ if(n==max_n) grow();
+ data[n++]=value;
+ }
+
+ reverse_iterator rbegin(void)
+ { return reverse_iterator(end()); }
+
+ const_reverse_iterator rbegin(void) const
+ { return const_reverse_iterator(end()); }
+
+ reverse_iterator rend(void)
+ { return reverse_iterator(begin()); }
+
+ const_reverse_iterator rend(void) const
+ { return const_reverse_iterator(begin()); }
+
+ void reserve(unsigned long r)
+ {
+ if(r>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+ T *new_data=(T*)std::realloc(data, r*sizeof(T));
+ if(!new_data) throw std::bad_alloc();
+ data=new_data;
+ max_n=r;
+ }
+
+ void resize(unsigned long n_)
+ {
+ if(n_>max_n) reserve(n_);
+ n=n_;
+ }
+
+ void resize(unsigned long n_, const T& value)
+ {
+ if(n_>max_n) reserve(n_);
+ if(n<n_) for(unsigned long i=n; i<n_; ++i) data[i]=value;
+ n=n_;
+ }
+
+ void set_zero(void)
+ { std::memset(data, 0, n*sizeof(T)); }
+
+ unsigned long size(void) const
+ { return n; }
+
+ void swap(Array1<T>& x)
+ {
+ std::swap(n, x.n);
+ std::swap(max_n, x.max_n);
+ std::swap(data, x.data);
+ }
+
+ // resize the array to avoid wasted space, without changing contents
+ // (Note: realloc, at least on some platforms, will not do the trick)
+ void trim(void)
+ {
+ if(n==max_n) return;
+ T *new_data=(T*)std::malloc(n*sizeof(T));
+ if(!new_data) return;
+ std::memcpy(new_data, data, n*sizeof(T));
+ std::free(data);
+ data=new_data;
+ max_n=n;
+ }
+};
+
+// some common arrays
+
+typedef Array1<double> Array1d;
+typedef Array1<float> Array1f;
+typedef Array1<long long> Array1ll;
+typedef Array1<unsigned long long> Array1ull;
+typedef Array1<int> Array1i;
+typedef Array1<unsigned int> Array1ui;
+typedef Array1<short> Array1s;
+typedef Array1<unsigned short> Array1us;
+typedef Array1<char> Array1c;
+typedef Array1<unsigned char> Array1uc;
+
+//============================================================================
+template<typename T>
+struct WrapArray1
+{
+ // STL-friendly typedefs
+
+ typedef T* iterator;
+ typedef const T* const_iterator;
+ typedef unsigned long size_type;
+ typedef long difference_type;
+ typedef T& reference;
+ typedef const T& const_reference;
+ typedef T value_type;
+ typedef T* pointer;
+ typedef const T* const_pointer;
+ typedef std::reverse_iterator<iterator> reverse_iterator;
+ typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
+
+ // the actual representation
+
+ unsigned long n;
+ unsigned long max_n;
+ T* data;
+
+ // most of STL vector's interface, with a few changes
+
+ WrapArray1(void)
+ : n(0), max_n(0), data(0)
+ {}
+
+ WrapArray1(unsigned long n_, T* data_)
+ : n(n_), max_n(n_), data(data_)
+ { assert(data || max_n==0); }
+
+ WrapArray1(unsigned long n_, T* data_, unsigned long max_n_)
+ : n(n_), max_n(max_n_), data(data_)
+ {
+ assert(n<=max_n);
+ assert(data || max_n==0);
+ }
+
+ // Allow for simple shallow copies of existing arrays
+ // Note that if the underlying arrays change where their data is, the WrapArray may be screwed up
+
+ WrapArray1(Array1<T>& a)
+ : n(a.n), max_n(a.max_n), data(a.data)
+ {}
+
+ WrapArray1(std::vector<T>& a)
+ : n(a.size()), max_n(a.capacity()), data(&a[0])
+ {}
+
+ void init(unsigned long n_, T* data_, unsigned long max_n_)
+ {
+ assert(n_<=max_n_);
+ assert(data_ || max_n_==0);
+ n=n_;
+ max_n=max_n_;
+ data=data_;
+ }
+
+ const T& operator[](unsigned long i) const
+ { return data[i]; }
+
+ T& operator[](unsigned long i)
+ { return data[i]; }
+
+ // these are range-checked (in debug mode) versions of operator[], like at()
+ const T& operator()(unsigned long i) const
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ T& operator()(unsigned long i)
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ bool operator==(const WrapArray1<T>& x) const
+ {
+ if(n!=x.n) return false;
+ for(unsigned long i=0; i<n; ++i) if(!(data[i]==x.data[i])) return false;
+ return true;
+ }
+
+ bool operator!=(const WrapArray1<T>& x) const
+ {
+ if(n!=x.n) return true;
+ for(unsigned long i=0; i<n; ++i) if(data[i]!=x.data[i]) return true;
+ return false;
+ }
+
+ bool operator<(const WrapArray1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]<x[i]) return true;
+ else if(x[i]<data[i]) return false;
+ }
+ return n<x.n;
+ }
+
+ bool operator>(const WrapArray1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]>x[i]) return true;
+ else if(x[i]>data[i]) return false;
+ }
+ return n>x.n;
+ }
+
+ bool operator<=(const WrapArray1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]<x[i]) return true;
+ else if(x[i]<data[i]) return false;
+ }
+ return n<=x.n;
+ }
+
+ bool operator>=(const WrapArray1<T>& x) const
+ {
+ for(unsigned long i=0; i<n && i<x.n; ++i){
+ if(data[i]>x[i]) return true;
+ else if(x[i]>data[i]) return false;
+ }
+ return n>=x.n;
+ }
+
+ void add_unique(const T& value)
+ {
+ for(unsigned long i=0; i<n; ++i) if(data[i]==value) return;
+ assert(n<max_n);
+ data[n++]=value;
+ }
+
+ void assign(const T& value)
+ { for(unsigned long i=0; i<n; ++i) data[i]=value; }
+
+ void assign(unsigned long num, const T& value)
+ { fill(num, value); }
+
+ // note: copydata may not alias this array's data, and this should not be
+ // used when T is a full object (which defines its own copying operation)
+ void assign(unsigned long num, const T* copydata)
+ {
+ assert(num==0 || copydata);
+ assert(num<=max_n);
+ n=num;
+ std::memcpy(data, copydata, n*sizeof(T));
+ }
+
+ template<typename InputIterator>
+ void assign(InputIterator first, InputIterator last)
+ { assign_(first, last, typename Array1IsIntegral<InputIterator>::type()); }
+
+ template<typename InputIterator>
+ void assign_(InputIterator first, InputIterator last, Array1True check)
+ { fill(first, last); }
+
+ template<typename InputIterator>
+ void assign_(InputIterator first, InputIterator last, Array1False check)
+ {
+ unsigned long i=0;
+ InputIterator p=first;
+ for(; p!=last; ++p, ++i){
+ assert(i<max_n);
+ data[i]=*p;
+ }
+ n=i;
+ }
+
+ const T& at(unsigned long i) const
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ T& at(unsigned long i)
+ {
+ assert(i<n);
+ return data[i];
+ }
+
+ const T& back(void) const
+ {
+ assert(data && n>0);
+ return data[n-1];
+ }
+
+ T& back(void)
+ {
+ assert(data && n>0);
+ return data[n-1];
+ }
+
+ const T* begin(void) const
+ { return data; }
+
+ T* begin(void)
+ { return data; }
+
+ unsigned long capacity(void) const
+ { return max_n; }
+
+ void clear(void)
+ { n=0; }
+
+ bool empty(void) const
+ { return n==0; }
+
+ const T* end(void) const
+ { return data+n; }
+
+ T* end(void)
+ { return data+n; }
+
+ void erase(unsigned long index)
+ {
+ assert(index<n);
+ for(unsigned long i=index; i<n-1; ++i)
+ data[i]=data[i-1];
+ pop_back();
+ }
+
+ void fill(unsigned long num, const T& value)
+ {
+ assert(num<=max_n);
+ n=num;
+ for(unsigned long i=0; i<n; ++i) data[i]=value;
+ }
+
+ const T& front(void) const
+ {
+ assert(n>0);
+ return *data;
+ }
+
+ T& front(void)
+ {
+ assert(n>0);
+ return *data;
+ }
+
+ void insert(unsigned long index, const T& entry)
+ {
+ assert(index<=n);
+ push_back(back());
+ for(unsigned long i=n-1; i>index; --i)
+ data[i]=data[i-1];
+ data[index]=entry;
+ }
+
+ unsigned long max_size(void) const
+ { return max_n; }
+
+ void pop_back(void)
+ {
+ assert(n>0);
+ --n;
+ }
+
+ void push_back(const T& value)
+ {
+ assert(n<max_n);
+ data[n++]=value;
+ }
+
+ reverse_iterator rbegin(void)
+ { return reverse_iterator(end()); }
+
+ const_reverse_iterator rbegin(void) const
+ { return const_reverse_iterator(end()); }
+
+ reverse_iterator rend(void)
+ { return reverse_iterator(begin()); }
+
+ const_reverse_iterator rend(void) const
+ { return const_reverse_iterator(begin()); }
+
+ void reserve(unsigned long r)
+ { assert(r<=max_n); }
+
+ void resize(unsigned long n_)
+ {
+ assert(n_<=max_n);
+ n=n_;
+ }
+
+ void resize(unsigned long n_, const T& value)
+ {
+ assert(n_<=max_n);
+ if(n<n_) for(unsigned long i=n; i<n_; ++i) data[i]=value;
+ n=n_;
+ }
+
+ // note: shouldn't be used when T is a full object (setting to zero may not make sense)
+ void set_zero(void)
+ { std::memset(data, 0, n*sizeof(T)); }
+
+ unsigned long size(void) const
+ { return n; }
+
+ void swap(WrapArray1<T>& x)
+ {
+ std::swap(n, x.n);
+ std::swap(max_n, x.max_n);
+ std::swap(data, x.data);
+ }
+};
+
+// some common arrays
+
+typedef WrapArray1<double> WrapArray1d;
+typedef WrapArray1<float> WrapArray1f;
+typedef WrapArray1<long long> WrapArray1ll;
+typedef WrapArray1<unsigned long long> WrapArray1ull;
+typedef WrapArray1<int> WrapArray1i;
+typedef WrapArray1<unsigned int> WrapArray1ui;
+typedef WrapArray1<short> WrapArray1s;
+typedef WrapArray1<unsigned short> WrapArray1us;
+typedef WrapArray1<char> WrapArray1c;
+typedef WrapArray1<unsigned char> WrapArray1uc;
+
+#endif
280 array2.h
@@ -0,0 +1,280 @@
+#ifndef ARRAY2_H
+#define ARRAY2_H
+
+#include "array1.h"
+#include <algorithm>
+#include <cassert>
+#include <vector>
+
+template<class T, class ArrayT=std::vector<T> >
+struct Array2
+{
+ // STL-friendly typedefs
+
+ typedef typename ArrayT::iterator iterator;
+ typedef typename ArrayT::const_iterator const_iterator;
+ typedef typename ArrayT::size_type size_type;
+ typedef long difference_type;
+ typedef T& reference;
+ typedef const T& const_reference;
+ typedef T value_type;
+ typedef T* pointer;
+ typedef const T* const_pointer;
+ typedef typename ArrayT::reverse_iterator reverse_iterator;
+ typedef typename ArrayT::const_reverse_iterator const_reverse_iterator;
+
+ // the actual representation
+
+ int ni, nj;
+ ArrayT a;
+
+ // the interface
+
+ Array2(void)
+ : ni(0), nj(0)
+ {}
+
+ Array2(int ni_, int nj_)
+ : ni(ni_), nj(nj_), a(ni_*nj_)
+ { assert(ni_>=0 && nj>=0); }
+
+ Array2(int ni_, int nj_, ArrayT& a_)
+ : ni(ni_), nj(nj_), a(a_)
+ { assert(ni_>=0 && nj>=0); }
+
+ Array2(int ni_, int nj_, const T& value)
+ : ni(ni_), nj(nj_), a(ni_*nj_, value)
+ { assert(ni_>=0 && nj>=0); }
+
+ Array2(int ni_, int nj_, const T& value, size_type max_n_)
+ : ni(ni_), nj(nj_), a(ni_*nj_, value, max_n_)
+ { assert(ni_>=0 && nj>=0); }
+
+ Array2(int ni_, int nj_, T* data_)
+ : ni(ni_), nj(nj_), a(ni_*nj_, data_)
+ { assert(ni_>=0 && nj>=0); }
+
+ Array2(int ni_, int nj_, T* data_, size_type max_n_)
+ : ni(ni_), nj(nj_), a(ni_*nj_, data_, max_n_)
+ { assert(ni_>=0 && nj>=0); }
+
+ template<class OtherArrayT>
+ Array2(Array2<T, OtherArrayT>& other)
+ : ni(other.ni), nj(other.nj), a(other.a)
+ {}
+
+ ~Array2(void)
+ {
+#ifndef NDEBUG
+ ni=nj=0;
+#endif
+ }
+
+ const T& operator()(int i, int j) const
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj);
+ return a[i+ni*j];
+ }
+
+ T& operator()(int i, int j)
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj);
+ return a[i+ni*j];
+ }
+
+ bool operator==(const Array2<T>& x) const
+ { return ni==x.ni && nj==x.nj && a==x.a; }
+
+ bool operator!=(const Array2<T>& x) const
+ { return ni!=x.ni || nj!=x.nj || a!=x.a; }
+
+ bool operator<(const Array2<T>& x) const
+ {
+ if(ni<x.ni) return true; else if(ni>x.ni) return false;
+ if(nj<x.nj) return true; else if(nj>x.nj) return false;
+ return a<x.a;
+ }
+
+ bool operator>(const Array2<T>& x) const
+ {
+ if(ni>x.ni) return true; else if(ni<x.ni) return false;
+ if(nj>x.nj) return true; else if(nj<x.nj) return false;
+ return a>x.a;
+ }
+
+ bool operator<=(const Array2<T>& x) const
+ {
+ if(ni<x.ni) return true; else if(ni>x.ni) return false;
+ if(nj<x.nj) return true; else if(nj>x.nj) return false;
+ return a<=x.a;
+ }
+
+ bool operator>=(const Array2<T>& x) const
+ {
+ if(ni>x.ni) return true; else if(ni<x.ni) return false;
+ if(nj>x.nj) return true; else if(nj<x.nj) return false;
+ return a>=x.a;
+ }
+
+ void assign(const T& value)
+ { a.assign(value); }
+
+ void assign(int ni_, int nj_, const T& value)
+ {
+ a.assign(ni_*nj_, value);
+ ni=ni_;
+ nj=nj_;
+ }
+
+ void assign(int ni_, int nj_, const T* copydata)
+ {
+ a.assign(ni_*nj_, copydata);
+ ni=ni_;
+ nj=nj_;
+ }
+
+ const T& at(int i, int j) const
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj);
+ return a[i+ni*j];
+ }
+
+ T& at(int i, int j)
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj);
+ return a[i+ni*j];
+ }
+
+ const T& back(void) const
+ {
+ assert(a.size());
+ return a.back();
+ }
+
+ T& back(void)
+ {
+ assert(a.size());
+ return a.back();
+ }
+
+ const_iterator begin(void) const
+ { return a.begin(); }
+
+ iterator begin(void)
+ { return a.begin(); }
+
+ size_type capacity(void) const
+ { return a.capacity(); }
+
+ void clear(void)
+ {
+ a.clear();
+ ni=nj=0;
+ }
+
+ bool empty(void) const
+ { return a.empty(); }
+
+ const_iterator end(void) const
+ { return a.end(); }
+
+ iterator end(void)
+ { return a.end(); }
+
+ void fill(int ni_, int nj_, const T& value)
+ {
+ a.fill(ni_*nj_, value);
+ ni=ni_;
+ nj=nj_;
+ }
+
+ const T& front(void) const
+ {
+ assert(a.size());
+ return a.front();
+ }
+
+ T& front(void)
+ {
+ assert(a.size());
+ return a.front();
+ }
+
+ size_type max_size(void) const
+ { return a.max_size(); }
+
+ reverse_iterator rbegin(void)
+ { return reverse_iterator(end()); }
+
+ const_reverse_iterator rbegin(void) const
+ { return const_reverse_iterator(end()); }
+
+ reverse_iterator rend(void)
+ { return reverse_iterator(begin()); }
+
+ const_reverse_iterator rend(void) const
+ { return const_reverse_iterator(begin()); }
+
+ void reserve(int reserve_ni, int reserve_nj)
+ { a.reserve(reserve_ni*reserve_nj); }
+
+ void resize(int ni_, int nj_)
+ {
+ assert(ni_>=0 && nj_>=0);
+ a.resize(ni_*nj_);
+ ni=ni_;
+ nj=nj_;
+ }
+
+ void resize(int ni_, int nj_, const T& value)
+ {
+ assert(ni_>=0 && nj_>=0);
+ a.resize(ni_*nj_, value);
+ ni=ni_;
+ nj=nj_;
+ }
+
+ void set_zero(void)
+ { a.set_zero(); }
+
+ size_type size(void) const
+ { return a.size(); }
+
+ void swap(Array2<T>& x)
+ {
+ std::swap(ni, x.ni);
+ std::swap(nj, x.nj);
+ a.swap(x.a);
+ }
+
+ void trim(void)
+ { a.trim(); }
+};
+
+// some common arrays
+
+typedef Array2<double, Array1<double> > Array2d;
+typedef Array2<float, Array1<float> > Array2f;
+typedef Array2<long long, Array1<long long> > Array2ll;
+typedef Array2<unsigned long long, Array1<unsigned long long> > Array2ull;
+typedef Array2<int, Array1<int> > Array2i;
+typedef Array2<unsigned int, Array1<unsigned int> > Array2ui;
+typedef Array2<short, Array1<short> > Array2s;
+typedef Array2<unsigned short, Array1<unsigned short> > Array2us;
+typedef Array2<char, Array1<char> > Array2c;
+typedef Array2<unsigned char, Array1<unsigned char> > Array2uc;
+
+// and wrapped versions
+
+typedef Array2<double, WrapArray1<double> > WrapArray2d;
+typedef Array2<float, WrapArray1<float> > WrapArray2f;
+typedef Array2<long long, WrapArray1<long long> > WrapArray2ll;
+typedef Array2<unsigned long long, WrapArray1<unsigned long long> > WrapArray2ull;
+typedef Array2<int, WrapArray1<int> > WrapArray2i;
+typedef Array2<unsigned int, WrapArray1<unsigned int> > WrapArray2ui;
+typedef Array2<short, WrapArray1<short> > WrapArray2s;
+typedef Array2<unsigned short, WrapArray1<unsigned short> > WrapArray2us;
+typedef Array2<char, WrapArray1<char> > WrapArray2c;
+typedef Array2<unsigned char, WrapArray1<unsigned char> > WrapArray2uc;
+
+#endif
63 array2_utils.h
@@ -0,0 +1,63 @@
+#ifndef ARRAY2_UTILS_H
+#define ARRAY2_UTILS_H
+
+#include "vec.h"
+#include "array2.h"
+#include "util.h"
+
+template<class S, class T>
+T interpolate_value(const Vec<2,S>& point, const Array2<T, Array1<T> >& grid) {
+ int i,j;
+ S fx,fy;
+
+ get_barycentric(point[0], i, fx, 0, grid.ni);
+ get_barycentric(point[1], j, fy, 0, grid.nj);
+
+ return bilerp(
+ grid(i,j), grid(i+1,j),
+ grid(i,j+1), grid(i+1,j+1),
+ fx, fy);
+}
+
+template<class S, class T>
+float interpolate_gradient(Vec<2,T>& gradient, const Vec<2,S>& point, const Array2<T, Array1<float> >& grid) {
+ int i,j;
+ S fx,fy;
+ get_barycentric(point[0], i, fx, 0, grid.ni);
+ get_barycentric(point[1], j, fy, 0, grid.nj);
+
+ T v00 = grid(i,j);
+ T v01 = grid(i,j+1);
+ T v10 = grid(i+1,j);
+ T v11 = grid(i+1,j+1);
+
+ T ddy0 = (v01 - v00);
+ T ddy1 = (v11 - v10);
+
+ T ddx0 = (v10 - v00);
+ T ddx1 = (v11 - v01);
+
+ gradient[0] = lerp(ddx0,ddx1,fy);
+ gradient[1] = lerp(ddy0,ddy1,fx);
+
+ //may as well return value too
+ return bilerp(v00, v10, v01, v11, fx, fy);
+}
+
+template<class T>
+void write_matlab_array(std::ostream &output, Array2<T, Array1<T> >&a, const char *variable_name, bool transpose=false)
+{
+ output<<variable_name<<"=[";
+ for(int j = 0; j < a.nj; ++j) {
+ for(int i = 0; i < a.ni; ++i) {
+ output<<a(i,j)<<" ";
+ }
+ output<<";";
+ }
+ output<<"]";
+ if(transpose)
+ output<<"'";
+ output<<";"<<std::endl;
+}
+
+#endif
272 array3.h
@@ -0,0 +1,272 @@
+#ifndef ARRAY3_H
+#define ARRAY3_H
+
+#include "array1.h"
+#include <algorithm>
+#include <cassert>
+#include <vector>
+
+template<class T, class ArrayT=std::vector<T> >
+struct Array3
+{
+ // STL-friendly typedefs
+
+ typedef typename ArrayT::iterator iterator;
+ typedef typename ArrayT::const_iterator const_iterator;
+ typedef typename ArrayT::size_type size_type;
+ typedef long difference_type;
+ typedef T& reference;
+ typedef const T& const_reference;
+ typedef T value_type;
+ typedef T* pointer;
+ typedef const T* const_pointer;
+ typedef typename ArrayT::reverse_iterator reverse_iterator;
+ typedef typename ArrayT::const_reverse_iterator const_reverse_iterator;
+
+ // the actual representation
+
+ int ni, nj, nk;
+ ArrayT a;
+
+ // the interface
+
+ Array3(void)
+ : ni(0), nj(0), nk(0)
+ {}
+
+ Array3(int ni_, int nj_, int nk_)
+ : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ Array3(int ni_, int nj_, int nk_, ArrayT& a_)
+ : ni(ni_), nj(nj_), nk(nk_), a(a_)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ Array3(int ni_, int nj_, int nk_, const T& value)
+ : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ Array3(int ni_, int nj_, int nk_, const T& value, size_type max_n_)
+ : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value, max_n_)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ Array3(int ni_, int nj_, int nk_, T* data_)
+ : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ Array3(int ni_, int nj_, int nk_, T* data_, size_type max_n_)
+ : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_, max_n_)
+ { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+ ~Array3(void)
+ {
+#ifndef NDEBUG
+ ni=nj=0;
+#endif
+ }
+
+ const T& operator()(int i, int j, int k) const
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+ return a[i+ni*(j+nj*k)];
+ }
+
+ T& operator()(int i, int j, int k)
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+ return a[i+ni*(j+nj*k)];
+ }
+
+ bool operator==(const Array3<T>& x) const
+ { return ni==x.ni && nj==x.nj && nk==x.nk && a==x.a; }
+
+ bool operator!=(const Array3<T>& x) const
+ { return ni!=x.ni || nj!=x.nj || nk!=x.nk || a!=x.a; }
+
+ bool operator<(const Array3<T>& x) const
+ {
+ if(ni<x.ni) return true; else if(ni>x.ni) return false;
+ if(nj<x.nj) return true; else if(nj>x.nj) return false;
+ if(nk<x.nk) return true; else if(nk>x.nk) return false;
+ return a<x.a;
+ }
+
+ bool operator>(const Array3<T>& x) const
+ {
+ if(ni>x.ni) return true; else if(ni<x.ni) return false;
+ if(nj>x.nj) return true; else if(nj<x.nj) return false;
+ if(nk>x.nk) return true; else if(nk<x.nk) return false;
+ return a>x.a;
+ }
+
+ bool operator<=(const Array3<T>& x) const
+ {
+ if(ni<x.ni) return true; else if(ni>x.ni) return false;
+ if(nj<x.nj) return true; else if(nj>x.nj) return false;
+ if(nk<x.nk) return true; else if(nk>x.nk) return false;
+ return a<=x.a;
+ }
+
+ bool operator>=(const Array3<T>& x) const
+ {
+ if(ni>x.ni) return true; else if(ni<x.ni) return false;
+ if(nj>x.nj) return true; else if(nj<x.nj) return false;
+ if(nk>x.nk) return true; else if(nk<x.nk) return false;
+ return a>=x.a;
+ }
+
+ void assign(const T& value)
+ { a.assign(value); }
+
+ void assign(int ni_, int nj_, int nk_, const T& value)
+ {
+ a.assign(ni_*nj_*nk_, value);
+ ni=ni_;
+ nj=nj_;
+ nk=nk_;
+ }
+
+ void assign(int ni_, int nj_, int nk_, const T* copydata)
+ {
+ a.assign(ni_*nj_*nk_, copydata);
+ ni=ni_;
+ nj=nj_;
+ nk=nk_;
+ }
+
+ const T& at(int i, int j, int k) const
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+ return a[i+ni*(j+nj*k)];
+ }
+
+ T& at(int i, int j, int k)
+ {
+ assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+ return a[i+ni*(j+nj*k)];
+ }
+
+ const T& back(void) const
+ {
+ assert(a.size());
+ return a.back();
+ }
+
+ T& back(void)
+ {
+ assert(a.size());
+ return a.back();
+ }
+
+ const_iterator begin(void) const
+ { return a.begin(); }
+
+ iterator begin(void)
+ { return a.begin(); }
+
+ size_type capacity(void) const
+ { return a.capacity(); }
+
+ void clear(void)
+ {
+ a.clear();
+ ni=nj=nk=0;
+ }
+
+ bool empty(void) const
+ { return a.empty(); }
+
+ const_iterator end(void) const
+ { return a.end(); }
+
+ iterator end(void)
+ { return a.end(); }
+
+ void fill(int ni_, int nj_, int nk_, const T& value)
+ {
+ a.fill(ni_*nj_*nk_, value);
+ ni=ni_;
+ nj=nj_;
+ nk=nk_;
+ }
+
+ const T& front(void) const
+ {
+ assert(a.size());
+ return a.front();
+ }
+
+ T& front(void)
+ {
+ assert(a.size());
+ return a.front();
+ }
+
+ size_type max_size(void) const
+ { return a.max_size(); }
+
+ reverse_iterator rbegin(void)
+ { return reverse_iterator(end()); }
+
+ const_reverse_iterator rbegin(void) const
+ { return const_reverse_iterator(end()); }
+
+ reverse_iterator rend(void)
+ { return reverse_iterator(begin()); }
+
+ const_reverse_iterator rend(void) const
+ { return const_reverse_iterator(begin()); }
+
+ void reserve(int reserve_ni, int reserve_nj, int reserve_nk)
+ { a.reserve(reserve_ni*reserve_nj*reserve_nk); }
+
+ void resize(int ni_, int nj_, int nk_)
+ {
+ assert(ni_>=0 && nj_>=0 && nk_>=0);
+ a.resize(ni_*nj_*nk_);
+ ni=ni_;
+ nj=nj_;
+ nk=nk_;
+ }
+
+ void resize(int ni_, int nj_, int nk_, const T& value)
+ {
+ assert(ni_>=0 && nj_>=0 && nk_>=0);
+ a.resize(ni_*nj_*nk_, value);
+ ni=ni_;
+ nj=nj_;
+ nk=nk_;
+ }
+
+ void set_zero(void)
+ { a.set_zero(); }
+
+ size_type size(void) const
+ { return a.size(); }
+
+ void swap(Array3<T>& x)
+ {
+ std::swap(ni, x.ni);
+ std::swap(nj, x.nj);
+ std::swap(nk, x.nk);
+ a.swap(x.a);
+ }
+
+ void trim(void)
+ { a.trim(); }
+};
+
+// some common arrays
+
+typedef Array3<double, Array1<double> > Array3d;
+typedef Array3<float, Array1<float> > Array3f;
+typedef Array3<long long, Array1<long long> > Array3ll;
+typedef Array3<unsigned long long, Array1<unsigned long long> > Array3ull;
+typedef Array3<int, Array1<int> > Array3i;
+typedef Array3<unsigned int, Array1<unsigned int> > Array3ui;
+typedef Array3<short, Array1<short> > Array3s;
+typedef Array3<unsigned short, Array1<unsigned short> > Array3us;
+typedef Array3<char, Array1<char> > Array3c;
+typedef Array3<unsigned char, Array1<unsigned char> > Array3uc;
+
+#endif
72 array3_utils.h
@@ -0,0 +1,72 @@
+#ifndef ARRAY3_UTILS_H
+#define ARRAY3_UTILS_H
+
+#include "vec.h"
+#include "array3.h"
+#include "util.h"
+
+template<class S, class T>
+T interpolate_value(const Vec<3,S>& point, const Array3<T, Array1<T> >& grid) {
+ int i,j,k;
+ S fi,fj,fk;
+
+ get_barycentric(point[0], i, fi, 0, grid.ni);
+ get_barycentric(point[1], j, fj, 0, grid.nj);
+ get_barycentric(point[2], k, fk, 0, grid.nk);
+
+ return trilerp(
+ grid(i,j,k), grid(i+1,j,k), grid(i,j+1,k), grid(i+1,j+1,k),
+ grid(i,j,k+1), grid(i+1,j,k+1), grid(i,j+1,k+1), grid(i+1,j+1,k+1),
+ fi,fj,fk);
+}
+
+template<class S,class T>
+T interpolate_gradient(Vec<3,T>& gradient, const Vec<3,S>& point, const Array3<T, Array1<T> >& grid) {
+ int i,j,k;
+ S fx,fy,fz;
+
+ get_barycentric(point[0], i, fx, 0, grid.ni);
+ get_barycentric(point[1], j, fy, 0, grid.nj);
+ get_barycentric(point[2], k, fz, 0, grid.nk);
+
+ T v000 = grid(i,j,k);
+ T v001 = grid(i,j,k+1);
+ T v010 = grid(i,j+1,k);
+ T v011 = grid(i,j+1,k+1);
+ T v100 = grid(i+1,j,k);
+ T v101 = grid(i+1,j,k+1);
+ T v110 = grid(i+1,j+1,k);
+ T v111 = grid(i+1,j+1,k+1);
+
+ T ddx00 = (v100 - v000);
+ T ddx10 = (v110 - v010);
+ T ddx01 = (v101 - v001);
+ T ddx11 = (v111 - v011);
+ T dv_dx = bilerp(ddx00,ddx10,ddx01,ddx11, fy,fz);
+
+ T ddy00 = (v010 - v000);
+ T ddy10 = (v110 - v100);
+ T ddy01 = (v011 - v001);
+ T ddy11 = (v111 - v101);
+ T dv_dy = bilerp(ddy00,ddy10,ddy01,ddy11, fx,fz);
+
+ T ddz00 = (v001 - v000);
+ T ddz10 = (v101 - v100);
+ T ddz01 = (v011 - v010);
+ T ddz11 = (v111 - v110);
+ T dv_dz = bilerp(ddz00,ddz10,ddz01,ddz11, fx,fy);
+
+ gradient[0] = dv_dx;
+ gradient[1] = dv_dy;
+ gradient[2] = dv_dz;
+
+ //return value for good measure.
+ return trilerp(
+ v000, v100,
+ v010, v110,
+ v001, v101,
+ v011, v111,
+ fx, fy, fz);
+}
+
+#endif
610 fluidsim.cpp
@@ -0,0 +1,610 @@
+#include "fluidsim.h"
+
+#include "array3_utils.h"
+#include "levelset_util.h"
+#include "pcgsolver/sparse_matrix.h"
+#include "pcgsolver/pcg_solver.h"
+
+void extrapolate(Array3f& grid, Array3c& valid);
+
+void FluidSim::initialize(float width, int ni_, int nj_, int nk_) {
+ ni = ni_;
+ nj = nj_;
+ nk = nk_;
+ dx = width / (float)ni;
+ u.resize(ni+1,nj,nk); temp_u.resize(ni+1,nj,nk); u_weights.resize(ni+1,nj,nk); u_valid.resize(ni+1,nj,nk);
+ v.resize(ni,nj+1,nk); temp_v.resize(ni,nj+1,nk); v_weights.resize(ni,nj+1,nk); v_valid.resize(ni,nj+1,nk);
+ w.resize(ni,nj,nk+1); temp_w.resize(ni,nj,nk+1); w_weights.resize(ni,nj,nk+1); w_valid.resize(ni,nj,nk+1);
+
+ particle_radius = (float)(dx * 1.01*sqrt(3.0)/2.0);
+ //make the particles large enough so they always appear on the grid
+
+ u.set_zero();
+ v.set_zero();
+ w.set_zero();
+ nodal_solid_phi.resize(ni+1,nj+1,nk+1);
+ valid.resize(ni+1, nj+1, nk+1);
+ old_valid.resize(ni+1, nj+1, nk+1);
+ liquid_phi.resize(ni,nj,nk);
+
+}
+
+//Initialize the grid-based signed distance field that dictates the position of the solid boundary
+void FluidSim::set_boundary(float (*phi)(const Vec3f&)) {
+
+ for(int k = 0; k < nk+1; ++k) for(int j = 0; j < nj+1; ++j) for(int i = 0; i < ni+1; ++i) {
+ Vec3f pos(i*dx,j*dx,k*dx);
+ nodal_solid_phi(i,j,k) = phi(pos);
+ }
+
+}
+
+void FluidSim::set_liquid(float (*phi)(const Vec3f&)) {
+ //surface.reset_phi(phi, dx, Vec3f(0.5f*dx,0.5f*dx,0.5f*dx), ni, nj, nk);
+
+ //initialize particles
+ int seed = 0;
+ for(int k = 0; k < nk; ++k) for(int j = 0; j < nj; ++j) for(int i = 0; i < ni; ++i) {
+ Vec3f pos(i*dx,j*dx,k*dx);
+ float a = randhashf(seed++); float b = randhashf(seed++); float c = randhashf(seed++);
+ pos += dx * Vec3f(a,b,c);
+
+ if(phi(pos) <= -particle_radius) {
+ float solid_phi = interpolate_value(pos/dx, nodal_solid_phi);
+ if(solid_phi >= 0)
+ particles.push_back(pos);
+ }
+ }
+}
+
+//The main fluid simulation step
+void FluidSim::advance(float dt) {
+ float t = 0;
+
+ while(t < dt) {
+ float substep = cfl();
+ if(t + substep > dt)
+ substep = dt - t;
+ printf("Taking substep of size %f (to %0.3f%% of the frame)\n", substep, 100 * (t+substep)/dt);
+
+ printf(" Surface (particle) advection\n");
+ advect_particles(substep);
+
+ printf(" Velocity advection\n");
+ //Advance the velocity
+ advect(substep);
+ add_force(substep);
+
+ printf(" Pressure projection\n");
+ project(substep);
+
+ //Pressure projection only produces valid velocities in faces with non-zero associated face area.
+ //Because the advection step may interpolate from these invalid faces,
+ //we must extrapolate velocities from the fluid domain into these invalid faces.
+ printf(" Extrapolation\n");
+ extrapolate(u, u_valid);
+ extrapolate(v, v_valid);
+ extrapolate(w, w_valid);
+
+ //For extrapolated velocities, replace the normal component with
+ //that of the object.
+ printf(" Constrain boundary velocities\n");
+ constrain_velocity();
+
+ t+=substep;
+ }
+}
+
+
+float FluidSim::cfl() {
+
+ float maxvel = 0;
+ for(unsigned int i = 0; i < u.a.size(); ++i)
+ maxvel = max(maxvel, fabs(u.a[i]));
+ for(unsigned int i = 0; i < v.a.size(); ++i)
+ maxvel = max(maxvel, fabs(v.a[i]));
+ for(unsigned int i = 0; i < w.a.size(); ++i)
+ maxvel = max(maxvel, fabs(w.a[i]));
+
+ return dx / maxvel;
+}
+
+void FluidSim::add_particle(const Vec3f& pos) {
+ particles.push_back(pos);
+}
+
+void FluidSim::add_force(float dt) {
+
+ //gravity
+ for(int k = 0;k < nk; ++k) for(int j = 0; j < nj+1; ++j) for(int i = 0; i < ni; ++i) {
+ v(i,j,k) -= 9.81f * dt;
+ }
+
+}
+
+
+
+//For extrapolated points, replace the normal component
+//of velocity with the object velocity (in this case zero).
+void FluidSim::constrain_velocity() {
+ temp_u = u;
+ temp_v = v;
+ temp_w = w;
+
+ //(At lower grid resolutions, the normal estimate from the signed
+ //distance function can be poor, so it doesn't work quite as well.
+ //An exact normal would do better if we had it for the geometry.)
+
+ //constrain u
+ for(int k = 0; k < u.nk;++k) for(int j = 0; j < u.nj; ++j) for(int i = 0; i < u.ni; ++i) {
+ if(u_weights(i,j,k) == 0) {
+ //apply constraint
+ Vec3f pos(i*dx, (j+0.5f)*dx, (k+0.5f)*dx);
+ Vec3f vel = get_velocity(pos);
+ Vec3f normal(0,0,0);
+ interpolate_gradient(normal, pos/dx, nodal_solid_phi);
+ normalize(normal);
+ float perp_component = dot(vel, normal);
+ vel -= perp_component*normal;
+ temp_u(i,j,k) = vel[0];
+ }
+ }
+
+ //constrain v
+ for(int k = 0; k < v.nk;++k) for(int j = 0; j < v.nj; ++j) for(int i = 0; i < v.ni; ++i) {
+ if(v_weights(i,j,k) == 0) {
+ //apply constraint
+ Vec3f pos((i+0.5f)*dx, j*dx, (k+0.5f)*dx);
+ Vec3f vel = get_velocity(pos);
+ Vec3f normal(0,0,0);
+ interpolate_gradient(normal, pos/dx, nodal_solid_phi);
+ normalize(normal);
+ float perp_component = dot(vel, normal);
+ vel -= perp_component*normal;
+ temp_v(i,j,k) = vel[1];
+ }
+ }
+
+ //constrain w
+ for(int k = 0; k < w.nk;++k) for(int j = 0; j < w.nj; ++j) for(int i = 0; i < w.ni; ++i) {
+ if(w_weights(i,j,k) == 0) {
+ //apply constraint
+ Vec3f pos((i+0.5f)*dx, (j+0.5f)*dx, k*dx);
+ Vec3f vel = get_velocity(pos);
+ Vec3f normal(0,0,0);
+ interpolate_gradient(normal, pos/dx, nodal_solid_phi);
+ normalize(normal);
+ float perp_component = dot(vel, normal);
+ vel -= perp_component*normal;
+ temp_w(i,j,k) = vel[2];
+ }
+ }
+
+ //update
+ u = temp_u;
+ v = temp_v;
+ w = temp_w;
+
+}
+
+void FluidSim::advect_particles(float dt) {
+ for(unsigned int p = 0; p < particles.size(); ++p) {
+ particles[p] = trace_rk2(particles[p], dt);
+
+ //check boundaries and project exterior particles back in
+ float phi_val = interpolate_value(particles[p]/dx, nodal_solid_phi);
+ if(phi_val < 0) {
+ Vec3f grad;
+ interpolate_gradient(grad, particles[p]/dx, nodal_solid_phi);
+ if(mag(grad) > 0)
+ normalize(grad);
+ particles[p] -= phi_val * grad;
+ }
+ }
+
+
+}
+
+//Basic first order semi-Lagrangian advection of velocities
+void FluidSim::advect(float dt) {
+
+ temp_u.assign(0);
+ temp_v.assign(0);
+ temp_w.assign(0);
+
+ //semi-Lagrangian advection on u-component of velocity
+ for(int k = 0; k < nk; ++k) for(int j = 0; j < nj; ++j) for(int i = 0; i < ni+1; ++i) {
+ Vec3f pos(i*dx, (j+0.5f)*dx, (k+0.5f)*dx);
+ pos = trace_rk2(pos, -dt);
+ temp_u(i,j,k) = get_velocity(pos)[0];
+ }
+
+ //semi-Lagrangian advection on v-component of velocity
+ for(int k = 0; k < nk; ++k) for(int j = 0; j < nj+1; ++j) for(int i = 0; i < ni; ++i) {
+ Vec3f pos((i+0.5f)*dx, j*dx, (k+0.5f)*dx);
+ pos = trace_rk2(pos, -dt);
+ temp_v(i,j,k) = get_velocity(pos)[1];
+ }
+
+ //semi-Lagrangian advection on w-component of velocity
+ for(int k = 0; k < nk+1; ++k) for(int j = 0; j < nj; ++j) for(int i = 0; i < ni; ++i) {
+ Vec3f pos((i+0.5f)*dx, (j+0.5f)*dx, k*dx);
+ pos = trace_rk2(pos, -dt);
+ temp_w(i,j,k) = get_velocity(pos)[2];
+ }
+
+ //move update velocities into u/v vectors
+ u = temp_u;
+ v = temp_v;
+ w = temp_w;
+}
+
+void FluidSim::compute_phi() {
+
+ //grab from particles
+ liquid_phi.assign(3*dx);
+ for(unsigned int p = 0; p < particles.size(); ++p) {
+ Vec3i cell_ind(particles[p] / dx);
+ for(int k = max(0,cell_ind[2] - 1); k <= min(cell_ind[2]+1,nk-1); ++k) {
+ for(int j = max(0,cell_ind[1] - 1); j <= min(cell_ind[1]+1,nj-1); ++j) {
+ for(int i = max(0,cell_ind[0] - 1); i <= min(cell_ind[0]+1,ni-1); ++i) {
+ Vec3f sample_pos((i+0.5f)*dx, (j+0.5f)*dx,(k+0.5f)*dx);
+ float test_val = dist(sample_pos, particles[p]) - particle_radius;
+ if(test_val < liquid_phi(i,j,k))
+ liquid_phi(i,j,k) = test_val;
+ }
+ }
+ }
+ }
+
+ //extend phi slightly into solids (this is a simple, naive approach, but works reasonably well)
+ Array3f phi_temp = liquid_phi;
+ for(int k = 0; k < nk; ++k) {
+ for(int j = 0; j < nj; ++j) {
+ for(int i = 0; i < ni; ++i) {
+ if(liquid_phi(i,j,k) < 0.5*dx) {
+ float solid_phi_val = 0.125f*(nodal_solid_phi(i,j,k) + nodal_solid_phi(i+1,j,k) + nodal_solid_phi(i,j+1,k) + nodal_solid_phi(i+1,j+1,k)
+ + nodal_solid_phi(i,j,k+1) + nodal_solid_phi(i+1,j,k+1) + nodal_solid_phi(i,j+1,k+1) + nodal_solid_phi(i+1,j+1,k+1));
+ if(solid_phi_val < 0)
+ phi_temp(i,j,k) = -0.5f*dx;
+ }
+ }
+ }
+ }
+ liquid_phi = phi_temp;
+
+
+}
+
+
+
+void FluidSim::project(float dt) {
+
+ //Estimate the liquid signed distance
+ compute_phi();
+
+ //Compute finite-volume type face area weight for each velocity sample.
+ compute_weights();
+
+ //Set up and solve the variational pressure solve.
+ solve_pressure(dt);
+
+}
+
+
+//Apply RK2 to advect a point in the domain.
+Vec3f FluidSim::trace_rk2(const Vec3f& position, float dt) {
+ Vec3f input = position;
+ Vec3f velocity = get_velocity(input);
+ velocity = get_velocity(input + 0.5f*dt*velocity);
+ input += dt*velocity;
+ return input;
+}
+
+//Interpolate velocity from the MAC grid.
+Vec3f FluidSim::get_velocity(const Vec3f& position) {
+
+ //Interpolate the velocity from the u and v grids
+ float u_value = interpolate_value(position / dx - Vec3f(0, 0.5f, 0.5f), u);
+ float v_value = interpolate_value(position / dx - Vec3f(0.5f, 0, 0.5f), v);
+ float w_value = interpolate_value(position / dx - Vec3f(0.5f, 0.5f, 0), w);
+
+ return Vec3f(u_value, v_value, w_value);
+}
+
+
+
+//Compute finite-volume style face-weights for fluid from nodal signed distances
+void FluidSim::compute_weights() {
+
+ //Compute face area fractions (using marching squares cases).
+ for(int k = 0; k < nk; ++k) for(int j = 0; j < nj; ++j) for(int i = 0; i < ni+1; ++i) {
+ u_weights(i,j,k) = 1 - fraction_inside(nodal_solid_phi(i,j, k),
+ nodal_solid_phi(i,j+1,k),
+ nodal_solid_phi(i,j, k+1),
+ nodal_solid_phi(i,j+1,k+1));
+ u_weights(i,j,k) = clamp(u_weights(i,j,k),0.0f,1.0f);
+ }
+ for(int k = 0; k < nk; ++k) for(int j = 0; j < nj+1; ++j) for(int i = 0; i < ni; ++i) {
+ v_weights(i,j,k) = 1 - fraction_inside(nodal_solid_phi(i, j,k),
+ nodal_solid_phi(i, j,k+1),
+ nodal_solid_phi(i+1,j,k),
+ nodal_solid_phi(i+1,j,k+1));
+ v_weights(i,j,k) = clamp(v_weights(i,j,k),0.0f,1.0f);
+ }
+ for(int k = 0; k < nk+1; ++k) for(int j = 0; j < nj; ++j) for(int i = 0; i < ni; ++i) {
+ w_weights(i,j,k) = 1 - fraction_inside(nodal_solid_phi(i, j, k),
+ nodal_solid_phi(i, j+1,k),
+ nodal_solid_phi(i+1,j, k),
+ nodal_solid_phi(i+1,j+1,k));
+ w_weights(i,j,k) = clamp(w_weights(i,j,k),0.0f,1.0f);
+ }
+
+
+}
+
+//An implementation of the variational pressure projection solve for static geometry
+void FluidSim::solve_pressure(float dt) {
+
+
+ int ni = v.ni;
+ int nj = u.nj;
+ int nk = u.nk;
+
+ int system_size = ni*nj*nk;
+ if(rhs.size() != system_size) {
+ rhs.resize(system_size);
+ pressure.resize(system_size);
+ matrix.resize(system_size);
+ }
+
+ matrix.zero();
+ rhs.assign(rhs.size(), 0);
+ pressure.assign(pressure.size(), 0);
+
+ //Build the linear system for pressure
+ for(int k = 1; k < nk-1; ++k) {
+ for(int j = 1; j < nj-1; ++j) {
+ for(int i = 1; i < ni-1; ++i) {
+ int index = i + ni*j + ni*nj*k;
+
+ rhs[index] = 0;
+ pressure[index] = 0;
+ float centre_phi = liquid_phi(i,j,k);
+ if(centre_phi < 0) {
+
+ //right neighbour
+ float term = u_weights(i+1,j,k) * dt / sqr(dx);
+ float right_phi = liquid_phi(i+1,j,k);
+ if(right_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index + 1, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, right_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] -= u_weights(i+1,j,k)*u(i+1,j,k) / dx;
+
+ //left neighbour
+ term = u_weights(i,j,k) * dt / sqr(dx);
+ float left_phi = liquid_phi(i-1,j,k);
+ if(left_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index - 1, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, left_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] += u_weights(i,j,k)*u(i,j,k) / dx;
+
+ //top neighbour
+ term = v_weights(i,j+1,k) * dt / sqr(dx);
+ float top_phi = liquid_phi(i,j+1,k);
+ if(top_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index + ni, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, top_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] -= v_weights(i,j+1,k)*v(i,j+1,k) / dx;
+
+ //bottom neighbour
+ term = v_weights(i,j,k) * dt / sqr(dx);
+ float bot_phi = liquid_phi(i,j-1,k);
+ if(bot_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index - ni, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, bot_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] += v_weights(i,j,k)*v(i,j,k) / dx;
+
+
+ //far neighbour
+ term = w_weights(i,j,k+1) * dt / sqr(dx);
+ float far_phi = liquid_phi(i,j,k+1);
+ if(far_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index + ni*nj, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, far_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] -= w_weights(i,j,k+1)*w(i,j,k+1) / dx;
+
+ //near neighbour
+ term = w_weights(i,j,k) * dt / sqr(dx);
+ float near_phi = liquid_phi(i,j,k-1);
+ if(near_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index - ni*nj, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, near_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] += w_weights(i,j,k)*w(i,j,k) / dx;
+
+ /*
+ //far neighbour
+ term = w_weights(i,j,k+1) * dt / sqr(dx);
+ float far_phi = liquid_phi(i,j,k+1);
+ if(far_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index + ni*nj, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, far_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] -= w_weights(i,j,k+1)*w(i,j,k+1) / dx;
+
+ //near neighbour
+ term = w_weights(i,j,k) * dt / sqr(dx);
+ float near_phi = liquid_phi(i,j,k-1);
+ if(near_phi < 0) {
+ matrix.add_to_element(index, index, term);
+ matrix.add_to_element(index, index - ni*nj, -term);
+ }
+ else {
+ float theta = fraction_inside(centre_phi, near_phi);
+ if(theta < 0.01f) theta = 0.01f;
+ matrix.add_to_element(index, index, term/theta);
+ }
+ rhs[index] += w_weights(i,j,k)*w(i,j,k) / dx;
+ */
+
+ }
+ }
+ }
+ }
+
+ //Solve the system using Robert Bridson's incomplete Cholesky PCG solver
+
+ double tolerance;
+ int iterations;
+ solver.set_solver_parameters(1e-18, 1000);
+ bool success = solver.solve(matrix, rhs, pressure, tolerance, iterations);
+ //printf("Solver took %d iterations and had residual %e\n", iterations, tolerance);
+ if(!success) {
+ printf("WARNING: Pressure solve failed!************************************************\n");
+ }
+
+ //Apply the velocity update
+ u_valid.assign(0);
+ for(int k = 0; k < u.nk; ++k) for(int j = 0; j < u.nj; ++j) for(int i = 1; i < u.ni-1; ++i) {
+ int index = i + j*ni + k*ni*nj;
+ if(u_weights(i,j,k) > 0 && (liquid_phi(i,j,k) < 0 || liquid_phi(i-1,j,k) < 0)) {
+ float theta = 1;
+ if(liquid_phi(i,j,k) >= 0 || liquid_phi(i-1,j,k) >= 0)
+ theta = fraction_inside(liquid_phi(i-1,j,k), liquid_phi(i,j,k));
+ if(theta < 0.01f) theta = 0.01f;
+ u(i,j,k) -= dt * (float)(pressure[index] - pressure[index-1]) / dx / theta;
+ u_valid(i,j,k) = 1;
+ }
+ }
+
+ v_valid.assign(0);
+ for(int k = 0; k < v.nk; ++k) for(int j = 1; j < v.nj-1; ++j) for(int i = 0; i < v.ni; ++i) {
+ int index = i + j*ni + k*ni*nj;
+ if(v_weights(i,j,k) > 0 && (liquid_phi(i,j,k) < 0 || liquid_phi(i,j-1,k) < 0)) {
+ float theta = 1;
+ if(liquid_phi(i,j,k) >= 0 || liquid_phi(i,j-1,k) >= 0)
+ theta = fraction_inside(liquid_phi(i,j-1,k), liquid_phi(i,j,k));
+ if(theta < 0.01f) theta = 0.01f;
+ v(i,j,k) -= dt * (float)(pressure[index] - pressure[index-ni]) / dx / theta;
+ v_valid(i,j,k) = 1;
+ }
+ }
+
+ w_valid.assign(0);
+ for(int k = 0; k < w.nk; ++k) for(int j = 0; j < w.nj; ++j) for(int i = 1; i < w.ni-1; ++i) {
+ int index = i + j*ni + k*ni*nj;
+ if(w_weights(i,j,k) > 0 && (liquid_phi(i,j,k) < 0 || liquid_phi(i,j,k-1) < 0)) {
+ float theta = 1;
+ if(liquid_phi(i,j,k) >= 0 || liquid_phi(i,j,k-1) >= 0)
+ theta = fraction_inside(liquid_phi(i,j,k-1), liquid_phi(i,j,k));
+ if(theta < 0.01f) theta = 0.01f;
+ w(i,j,k) -= dt * (float)(pressure[index] - pressure[index-ni*nj]) / dx / theta;
+ w_valid(i,j,k) = 1;
+ }
+ }
+
+ for(unsigned int i = 0; i < u_valid.a.size(); ++i)
+ if(u_valid.a[i] == 0)
+ u.a[i] = 0;
+ for(unsigned int i = 0; i < v_valid.a.size(); ++i)
+ if(v_valid.a[i] == 0)
+ v.a[i] = 0;
+ for(unsigned int i = 0; i < w_valid.a.size(); ++i)
+ if(w_valid.a[i] == 0)
+ w.a[i] = 0;
+}
+
+
+//Apply several iterations of a very simple propagation of valid velocity data in all directions
+void extrapolate(Array3f& grid, Array3c& valid) {
+
+ Array3f temp_grid = grid;
+ Array3c old_valid(valid.ni,valid.nj,valid.nk);
+ for(int layers = 0; layers < 10; ++layers) {
+ old_valid = valid;
+ for(int k = 1; k < grid.nk-1; ++k) for(int j = 1; j < grid.nj-1; ++j) for(int i = 1; i < grid.ni-1; ++i) {
+ float sum = 0;
+ int count = 0;
+
+ if(!old_valid(i,j,k)) {
+
+ if(old_valid(i+1,j,k)) {
+ sum += grid(i+1,j,k);
+ ++count;
+ }
+ if(old_valid(i-1,j,k)) {
+ sum += grid(i-1,j,k);
+ ++count;
+ }
+ if(old_valid(i,j+1,k)) {
+ sum += grid(i,j+1,k);
+ ++count;
+ }
+ if(old_valid(i,j-1,k)) {
+ sum += grid(i,j-1,k);
+ ++count;
+ }
+ if(old_valid(i,j,k+1)) {
+ sum += grid(i,j,k+1);
+ ++count;
+ }
+ if(old_valid(i,j,k-1)) {
+ sum += grid(i,j,k-1);
+ ++count;
+ }
+
+ //If any of neighbour cells were valid,
+ //assign the cell their average value and tag it as valid
+ if(count > 0) {
+ temp_grid(i,j,k) = sum /(float)count;
+ valid(i,j,k) = 1;
+ }
+
+ }
+ }
+ grid = temp_grid;
+
+ }
+
+}
72 fluidsim.h
@@ -0,0 +1,72 @@
+#ifndef FLUIDSIM_H
+#define FLUIDSIM_H
+
+#include "array3.h"
+#include "vec.h"
+#include "pcgsolver/sparse_matrix.h"
+#include "pcgsolver/pcg_solver.h"
+
+#include <vector>
+
+class FluidSim {
+
+public:
+ void initialize(float width, int ni_, int nj_, int nk_);
+ void set_boundary(float (*phi)(const Vec3f&));
+ void set_liquid(float (*phi)(const Vec3f&));
+ void add_particle(const Vec3f& pos);
+
+ void advance(float dt);
+
+ //Grid dimensions
+ int ni,nj,nk;
+ float dx;
+
+ //Fluid velocity
+ Array3f u, v, w;
+ Array3f temp_u, temp_v, temp_w;
+
+ //Static geometry representation
+ Array3f nodal_solid_phi;
+ Array3f u_weights, v_weights, w_weights;
+ Array3c u_valid, v_valid, w_valid;
+
+ std::vector<Vec3f> particles;
+ float particle_radius;
+
+ Array3f liquid_phi;
+
+ //Data arrays for extrapolation
+ Array3c valid, old_valid;
+
+ //Solver data
+ PCGSolver<double> solver;
+ SparseMatrixd matrix;
+ std::vector<double> rhs;
+ std::vector<double> pressure;
+
+ Vec3f get_velocity(const Vec3f& position);
+
+private:
+
+ Vec3f trace_rk2(const Vec3f& position, float dt);
+
+ float cfl();
+
+ void advect_particles(float dt);
+ void advect(float dt);
+ void add_force(float dt);
+ void project(float dt);
+ void constrain_velocity();
+
+ //helpers for pressure projection
+ void compute_weights();
+ void solve_pressure(float dt);
+ void compute_phi();
+
+
+
+};
+
+
+#endif
103 levelset_util.cpp
@@ -0,0 +1,103 @@
+#include "levelset_util.h"
+
+//Given two signed distance values (line endpoints), determine what fraction of a connecting segment is "inside"
+float fraction_inside(float phi_left, float phi_right) {
+ if(phi_left < 0 && phi_right < 0)
+ return 1;
+ if (phi_left < 0 && phi_right >= 0)
+ return phi_left / (phi_left - phi_right);
+ if(phi_left >= 0 && phi_right < 0)
+ return phi_right / (phi_right - phi_left);
+ else
+ return 0;
+}
+
+static void cycle_array(float* arr, int size) {
+ float t = arr[0];
+ for(int i = 0; i < size-1; ++i)
+ arr[i] = arr[i+1];
+ arr[size-1] = t;
+}
+
+//Given four signed distance values (square corners), determine what fraction of the square is "inside"
+float fraction_inside(float phi_bl, float phi_br, float phi_tl, float phi_tr) {
+
+ int inside_count = (phi_bl<0?1:0) + (phi_tl<0?1:0) + (phi_br<0?1:0) + (phi_tr<0?1