From 5908d5171c270c1fc99faa46cc7a115a9c740496 Mon Sep 17 00:00:00 2001 From: Christopher Batty Date: Fri, 29 Mar 2013 14:52:59 -0400 Subject: [PATCH] Create initial project. --- Fluid3D.sln | 26 + Fluid3D.vcproj | 268 ++++++++ README | 6 + array1.h | 793 ++++++++++++++++++++++++ array2.h | 280 +++++++++ array2_utils.h | 63 ++ array3.h | 272 +++++++++ array3_utils.h | 72 +++ fluidsim.cpp | 676 +++++++++++++++++++++ fluidsim.h | 83 +++ levelset_util.cpp | 103 ++++ levelset_util.h | 7 + main.cpp | 93 +++ pcgsolver/blas_wrapper.h | 54 ++ pcgsolver/pcg_solver.h | 312 ++++++++++ pcgsolver/sparse_matrix.h | 290 +++++++++ util.h | 472 ++++++++++++++ vec.h | 472 ++++++++++++++ viewpls3D/ViewFluid3D.vcproj | 211 +++++++ viewpls3D/gluvi.cpp | 1114 ++++++++++++++++++++++++++++++++++ viewpls3D/gluvi.h | 200 ++++++ viewpls3D/main.cpp | 370 +++++++++++ viewpls3D/util.h | 470 ++++++++++++++ viewpls3D/vec.h | 472 ++++++++++++++ viscosity3d.cpp | 487 +++++++++++++++ viscosity3d.h | 12 + volume_fractions.cpp | 139 +++++ volume_fractions.h | 30 + 28 files changed, 7847 insertions(+) create mode 100644 Fluid3D.sln create mode 100644 Fluid3D.vcproj create mode 100644 README create mode 100644 array1.h create mode 100644 array2.h create mode 100644 array2_utils.h create mode 100644 array3.h create mode 100644 array3_utils.h create mode 100644 fluidsim.cpp create mode 100644 fluidsim.h create mode 100644 levelset_util.cpp create mode 100644 levelset_util.h create mode 100644 main.cpp create mode 100644 pcgsolver/blas_wrapper.h create mode 100644 pcgsolver/pcg_solver.h create mode 100644 pcgsolver/sparse_matrix.h create mode 100644 util.h create mode 100644 vec.h create mode 100644 viewpls3D/ViewFluid3D.vcproj create mode 100644 viewpls3D/gluvi.cpp create mode 100644 viewpls3D/gluvi.h create mode 100644 viewpls3D/main.cpp create mode 100644 viewpls3D/util.h create mode 100644 viewpls3D/vec.h create mode 100644 viscosity3d.cpp create mode 100644 viscosity3d.h create mode 100644 volume_fractions.cpp create mode 100644 volume_fractions.h diff --git a/Fluid3D.sln b/Fluid3D.sln new file mode 100644 index 0000000..8980210 --- /dev/null +++ b/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 diff --git a/Fluid3D.vcproj b/Fluid3D.vcproj new file mode 100644 index 0000000..fba477f --- /dev/null +++ b/Fluid3D.vcproj @@ -0,0 +1,268 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/README b/README new file mode 100644 index 0000000..5b59c52 --- /dev/null +++ b/README @@ -0,0 +1,6 @@ +This is a barebones 3D fluid simulator for animating inviscid free surface liquids. +The "viewpls3D" sub-project provides a simple viewer for the resulting data. + +Questions or comments, contact: +Christopher Batty +christopherbatty@yahoo.com \ No newline at end of file diff --git a/array1.h b/array1.h new file mode 100644 index 0000000..cf3fc4c --- /dev/null +++ b/array1.h @@ -0,0 +1,793 @@ +#ifndef ARRAY1_H +#define ARRAY1_H + +#include +#include +#include +#include +#include +#include +#include + +// In this file: +// Array1: a dynamic 1D array for plain-old-data (not objects) +// WrapArray1: 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 struct Array1IsIntegral{ typedef Array1False type; }; // default: no (specializations to yes follow) +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; +template<> struct Array1IsIntegral{ typedef Array1True type; }; + +//============================================================================ +template +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 reverse_iterator; + typedef std::reverse_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; iULONG_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; iULONG_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 &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& operator=(const Array1& x) + { + if(max_n& x) const + { + if(n!=x.n) return false; + for(unsigned long i=0; i& x) const + { + if(n!=x.n) return true; + for(unsigned long i=0; i& x) const + { + for(unsigned long i=0; i(const Array1& x) const + { + for(unsigned long i=0; ix[i]) return true; + else if(x[i]>data[i]) return false; + } + return n>x.n; + } + + bool operator<=(const Array1& x) const + { + for(unsigned long i=0; i=(const Array1& x) const + { + for(unsigned long i=0; ix[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; imax_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 + void assign(InputIterator first, InputIterator last) + { assign_(first, last, typename Array1IsIntegral::type()); } + + template + void assign_(InputIterator first, InputIterator last, Array1True check) + { fill(first, last); } + + template + 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(i0); + 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(indexmax_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; i0); + return *data; + } + + T& front(void) + { + assert(n>0); + return *data; + } + + void grow(void) + { + unsigned long new_size=(max_n*sizeof(T)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& 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 Array1d; +typedef Array1 Array1f; +typedef Array1 Array1ll; +typedef Array1 Array1ull; +typedef Array1 Array1i; +typedef Array1 Array1ui; +typedef Array1 Array1s; +typedef Array1 Array1us; +typedef Array1 Array1c; +typedef Array1 Array1uc; + +//============================================================================ +template +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 reverse_iterator; + typedef std::reverse_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& a) + : n(a.n), max_n(a.max_n), data(a.data) + {} + + WrapArray1(std::vector& 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& x) const + { + if(n!=x.n) return false; + for(unsigned long i=0; i& x) const + { + if(n!=x.n) return true; + for(unsigned long i=0; i& x) const + { + for(unsigned long i=0; i(const WrapArray1& x) const + { + for(unsigned long i=0; ix[i]) return true; + else if(x[i]>data[i]) return false; + } + return n>x.n; + } + + bool operator<=(const WrapArray1& x) const + { + for(unsigned long i=0; i=(const WrapArray1& x) const + { + for(unsigned long i=0; ix[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 + void assign(InputIterator first, InputIterator last) + { assign_(first, last, typename Array1IsIntegral::type()); } + + template + void assign_(InputIterator first, InputIterator last, Array1True check) + { fill(first, last); } + + template + void assign_(InputIterator first, InputIterator last, Array1False check) + { + unsigned long i=0; + InputIterator p=first; + for(; p!=last; ++p, ++i){ + assert(i0); + 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(index0); + 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& x) + { + std::swap(n, x.n); + std::swap(max_n, x.max_n); + std::swap(data, x.data); + } +}; + +// some common arrays + +typedef WrapArray1 WrapArray1d; +typedef WrapArray1 WrapArray1f; +typedef WrapArray1 WrapArray1ll; +typedef WrapArray1 WrapArray1ull; +typedef WrapArray1 WrapArray1i; +typedef WrapArray1 WrapArray1ui; +typedef WrapArray1 WrapArray1s; +typedef WrapArray1 WrapArray1us; +typedef WrapArray1 WrapArray1c; +typedef WrapArray1 WrapArray1uc; + +#endif diff --git a/array2.h b/array2.h new file mode 100644 index 0000000..586fafc --- /dev/null +++ b/array2.h @@ -0,0 +1,280 @@ +#ifndef ARRAY2_H +#define ARRAY2_H + +#include "array1.h" +#include +#include +#include + +template > +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 + Array2(Array2& 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=0 && j=0 && i=0 && j& x) const + { return ni==x.ni && nj==x.nj && a==x.a; } + + bool operator!=(const Array2& x) const + { return ni!=x.ni || nj!=x.nj || a!=x.a; } + + bool operator<(const Array2& x) const + { + if(nix.ni) return false; + if(njx.nj) return false; + return a(const Array2& x) const + { + if(ni>x.ni) return true; else if(nix.nj) return true; else if(njx.a; + } + + bool operator<=(const Array2& x) const + { + if(nix.ni) return false; + if(njx.nj) return false; + return a<=x.a; + } + + bool operator>=(const Array2& x) const + { + if(ni>x.ni) return true; else if(nix.nj) return true; else if(nj=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=0 && j=0 && i=0 && j=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& 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 > Array2d; +typedef Array2 > Array2f; +typedef Array2 > Array2ll; +typedef Array2 > Array2ull; +typedef Array2 > Array2i; +typedef Array2 > Array2ui; +typedef Array2 > Array2s; +typedef Array2 > Array2us; +typedef Array2 > Array2c; +typedef Array2 > Array2uc; + +// and wrapped versions + +typedef Array2 > WrapArray2d; +typedef Array2 > WrapArray2f; +typedef Array2 > WrapArray2ll; +typedef Array2 > WrapArray2ull; +typedef Array2 > WrapArray2i; +typedef Array2 > WrapArray2ui; +typedef Array2 > WrapArray2s; +typedef Array2 > WrapArray2us; +typedef Array2 > WrapArray2c; +typedef Array2 > WrapArray2uc; + +#endif diff --git a/array2_utils.h b/array2_utils.h new file mode 100644 index 0000000..127825d --- /dev/null +++ b/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 +T interpolate_value(const Vec<2,S>& point, const Array2 >& 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 +float interpolate_gradient(Vec<2,T>& gradient, const Vec<2,S>& point, const Array2 >& 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 +void write_matlab_array(std::ostream &output, Array2 >&a, const char *variable_name, bool transpose=false) +{ + output< +#include +#include + +template > +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=0 && j=0 && k=0 && i=0 && j=0 && k& x) const + { return ni==x.ni && nj==x.nj && nk==x.nk && a==x.a; } + + bool operator!=(const Array3& x) const + { return ni!=x.ni || nj!=x.nj || nk!=x.nk || a!=x.a; } + + bool operator<(const Array3& x) const + { + if(nix.ni) return false; + if(njx.nj) return false; + if(nkx.nk) return false; + return a(const Array3& x) const + { + if(ni>x.ni) return true; else if(nix.nj) return true; else if(njx.nk) return true; else if(nkx.a; + } + + bool operator<=(const Array3& x) const + { + if(nix.ni) return false; + if(njx.nj) return false; + if(nkx.nk) return false; + return a<=x.a; + } + + bool operator>=(const Array3& x) const + { + if(ni>x.ni) return true; else if(nix.nj) return true; else if(njx.nk) return true; else if(nk=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=0 && j=0 && k=0 && i=0 && j=0 && k=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& 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 > Array3d; +typedef Array3 > Array3f; +typedef Array3 > Array3ll; +typedef Array3 > Array3ull; +typedef Array3 > Array3i; +typedef Array3 > Array3ui; +typedef Array3 > Array3s; +typedef Array3 > Array3us; +typedef Array3 > Array3c; +typedef Array3 > Array3uc; + +#endif diff --git a/array3_utils.h b/array3_utils.h new file mode 100644 index 0000000..8360258 --- /dev/null +++ b/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 +T interpolate_value(const Vec<3,S>& point, const Array3 >& 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 +T interpolate_gradient(Vec<3,T>& gradient, const Vec<3,S>& point, const Array3 >& 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 \ No newline at end of file diff --git a/fluidsim.cpp b/fluidsim.cpp new file mode 100644 index 0000000..1a5e558 --- /dev/null +++ b/fluidsim.cpp @@ -0,0 +1,676 @@ +#include "fluidsim.h" + +#include "array3_utils.h" +#include "levelset_util.h" +#include "pcgsolver/sparse_matrix.h" +#include "pcgsolver/pcg_solver.h" + +#include "volume_fractions.h" +#include "viscosity3d.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); + cell_solid_phi.resize(ni,nj,nk); + valid.resize(ni+1, nj+1, nk+1); + old_valid.resize(ni+1, nj+1, nk+1); + liquid_phi.resize(ni,nj,nk); + + //set viscosity to 1 + viscosity.resize(ni+1, nj+1, nk+1, 1.0); + + c_vol_liquid.resize(ni,nj,nk, 0); + u_vol_liquid.resize(ni+1,nj,nk, 0); + v_vol_liquid.resize(ni,nj+1,nk, 0); + w_vol_liquid.resize(ni,nj,nk+1, 0); + ex_vol_liquid.resize(ni,nj+1,nk+1, 0); + ey_vol_liquid.resize(ni+1,nj,nk+1, 0); + ez_vol_liquid.resize(ni+1,nj+1,nk, 0); + +} + +//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); + } + + for(int k = 0; k < nk; ++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+0.5f)*dx); + cell_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(" Solve viscosity"); + apply_viscosity(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; + } + +} + + +//Perform the viscosity solve + +void FluidSim::apply_viscosity(float dt) { + + printf("Computing weights\n"); + //Estimate weights at velocity and stress positions + compute_viscosity_weights(); + + printf("Setting up solve\n"); + //Set up and solve the linear system + solve_viscosity(dt); +} + +void FluidSim::solve_viscosity(float dt) { + + advance_viscosity_implicit_weighted(u, v, w, + u_vol_liquid, v_vol_liquid, w_vol_liquid, + c_vol_liquid, ex_vol_liquid, ey_vol_liquid, ez_vol_liquid, cell_solid_phi, viscosity, dt, dx); + +} + +float interpolate_phi(const Vec3f& point, const Array3f& grid, const Vec3f& origin, const float dx) { + float inv_dx = 1/dx; + Vec3f temp = (point-origin)*inv_dx; + return interpolate_value(temp, grid); +} + +void estimate_volume_fractions(Array3f& volumes, + const Vec3f& start_centre, const float dx, + const Array3f& phi, const Vec3f& phi_origin, const float phi_dx) +{ + + for(int k = 0; k < volumes.nk; ++k) for(int j = 0; j < volumes.nj; ++j) for(int i = 0; i < volumes.ni; ++i) { + Vec3f centre = start_centre + Vec3f(i*dx, j*dx, k*dx); + + float offset = 0.5f*dx; + + float phi000 = interpolate_phi(centre + Vec3f(-offset,-offset,-offset), phi, phi_origin, phi_dx); + float phi001 = interpolate_phi(centre + Vec3f(-offset,-offset,+offset), phi, phi_origin, phi_dx); + float phi010 = interpolate_phi(centre + Vec3f(-offset,+offset,-offset), phi, phi_origin, phi_dx); + float phi011 = interpolate_phi(centre + Vec3f(-offset,+offset,+offset), phi, phi_origin, phi_dx); + float phi100 = interpolate_phi(centre + Vec3f(+offset,-offset,-offset), phi, phi_origin, phi_dx); + float phi101 = interpolate_phi(centre + Vec3f(+offset,-offset,+offset), phi, phi_origin, phi_dx); + float phi110 = interpolate_phi(centre + Vec3f(+offset,+offset,-offset), phi, phi_origin, phi_dx); + float phi111 = interpolate_phi(centre + Vec3f(+offset,+offset,+offset), phi, phi_origin, phi_dx); + + volumes(i,j,k) = volume_fraction(phi000, phi100, phi010, phi110, phi001, phi101, phi011, phi111); + + } + +} + +void FluidSim::compute_viscosity_weights() { + + estimate_volume_fractions(c_vol_liquid, Vec3f(0.5f*dx, 0.5f*dx, 0.5f*dx), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(u_vol_liquid, Vec3f(0, 0.5f*dx, 0.5f*dx), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(v_vol_liquid, Vec3f(0.5f*dx, 0, 0.5f*dx), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(w_vol_liquid, Vec3f(0.5f*dx, 0.5f*dx, 0), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(ex_vol_liquid, Vec3f(0.5f*dx, 0, 0), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(ey_vol_liquid, Vec3f(0, 0.5f*dx, 0), dx, liquid_phi, Vec3f(0,0,0), dx); + estimate_volume_fractions(ez_vol_liquid, Vec3f(0, 0, 0.5f*dx), dx, liquid_phi, Vec3f(0,0,0), dx); + +} + + +//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 + temp_u(i,j,k) = 0;//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 + temp_v(i,j,k) = 0; //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 + temp_w(i,j,k) = 0; //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 = 1; k < w.nk-1; ++k) for(int j = 0; j < w.nj; ++j) for(int i = 0; i < w.ni; ++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; + + } + +} diff --git a/fluidsim.h b/fluidsim.h new file mode 100644 index 0000000..d90b3bf --- /dev/null +++ b/fluidsim.h @@ -0,0 +1,83 @@ +#ifndef FLUIDSIM_H +#define FLUIDSIM_H + +#include "array3.h" +#include "vec.h" +#include "pcgsolver/sparse_matrix.h" +#include "pcgsolver/pcg_solver.h" + +#include + +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 cell_solid_phi; + Array3f u_weights, v_weights, w_weights; + Array3c u_valid, v_valid, w_valid; + + + //Data for viscosity solve + Array3f u_vol_liquid, v_vol_liquid, w_vol_liquid, + ex_vol_liquid, ey_vol_liquid, ez_vol_liquid, c_vol_liquid; + + Array3f viscosity; + + std::vector particles; + float particle_radius; + + Array3f liquid_phi; + + //Data arrays for extrapolation + Array3c valid, old_valid; + + //Solver data + PCGSolver solver; + SparseMatrixd matrix; + std::vector rhs; + std::vector 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 apply_viscosity(float dt); + void project(float dt); + void constrain_velocity(); + + //helpers for pressure projection + void compute_weights(); + void solve_pressure(float dt); + void compute_phi(); + + //helpers for viscosity + void compute_viscosity_weights(); + void solve_viscosity(float dt); + +}; + + +#endif \ No newline at end of file diff --git a/levelset_util.cpp b/levelset_util.cpp new file mode 100644 index 0000000..20e1ccc --- /dev/null +++ b/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:0); + float list[] = { phi_bl, phi_br, phi_tr, phi_tl }; + + if(inside_count == 4) + return 1; + else if (inside_count == 3) { + //rotate until the positive value is in the first position + while(list[0] < 0) { + cycle_array(list,4); + } + + //Work out the area of the exterior triangle + float side0 = 1-fraction_inside(list[0], list[3]); + float side1 = 1-fraction_inside(list[0], list[1]); + return 1 - 0.5f*side0*side1; + } + else if(inside_count == 2) { + + //rotate until a negative value is in the first position, and the next negative is in either slot 1 or 2. + while(list[0] >= 0 || !(list[1] < 0 || list[2] < 0)) { + cycle_array(list,4); + } + + if(list[1] < 0) { //the matching signs are adjacent + float side_left = fraction_inside(list[0], list[3]); + float side_right = fraction_inside(list[1], list[2]); + return 0.5f*(side_left + side_right); + } + else { //matching signs are diagonally opposite + //determine the centre point's sign to disambiguate this case + float middle_point = 0.25f*(list[0] + list[1] + list[2] + list[3]); + if(middle_point < 0) { + float area = 0; + + //first triangle (top left) + float side1 = 1-fraction_inside(list[0], list[3]); + float side3 = 1-fraction_inside(list[2], list[3]); + + area += 0.5f*side1*side3; + + //second triangle (top right) + float side2 = 1-fraction_inside(list[2], list[1]); + float side0 = 1-fraction_inside(list[0], list[1]); + area += 0.5f*side0*side2; + + return 1-area; + } + else { + float area = 0; + + //first triangle (bottom left) + float side0 = fraction_inside(list[0], list[1]); + float side1 = fraction_inside(list[0], list[3]); + area += 0.5f*side0*side1; + + //second triangle (top right) + float side2 = fraction_inside(list[2], list[1]); + float side3 = fraction_inside(list[2], list[3]); + area += 0.5f*side2*side3; + return area; + } + + } + } + else if(inside_count == 1) { + //rotate until the negative value is in the first position + while(list[0] >= 0) { + cycle_array(list,4); + } + + //Work out the area of the interior triangle, and subtract from 1. + float side0 = fraction_inside(list[0], list[3]); + float side1 = fraction_inside(list[0], list[1]); + return 0.5f*side0*side1; + } + else + return 0; + +} diff --git a/levelset_util.h b/levelset_util.h new file mode 100644 index 0000000..58e84fb --- /dev/null +++ b/levelset_util.h @@ -0,0 +1,7 @@ +#ifndef LEVELSET_UTIL_H +#define LEVELSET_UTIL_H + +float fraction_inside(float phi_left, float phi_right); +float fraction_inside(float phi_bl, float phi_br, float phi_tl, float phi_tr); + +#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..073b51a --- /dev/null +++ b/main.cpp @@ -0,0 +1,93 @@ +#include +#include +#include +#include +#include +#include + +#include "fluidsim.h" +#include "array3_utils.h" + +using namespace std; + +int grid_resolution = 60; +float timestep = 0.01f; +int frame = 0; + +float grid_width = 1; + +FluidSim sim; + +float sphere_phi(const Vec3f& position, const Vec3f& centre, float radius) { + return (dist(position,centre) - radius); +} + +Vec3f c0(0.5f,0.5f,0.5f); +float rad0 = 0.35f; + +float boundary_phi(const Vec3f& position) { + return -sphere_phi(position, c0, rad0); +} + +float liquid_phi(const Vec3f& position) { + return sphere_phi(position, Vec3f(0.55f, 0.55f, 0.4f), 0.23f); +} + +void export_particles(string path, int frame, const std::vector& particles, float radius); + +//Main testing code +//------------- +int main(int argc, char **argv) +{ + + if(argc!=2){ + cerr << "The first parameter should be the folder to write the output liquid meshes into. (eg. c:\\output\\)" << endl; + return 1; + } + + string outpath(argv[1]); + + printf("Initializing data\n"); + sim.initialize(grid_width, grid_resolution, grid_resolution, grid_resolution); + + printf("Initializing boundary\n"); + sim.set_boundary(boundary_phi); + + printf("Initializing liquid\n"); + sim.set_liquid(liquid_phi); + + printf("Exporting initial data\n"); + export_particles(outpath, 0, sim.particles, sim.particle_radius); + + for(frame = 1; frame <1000; ++frame) { + printf("--------------------\nFrame %d\n", frame); + + //Simulate + printf("Simulating liquid\n"); + sim.advance(timestep); + + printf("Exporting particle data\n"); + export_particles(outpath, frame, sim.particles, sim.particle_radius); + } + + return 0; +} + + +void export_particles(string path, int frame, const std::vector& particles, float radius) { + //Write the output + + std::stringstream strout; + strout << path << "particles_" << frame << ".txt"; + string filepath = strout.str(); + + ofstream outfile(filepath.c_str()); + //write vertex count and particle radius + outfile << particles.size() << " " << radius << std::endl; + //write vertices + for(unsigned int i = 0; i < particles.size(); ++i) + outfile << particles[i][0] << " " << particles[i][1] << " " << particles[i][2] << std::endl; + outfile.close(); +} + + diff --git a/pcgsolver/blas_wrapper.h b/pcgsolver/blas_wrapper.h new file mode 100644 index 0000000..699080b --- /dev/null +++ b/pcgsolver/blas_wrapper.h @@ -0,0 +1,54 @@ +#ifndef BLAS_WRAPPER_H +#define BLAS_WRAPPER_H + +// Simple placeholder code for BLAS calls - replace with calls to a real BLAS library + +#include + +namespace BLAS{ + +// dot products ============================================================== + +inline double dot(const std::vector &x, const std::vector &y) +{ + //return cblas_ddot((int)x.size(), &x[0], 1, &y[0], 1); + + double sum = 0; + for(unsigned int i = 0; i < x.size(); ++i) + sum += x[i]*y[i]; + return sum; +} + +// inf-norm (maximum absolute value: index of max returned) ================== + +inline int index_abs_max(const std::vector &x) +{ + //return cblas_idamax((int)x.size(), &x[0], 1); + int maxind = 0; + double maxvalue = 0; + for(unsigned int i = 0; i < x.size(); ++i) { + if(fabs(x[i]) > maxvalue) { + maxvalue = fabs(x[i]); + maxind = i; + } + } + return maxind; +} + +// inf-norm (maximum absolute value) ========================================= +// technically not part of BLAS, but useful + +inline double abs_max(const std::vector &x) +{ return std::fabs(x[index_abs_max(x)]); } + +// saxpy (y=alpha*x+y) ======================================================= + +inline void add_scaled(double alpha, const std::vector &x, std::vector &y) +{ + //cblas_daxpy((int)x.size(), alpha, &x[0], 1, &y[0], 1); + for(unsigned int i = 0; i < x.size(); ++i) + y[i] += alpha*x[i]; +} + +} +#endif diff --git a/pcgsolver/pcg_solver.h b/pcgsolver/pcg_solver.h new file mode 100644 index 0000000..f81cb22 --- /dev/null +++ b/pcgsolver/pcg_solver.h @@ -0,0 +1,312 @@ +#ifndef PCG_SOLVER_H +#define PCG_SOLVER_H + +// Implements PCG with Modified Incomplete Cholesky (0) preconditioner. +// PCGSolver is the main class for setting up and solving a linear system. +// Note that this only handles symmetric positive (semi-)definite matrices, +// with guarantees made only for M-matrices (where off-diagonal entries are all +// non-positive, and row sums are non-negative). + +#include +#include "sparse_matrix.h" +#include "blas_wrapper.h" + +//============================================================================ +// A simple compressed sparse column data structure (with separate diagonal) +// for lower triangular matrices + +template +struct SparseColumnLowerFactor +{ + unsigned int n; + std::vector invdiag; // reciprocals of diagonal elements + std::vector value; // values below the diagonal, listed column by column + std::vector rowindex; // a list of all row indices, for each column in turn + std::vector colstart; // where each column begins in rowindex (plus an extra entry at the end, of #nonzeros) + std::vector adiag; // just used in factorization: minimum "safe" diagonal entry allowed + + explicit SparseColumnLowerFactor(unsigned int n_=0) + : n(n_), invdiag(n_), colstart(n_+1), adiag(n_) + {} + + void clear(void) + { + n=0; + invdiag.clear(); + value.clear(); + rowindex.clear(); + colstart.clear(); + adiag.clear(); + } + + void resize(unsigned int n_) + { + n=n_; + invdiag.resize(n); + colstart.resize(n+1); + adiag.resize(n); + } + + void write_matlab(std::ostream &output, const char *variable_name) + { + output< +void factor_modified_incomplete_cholesky0(const SparseMatrix &matrix, SparseColumnLowerFactor &factor, + T modification_parameter=0.97, T min_diagonal_ratio=0.25) +{ + // first copy lower triangle of matrix into factor (Note: assuming A is symmetric of course!) + factor.resize(matrix.n); + zero(factor.invdiag); // important: eliminate old values from previous solves! + factor.value.resize(0); + factor.rowindex.resize(0); + zero(factor.adiag); + for(unsigned int i=0; ii){ + factor.rowindex.push_back(matrix.index[i][j]); + factor.value.push_back(matrix.value[i][j]); + }else if(matrix.index[i][j]==i){ + factor.invdiag[i]=factor.adiag[i]=matrix.value[i][j]; + } + } + } + factor.colstart[matrix.n]=(unsigned int)factor.rowindex.size(); + // now do the incomplete factorization (figure out numerical values) + + // MATLAB code: + // L=tril(A); + // for k=1:size(L,2) + // L(k,k)=sqrt(L(k,k)); + // L(k+1:end,k)=L(k+1:end,k)/L(k,k); + // for j=find(L(:,k))' + // if j>k + // fullupdate=L(:,k)*L(j,k); + // incompleteupdate=fullupdate.*(A(:,j)~=0); + // missing=sum(fullupdate-incompleteupdate); + // L(j:end,j)=L(j:end,j)-incompleteupdate(j:end); + // L(j,j)=L(j,j)-omega*missing; + // end + // end + // end + + for(unsigned int k=0; k +void solve_lower(const SparseColumnLowerFactor &factor, const std::vector &rhs, std::vector &result) +{ + assert(factor.n==rhs.size()); + assert(factor.n==result.size()); + result=rhs; + for(unsigned int i=0; i +void solve_lower_transpose_in_place(const SparseColumnLowerFactor &factor, std::vector &x) +{ + assert(factor.n==x.size()); + assert(factor.n>0); + unsigned int i=factor.n; + do{ + --i; + for(unsigned int j=factor.colstart[i]; j +struct PCGSolver +{ + PCGSolver(void) + { + set_solver_parameters(1e-12, 100, 0.97, 0.25); + } + + void set_solver_parameters(T tolerance_factor_, int max_iterations_, T modified_incomplete_cholesky_parameter_=0.97, T min_diagonal_ratio_=0.25) + { + tolerance_factor=tolerance_factor_; + if(tolerance_factor<1e-30) tolerance_factor=1e-30; + max_iterations=max_iterations_; + modified_incomplete_cholesky_parameter=modified_incomplete_cholesky_parameter_; + min_diagonal_ratio=min_diagonal_ratio_; + } + + bool solve(const SparseMatrix &matrix, const std::vector &rhs, std::vector &result, T &residual_out, int &iterations_out) + { + unsigned int n=matrix.n; + if(m.size()!=n){ m.resize(n); s.resize(n); z.resize(n); r.resize(n); } + zero(result); + r=rhs; + residual_out=BLAS::abs_max(r); + if(residual_out==0) { + iterations_out=0; + return true; + } + double tol=tolerance_factor*residual_out; + + form_preconditioner(matrix); + apply_preconditioner(r, z); + double rho=BLAS::dot(z, r); + if(rho==0 || rho!=rho) { + iterations_out=0; + return false; + } + + s=z; + fixed_matrix.construct_from_matrix(matrix); + int iteration; + for(iteration=0; iteration ic_factor; // modified incomplete cholesky factor + std::vector m, z, s, r; // temporary vectors for PCG + FixedSparseMatrix fixed_matrix; // used within loop + + // parameters + T tolerance_factor; + int max_iterations; + T modified_incomplete_cholesky_parameter; + T min_diagonal_ratio; + + void form_preconditioner(const SparseMatrix& matrix) + { + factor_modified_incomplete_cholesky0(matrix, ic_factor); + } + + void apply_preconditioner(const std::vector &x, std::vector &result) + { + solve_lower(ic_factor, x, result); + solve_lower_transpose_in_place(ic_factor,result); + } +}; + +#endif diff --git a/pcgsolver/sparse_matrix.h b/pcgsolver/sparse_matrix.h new file mode 100644 index 0000000..9b54cfa --- /dev/null +++ b/pcgsolver/sparse_matrix.h @@ -0,0 +1,290 @@ +#ifndef SPARSE_MATRIX_H +#define SPARSE_MATRIX_H + +#include +#include +#include "util.h" + +//============================================================================ +// Dynamic compressed sparse row matrix. + +template +struct SparseMatrix +{ + unsigned int n; // dimension + std::vector > index; // for each row, a list of all column indices (sorted) + std::vector > value; // values corresponding to index + + explicit SparseMatrix(unsigned int n_=0, unsigned int expected_nonzeros_per_row=7) + : n(n_), index(n_), value(n_) + { + for(unsigned int i=0; ij) return 0; + } + return 0; + } + + void set_element(unsigned int i, unsigned int j, T new_value) + { + unsigned int k=0; + for(; kj){ + insert(index[i], k, j); + insert(value[i], k, new_value); + return; + } + } + index[i].push_back(j); + value[i].push_back(new_value); + } + + void add_to_element(unsigned int i, unsigned int j, T increment_value) + { + unsigned int k=0; + for(; kj){ + insert(index[i], k, j); + insert(value[i], k, increment_value); + return; + } + } + index[i].push_back(j); + value[i].push_back(increment_value); + } + + // assumes indices is already sorted + void add_sparse_row(unsigned int i, const std::vector &indices, const std::vector &values) + { + unsigned int j=0, k=0; + while(jindices[j]){ + insert(index[i], k, indices[j]); + insert(value[i], k, values[j]); + ++j; + }else{ + value[i][k]+=values[j]; + ++j; + ++k; + } + } + for(;j SparseMatrixf; +typedef SparseMatrix SparseMatrixd; + +// perform result=matrix*x +template +void multiply(const SparseMatrix &matrix, const std::vector &x, std::vector &result) +{ + assert(matrix.n==x.size()); + result.resize(matrix.n); + for(unsigned int i=0; i +void multiply_and_subtract(const SparseMatrix &matrix, const std::vector &x, std::vector &result) +{ + assert(matrix.n==x.size()); + result.resize(matrix.n); + for(unsigned int i=0; i +struct FixedSparseMatrix +{ + unsigned int n; // dimension + std::vector value; // nonzero values row by row + std::vector colindex; // corresponding column indices + std::vector rowstart; // where each row starts in value and colindex (and last entry is one past the end, the number of nonzeros) + + explicit FixedSparseMatrix(unsigned int n_=0) + : n(n_), value(0), colindex(0), rowstart(n_+1) + {} + + void clear(void) + { + n=0; + value.clear(); + colindex.clear(); + rowstart.clear(); + } + + void resize(int n_) + { + n=n_; + rowstart.resize(n+1); + } + + void construct_from_matrix(const SparseMatrix &matrix) + { + resize(matrix.n); + rowstart[0]=0; + for(unsigned int i=0; i FixedSparseMatrixf; +typedef FixedSparseMatrix FixedSparseMatrixd; + +// perform result=matrix*x +template +void multiply(const FixedSparseMatrix &matrix, const std::vector &x, std::vector &result) +{ + assert(matrix.n==x.size()); + result.resize(matrix.n); + for(unsigned int i=0; i +void multiply_and_subtract(const FixedSparseMatrix &matrix, const std::vector &x, std::vector &result) +{ + assert(matrix.n==x.size()); + result.resize(matrix.n); + for(unsigned int i=0; i +#include +#include +#include + +#ifndef M_PI +const double M_PI = 3.1415926535897932384626433832795; +#endif + +#ifdef WIN32 +#undef min +#undef max +#endif + +using std::min; +using std::max; +using std::swap; + +template +inline T sqr(const T& x) +{ return x*x; } + +template +inline T cube(const T& x) +{ return x*x*x; } + +template +inline T min(T a1, T a2, T a3) +{ return min(a1, min(a2, a3)); } + +template +inline T min(T a1, T a2, T a3, T a4) +{ return min(min(a1, a2), min(a3, a4)); } + +template +inline T min(T a1, T a2, T a3, T a4, T a5) +{ return min(min(a1, a2), min(a3, a4), a5); } + +template +inline T min(T a1, T a2, T a3, T a4, T a5, T a6) +{ return min(min(a1, a2), min(a3, a4), min(a5, a6)); } + +template +inline T max(T a1, T a2, T a3) +{ return max(a1, max(a2, a3)); } + +template +inline T max(T a1, T a2, T a3, T a4) +{ return max(max(a1, a2), max(a3, a4)); } + +template +inline T max(T a1, T a2, T a3, T a4, T a5) +{ return max(max(a1, a2), max(a3, a4), a5); } + +template +inline T max(T a1, T a2, T a3, T a4, T a5, T a6) +{ return max(max(a1, a2), max(a3, a4), max(a5, a6)); } + +template +inline void minmax(T a1, T a2, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T a4, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T a4, T a5, T& amin, T& amax) +{ + //@@@ the logic could be shortcircuited a lot! + amin=min(a1,a2,a3,a4,a5); + amax=max(a1,a2,a3,a4,a5); +} + +template +inline void minmax(T a1, T a2, T a3, T a4, T a5, T a6, T& amin, T& amax) +{ + //@@@ the logic could be shortcircuited a lot! + amin=min(a1,a2,a3,a4,a5,a6); + amax=max(a1,a2,a3,a4,a5,a6); +} + +template +inline void update_minmax(T a1, T& amin, T& amax) +{ + if(a1amax) amax=a1; +} + +// swap so that a +inline void sort(T &a, T &b) +{ + if(a>b) swap(a,b); +} + +// swap so that a +inline void sort(T &a, T &b, T &c) +{ + if(a>b) swap(a,b); + if(a>c) swap(a,c); + if(b>c) swap(b,c); +} + +// swap so that a +inline void sort(T &a, T &b, T &c, T &d) +{ + if(a>b) swap(a,b); + if(c>d) swap(c,d); + if(a>c) swap(a,c); + if(b>d) swap(b,d); + if(b>c) swap(b,c); +} + +template +inline T clamp(T a, T lower, T upper) +{ + if(aupper) return upper; + else return a; +} + +// only makes sense with T=float or double +template +inline T smooth_step(T r) +{ + if(r<0) return 0; + else if(r>1) return 1; + return r*r*r*(10+r*(-15+r*6)); +} + +// only makes sense with T=float or double +template +inline T smooth_step(T r, T r_lower, T r_upper, T value_lower, T value_upper) +{ return value_lower + smooth_step((r-r_lower)/(r_upper-r_lower)) * (value_upper-value_lower); } + +// only makes sense with T=float or double +template +inline T ramp(T r) +{ return smooth_step((r+1)/2)*2-1; } + +#ifdef WIN32 +inline int lround(double x) +{ + if(x>0) + return (x-floor(x)<0.5) ? (int)floor(x) : (int)ceil(x); + else + return (x-floor(x)<=0.5) ? (int)floor(x) : (int)ceil(x); +} + +inline double remainder(double x, double y) +{ + return x-std::floor(x/y+0.5)*y; +} +#endif + +inline unsigned int round_up_to_power_of_two(unsigned int n) +{ + int exponent=0; + --n; + while(n){ + ++exponent; + n>>=1; + } + return 1<1){ + ++exponent; + n>>=1; + } + return 1<>16); + i*=2654435769u; + i^=(i>>16); + i*=2654435769u; + return i; +} + +// the inverse of randhash +inline unsigned int unhash(unsigned int h) +{ + h*=340573321u; + h^=(h>>16); + h*=340573321u; + h^=(h>>16); + h*=340573321u; + h^=0xA3C59AC3u; + return h; +} + +// returns repeatable stateless pseudo-random number in [0,1] +inline double randhashd(unsigned int seed) +{ return randhash(seed)/(double)UINT_MAX; } +inline float randhashf(unsigned int seed) +{ return randhash(seed)/(float)UINT_MAX; } + +// returns repeatable stateless pseudo-random number in [a,b] +inline double randhashd(unsigned int seed, double a, double b) +{ return (b-a)*randhash(seed)/(double)UINT_MAX + a; } +inline float randhashf(unsigned int seed, float a, float b) +{ return ( (b-a)*randhash(seed)/(float)UINT_MAX + a); } + +inline int intlog2(int x) +{ + int exp=-1; + while(x){ + x>>=1; + ++exp; + } + return exp; +} + +template +inline void get_barycentric(T x, int& i, T& f, int i_low, int i_high) +{ + T s=std::floor(x); + i=(int)s; + if(ii_high-2){ + i=i_high-2; + f=1; + }else + f=(T)(x-s); +} + +template +inline S lerp(const S& value0, const S& value1, T f) +{ return (1-f)*value0 + f*value1; } + +template +inline S bilerp(const S& v00, const S& v10, + const S& v01, const S& v11, + T fx, T fy) +{ + return lerp(lerp(v00, v10, fx), + lerp(v01, v11, fx), + fy); +} + +template +inline S trilerp(const S& v000, const S& v100, + const S& v010, const S& v110, + const S& v001, const S& v101, + const S& v011, const S& v111, + T fx, T fy, T fz) +{ + return lerp(bilerp(v000, v100, v010, v110, fx, fy), + bilerp(v001, v101, v011, v111, fx, fy), + fz); +} + +template +inline S quadlerp(const S& v0000, const S& v1000, + const S& v0100, const S& v1100, + const S& v0010, const S& v1010, + const S& v0110, const S& v1110, + const S& v0001, const S& v1001, + const S& v0101, const S& v1101, + const S& v0011, const S& v1011, + const S& v0111, const S& v1111, + T fx, T fy, T fz, T ft) +{ + return lerp(trilerp(v0000, v1000, v0100, v1100, v0010, v1010, v0110, v1110, fx, fy, fz), + trilerp(v0001, v1001, v0101, v1101, v0011, v1011, v0111, v1111, fx, fy, fz), + ft); +} + +// f should be between 0 and 1, with f=0.5 corresponding to balanced weighting between w0 and w2 +template +inline void quadratic_bspline_weights(T f, T& w0, T& w1, T& w2) +{ + w0=T(0.5)*sqr(f-1); + w1=T(0.75)-sqr(f-T(0.5));; + w2=T(0.5)*sqr(f); +} + +// f should be between 0 and 1 +template +inline void cubic_interp_weights(T f, T& wneg1, T& w0, T& w1, T& w2) +{ + T f2(f*f), f3(f2*f); + wneg1=-T(1./3)*f+T(1./2)*f2-T(1./6)*f3; + w0=1-f2+T(1./2)*(f3-f); + w1=f+T(1./2)*(f2-f3); + w2=T(1./6)*(f3-f); +} + +template +inline S cubic_interp(const S& value_neg1, const S& value0, const S& value1, const S& value2, T f) +{ + T wneg1, w0, w1, w2; + cubic_interp_weights(f, wneg1, w0, w1, w2); + return wneg1*value_neg1 + w0*value0 + w1*value1 + w2*value2; +} + +template +void zero(std::vector& v) +{ for(int i=(int)v.size()-1; i>=0; --i) v[i]=0; } + +template +T abs_max(const std::vector& v) +{ + T m=0; + for(int i=(int)v.size()-1; i>=0; --i){ + if(std::fabs(v[i])>m) + m=std::fabs(v[i]); + } + return m; +} + +template +bool contains(const std::vector& a, T e) +{ + for(unsigned int i=0; i +void add_unique(std::vector& a, T e) +{ + for(unsigned int i=0; i +void insert(std::vector& a, unsigned int index, T e) +{ + a.push_back(a.back()); + for(unsigned int i=(unsigned int)a.size()-1; i>index; --i) + a[i]=a[i-1]; + a[index]=e; +} + +template +void erase(std::vector& a, unsigned int index) +{ + for(unsigned int i=index; i +void erase_swap(std::vector& a, unsigned int index) +{ + for(unsigned int i=index; i +void erase_unordered(std::vector& a, unsigned int index) +{ + a[index]=a.back(); + a.pop_back(); +} + +template +void erase_unordered_swap(std::vector& a, unsigned int index) +{ + swap(a[index], a.back()); + a.pop_back(); +} + +template +void find_and_erase_unordered(std::vector& a, const T& doomed_element) +{ + for(unsigned int i=0; i +void replace_once(std::vector& a, const T& old_element, const T& new_element) +{ + for(unsigned int i=0; i +void write_matlab(std::ostream& output, const std::vector& a, const char *variable_name, bool column_vector=true, int significant_digits=18) +{ + output< +#include +#include +#include "util.h" + +// Defines a thin wrapper around fixed size C-style arrays, using template parameters, +// which is useful for dealing with vectors of different dimensions. +// For example, float[3] is equivalent to Vec<3,float>. +// Entries in the vector are accessed with the overloaded [] operator, so +// for example if x is a Vec<3,float>, then the middle entry is x[1]. +// For convenience, there are a number of typedefs for abbreviation: +// Vec<3,float> -> Vec3f +// Vec<2,int> -> Vec2i +// and so on. +// Arithmetic operators are appropriately overloaded, and functions are defined +// for additional operations (such as dot-products, norms, cross-products, etc.) + +template +struct Vec +{ + T v[N]; + + Vec(void) + {} + + Vec(T value_for_all) + { for(unsigned int i=0; i + Vec(const S *source) + { for(unsigned int i=0; i + explicit Vec(const Vec& source) + { for(unsigned int i=0; i(T v0, T v1) + { + assert(N==2); + v[0]=v0; v[1]=v1; + } + + Vec(T v0, T v1, T v2) + { + assert(N==3); + v[0]=v0; v[1]=v1; v[2]=v2; + } + + Vec(T v0, T v1, T v2, T v3) + { + assert(N==4); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; + } + + Vec(T v0, T v1, T v2, T v3, T v4) + { + assert(N==5); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; + } + + Vec(T v0, T v1, T v2, T v3, T v4, T v5) + { + assert(N==6); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; v[5]=v5; + } + + T &operator[](int index) + { + assert(0<=index && (unsigned int)index operator+=(const Vec &w) + { + for(unsigned int i=0; i operator+(const Vec &w) const + { + Vec sum(*this); + sum+=w; + return sum; + } + + Vec operator-=(const Vec &w) + { + for(unsigned int i=0; i operator-(void) const // unary minus + { + Vec negative; + for(unsigned int i=0; i operator-(const Vec &w) const // (binary) subtraction + { + Vec diff(*this); + diff-=w; + return diff; + } + + Vec operator*=(T a) + { + for(unsigned int i=0; i operator*(T a) const + { + Vec w(*this); + w*=a; + return w; + } + + Vec operator*(const Vec &w) const + { + Vec componentwise_product; + for(unsigned int i=0; i operator/=(T a) + { + for(unsigned int i=0; i operator/(T a) const + { + Vec w(*this); + w/=a; + return w; + } +}; + +typedef Vec<2,double> Vec2d; +typedef Vec<2,float> Vec2f; +typedef Vec<2,int> Vec2i; +typedef Vec<2,unsigned int> Vec2ui; +typedef Vec<2,short> Vec2s; +typedef Vec<2,unsigned short> Vec2us; +typedef Vec<2,char> Vec2c; +typedef Vec<2,unsigned char> Vec2uc; + +typedef Vec<3,double> Vec3d; +typedef Vec<3,float> Vec3f; +typedef Vec<3,int> Vec3i; +typedef Vec<3,unsigned int> Vec3ui; +typedef Vec<3,short> Vec3s; +typedef Vec<3,unsigned short> Vec3us; +typedef Vec<3,char> Vec3c; +typedef Vec<3,unsigned char> Vec3uc; + +typedef Vec<4,double> Vec4d; +typedef Vec<4,float> Vec4f; +typedef Vec<4,int> Vec4i; +typedef Vec<4,unsigned int> Vec4ui; +typedef Vec<4,short> Vec4s; +typedef Vec<4,unsigned short> Vec4us; +typedef Vec<4,char> Vec4c; +typedef Vec<4,unsigned char> Vec4uc; + +typedef Vec<6,double> Vec6d; +typedef Vec<6,float> Vec6f; +typedef Vec<6,unsigned int> Vec6ui; +typedef Vec<6,int> Vec6i; +typedef Vec<6,short> Vec6s; +typedef Vec<6,unsigned short> Vec6us; +typedef Vec<6,char> Vec6c; +typedef Vec<6,unsigned char> Vec6uc; + + +template +T mag2(const Vec &a) +{ + T l=sqr(a.v[0]); + for(unsigned int i=1; i +T mag(const Vec &a) +{ return sqrt(mag2(a)); } + +template +inline T dist2(const Vec &a, const Vec &b) +{ + T d=sqr(a.v[0]-b.v[0]); + for(unsigned int i=1; i +inline T dist(const Vec &a, const Vec &b) +{ return std::sqrt(dist2(a,b)); } + +template +inline void normalize(Vec &a) +{ a/=mag(a); } + +template +inline Vec normalized(const Vec &a) +{ return a/mag(a); } + +template +inline T infnorm(const Vec &a) +{ + T d=std::fabs(a.v[0]); + for(unsigned int i=1; i +void zero(Vec &a) +{ + for(unsigned int i=0; i +std::ostream &operator<<(std::ostream &out, const Vec &v) +{ + out< +std::istream &operator>>(std::istream &in, Vec &v) +{ + in>>v.v[0]; + for(unsigned int i=1; i>v.v[i]; + return in; +} + +template +inline bool operator==(const Vec &a, const Vec &b) +{ + bool t = (a.v[0] == b.v[0]); + unsigned int i=1; + while(i +inline bool operator!=(const Vec &a, const Vec &b) +{ + bool t = (a.v[0] != b.v[0]); + unsigned int i=1; + while(i +inline Vec operator*(T a, const Vec &v) +{ + Vec w(v); + w*=a; + return w; +} + +template +inline T min(const Vec &a) +{ + T m=a.v[0]; + for(unsigned int i=1; i +inline Vec min_union(const Vec &a, const Vec &b) +{ + Vec m; + for(unsigned int i=0; i +inline Vec max_union(const Vec &a, const Vec &b) +{ + Vec m; + for(unsigned int i=0; i b.v[i]) ? m.v[i]=a.v[i] : m.v[i]=b.v[i]; + return m; +} + +template +inline T max(const Vec &a) +{ + T m=a.v[0]; + for(unsigned int i=1; im) m=a.v[i]; + return m; +} + +template +inline T dot(const Vec &a, const Vec &b) +{ + T d=a.v[0]*b.v[0]; + for(unsigned int i=1; i +inline Vec<2,T> rotate(const Vec<2,T>& a, float angle) +{ + T c = cos(angle); + T s = sin(angle); + return Vec<2,T>(c*a[0] - s*a[1],s*a[0] + c*a[1]); // counter-clockwise rotation +} + +template +inline Vec<2,T> perp(const Vec<2,T> &a) +{ return Vec<2,T>(-a.v[1], a.v[0]); } // counter-clockwise rotation by 90 degrees + +template +inline T cross(const Vec<2,T> &a, const Vec<2,T> &b) +{ return a.v[0]*b.v[1]-a.v[1]*b.v[0]; } + +template +inline Vec<3,T> cross(const Vec<3,T> &a, const Vec<3,T> &b) +{ return Vec<3,T>(a.v[1]*b.v[2]-a.v[2]*b.v[1], a.v[2]*b.v[0]-a.v[0]*b.v[2], a.v[0]*b.v[1]-a.v[1]*b.v[0]); } + +template +inline T triple(const Vec<3,T> &a, const Vec<3,T> &b, const Vec<3,T> &c) +{ return a.v[0]*(b.v[1]*c.v[2]-b.v[2]*c.v[1]) + +a.v[1]*(b.v[2]*c.v[0]-b.v[0]*c.v[2]) + +a.v[2]*(b.v[0]*c.v[1]-b.v[1]*c.v[0]); } + +template +inline unsigned int hash(const Vec &a) +{ + unsigned int h=a.v[0]; + for(unsigned int i=1; i +inline void assign(const Vec &a, T &a0, T &a1) +{ + assert(N==2); + a0=a.v[0]; a1=a.v[1]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2) +{ + assert(N==3); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2, T &a3) +{ + assert(N==4); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2, T &a3, T &a4, T &a5) +{ + assert(N==6); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3]; a4=a.v[4]; a5=a.v[5]; +} + +template +inline Vec round(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec floor(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec ceil(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec fabs(const Vec &a) +{ + Vec result; + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, + Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, const Vec &x4, + Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, const Vec &x4, + const Vec &x5, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void update_minmax(const Vec &x, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/viewpls3D/gluvi.cpp b/viewpls3D/gluvi.cpp new file mode 100644 index 0000000..72a12fe --- /dev/null +++ b/viewpls3D/gluvi.cpp @@ -0,0 +1,1114 @@ +#include +#include +#include +#include +#include "gluvi.h" +#include "vec.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846f +#endif + +using namespace std; + +namespace Gluvi{ + +Target3D:: +Target3D(float target_[3], float dist_, float heading_, float pitch_, float fovy_, float near_clip_factor_, float far_clip_factor_) + : dist(dist_), heading(heading_), pitch(pitch_), fovy(fovy_), + near_clip_factor(near_clip_factor_), far_clip_factor(far_clip_factor_), action_mode(INACTIVE) +{ + if(target_){ + target[0]=target_[0]; + target[1]=target_[1]; + target[2]=target_[2]; + }else{ + target[0]=0; + target[1]=0; + target[2]=0; + } + default_target[0]=target[0]; + default_target[1]=target[1]; + default_target[2]=target[2]; + default_dist=dist; + default_heading=heading; + default_pitch=pitch; +} + +void Target3D:: +click(int button, int state, int x, int y) +{ + if(state==GLUT_UP) + action_mode=INACTIVE; + else if(button==GLUT_LEFT_BUTTON) + action_mode=ROTATE; + else if(button==GLUT_MIDDLE_BUTTON) + action_mode=TRUCK; + else if(button==GLUT_RIGHT_BUTTON) + action_mode=DOLLY; + oldmousex=x; + oldmousey=y; +} + +void Target3D:: +drag(int x, int y) +{ + switch(action_mode){ + case INACTIVE: + return; // nothing to do + case ROTATE: + heading+=0.007f*(oldmousex-x); + if(heading<-M_PI) heading+=2*M_PI; + else if(heading>M_PI) heading-=2*M_PI; + pitch+=0.007f*(oldmousey-y); + if(pitch<-0.5f*M_PI) pitch=-0.5f*M_PI; + else if(pitch>0.5f*M_PI) pitch=0.5f*M_PI; + break; + case TRUCK: + target[0]+=(0.002f*dist)*cos(heading)*(oldmousex-x); + target[1]-=(0.002f*dist)*(oldmousey-y); + target[2]-=(0.002f*dist)*sin(heading)*(oldmousex-x); + break; + case DOLLY: + dist*=pow(1.01f, oldmousey-y + x-oldmousex); + break; + } + oldmousex=x; + oldmousey=y; + glutPostRedisplay(); +} + +void Target3D:: +return_to_default(void) +{ + target[0]=default_target[0]; + target[1]=default_target[1]; + target[2]=default_target[2]; + dist=default_dist; + heading=default_heading; + pitch=default_pitch; +} + +void Target3D:: +transform_mouse(int x, int y, float ray_origin[3], float ray_direction[3]) +{ + float ch=cos(heading), sh=sin(heading); + float cp=cos(pitch), sp=sin(pitch); + + ray_origin[0]=target[0]+dist*sh*cp; + ray_origin[1]=target[1]-dist*sp; + ray_origin[2]=target[2]+dist*ch*cp; + + float scale=0.5f*tan(fovy)/winheight; + float camx=(x-0.5f*winwidth)*scale, camy=(0.5f*winheight-y)*scale, camz=-1.0; // in camera coordinates, this is ray_direction (but not normalized) + // now need to rotate into world space from camera space + float px=camx, py=camy*cp-camz*sp, pz=camy*sp+camz*cp; + ray_direction[0]=px*ch+pz*sh; + ray_direction[1]=py; + ray_direction[2]=-px*sh+pz*ch; + + //@@@ is this right to get us to the near clipping plane? + ray_origin[0]+=near_clip_factor*dist*ray_direction[0]; + ray_origin[1]+=near_clip_factor*dist*ray_direction[1]; + ray_origin[2]+=near_clip_factor*dist*ray_direction[2]; + + // normalize direction vector + float mag=sqrt(ray_direction[0]*ray_direction[0] + + ray_direction[1]*ray_direction[1] + + ray_direction[2]*ray_direction[2]); + ray_direction[0]/=mag; + ray_direction[1]/=mag; + ray_direction[2]/=mag; +} + +void Target3D:: +get_viewing_direction(float direction[3]) +{ + float ch=cos(heading), sh=sin(heading); + float cp=cos(pitch), sp=sin(pitch); + direction[0]=-sh*cp; + direction[1]=sp; + direction[2]=-ch*cp; +} + +void Target3D:: +gl_transform(void) +{ + glViewport(0, 0, (GLsizei)winwidth, (GLsizei)winheight); + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluPerspective(fovy, winwidth/(float)winheight, near_clip_factor*dist, far_clip_factor*dist); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + GLfloat pos[3]; + pos[0]=target[0]-dist*sin(heading)*cos(pitch); + pos[1]=target[1]-dist*sin(pitch); + pos[2]=target[2]-dist*cos(heading)*cos(pitch); + glTranslatef(0, 0, -dist); // translate target dist away in the z direction + glRotatef(-180/M_PI*pitch, 1, 0, 0); // rotate pitch in the yz plane + glRotatef(-180/M_PI*heading, 0, 1, 0); // rotate heading in the xz plane + glTranslatef(-target[0], -target[1], -target[2]); // translate target to origin +} + +void Target3D:: +export_rib(ostream &output) +{ + output<<"Clipping "<M_PI) heading-=2*M_PI; + pitch+=0.007f*(oldmousey-y); + if(pitch<-0.5f*M_PI) pitch=-0.5f*M_PI; + else if(pitch>0.5f*M_PI) pitch=0.5f*M_PI; + break; + case TRUCK: + target[0]+=(0.002f*dist)*cos(heading)*(oldmousex-x); + target[1]-=(0.002f*dist)*(oldmousey-y); + target[2]-=(0.002f*dist)*sin(heading)*(oldmousex-x); + break; + case DOLLY: + dist*=pow(1.01f, oldmousey-y + x-oldmousex); + break; + } + oldmousex=x; + oldmousey=y; + glutPostRedisplay(); +} + +void TargetOrtho3D:: +return_to_default(void) +{ + target[0]=default_target[0]; + target[1]=default_target[1]; + target[2]=default_target[2]; + dist=default_dist; + heading=default_heading; + pitch=default_pitch; +} + +void TargetOrtho3D:: +transform_mouse(int x, int y, float ray_origin[3], float ray_direction[3]) +{ + // @@@ unimplemented +} + +void TargetOrtho3D:: +get_viewing_direction(float direction[3]) +{ + float ch=cos(heading), sh=sin(heading); + float cp=cos(pitch), sp=sin(pitch); + direction[0]=-sh*cp; + direction[1]=sp; + direction[2]=-ch*cp; +} + +void TargetOrtho3D:: +gl_transform(void) +{ + glViewport(0, 0, (GLsizei)winwidth, (GLsizei)winheight); + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + float halfheight=0.5f*height_factor*dist, halfwidth=halfheight*winwidth/(float)winheight; + glOrtho(-halfwidth, halfwidth, -halfheight, halfheight, near_clip_factor*dist, far_clip_factor*dist); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + GLfloat pos[3]; + pos[0]=target[0]-dist*sin(heading)*cos(pitch); + pos[1]=target[1]-dist*sin(pitch); + pos[2]=target[2]-dist*cos(heading)*cos(pitch); + glTranslatef(0, 0, -dist); // translate target dist away in the z direction + glRotatef(-180/M_PI*pitch, 1, 0, 0); // rotate pitch in the yz plane + glRotatef(-180/M_PI*heading, 0, 1, 0); // rotate heading in the xz plane + glTranslatef(-target[0], -target[1], -target[2]); // translate target to origin +} + +void TargetOrtho3D:: +export_rib(ostream &output) +{ + output<<"Clipping "< desired_height*winwidth) + desired_height=winheight*desired_width/winwidth; + else + desired_width=winwidth*desired_height/winheight; + left+=0.5f*(x+clickx)*height/winheight-0.5f*desired_width; + bottom+=(winheight-0.5f*(y+clicky))*height/winheight-0.5f*desired_height; + height=desired_height; + }else{ + // zoom in by some constant factor on the mouse click + float factor=0.70710678118654752440084f; + left+=(1-factor)*height*(x/(float)winheight); + bottom+=(1-factor)*height*(1-y/(float)winheight); + height*=factor; + } + glutPostRedisplay(); + break; + case ZOOM_OUT: + // zoom out by some constant factor + { + float factor=1.41421356237309504880168f; + left-=0.5f*(factor-1)*winwidth*height/winheight; + bottom-=0.5f*(factor-1)*height; + height*=factor; + } + glutPostRedisplay(); + break; + default: + ;// nothing to do + } + action_mode=INACTIVE; + + }else if(button==GLUT_LEFT_BUTTON) + action_mode=PAN; + else if(button==GLUT_MIDDLE_BUTTON){ + clickx=x; + clicky=y; + action_mode=ZOOM_IN; + }else if(button==GLUT_RIGHT_BUTTON) + action_mode=ZOOM_OUT; + moved_since_mouse_down=false; + oldmousex=x; + oldmousey=y; +} + +void PanZoom2D:: +drag(int x, int y) +{ + if(x!=oldmousex || y!=oldmousey){ + moved_since_mouse_down=true; + if(action_mode==PAN){ + float r=height/winheight; + left-=r*(x-oldmousex); + bottom+=r*(y-oldmousey); + glutPostRedisplay(); + }else if(action_mode==ZOOM_IN) + glutPostRedisplay(); + oldmousex=x; + oldmousey=y; + } +} + +void PanZoom2D:: +return_to_default(void) +{ + bottom=default_bottom; + left=default_left; + height=default_height; +} + +void PanZoom2D:: +transform_mouse(int x, int y, float coords[2]) +{ + float r=height/winheight; + coords[0]=x*r+left; + coords[1]=(winheight-y)*r+bottom; +} + +void PanZoom2D:: +gl_transform(void) +{ + glViewport(0, 0, (GLsizei)winwidth, (GLsizei)winheight); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + glOrtho(left, left+(height*winwidth)/winheight, bottom, bottom+height, 0, 1); + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); +} + +void PanZoom2D:: +export_rib(ostream &output) +{ + // no projection matrix + output<<"Clipping 1 2000"<winheight) scalefactor=2.0f/height; + else scalefactor=2.0f/(winwidth*height/winheight); + output<<"Scale "<dispx && x<=dispx+width && y=dispy-height){ + status=HIGHLIGHTED; + glutPostRedisplay(); + return true; + }else if(state==GLUT_UP && status!=UNINVOLVED){ + status=UNINVOLVED; + glutPostRedisplay(); + if(x>=dispx && x=dispy-height) + action(); + return true; + }else + return false; +} + +void Button:: +drag(int x, int y) +{ + // needs to control highlighting (SELECTED vs. HIGHLIGHTED) + if(status==SELECTED && x>=dispx && x=dispy-height){ + status=HIGHLIGHTED; + glutPostRedisplay(); + }else if(status==HIGHLIGHTED && !(x>=dispx && x=dispy-height)){ + status=SELECTED; + glutPostRedisplay(); + } +} + +//================================================================================= + +Slider:: +Slider(const char *text_, int length_, int position_, int justify_) + : status(UNINVOLVED), text(text_), length(length_), justify(justify_), position(position_) +{} + +void Slider:: +display(int x, int y) +{ + dispx=x; + dispy=y; + width=glutBitmapLength(GLUT_BITMAP_HELVETICA_12, (const unsigned char*)text); + if(widthscrollxmin+position+2 && x<=scrollxmin+position+11 && y=scrollymin+2){ + status=SELECTED; + clickx=x; + glutPostRedisplay(); + return true; + }else if(status!=UNINVOLVED && state==GLUT_UP){ + status=UNINVOLVED; + glutPostRedisplay(); + return true; + }else + return false; +} + +void Slider:: +drag(int x, int y) +{ + if(status==SELECTED){ + glutPostRedisplay(); + int newposition=position+(x-clickx); + clickx=x; + if(newposition<0){ + clickx+=(0-newposition); + newposition=0; + }else if(newposition>length){ + clickx+=(length-newposition); + newposition=length; + } + if(newposition!=position){ + position=newposition; + action(); + glutPostRedisplay(); + } + } +} + +//================================================================================= + +WidgetList:: +WidgetList(int indent_, bool hidden_) + : indent(indent_), hidden(hidden_), downclicked_member(-1) +{ +} + +void WidgetList:: +display(int x, int y) +{ + dispx=x; + dispy=y; + if(hidden){ + width=height=0; + }else{ + height=0; + for(unsigned int i=0; idisplay(x+indent, y-height); + height+=list[i]->height; + width=(widthwidth) ? indent+list[i]->width : width; + } + } +} + +bool WidgetList:: +click(int state, int x, int y) +{ + //if(hidden || x=dispx+width || y>=dispy || yclick(state, x, y)){ + downclicked_member=i; + return true; + } + } + }else if(state==GLUT_UP && downclicked_member>=0){ + list[downclicked_member]->click(state, x, y); + downclicked_member=-1; + } + return false; +} + +void WidgetList:: +drag(int x, int y) +{ + if(downclicked_member>=0) + list[downclicked_member]->drag(x, y); +} + +//================================================================================= + +static void gluviReshape(int w, int h) +{ + winwidth=w; + winheight=h; + glutPostRedisplay(); // triggers the camera to adjust itself to the new dimensions +} + +//================================================================================= + +static void gluviDisplay() +{ + glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT); + + // draw the scene + if(camera) camera->gl_transform(); + if(userDisplayFunc) userDisplayFunc(); + + // now draw widgets on top + glPushAttrib(GL_CURRENT_BIT|GL_ENABLE_BIT|GL_LINE_BIT); + glDisable(GL_DEPTH_TEST); + glDisable(GL_LIGHTING); + glLineWidth(1); + // and probably more needs setting before widgets + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0, winwidth, 0, winheight); + + root.display(0, winheight); + + // and allow the camera to draw something on screen (e.g. for zooming extent) + if(camera) camera->display_screen(); + + glPopAttrib(); + + glutSwapBuffers(); +} + +//================================================================================= + +static enum {NOBODY, CAMERA, WIDGETS, USER} mouse_owner=NOBODY; + +static void gluviMouse(int button, int state, int x, int y) +{ + if(state==GLUT_DOWN){ + int mods=glutGetModifiers(); + if(camera && mods==GLUT_ACTIVE_SHIFT){ + camera->click(button, state, x, y); + mouse_owner=CAMERA; + }else if(button==GLUT_LEFT_BUTTON && root.click(state, x, winheight-y)){ + mouse_owner=WIDGETS; + }else if(userMouseFunc){ + userMouseFunc(button, state, x, y); + mouse_owner=USER; + } + }else{ // mouse up - send event to whoever got the mouse down + switch(mouse_owner){ + case CAMERA: + camera->click(button, state, x, y); + break; + case WIDGETS: + root.click(state, x, winheight-y); + break; + case USER: + if(userMouseFunc) userMouseFunc(button, state, x, y); + break; + default: + ;// nothing to do + } + mouse_owner=NOBODY; + } +} + +//================================================================================= + +static void gluviDrag(int x, int y) +{ + switch(mouse_owner){ + case CAMERA: + camera->drag(x, y); + break; + case WIDGETS: + root.drag(x, winheight-y); + break; + case USER: + if(userDragFunc) userDragFunc(x, y); + break; + default: + ;// nothing to do + } +} + +//================================================================================= + +void ppm_screenshot(const char *filename_format, ...) +{ + va_list ap; + va_start(ap, filename_format); +#ifdef _MSC_VER +#define FILENAMELENGTH 256 + char filename[FILENAMELENGTH]; + _vsnprintf(filename, FILENAMELENGTH, filename_format, ap); + ofstream out(filename, ofstream::binary); +#else + char *filename; + vasprintf(&filename, filename_format, ap); + ofstream out(filename, ofstream::binary); + free(filename); +#endif + if(!out) return; + GLubyte *image_buffer=new GLubyte[3*winwidth*winheight]; + glReadBuffer(GL_FRONT); + glReadPixels(0, 0, winwidth, winheight, GL_RGB, GL_UNSIGNED_BYTE, image_buffer); + out<<"P6\n"<>8)%256); + output.put(v%256); +} + +static void write_big_endian_uint(std::ostream &output, unsigned int v) +{ + output.put((v>>24)%256); + output.put((v>>16)%256); + output.put((v>>8)%256); + output.put(v%256); +} + +void sgi_screenshot(const char *filename_format, ...) +{ + va_list ap; + va_start(ap, filename_format); +#ifdef _MSC_VER +#define FILENAMELENGTH 256 + char filename[FILENAMELENGTH]; + _vsnprintf(filename, FILENAMELENGTH, filename_format, ap); + ofstream output(filename, ofstream::binary); +#else + char *filename; + vasprintf(&filename, filename_format, ap); + ofstream output(filename, ofstream::binary); +#endif + if(!output) return; + // first write the SGI header + write_big_endian_ushort(output, 474); // magic number to identify this as an SGI image file + output.put(0); // uncompressed + output.put(1); // use 8-bit colour depth + write_big_endian_ushort(output, 3); // number of dimensions + write_big_endian_ushort(output, winwidth); // x size + write_big_endian_ushort(output, winheight); // y size + write_big_endian_ushort(output, 3); // three colour channels (z size) + write_big_endian_uint(output, 0); // minimum pixel value + write_big_endian_uint(output, 255); // maximum pixel value + write_big_endian_uint(output, 0); // dummy spacing + // image name + int i; + for(i=0; i<80 && filename[i]; ++i) + output.put(filename[i]); + for(; i<80; ++i) + output.put(0); + write_big_endian_uint(output, 0); // colormap is normal + for(i=0; i<404; ++i) output.put(0); // filler to complete header + // now write the SGI image data + GLubyte *image_buffer=new GLubyte[winwidth*winheight]; + glReadBuffer(GL_FRONT); + glReadPixels(0, 0, winwidth, winheight, GL_RED, GL_UNSIGNED_BYTE, image_buffer); + output.write((const char*)image_buffer, winwidth*winheight); + glReadPixels(0, 0, winwidth, winheight, GL_GREEN, GL_UNSIGNED_BYTE, image_buffer); + output.write((const char*)image_buffer, winwidth*winheight); + glReadPixels(0, 0, winwidth, winheight, GL_BLUE, GL_UNSIGNED_BYTE, image_buffer); + output.write((const char*)image_buffer, winwidth*winheight); + delete[] image_buffer; +#ifndef _MSC_VER + free(filename); +#endif +} + +void set_generic_lights(void) +{ + glEnable(GL_LIGHTING); + { + GLfloat ambient[4] = {.3f, .3f, .3f, 1}; + glLightModelfv (GL_LIGHT_MODEL_AMBIENT,ambient); + } + { + GLfloat color[4] = {.8f, .8f, .8f, 1}; + glLightfv (GL_LIGHT0, GL_DIFFUSE, color); + glLightfv (GL_LIGHT0, GL_SPECULAR, color); + glEnable (GL_LIGHT0); + } + { + GLfloat color[4] = {.4f, .4f, .4f, 1}; + glLightfv (GL_LIGHT1, GL_DIFFUSE, color); + glLightfv (GL_LIGHT1, GL_SPECULAR, color); + glEnable (GL_LIGHT1); + } + { + GLfloat color[4] = {.2f, .2f, .2f, 1}; + glLightfv (GL_LIGHT2, GL_DIFFUSE, color); + glLightfv (GL_LIGHT2, GL_SPECULAR, color); + glEnable (GL_LIGHT2); + } +} + +void set_generic_material(float r, float g, float b, GLenum face) +{ + GLfloat ambient[4], diffuse[4], specular[4]; + ambient[0]=0.1f*r+0.03f; ambient[1]=0.1f*g+0.03f; ambient[2]=0.1f*b+0.03f; ambient[3]=1; + diffuse[0]=0.7f*r; diffuse[1]=0.7f*g; diffuse[2]=0.7f*b; diffuse[3]=1; + specular[0]=0.1f*r+0.1f; specular[1]=0.1f*g+0.1f; specular[2]=0.1f*b+0.1f; specular[3]=1; + glMaterialfv(face, GL_AMBIENT, ambient); + glMaterialfv(face, GL_DIFFUSE, diffuse); + glMaterialfv(face, GL_SPECULAR, specular); + glMaterialf(face, GL_SHININESS, 32); +} + +void set_matte_material(float r, float g, float b, GLenum face) +{ + GLfloat ambient[4], diffuse[4], specular[4]; + ambient[0]=0.1f*r+0.03f; ambient[1]=0.1f*g+0.03f; ambient[2]=0.1f*b+0.03f; ambient[3]=1; + diffuse[0]=0.9f*r; diffuse[1]=0.9f*g; diffuse[2]=0.9f*b; diffuse[3]=1; + specular[0]=0; specular[1]=0; specular[2]=0; specular[3]=1; + glMaterialfv(face, GL_AMBIENT, ambient); + glMaterialfv(face, GL_DIFFUSE, diffuse); + glMaterialfv(face, GL_SPECULAR, specular); +} + +/** + * Draw a vector in 3D. If arrow_head_length not specified, set it to 10% of vector length. + * Mar 29, 2006 + */ +void draw_3d_arrow(const float base[3], const float point[3], float arrow_head_length) +{ + //glPushAttrib(GL_CURRENT_BIT|GL_ENABLE_BIT|GL_LINE_BIT); + //glDisable(GL_LIGHTING); + glLineWidth(1); + glColor3f(0.5,0.5,0.5); + + Vec3f w(point[0]-base[0], point[1]-base[1], point[2]-base[2]); + double len = mag(w); + w = w / (float)len; // normalize to build coordinate system + + // create a coordinate system from the vector: + // get a vector perp to w and y-axis + Vec3f u = cross(w, Vec3f(0,1,0)); + + // If vector is parallel to the y-axis, use the x-axis + if (mag(u) == 0) + u = cross(w, Vec3f(1,0,0)); + + // get a vector perp to w and u + Vec3f v = cross(w, u/mag(u)); + v = v/mag(v); + + if (!arrow_head_length) + arrow_head_length = 0.1f * (float)len; + + // arrow head points + //@@@@@@@ POSSIBILITY: CREATE FOUR ARROW HEAD SEGMENTS INSTEAD OF TWO + Vec3f arrow1, arrow2; + arrow1[0] = point[0] + arrow_head_length * (v[0] - w[0]); + arrow1[1] = point[1] + arrow_head_length * (v[1] - w[1]); + arrow1[2] = point[2] + arrow_head_length * (v[2] - w[2]); + arrow2[0] = point[0] + arrow_head_length * (-v[0] - w[0]); + arrow2[1] = point[1] + arrow_head_length * (-v[1] - w[1]); + arrow2[2] = point[2] + arrow_head_length * (-v[2] - w[2]); + + glBegin(GL_LINES); + glVertex3f(base[0], base[1], base[2]); + glVertex3f(point[0], point[1], point[2]); + glVertex3f(point[0], point[1], point[2]); + glVertex3f(arrow1[0], arrow1[1], arrow1[2]); + glVertex3f(point[0], point[1], point[2]); + glVertex3f(arrow2[0], arrow2[1], arrow2[2]); + glEnd(); + + //glPopAttrib(); +} + +// draw_2d_arrow assumptions: +// line width, point size, and color are set by the user prior to calling the routine +void draw_2d_arrow(const Vec2f base, const Vec2f point, float arrow_head_length) +{ + //glPushAttrib(GL_CURRENT_BIT|GL_ENABLE_BIT|GL_LINE_BIT); + //glDisable(GL_LIGHTING); + + Vec2f w = point-base; + double len = mag(w); + + if (len==0) { + glBegin(GL_POINTS); + glVertex2f(base[0], base[1]); + glEnd(); + return; + } + + w = w / (float)len; // normalize to build coordinate system + + // u = w + 90 + // using rotation matrix 0 1 + // -1 0 + Vec2f u = Vec2f(1*w[1], -1*w[0]); + u = u/mag(u); + + // v = w - 90 (in fact v=-u) + Vec2f v = Vec2f(-1*w[1], 1*w[0]); + v = v/mag(v); + + if (!arrow_head_length) { + arrow_head_length = 0.1f * (float)len; + } + + // arrow head points + Vec2f arrow1, arrow2; + arrow1 = point + arrow_head_length * (v-w); + arrow2 = point + arrow_head_length * (u-w); + + glBegin(GL_LINES); + glVertex2f(base[0], base[1]); + glVertex2f(point[0], point[1]); + glVertex2f(point[0], point[1]); + glVertex2f(arrow1[0], arrow1[1]); + glVertex2f(point[0], point[1]); + glVertex2f(arrow2[0], arrow2[1]); + glEnd(); + + //glPopAttrib(); +} + + +void draw_coordinate_grid(float size, int spacing) +{ + glPushAttrib(GL_CURRENT_BIT|GL_ENABLE_BIT|GL_LINE_BIT); + glDisable(GL_LIGHTING); + glLineWidth(1); + + glBegin(GL_LINES); + glColor3f(0.5,0.5,0.5); + for(int i=-spacing; i<=spacing; ++i){ + glVertex3f(-size,0,i*size/spacing); + glVertex3f((i!=0)*size,0,i*size/spacing); + glVertex3f(i*size/spacing,0,-size); + glVertex3f(i*size/spacing,0,(i!=0)*size); + } + glColor3f(1,0,0); + glVertex3f(0,0,0); + glVertex3f(size,0,0); + glColor3f(0,1,0); + glVertex3f(0,0,0); + glVertex3f(0,size,0); + glColor3f(0,0,1); + glVertex3f(0,0,0); + glVertex3f(0,0,size); + glEnd(); + + glPopAttrib(); +} + +void draw_text(const float point[3], const char *text, int fontsize) +{ + // please implement me! +} + +//================================================================================= + +void init(const char *windowtitle, int *argc, char **argv) +{ + glutInit(argc, argv); + glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGBA|GLUT_DEPTH); + glutInitWindowSize(winwidth, winheight); + glutCreateWindow(windowtitle); + glutReshapeFunc(gluviReshape); + glutDisplayFunc(gluviDisplay); + glutMouseFunc(gluviMouse); + glutMotionFunc(gluviDrag); + glEnable(GL_DEPTH_TEST); + glClearColor(0, 0, 0, 0); + glClearDepth(1); + glPixelStorei(GL_PACK_ALIGNMENT, 1); + glPixelStorei(GL_UNPACK_ALIGNMENT, 1); +} + +//================================================================================= + +void (*userDisplayFunc)(void)=0; +void (*userMouseFunc)(int button, int state, int x, int y)=0; +void (*userDragFunc)(int x, int y)=0; +Camera *camera=0; +WidgetList root(0); +int winwidth=720, winheight=480; + +//================================================================================= + +void run(void) +{ + glutMainLoop(); +} + +}; + diff --git a/viewpls3D/gluvi.h b/viewpls3D/gluvi.h new file mode 100644 index 0000000..c83cb9e --- /dev/null +++ b/viewpls3D/gluvi.h @@ -0,0 +1,200 @@ +#ifndef GLUVI_H +#define GLUVI_H + +#include +#include + +#ifdef __APPLE__ +#include // why does Apple have to put glut.h here... +#else +#include // ...when everyone else puts it here? +#endif + +#include "vec.h" + +namespace Gluvi{ + +/* +For portions of the code which export RIB commands (e.g. from cameras): + +The RenderMan (R) Interface Procedures and Protocol are: +Copyright 1988, 1989, Pixar +All Rights Reserved +*/ + +// the camera capability for Gluvi +struct Camera +{ + virtual ~Camera(void) {} + virtual void click(int button, int state, int x, int y) = 0; + virtual void drag(int x, int y) = 0; + virtual void return_to_default(void) = 0; + //@@@ add these to be called by a user glutKeyboardFunc() thing so that + //@@@ cameras can easily add keyboard shortcuts for e.g. return_to_default, transformation, etc. + //virtual void navigation_keyboard_handler(unsigned char key, int x, int y) = 0; + //virtual void navigation_special_key_handler(int key, int x, int y) = 0; + virtual void gl_transform(void) = 0; + virtual void export_rib(std::ostream &output) = 0; + virtual void display_screen(void) {} // in case the camera needs to show anything on screen +}; + +struct Target3D : public Camera +{ + float target[3], dist; + float heading, pitch; + float default_target[3], default_dist; + float default_heading, default_pitch; + float fovy; + float near_clip_factor, far_clip_factor; + enum {INACTIVE, ROTATE, TRUCK, DOLLY} action_mode; + int oldmousex, oldmousey; + + Target3D(float target_[3]=0, float dist_=1, float heading_=0, float pitch_=0, float fovy_=45, + float near_clip_factor_=0.01, float far_clip_factor_=100); + void click(int button, int state, int x, int y); + void drag(int x, int y); + void return_to_default(void); + void transform_mouse(int x, int y, float ray_origin[3], float ray_direction[3]); + void get_viewing_direction(float direction[3]); + void gl_transform(void); + void export_rib(std::ostream &output); +}; + +// same as above, but with orthographic projection +struct TargetOrtho3D : public Camera +{ + float target[3], dist; + float heading, pitch; + float default_target[3], default_dist; + float default_heading, default_pitch; + float height_factor; + float near_clip_factor, far_clip_factor; + enum {INACTIVE, ROTATE, TRUCK, DOLLY} action_mode; // @@@@ WHAT ABOUT ZOOMING??? IS WIDTH ALWAYS A FUNCTION OF DIST? + int oldmousex, oldmousey; + + TargetOrtho3D(float target_[3]=0, float dist_=1, float heading_=0, float pitch_=0, float height_factor_=1, + float near_clip_factor_=0.01, float far_clip_factor_=100); + void click(int button, int state, int x, int y); + void drag(int x, int y); + void return_to_default(void); + void transform_mouse(int x, int y, float ray_origin[3], float ray_direction[3]); + void get_viewing_direction(float direction[3]); + void gl_transform(void); + void export_rib(std::ostream &output); +}; + +struct PanZoom2D : public Camera +{ + float bottom, left, height; + float default_bottom, default_left, default_height; + enum {INACTIVE, PAN, ZOOM_IN, ZOOM_OUT} action_mode; + int oldmousex, oldmousey; + bool moved_since_mouse_down; // to distinuish simple clicks from drags + int clickx, clicky; + + PanZoom2D(float bottom_=0, float left_=0, float height_=1); + void click(int button, int state, int x, int y); + void drag(int x, int y); + void return_to_default(void); + void transform_mouse(int x, int y, float coords[2]); + void gl_transform(void); + void export_rib(std::ostream &output); + void display_screen(void); +}; + +// overlaid user-interface widgets +struct Widget +{ + int dispx, dispy, width, height; // set in display() + + virtual ~Widget() {} + virtual void display(int x, int y) = 0; + virtual bool click(int state, int x, int y) { return false; } // returns true if click handled by widget + virtual void drag(int x, int y) {} +}; + +struct StaticText : public Widget +{ + const char *text; + + StaticText(const char *text_); + void display(int x, int y); +}; + +struct Button : public Widget +{ + enum {UNINVOLVED, SELECTED, HIGHLIGHTED} status; + const char *text; + int minwidth; + + Button(const char *text_, int minwidth_=0); + void display(int x, int y); + bool click(int state, int x, int y); + void drag(int x, int y); + virtual void action() {} +}; + +struct Slider : public Widget +{ + enum {UNINVOLVED, SELECTED} status; + const char *text; + int length, justify; + int position; + int scrollxmin, scrollxmax, scrollymin, scrollymax; + int clickx; + + Slider(const char *text_, int length_=100, int position_=0, int justify_=0); + void display(int x, int y); + bool click(int state, int x, int y); + void drag(int x, int y); + virtual void action() {} +}; + +struct WidgetList : public Widget +{ + int indent; + bool hidden; + std::vector list; + int downclicked_member; + + WidgetList(int listindent_=12, bool hidden_=false); + void display(int x, int y); + bool click(int state, int x, int y); + void drag(int x, int y); +}; + +// display callback +extern void (*userDisplayFunc)(void); + +// mouse callbacks for events that Gluvi ignores (control not pressed, or mouse not on an active widget) +extern void (*userMouseFunc)(int button, int state, int x, int y); +extern void (*userDragFunc)(int x, int y); + +// user is free to do their own callbacks for everything else except glutReshape() + +// additional helpful functions +void ppm_screenshot(const char *filename_format, ...); +void sgi_screenshot(const char *filename_format, ...); +void set_generic_lights(void); +void set_generic_material(float r, float g, float b, GLenum face=GL_FRONT_AND_BACK); +void set_matte_material(float r, float g, float b, GLenum face=GL_FRONT_AND_BACK); +//@@@@@@@ USEFUL FUNCTIONALITY: +void draw_3d_arrow(const float base[3], const float point[3], float arrow_head_length=0); +void draw_2d_arrow(const Vec2f base, const Vec2f point, float arrow_head_length); +void draw_coordinate_grid(float size=1, int spacing=10); +void draw_text(const float point[3], const char *text, int fontsize=12); + +// call init first thing +void init(const char *windowtitle, int *argc, char **argv); + +// the Gluvi state +extern Camera *camera; +extern WidgetList root; +extern int winwidth, winheight; + +// then after setting the Gluvi state and doing any of your own set-up, call run() +void run(void); + +}; // end of namespace + +#endif diff --git a/viewpls3D/main.cpp b/viewpls3D/main.cpp new file mode 100644 index 0000000..8323821 --- /dev/null +++ b/viewpls3D/main.cpp @@ -0,0 +1,370 @@ +#include +#include +#include +#include +#include + +#include "gluvi.h" +#include "gl/glu.h" + +//Simple viewer for liquid simulator data +//Hold shift and use the mouse buttons to manipulate the camera + +using namespace std; + +string frame_number="frame 0"; + +string file_path; + +unsigned int frame= 0; +unsigned int highest_frame = 0; + +bool filming = false; +bool running = false; + +char * ppmfileformat; + +std::vector particles; +float particle_radius; + +bool read_frame(int newframe) +{ + if(newframe<0) newframe = highest_frame; + + std::ostringstream strout; + + strout << file_path << "particles_" << newframe << ".txt"; + printf("File path %s\n", strout.str().c_str()); + std::ifstream particles_in(strout.str().c_str()); + if(!particles_in.good()) { + printf("Failed to open particles!\n"); + return false; + } + unsigned int p_count; + particles_in >> p_count >> particle_radius; + particles.resize(p_count); + for(unsigned int p = 0; p < p_count; ++p) { + Vec3f pos; + particles_in >> pos[0] >> pos[1] >> pos[2]; + particles[p] = pos; + } + printf("Particle count: %d\n", p_count); + printf("Particle radius: %f\n", particle_radius); + + if(newframe > (int)highest_frame) + highest_frame = newframe; + + strout.str(""); + + frame=newframe; + + strout.str(""); + strout << "frame " << frame; + frame_number = strout.str(); + + return true; +} + +void set_view(Gluvi::Target3D &cam) +{ + cam.dist=0.5; +} + +void set_lights_and_material(int object) +{ + glEnable(GL_LIGHTING); + GLfloat global_ambient[4] = {0.1f, 0.1f, 0.1f, 1.0f}; + glLightModelfv(GL_LIGHT_MODEL_AMBIENT, global_ambient); + glShadeModel(GL_SMOOTH); + + //Light #1 + GLfloat color[4] = {1.0f, 1.0f, 1.0f, 1.0f}; + GLfloat position[3] = {1.0f, 1.0f, 1.0f}; + glLightfv(GL_LIGHT0, GL_SPECULAR, color); + glLightfv(GL_LIGHT0, GL_DIFFUSE, color); + glLightfv(GL_LIGHT0, GL_POSITION, position); + + //Light #2 + GLfloat color2[4] = {1.0f, 1.0f, 1.0f, 1.0f}; + GLfloat position2[3] = {-1.0f, -1.0f, 1.0f}; + glLightfv(GL_LIGHT1, GL_SPECULAR, color2); + glLightfv(GL_LIGHT1, GL_DIFFUSE, color2); + glLightfv(GL_LIGHT1, GL_POSITION, position2); + + GLfloat obj_color[4] = {.2f, .3f, .7f}; + glMaterialfv (GL_FRONT, GL_AMBIENT, obj_color); + glMaterialfv (GL_FRONT, GL_DIFFUSE, obj_color); + + GLfloat specular[4] = {.4f, .2f, .8f}; + glMaterialf (GL_FRONT, GL_SHININESS, 32); + glMaterialfv (GL_FRONT, GL_SPECULAR, specular); + glEnable(GL_LIGHT0); + glEnable(GL_LIGHT1); + +} + +void timer(int value) +{ + if(filming) { + Gluvi::ppm_screenshot(ppmfileformat, frame); + if(read_frame(frame+1)) { + if(frame == 0) { + filming = false; + } + } + glutPostRedisplay(); + glutTimerFunc(200, timer, 0); + } + + + if(running) { + read_frame(frame+1); + glutTimerFunc(1, timer, 0); + glutPostRedisplay(); + } + +} + + +void display(void) +{ + glClearColor(0.6f, 0.7f, 0.9f, 1); + + //Coordinate system + glDisable(GL_LIGHTING); + glBegin(GL_LINES); + glColor3f(1.0f,0,0); glVertex3f(0,0,0); glVertex3f(0.1f,0,0); + glColor3f(0.0f,1.0f,0); glVertex3f(0,0,0); glVertex3f(0,0.1f,0); + glColor3f(0.0f,0,1.0f); glVertex3f(0,0,0); glVertex3f(0,0,0.1f); + glEnd(); + + glEnable(GL_LIGHTING); + glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); + + set_lights_and_material(1); + + //Draw the liquid particles as simple spheres for now. + glPolygonMode(GL_FRONT_AND_BACK, GL_LINES); + GLUquadric* particle_sphere; + particle_sphere = gluNewQuadric(); + gluQuadricDrawStyle(particle_sphere, GLU_FILL ); + for(unsigned int p = 0; p < particles.size(); ++p) { + glPushMatrix(); + Vec3f pos = particles[p].v; + glTranslatef(pos[0], pos[1], pos[2]); + gluSphere(particle_sphere, particle_radius, 20, 20); + glPopMatrix(); + } + + //Draw the bound box for good measure + glDisable(GL_LIGHTING); + glColor3f(0,0,0); + glBegin(GL_LINES); + glVertex3f(0,0,0); + glVertex3f(0,0,1); + + glVertex3f(0,0,0); + glVertex3f(0,1,0); + + glVertex3f(0,0,0); + glVertex3f(1,0,0); + + glVertex3f(0,1,0); + glVertex3f(1,1,0); + + glVertex3f(1,1,0); + glVertex3f(1,0,0); + + glVertex3f(1,0,0); + glVertex3f(1,0,1); + + glVertex3f(0,1,0); + glVertex3f(0,1,1); + + glVertex3f(1,1,0); + glVertex3f(1,1,1); + + glVertex3f(0,1,1); + glVertex3f(1,1,1); + + glVertex3f(1,0,1); + glVertex3f(1,1,1); + + glVertex3f(0,0,1); + glVertex3f(1,0,1); + + glVertex3f(0,0,1); + glVertex3f(0,1,1); + + glEnd(); + + //Draw wireframe sphere geometry (specific to this scene). + glColor3f(0,0,0); + glPolygonMode(GL_FRONT_AND_BACK, GL_LINES); + GLUquadric* sphere; + sphere = gluNewQuadric(); + gluQuadricDrawStyle(sphere, GLU_LINE ); + glPushMatrix(); + glTranslatef(0.5f, 0.5f,0.5f); + gluSphere(sphere, 0.35, 20, 20); + glPopMatrix(); + + +} + +struct ScreenShotButton : public Gluvi::Button{ + const char *filename_format; + ScreenShotButton(const char *label, const char *filename_format_) : Gluvi::Button(label), filename_format(filename_format_) {} + void action() + { + Gluvi::ppm_screenshot(filename_format, frame); + } +}; + + +struct MovieButton : public Gluvi::Button{ + const char *filename_format; + MovieButton(const char *label, const char *filename_format_) : Gluvi::Button(label), filename_format(filename_format_) {} + void action() + { + if(!running) { + if(!filming) { + filming = true; + glutTimerFunc(1000, timer, 0); + } + else { + filming = false; + } + } + } +}; + +struct RunButton : public Gluvi::Button{ + RunButton(const char *label) : Gluvi::Button(label){} + void action() + { + if(!filming) { + if(!running) { + running = true; + } + else { + running = false; + } + } + } +}; + +void keyPress(unsigned char key, int x, int y) { + + if(key == 'r') { + if(!filming) { + if(!running) { + running = true; + glutTimerFunc(1000, timer, 0); + } + else { + running = false; + } + } + } + else if(key == 'f') { + if(!running) { + if(!filming) { + filming = true; + } + else { + filming = false; + } + } + } + glutPostRedisplay(); + + + +} + +void special_key_handler(int key, int x, int y) +{ + int mods=glutGetModifiers(); + switch(key){ +case GLUT_KEY_LEFT: + if(mods == GLUT_ACTIVE_SHIFT) { + if(read_frame(0)) + glutPostRedisplay(); + } + else { + if(read_frame(frame-1)) + glutPostRedisplay(); + } + break; +case GLUT_KEY_RIGHT: + if(mods == GLUT_ACTIVE_SHIFT) { + if(read_frame(highest_frame)) + glutPostRedisplay(); + } + else { + if(read_frame(frame+1)) + glutPostRedisplay(); + } + break; +default: + ; + } +} + + +Gluvi::Target3D* cam_local; +int main(int argc, char **argv) +{ + + Gluvi::init("Liquid Data Viewer", &argc, argv); + + if(argc!=2){ + cerr<<"Expecting path to simulation data folder"< +#include +#include +#include + +#ifndef M_PI +const double M_PI = 3.1415926535897932384626433832795; +#endif + +#ifdef WIN32 +#undef min +#undef max +#endif + +using std::min; +using std::max; +using std::swap; + +template +inline T sqr(const T& x) +{ return x*x; } + +template +inline T cube(const T& x) +{ return x*x*x; } + +template +inline T min(T a1, T a2, T a3) +{ return min(a1, min(a2, a3)); } + +template +inline T min(T a1, T a2, T a3, T a4) +{ return min(min(a1, a2), min(a3, a4)); } + +template +inline T min(T a1, T a2, T a3, T a4, T a5) +{ return min(min(a1, a2), min(a3, a4), a5); } + +template +inline T min(T a1, T a2, T a3, T a4, T a5, T a6) +{ return min(min(a1, a2), min(a3, a4), min(a5, a6)); } + +template +inline T max(T a1, T a2, T a3) +{ return max(a1, max(a2, a3)); } + +template +inline T max(T a1, T a2, T a3, T a4) +{ return max(max(a1, a2), max(a3, a4)); } + +template +inline T max(T a1, T a2, T a3, T a4, T a5) +{ return max(max(a1, a2), max(a3, a4), a5); } + +template +inline T max(T a1, T a2, T a3, T a4, T a5, T a6) +{ return max(max(a1, a2), max(a3, a4), max(a5, a6)); } + +template +inline void minmax(T a1, T a2, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T a4, T& amin, T& amax) +{ + if(a1 +inline void minmax(T a1, T a2, T a3, T a4, T a5, T& amin, T& amax) +{ + //@@@ the logic could be shortcircuited a lot! + amin=min(a1,a2,a3,a4,a5); + amax=max(a1,a2,a3,a4,a5); +} + +template +inline void minmax(T a1, T a2, T a3, T a4, T a5, T a6, T& amin, T& amax) +{ + //@@@ the logic could be shortcircuited a lot! + amin=min(a1,a2,a3,a4,a5,a6); + amax=max(a1,a2,a3,a4,a5,a6); +} + +template +inline void update_minmax(T a1, T& amin, T& amax) +{ + if(a1amax) amax=a1; +} + +template +inline void sort(T &a, T &b, T &c) +{ + T temp; + if(a +inline T clamp(T a, T lower, T upper) +{ + if(aupper) return upper; + else return a; +} + +// only makes sense with T=float or double +template +inline T smooth_step(T r) +{ + if(r<0) return 0; + else if(r>1) return 1; + return r*r*r*(10+r*(-15+r*6)); +} + +// only makes sense with T=float or double +template +inline T smooth_step(T r, T r_lower, T r_upper, T value_lower, T value_upper) +{ return value_lower + smooth_step((r-r_lower)/(r_upper-r_lower)) * (value_upper-value_lower); } + +// only makes sense with T=float or double +template +inline T ramp(T r) +{ return smooth_step((r+1)/2)*2-1; } + +#ifdef WIN32 +inline int lround(double x) +{ + if(x>0) + return (x-floor(x)<0.5) ? (int)floor(x) : (int)ceil(x); + else + return (x-floor(x)<=0.5) ? (int)floor(x) : (int)ceil(x); +} + +inline double remainder(double x, double y) +{ + return x-std::floor(x/y+0.5)*y; +} +#endif + +inline unsigned int round_up_to_power_of_two(unsigned int n) +{ + int exponent=0; + --n; + while(n){ + ++exponent; + n>>=1; + } + return 1<1){ + ++exponent; + n>>=1; + } + return 1<>16); + i*=2654435769u; + i^=(i>>16); + i*=2654435769u; + return i; +} + +// the inverse of randhash +inline unsigned int unhash(unsigned int h) +{ + h*=340573321u; + h^=(h>>16); + h*=340573321u; + h^=(h>>16); + h*=340573321u; + h^=0xA3C59AC3u; + return h; +} + +// returns repeatable stateless pseudo-random number in [0,1] +inline double randhashd(unsigned int seed) +{ return randhash(seed)/(double)UINT_MAX; } +inline float randhashf(unsigned int seed) +{ return randhash(seed)/(float)UINT_MAX; } + +// returns repeatable stateless pseudo-random number in [a,b] +inline double randhashd(unsigned int seed, double a, double b) +{ return (b-a)*randhash(seed)/(double)UINT_MAX + a; } +inline float randhashf(unsigned int seed, float a, float b) +{ return ( (b-a)*randhash(seed)/(float)UINT_MAX + a); } + +inline int intlog2(int x) +{ + int exp=-1; + while(x){ + x>>=1; + ++exp; + } + return exp; +} + +template +inline void get_barycentric(T x, int& i, T& f, int i_low, int i_high) +{ + T s=std::floor(x); + i=(int)s; + if(ii_high-2){ + i=i_high-2; + f=1; + }else + f=(T)(x-s); +} + +template +inline S lerp(const S& value0, const S& value1, T f) +{ return (1-f)*value0 + f*value1; } + +template +inline S bilerp(const S& v00, const S& v10, + const S& v01, const S& v11, + T fx, T fy) +{ + return lerp(lerp(v00, v10, fx), + lerp(v01, v11, fx), + fy); +} + +template +inline S trilerp(const S& v000, const S& v100, + const S& v010, const S& v110, + const S& v001, const S& v101, + const S& v011, const S& v111, + T fx, T fy, T fz) +{ + return lerp(bilerp(v000, v100, v010, v110, fx, fy), + bilerp(v001, v101, v011, v111, fx, fy), + fz); +} + +template +inline S quadlerp(const S& v0000, const S& v1000, + const S& v0100, const S& v1100, + const S& v0010, const S& v1010, + const S& v0110, const S& v1110, + const S& v0001, const S& v1001, + const S& v0101, const S& v1101, + const S& v0011, const S& v1011, + const S& v0111, const S& v1111, + T fx, T fy, T fz, T ft) +{ + return lerp(trilerp(v0000, v1000, v0100, v1100, v0010, v1010, v0110, v1110, fx, fy, fz), + trilerp(v0001, v1001, v0101, v1101, v0011, v1011, v0111, v1111, fx, fy, fz), + ft); +} + +// f should be between 0 and 1, with f=0.5 corresponding to balanced weighting between w0 and w2 +template +inline void quadratic_bspline_weights(T f, T& w0, T& w1, T& w2) +{ + w0=T(0.5)*sqr(f-1); + w1=T(0.75)-sqr(f-T(0.5));; + w2=T(0.5)*sqr(f); +} + +// f should be between 0 and 1 +template +inline void cubic_interp_weights(T f, T& wneg1, T& w0, T& w1, T& w2) +{ + T f2(f*f), f3(f2*f); + wneg1=-T(1./3)*f+T(1./2)*f2-T(1./6)*f3; + w0=1-f2+T(1./2)*(f3-f); + w1=f+T(1./2)*(f2-f3); + w2=T(1./6)*(f3-f); +} + +template +inline S cubic_interp(const S& value_neg1, const S& value0, const S& value1, const S& value2, T f) +{ + T wneg1, w0, w1, w2; + cubic_interp_weights(f, wneg1, w0, w1, w2); + return wneg1*value_neg1 + w0*value0 + w1*value1 + w2*value2; +} + +template +void zero(std::vector& v) +{ for(int i=(int)v.size()-1; i>=0; --i) v[i]=0; } + +template +T abs_max(const std::vector& v) +{ + T m=0; + for(int i=(int)v.size()-1; i>=0; --i){ + if(std::fabs(v[i])>m) + m=std::fabs(v[i]); + } + return m; +} + +template +bool contains(const std::vector& a, T e) +{ + for(unsigned int i=0; i +void add_unique(std::vector& a, T e) +{ + for(unsigned int i=0; i +void insert(std::vector& a, unsigned int index, T e) +{ + a.push_back(a.back()); + for(unsigned int i=(unsigned int)a.size()-1; i>index; --i) + a[i]=a[i-1]; + a[index]=e; +} + +template +void erase(std::vector& a, unsigned int index) +{ + for(unsigned int i=index; i +void erase_swap(std::vector& a, unsigned int index) +{ + for(unsigned int i=index; i +void erase_unordered(std::vector& a, unsigned int index) +{ + a[index]=a.back(); + a.pop_back(); +} + +template +void erase_unordered_swap(std::vector& a, unsigned int index) +{ + swap(a[index], a.back()); + a.pop_back(); +} + +template +void find_and_erase_unordered(std::vector& a, const T& doomed_element) +{ + for(unsigned int i=0; i +void replace_once(std::vector& a, const T& old_element, const T& new_element) +{ + for(unsigned int i=0; i +void write_matlab(std::ostream& output, const std::vector& a, const char *variable_name, bool column_vector=true, int significant_digits=18) +{ + output< +#include +#include +#include "util.h" + +// Defines a thin wrapper around fixed size C-style arrays, using template parameters, +// which is useful for dealing with vectors of different dimensions. +// For example, float[3] is equivalent to Vec<3,float>. +// Entries in the vector are accessed with the overloaded [] operator, so +// for example if x is a Vec<3,float>, then the middle entry is x[1]. +// For convenience, there are a number of typedefs for abbreviation: +// Vec<3,float> -> Vec3f +// Vec<2,int> -> Vec2i +// and so on. +// Arithmetic operators are appropriately overloaded, and functions are defined +// for additional operations (such as dot-products, norms, cross-products, etc.) + +template +struct Vec +{ + T v[N]; + + Vec(void) + {} + + Vec(T value_for_all) + { for(unsigned int i=0; i + Vec(const S *source) + { for(unsigned int i=0; i + explicit Vec(const Vec& source) + { for(unsigned int i=0; i(T v0, T v1) + { + assert(N==2); + v[0]=v0; v[1]=v1; + } + + Vec(T v0, T v1, T v2) + { + assert(N==3); + v[0]=v0; v[1]=v1; v[2]=v2; + } + + Vec(T v0, T v1, T v2, T v3) + { + assert(N==4); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; + } + + Vec(T v0, T v1, T v2, T v3, T v4) + { + assert(N==5); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; + } + + Vec(T v0, T v1, T v2, T v3, T v4, T v5) + { + assert(N==6); + v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; v[5]=v5; + } + + T &operator[](int index) + { + assert(0<=index && (unsigned int)index operator+=(const Vec &w) + { + for(unsigned int i=0; i operator+(const Vec &w) const + { + Vec sum(*this); + sum+=w; + return sum; + } + + Vec operator-=(const Vec &w) + { + for(unsigned int i=0; i operator-(void) const // unary minus + { + Vec negative; + for(unsigned int i=0; i operator-(const Vec &w) const // (binary) subtraction + { + Vec diff(*this); + diff-=w; + return diff; + } + + Vec operator*=(T a) + { + for(unsigned int i=0; i operator*(T a) const + { + Vec w(*this); + w*=a; + return w; + } + + Vec operator*(const Vec &w) const + { + Vec componentwise_product; + for(unsigned int i=0; i operator/=(T a) + { + for(unsigned int i=0; i operator/(T a) const + { + Vec w(*this); + w/=a; + return w; + } +}; + +typedef Vec<2,double> Vec2d; +typedef Vec<2,float> Vec2f; +typedef Vec<2,int> Vec2i; +typedef Vec<2,unsigned int> Vec2ui; +typedef Vec<2,short> Vec2s; +typedef Vec<2,unsigned short> Vec2us; +typedef Vec<2,char> Vec2c; +typedef Vec<2,unsigned char> Vec2uc; + +typedef Vec<3,double> Vec3d; +typedef Vec<3,float> Vec3f; +typedef Vec<3,int> Vec3i; +typedef Vec<3,unsigned int> Vec3ui; +typedef Vec<3,short> Vec3s; +typedef Vec<3,unsigned short> Vec3us; +typedef Vec<3,char> Vec3c; +typedef Vec<3,unsigned char> Vec3uc; + +typedef Vec<4,double> Vec4d; +typedef Vec<4,float> Vec4f; +typedef Vec<4,int> Vec4i; +typedef Vec<4,unsigned int> Vec4ui; +typedef Vec<4,short> Vec4s; +typedef Vec<4,unsigned short> Vec4us; +typedef Vec<4,char> Vec4c; +typedef Vec<4,unsigned char> Vec4uc; + +typedef Vec<6,double> Vec6d; +typedef Vec<6,float> Vec6f; +typedef Vec<6,unsigned int> Vec6ui; +typedef Vec<6,int> Vec6i; +typedef Vec<6,short> Vec6s; +typedef Vec<6,unsigned short> Vec6us; +typedef Vec<6,char> Vec6c; +typedef Vec<6,unsigned char> Vec6uc; + + +template +T mag2(const Vec &a) +{ + T l=sqr(a.v[0]); + for(unsigned int i=1; i +T mag(const Vec &a) +{ return sqrt(mag2(a)); } + +template +inline T dist2(const Vec &a, const Vec &b) +{ + T d=sqr(a.v[0]-b.v[0]); + for(unsigned int i=1; i +inline T dist(const Vec &a, const Vec &b) +{ return std::sqrt(dist2(a,b)); } + +template +inline void normalize(Vec &a) +{ a/=mag(a); } + +template +inline Vec normalized(const Vec &a) +{ return a/mag(a); } + +template +inline T infnorm(const Vec &a) +{ + T d=std::fabs(a.v[0]); + for(unsigned int i=1; i +void zero(Vec &a) +{ + for(unsigned int i=0; i +std::ostream &operator<<(std::ostream &out, const Vec &v) +{ + out< +std::istream &operator>>(std::istream &in, Vec &v) +{ + in>>v.v[0]; + for(unsigned int i=1; i>v.v[i]; + return in; +} + +template +inline bool operator==(const Vec &a, const Vec &b) +{ + bool t = (a.v[0] == b.v[0]); + unsigned int i=1; + while(i +inline bool operator!=(const Vec &a, const Vec &b) +{ + bool t = (a.v[0] != b.v[0]); + unsigned int i=1; + while(i +inline Vec operator*(T a, const Vec &v) +{ + Vec w(v); + w*=a; + return w; +} + +template +inline T min(const Vec &a) +{ + T m=a.v[0]; + for(unsigned int i=1; i +inline Vec min_union(const Vec &a, const Vec &b) +{ + Vec m; + for(unsigned int i=0; i +inline Vec max_union(const Vec &a, const Vec &b) +{ + Vec m; + for(unsigned int i=0; i b.v[i]) ? m.v[i]=a.v[i] : m.v[i]=b.v[i]; + return m; +} + +template +inline T max(const Vec &a) +{ + T m=a.v[0]; + for(unsigned int i=1; im) m=a.v[i]; + return m; +} + +template +inline T dot(const Vec &a, const Vec &b) +{ + T d=a.v[0]*b.v[0]; + for(unsigned int i=1; i +inline Vec<2,T> rotate(const Vec<2,T>& a, float angle) +{ + T c = cos(angle); + T s = sin(angle); + return Vec<2,T>(c*a[0] - s*a[1],s*a[0] + c*a[1]); // counter-clockwise rotation +} + +template +inline Vec<2,T> perp(const Vec<2,T> &a) +{ return Vec<2,T>(-a.v[1], a.v[0]); } // counter-clockwise rotation by 90 degrees + +template +inline T cross(const Vec<2,T> &a, const Vec<2,T> &b) +{ return a.v[0]*b.v[1]-a.v[1]*b.v[0]; } + +template +inline Vec<3,T> cross(const Vec<3,T> &a, const Vec<3,T> &b) +{ return Vec<3,T>(a.v[1]*b.v[2]-a.v[2]*b.v[1], a.v[2]*b.v[0]-a.v[0]*b.v[2], a.v[0]*b.v[1]-a.v[1]*b.v[0]); } + +template +inline T triple(const Vec<3,T> &a, const Vec<3,T> &b, const Vec<3,T> &c) +{ return a.v[0]*(b.v[1]*c.v[2]-b.v[2]*c.v[1]) + +a.v[1]*(b.v[2]*c.v[0]-b.v[0]*c.v[2]) + +a.v[2]*(b.v[0]*c.v[1]-b.v[1]*c.v[0]); } + +template +inline unsigned int hash(const Vec &a) +{ + unsigned int h=a.v[0]; + for(unsigned int i=1; i +inline void assign(const Vec &a, T &a0, T &a1) +{ + assert(N==2); + a0=a.v[0]; a1=a.v[1]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2) +{ + assert(N==3); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2, T &a3) +{ + assert(N==4); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3]; +} + +template +inline void assign(const Vec &a, T &a0, T &a1, T &a2, T &a3, T &a4, T &a5) +{ + assert(N==6); + a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3]; a4=a.v[4]; a5=a.v[5]; +} + +template +inline Vec round(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec floor(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec ceil(const Vec &a) +{ + Vec rounded; + for(unsigned int i=0; i +inline Vec fabs(const Vec &a) +{ + Vec result; + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, + Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, const Vec &x4, + Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void minmax(const Vec &x0, const Vec &x1, const Vec &x2, const Vec &x3, const Vec &x4, + const Vec &x5, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +inline void update_minmax(const Vec &x, Vec &xmin, Vec &xmax) +{ + for(unsigned int i=0; i +#include + +int u_ind(int i, int j, int k, int nx, int ny) { + return i + j*(nx+1) + k*(nx+1)*ny; +} + +int v_ind(int i, int j, int k, int nx, int ny, int nz){ + return i + j*nx + k*nx*(ny+1) + (nx+1)*ny*nz ; +} + +int w_ind(int i, int j, int k, int nx, int ny, int nz){ + return i + j*nx + k*nx*ny + (nx+1)*ny*nz + nx*(ny+1)*nz; +} + +Array3c u_state;//(nx+1,ny,nz,0); +Array3c v_state;//(nx,ny+1,nz,0); +Array3c w_state;//(nx,ny,nz+1,0); + + +SparseMatrixd matrix;//(dim,dim,15); +SparseMatrixd matrix2;//(dim,dim +std::vector rhs;//(dim); +std::vector soln;//(dim); + +void advance_viscosity_implicit_weighted(Array3f& u, Array3f& v, Array3f& w, + const Array3f& vol_u, const Array3f& vol_v, const Array3f& vol_w, + const Array3f& vol_c, const Array3f& vol_ex, const Array3f& vol_ey, const Array3f& vol_ez, + const Array3f& solid_phi, + const Array3f& viscosity, float dt, float dx) { + float over_dx = 1.0f/dx; + int nx = solid_phi.ni; + int ny = solid_phi.nj; + int nz = solid_phi.nk; + printf("Creating state arrays.\n"); + std::cout << "Phi-size:" << nx << " " << ny << " " << nz << std::endl; + int dim = (nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1); + if(u_state.ni != u.ni) { + printf("Creating matrices and vectors.\n"); + u_state.resize(nx+1,ny,nz); + v_state.resize(nx,ny+1,nz); + w_state.resize(nx,ny,nz+1); + matrix.resize(dim); + matrix2.resize(dim); + rhs.resize(dim); + soln.resize(dim); + printf("Done that for good."); + } + + u_state.assign(0); + v_state.assign(0); + w_state.assign(0); + matrix.zero(); + //matrix2.zero(); + rhs.assign(dim, 0); + soln.assign(dim, 0); + + const int SOLID = 3; + const int FLUID = 2; + //const int AIR = 1; + + //check if interpolated velocity positions are inside solid + for(int k = 0; k < nz; ++k) for(int j = 0; j < ny; ++j) for(int i = 0; i < nx+1; ++i) { + if(i - 1 < 0 || i >= nx || solid_phi(i-1,j,k) + solid_phi(i,j,k) <= 0) + u_state(i,j,k) = SOLID; + else + u_state(i,j,k) = FLUID; + } + + for(int k = 0; k < nz; ++k) for(int j = 0; j < ny+1; ++j) for(int i = 0; i < nx; ++i) { + if(j - 1 < 0 || j >= ny || solid_phi(i,j-1,k) + solid_phi(i,j,k) <= 0) + v_state(i,j,k) = SOLID; + else + v_state(i,j,k) = FLUID; + } + + for(int k = 0; k < nz+1; ++k) for(int j = 0; j < ny; ++j) for(int i = 0; i < nx; ++i) { + if(k - 1 < 0 || k >= nz || solid_phi(i,j,k-1) + solid_phi(i,j,k) <= 0) + w_state(i,j,k) = SOLID; + else + w_state(i,j,k) = FLUID; + } + + float factor = dt*sqr(over_dx); + //u-terms + //2u_xx+ v_xy +uyy + u_zz + w_xz + printf("Building u-components.\n"); + for(int k = 1; k < nz; ++k) for(int j = 1; j < ny; ++j) for(int i = 1; i < nx; ++i) { + + if(u_state(i,j,k) == FLUID) { + int index = u_ind(i,j,k,nx,ny); + + rhs[index] = vol_u(i,j,k)*u(i,j,k); + matrix.set_element(index,index,vol_u(i,j,k)); + + float visc_right = viscosity(i,j,k); + float visc_left = viscosity(i-1,j,k); + float vol_right = vol_c(i,j,k); + float vol_left = vol_c(i-1,j,k); + + float visc_top = 0.25f*(viscosity(i-1,j+1,k) + viscosity(i-1,j,k) + viscosity(i,j+1,k) + viscosity(i,j,k)); + float visc_bottom = 0.25f*(viscosity(i-1,j,k) + viscosity(i-1,j-1,k) + viscosity(i,j,k) + viscosity(i,j-1,k)); + float vol_top = vol_ez(i,j+1,k); + float vol_bottom = vol_ez(i,j,k); + + float visc_front = 0.25f*(viscosity(i-1,j,k+1) + viscosity(i-1,j,k) + viscosity(i,j,k+1) + viscosity(i,j,k)); + float visc_back = 0.25f*(viscosity(i-1,j,k) + viscosity(i-1,j,k-1) + viscosity(i,j,k) + viscosity(i,j,k-1)); + float vol_front = vol_ey(i,j,k+1); + float vol_back = vol_ey(i,j,k); + + //u_x_right + matrix.add_to_element(index,index, 2*factor*visc_right*vol_right); + if(u_state(i+1,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i+1,j,k,nx,ny), -2*factor*visc_right*vol_right); + else if(u_state(i+1,j,k) == SOLID) + rhs[index] -= -2*factor*visc_right*vol_right*u(i+1,j,k); + + //u_x_left + matrix.add_to_element(index,index, 2*factor*visc_left*vol_left); + if(u_state(i-1,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i-1,j,k,nx,ny), -2*factor*visc_left*vol_left); + else if(u_state(i-1,j,k) == SOLID) + rhs[index] -= -2*factor*visc_left*vol_left*u(i-1,j,k); + + //u_y_top + matrix.add_to_element(index,index, +factor*visc_top*vol_top); + if(u_state(i,j+1,k) == FLUID) + matrix.add_to_element(index,u_ind(i,j+1,k,nx,ny), -factor*visc_top*vol_top); + else if(u_state(i,j+1,k) == SOLID) + rhs[index] -= -u(i,j+1,k)*factor*visc_top*vol_top; + + //u_y_bottom + matrix.add_to_element(index,index, +factor*visc_bottom*vol_bottom); + if(u_state(i,j-1,k) == FLUID) + matrix.add_to_element(index,u_ind(i,j-1,k,nx,ny), -factor*visc_bottom*vol_bottom); + else if(u_state(i,j-1,k) == SOLID) + rhs[index] -= -u(i,j-1,k)*factor*visc_bottom*vol_bottom; + + //u_z_front + matrix.add_to_element(index,index, +factor*visc_front*vol_front); + if(u_state(i,j,k+1) == FLUID) + matrix.add_to_element(index,u_ind(i,j,k+1,nx,ny), -factor*visc_front*vol_front); + else if(u_state(i,j,k+1) == SOLID) + rhs[index] -= -u(i,j,k+1)*factor*visc_front*vol_front; + + //u_z_back + matrix.add_to_element(index,index, +factor*visc_back*vol_back); + if(u_state(i,j,k-1) == FLUID) + matrix.add_to_element(index,u_ind(i,j,k-1,nx,ny), -factor*visc_back*vol_back); + else if(u_state(i,j,k-1) == SOLID) + rhs[index] -= -u(i,j,k-1)*factor*visc_back*vol_back; + + //v_x_top + if(v_state(i,j+1,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j+1,k,nx,ny,nz), -factor*visc_top*vol_top); + else if(v_state(i,j+1,k) == SOLID) + rhs[index] -= -v(i,j+1,k)*factor*visc_top*vol_top; + + if(v_state(i-1,j+1,k) == FLUID) + matrix.add_to_element(index,v_ind(i-1,j+1,k,nx,ny,nz), factor*visc_top*vol_top); + else if(v_state(i-1,j+1,k) == SOLID) + rhs[index] -= v(i-1,j+1,k)*factor*visc_top*vol_top; + + //v_x_bottom + if(v_state(i,j,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j,k,nx,ny,nz), +factor*visc_bottom*vol_bottom); + else if(v_state(i,j,k) == SOLID) + rhs[index] -= v(i,j,k)*factor*visc_bottom*vol_bottom; + + if(v_state(i-1,j,k) == FLUID) + matrix.add_to_element(index,v_ind(i-1,j,k,nx,ny,nz), -factor*visc_bottom*vol_bottom); + else if(v_state(i-1,j,k) == SOLID) + rhs[index] -= -v(i-1,j,k)*factor*visc_bottom*vol_bottom; + + //w_x_front + if(w_state(i,j,k+1) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k+1,nx,ny,nz), -factor*visc_front*vol_front); + else if(w_state(i,j,k+1) == SOLID) + rhs[index] -= -w(i,j,k+1)*factor*visc_front*vol_front; + + if(w_state(i-1,j,k+1) == FLUID) + matrix.add_to_element(index,w_ind(i-1,j,k+1,nx,ny,nz), factor*visc_front*vol_front); + else if(w_state(i-1,j,k+1) == SOLID) + rhs[index] -= w(i-1,j,k+1)*factor*visc_front*vol_front; + + //w_x_back + if(w_state(i,j,k) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k,nx,ny,nz), +factor*visc_back*vol_back); + else if(w_state(i,j,k) == SOLID) + rhs[index] -= w(i,j,k)*factor*visc_back*vol_back; + + if(w_state(i-1,j,k) == FLUID) + matrix.add_to_element(index,w_ind(i-1,j,k,nx,ny,nz), -factor*visc_back*vol_back); + else if(w_state(i-1,j,k) == SOLID) + rhs[index] -= -w(i-1,j,k)*factor*visc_back*vol_back; + } + } + + //v-terms + //vxx + 2vyy + vzz + u_yx + w_yz + printf("Building v-components.\n"); + for(int k = 1; k < nz; ++k) for(int j = 1; j < ny; ++j) for(int i = 1; i < nx; ++i) { + if(v_state(i,j,k) == FLUID) { + int index = v_ind(i,j,k,nx,ny,nz); + + rhs[index] = vol_v(i,j,k)*v(i,j,k); + matrix.set_element(index, index, vol_v(i,j,k)); + + float visc_right = 0.25f*(viscosity(i,j-1,k) + viscosity(i+1,j-1,k) + viscosity(i,j,k) + viscosity(i+1,j,k)); + float visc_left = 0.25f*(viscosity(i,j-1,k) + viscosity(i-1,j-1,k) + viscosity(i,j,k) + viscosity(i-1,j,k)); + float vol_right = vol_ez(i+1,j,k); + float vol_left = vol_ez(i,j,k); + + float visc_top = viscosity(i,j,k); + float visc_bottom = viscosity(i,j-1,k); + float vol_top = vol_c(i,j,k); + float vol_bottom = vol_c(i,j-1,k); + + float visc_front = 0.25f*(viscosity(i,j-1,k) + viscosity(i,j-1,k+1) + viscosity(i,j,k) + viscosity(i,j,k+1)); + float visc_back = 0.25f*(viscosity(i,j-1,k) + viscosity(i,j-1,k-1) + viscosity(i,j,k) + viscosity(i,j,k-1)); + float vol_front = vol_ex(i,j,k+1); + float vol_back = vol_ex(i,j,k); + + //v_x_right + matrix.add_to_element(index,index, +factor*visc_right*vol_right); + if(v_state(i+1,j,k) == FLUID) + matrix.add_to_element(index,v_ind(i+1,j,k,nx,ny,nz), -factor*visc_right*vol_right); + else if(v_state(i+1,j,k) == SOLID) + rhs[index] -= -v(i+1,j,k)*factor*visc_right*vol_right; + + //v_x_left + matrix.add_to_element(index,index, +factor*visc_left*vol_left); + if(v_state(i-1,j,k) == FLUID) + matrix.add_to_element(index,v_ind(i-1,j,k,nx,ny,nz), -factor*visc_left*vol_left); + else if(v_state(i-1,j,k) == SOLID) + rhs[index] -= -v(i-1,j,k)*factor*visc_left*vol_left; + + //vy_top + matrix.add_to_element(index,index, +2*factor*visc_top*vol_top); + if(v_state(i,j+1,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j+1,k,nx,ny,nz), -2*factor*visc_top*vol_top); + else if (v_state(i,j+1,k) == SOLID) + rhs[index] -= -2*factor*visc_top*vol_top*v(i,j+1,k); + + //vy_bottom + matrix.add_to_element(index,index, +2*factor*visc_bottom*vol_bottom); + if(v_state(i,j-1,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j-1,k,nx,ny,nz), -2*factor*visc_bottom*vol_bottom); + else if(v_state(i,j-1,k) == SOLID) + rhs[index] -= -2*factor*visc_bottom*vol_bottom*v(i,j-1,k); + + //v_z_front + matrix.add_to_element(index,index, +factor*visc_front*vol_front); + if(v_state(i,j,k+1) == FLUID) + matrix.add_to_element(index,v_ind(i,j,k+1,nx,ny,nz), -factor*visc_front*vol_front); + else if(v_state(i+1,j,k) == SOLID) + rhs[index] -= -v(i,j,k+1)*factor*visc_front*vol_front; + + //v_z_back + matrix.add_to_element(index,index, +factor*visc_back*vol_back); + if(v_state(i,j,k-1) == FLUID) + matrix.add_to_element(index,v_ind(i,j,k-1,nx,ny,nz), -factor*visc_back*vol_back); + else if(v_state(i,j,k-1) == SOLID) + rhs[index] -= -v(i,j,k-1)*factor*visc_back*vol_back; + + //u_y_right + if(u_state(i+1,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i+1,j,k,nx,ny), -factor*visc_right*vol_right); + else if(u_state(i+1,j,k) == SOLID) + rhs[index] -= -u(i+1,j,k)*factor*visc_right*vol_right; + + if(u_state(i+1,j-1,k) == FLUID) + matrix.add_to_element(index,u_ind(i+1,j-1,k,nx,ny), factor*visc_right*vol_right); + else if(u_state(i+1,j-1,k) == SOLID) + rhs[index] -= u(i+1,j-1,k)*factor*visc_right*vol_right; + + //u_y_left + if(u_state(i,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i,j,k,nx,ny), factor*visc_left*vol_left); + else if(u_state(i,j,k) == SOLID) + rhs[index] -= u(i,j,k)*factor*visc_left*vol_left; + + if(u_state(i,j-1,k) == FLUID) + matrix.add_to_element(index,u_ind(i,j-1,k,nx,ny), -factor*visc_left*vol_left); + else if(u_state(i,j-1,k) == SOLID) + rhs[index] -= -u(i,j-1,k)*factor*visc_left*vol_left; + + //w_y_front + if(w_state(i,j,k+1) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k+1,nx,ny,nz), -factor*visc_front*vol_front); + else if(w_state(i,j,k+1) == SOLID) + rhs[index] -= -w(i,j,k+1)*factor*visc_front*vol_front; + + if(w_state(i,j-1,k+1) == FLUID) + matrix.add_to_element(index,w_ind(i,j-1,k+1,nx,ny,nz), factor*visc_front*vol_front); + else if(w_state(i,j-1,k+1) == SOLID) + rhs[index] -= w(i,j-1,k+1)*factor*visc_front*vol_front; + + //w_y_back + if(w_state(i,j,k) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k,nx,ny,nz), factor*visc_back*vol_back); + else if(w_state(i,j,k) == SOLID) + rhs[index] -= w(i,j,k)*factor*visc_back*vol_back; + + if(w_state(i,j-1,k) == FLUID) + matrix.add_to_element(index,w_ind(i,j-1,k,nx,ny,nz), -factor*visc_back*vol_back); + else if(w_state(i,j-1,k) == SOLID) + rhs[index] -= -w(i,j-1,k)*factor*visc_back*vol_back; + } + } + + //w-terms + //wxx+ wyy+ 2wzz + u_zx + v_zy + printf("Building w-components.\n"); + for(int k = 1; k < nz; ++k) for(int j = 1; j < ny; ++j) for(int i = 1; i < nx; ++i) { + if(w_state(i,j,k) == FLUID) { + int index = w_ind(i,j,k,nx,ny,nz); + rhs[index] = vol_w(i,j,k)*w(i,j,k); + matrix.set_element(index,index, vol_w(i,j,k)); + + float visc_right = 0.25f*(viscosity(i,j,k) + viscosity(i,j,k-1) + viscosity(i+1,j,k) + viscosity(i+1,j,k-1)); + float visc_left = 0.25f*(viscosity(i,j,k) + viscosity(i,j,k-1) + viscosity(i-1,j,k) + viscosity(i-1,j,k-1)); + float vol_right = vol_ey(i+1,j,k); + float vol_left = vol_ey(i,j,k); + + float visc_top = 0.25f*(viscosity(i,j,k) + viscosity(i,j,k-1) + viscosity(i,j+1,k) + viscosity(i,j+1,k-1));; + float visc_bottom = 0.25f*(viscosity(i,j,k) + viscosity(i,j,k-1) + viscosity(i,j-1,k) + viscosity(i,j-1,k-1));; + float vol_top = vol_ex(i,j+1,k); + float vol_bottom = vol_ex(i,j,k); + + float visc_front = viscosity(i,j,k); + float visc_back = viscosity(i,j,k-1); + float vol_front = vol_c(i,j,k); + float vol_back = vol_c(i,j,k-1); + + //w_x_right + matrix.add_to_element(index,index, +factor*visc_right*vol_right); + if(w_state(i+1,j,k) == FLUID) + matrix.add_to_element(index,w_ind(i+1,j,k,nx,ny,nz), -factor*visc_right*vol_right); + else if (w_state(i+1,j,k) == SOLID) + rhs[index] -= -factor*visc_right*vol_right*w(i+1,j,k); + + //w_x_left + matrix.add_to_element(index,index, factor*visc_left*vol_left); + if(w_state(i-1,j,k) == FLUID) + matrix.add_to_element(index,w_ind(i-1,j,k,nx,ny,nz), -factor*visc_left*vol_left); + else if (w_state(i-1,j,k) == SOLID) + rhs[index] -= -factor*visc_left*vol_left*w(i-1,j,k); + + //w_y_top + matrix.add_to_element(index,index, +factor*visc_top*vol_top); + if(w_state(i,j+1,k) == FLUID) + matrix.add_to_element(index,w_ind(i,j+1,k,nx,ny,nz), -factor*visc_top*vol_top); + else if (w_state(i,j+1,k) == SOLID) + rhs[index] -= -factor*visc_top*vol_top*w(i,j+1,k); + + //w_y_bottom + matrix.add_to_element(index,index, factor*visc_bottom*vol_bottom); + if(w_state(i,j-1,k) == FLUID) + matrix.add_to_element(index,w_ind(i,j-1,k,nx,ny,nz), -factor*visc_bottom*vol_bottom); + else if (w_state(i,j-1,k) == SOLID) + rhs[index] -= -factor*visc_bottom*vol_bottom*w(i,j-1,k); + + //w_z_front + matrix.add_to_element(index,index, +2*factor*visc_front*vol_front); + if(w_state(i,j,k+1) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k+1,nx,ny,nz), -2*factor*visc_front*vol_front); + else if (w_state(i,j,k+1) == SOLID) + rhs[index] -= -2*factor*visc_front*vol_front*w(i,j,k+1); + + //w_z_back + matrix.add_to_element(index,index, +2*factor*visc_back*vol_back); + if(w_state(i,j,k-1) == FLUID) + matrix.add_to_element(index,w_ind(i,j,k-1,nx,ny,nz), -2*factor*visc_back*vol_back); + else if (w_state(i,j,k-1) == SOLID) + rhs[index] -= -2*factor*visc_back*vol_back*w(i,j,k-1); + + //u_z_right + if(u_state(i+1,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i+1,j,k,nx,ny), -factor*visc_right*vol_right); + else if(u_state(i+1,j,k) == SOLID) + rhs[index] -= -u(i+1,j,k)*factor*visc_right*vol_right; + + if(u_state(i+1,j,k-1) == FLUID) + matrix.add_to_element(index,u_ind(i+1,j,k-1,nx,ny), factor*visc_right*vol_right); + else if(u_state(i+1,j,k-1) == SOLID) + rhs[index] -= u(i+1,j,k-1)*factor*visc_right*vol_right; + + //u_z_left + if(u_state(i,j,k) == FLUID) + matrix.add_to_element(index,u_ind(i,j,k,nx,ny), factor*visc_left*vol_left); + else if(u_state(i,j,k) == SOLID) + rhs[index] -= u(i,j,k)*factor*visc_left*vol_left; + + if(u_state(i,j,k-1) == FLUID) + matrix.add_to_element(index,u_ind(i,j,k-1,nx,ny), -factor*visc_left*vol_left); + else if(u_state(i,j,k-1) == SOLID) + rhs[index] -= -u(i,j,k-1)*factor*visc_left*vol_left; + + //v_z_top + if(v_state(i,j+1,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j+1,k,nx,ny,nz), -factor*visc_top*vol_top); + else if(v_state(i,j+1,k) == SOLID) + rhs[index] -= -v(i,j+1,k)*factor*visc_top*vol_top; + + if(v_state(i,j+1,k-1) == FLUID) + matrix.add_to_element(index,v_ind(i,j+1,k-1,nx,ny,nz), factor*visc_top*vol_top); + else if(v_state(i,j+1,k-1) == SOLID) + rhs[index] -= v(i,j+1,k-1)*factor*visc_top*vol_top; + + //v_z_bottom + if(v_state(i,j,k) == FLUID) + matrix.add_to_element(index,v_ind(i,j,k,nx,ny,nz), +factor*visc_bottom*vol_bottom); + else if(v_state(i,j,k) == SOLID) + rhs[index] -= v(i,j,k)*factor*visc_bottom*vol_bottom; + + if(v_state(i,j,k-1) == FLUID) + matrix.add_to_element(index,v_ind(i,j,k-1,nx,ny,nz), -factor*visc_bottom*vol_bottom); + else if(v_state(i,j,k-1) == SOLID) + rhs[index] -= -v(i,j,k-1)*factor*visc_bottom*vol_bottom; + + } + } + + ////strip out near zero entries to speed this thing up! + //printf("Stripping out near-zeros.\n"); + ////SparseMatrixd matrix2(matrix.m,matrix.n,15); + //for(unsigned int row = 0; row < matrix.m; ++row){ + // for(unsigned int col = 0; col < matrix.index[row].size(); ++col) { + // int index = matrix.index[row][col]; + // double val = matrix(row,index); + // if(std::abs(val) > 1e-10) + // matrix2.set_element(row,index, val); + // } + //} + printf("Solving sparse system.\n"); + PCGSolver solver; + double res_out; + int iter_out; + solver.set_solver_parameters(1e-7, 500); + + printf("Launching CG\n"); + solver.solve(matrix, rhs, soln, res_out, iter_out); + + if(iter_out >= 500) + printf("\n\n\n**********FAILED**************\n\n\n"); + + printf("Copying back.\n"); + for(int k = 0; k < nz; ++k) + for(int j = 0; j < ny; ++j) + for(int i = 0; i < nx+1; ++i) { + if(u_state(i,j,k) == FLUID) { + u(i,j,k) = (float)soln[u_ind(i,j,k,nx,ny)]; + } + else { + u(i,j,k) = 0; + } + } + + for(int k = 0; k < nz; ++k) + for(int j = 0; j < ny+1; ++j) + for(int i = 0; i < nx; ++i) { + if(v_state(i,j,k) == FLUID) { + v(i,j,k) = (float)soln[v_ind(i,j,k,nx,ny,nz)]; + } + else + v(i,j,k) = 0; + } + + for(int k = 0; k < nz+1; ++k) + for(int j = 0; j < ny; ++j) + for(int i = 0; i < nx; ++i) { + if(w_state(i,j,k)== FLUID) { + w(i,j,k) = (float)soln[w_ind(i,j,k,nx,ny,nz)]; + } + else + w(i,j,k) = 0; + } + printf("Done copying back.\n"); + +} + diff --git a/viscosity3d.h b/viscosity3d.h new file mode 100644 index 0000000..d6309b8 --- /dev/null +++ b/viscosity3d.h @@ -0,0 +1,12 @@ +#ifndef VISCOSITY3D_H +#define VISCOSITY3D_H + +#include "array3.h" + +void advance_viscosity_implicit_weighted(Array3f& u, Array3f& v, Array3f& w, + const Array3f& vol_u, const Array3f& vol_v, const Array3f& vol_w, + const Array3f& vol_c, const Array3f& vol_ex, const Array3f& vol_ey, const Array3f& vol_ez, + const Array3f& solid_phi, + const Array3f& viscosity, float dt, float dx); +#endif + diff --git a/volume_fractions.cpp b/volume_fractions.cpp new file mode 100644 index 0000000..e03168b --- /dev/null +++ b/volume_fractions.cpp @@ -0,0 +1,139 @@ +#include "util.h" +#include "volume_fractions.h" + +// Assumes phi0<0 and phi1>=0, phi2>=0, or vice versa. +// In particular, phi0 must not equal either of phi1 or phi2. +template +static T sorted_triangle_fraction(T phi0, T phi1, T phi2) +{ + return phi0*phi0/(2*(phi0-phi1)*(phi0-phi2)); +} + +float area_fraction(float phi0, float phi1, float phi2) +{ + if(phi0<0){ + if(phi1<0){ + if(phi2<0) return 0; + else return 1-sorted_triangle_fraction(phi2, phi0, phi1); + }else if(phi2<0) return 1-sorted_triangle_fraction(phi1, phi2, phi0); + else return sorted_triangle_fraction(phi0, phi1, phi2); + }else if(phi1<0){ + if(phi2<0) return 1-sorted_triangle_fraction(phi0, phi1, phi2); + else return sorted_triangle_fraction(phi1, phi2, phi0); + }else if(phi2<0) return sorted_triangle_fraction(phi2, phi0, phi1); + else return 0; +} + +double area_fraction(double phi0, double phi1, double phi2) +{ + if(phi0<0){ + if(phi1<0){ + if(phi2<0) return 0; + else return 1-sorted_triangle_fraction(phi2, phi0, phi1); + }else if(phi2<0) return 1-sorted_triangle_fraction(phi1, phi2, phi0); + else return sorted_triangle_fraction(phi0, phi1, phi2); + }else if(phi1<0){ + if(phi2<0) return 1-sorted_triangle_fraction(phi0, phi1, phi2); + else return sorted_triangle_fraction(phi1, phi2, phi0); + }else if(phi2<0) return sorted_triangle_fraction(phi2, phi0, phi1); + else return 0; +} + +float area_fraction(float phi00, float phi10, float phi01, float phi11) +{ + float phimid=(phi00+phi10+phi01+phi11)/4; + return (area_fraction(phi00, phi10, phimid) + +area_fraction(phi10, phi11, phimid) + +area_fraction(phi11, phi01, phimid) + +area_fraction(phi01, phi00, phimid))/4; +} + +double area_fraction(double phi00, double phi10, double phi01, double phi11) +{ + double phimid=(phi00+phi10+phi01+phi11)/4; + return (area_fraction(phi00, phi10, phimid) + +area_fraction(phi10, phi11, phimid) + +area_fraction(phi11, phi01, phimid) + +area_fraction(phi01, phi00, phimid))/4; +} + +//============================================================================ + +// Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa. +// In particular, phi0 must not equal any of phi1, phi2 or phi3. +template +static T sorted_tet_fraction(T phi0, T phi1, T phi2, T phi3) +{ + return phi0*phi0*phi0/((phi0-phi1)*(phi0-phi2)*(phi0-phi3)); +} + +// Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa. +// In particular, phi0 and phi1 must not equal any of phi2 and phi3. +template +static T sorted_prism_fraction(T phi0, T phi1, T phi2, T phi3) +{ + T a=phi0/(phi0-phi2), + b=phi0/(phi0-phi3), + c=phi1/(phi1-phi3), + d=phi1/(phi1-phi2); + return a*b*(1-d)+b*(1-c)*d+c*d; +} + +float volume_fraction(float phi0, float phi1, float phi2, float phi3) +{ + sort(phi0, phi1, phi2, phi3); + if(phi3<=0) return 1; + else if(phi2<=0) return 1-sorted_tet_fraction(phi3, phi2, phi1, phi0); + else if(phi1<=0) return sorted_prism_fraction(phi0, phi1, phi2, phi3); + else if(phi0<=0) return sorted_tet_fraction(phi0, phi1, phi2, phi3); + else return 0; +} + +double volume_fraction(double phi0, double phi1, double phi2, double phi3) +{ + sort(phi0, phi1, phi2, phi3); + if(phi3<=0) return 1; + else if(phi2<=0) return 1-sorted_tet_fraction(phi3, phi2, phi1, phi0); + else if(phi1<=0) return sorted_prism_fraction(phi0, phi1, phi2, phi3); + else if(phi0<=0) return sorted_tet_fraction(phi0, phi1, phi2, phi3); + else return 0; +} + +float volume_fraction(float phi000, float phi100, + float phi010, float phi110, + float phi001, float phi101, + float phi011, float phi111) +{ + // This is the average of the two possible decompositions of the cube into + // five tetrahedra. + return ( volume_fraction(phi000, phi001, phi101, phi011) + +volume_fraction(phi000, phi101, phi100, phi110) + +volume_fraction(phi000, phi010, phi011, phi110) + +volume_fraction(phi101, phi011, phi111, phi110) + +2*volume_fraction(phi000, phi011, phi101, phi110) + +volume_fraction(phi100, phi101, phi001, phi111) + +volume_fraction(phi100, phi001, phi000, phi010) + +volume_fraction(phi100, phi110, phi111, phi010) + +volume_fraction(phi001, phi111, phi011, phi010) + +2*volume_fraction(phi100, phi111, phi001, phi010))/12; +} + +double volume_fraction(double phi000, double phi100, + double phi010, double phi110, + double phi001, double phi101, + double phi011, double phi111) +{ + // This is the average of the two possible decompositions of the cube into + // five tetrahedra. + return ( volume_fraction(phi000, phi001, phi101, phi011) + +volume_fraction(phi000, phi101, phi100, phi110) + +volume_fraction(phi000, phi010, phi011, phi110) + +volume_fraction(phi101, phi011, phi111, phi110) + +2*volume_fraction(phi000, phi011, phi101, phi110) + +volume_fraction(phi100, phi101, phi001, phi111) + +volume_fraction(phi100, phi001, phi000, phi010) + +volume_fraction(phi100, phi110, phi111, phi010) + +volume_fraction(phi001, phi111, phi011, phi010) + +2*volume_fraction(phi100, phi111, phi001, phi010))/12; +} + diff --git a/volume_fractions.h b/volume_fractions.h new file mode 100644 index 0000000..e5d6349 --- /dev/null +++ b/volume_fractions.h @@ -0,0 +1,30 @@ +#ifndef VOLUME_FRACTIONS_H +#define VOLUME_FRACTIONS_H + +// Given a triangle with level set values, use linear interpolation to +// estimate the fraction of the triangle occupied by the phi<0 part +float area_fraction(float phi0, float phi1, float phi2); +double area_fraction(double phi0, double phi1, double phi2); + +// Given a rectangle with level set values, estimate fraction occupied +// by the phi<0 part +float area_fraction(float phi00, float phi10, float phi01, float phi11); +double area_fraction(double phi00, double phi10, double phi01, double phi11); + +// Given a tetrahedron with level set values, use linear interpolation to +// estimate the fraction of the tetrahedron occupied by the phi<0 part +float volume_fraction(float phi0, float phi1, float phi2, float phi3); +double volume_fraction(double phi0, double phi1, double phi2, double phi3); + +// Given a parallelepiped (e.g. cube) with level set values, estimate +// fraction occupied by the phi<0 part +float volume_fraction(float phi000, float phi100, + float phi010, float phi110, + float phi001, float phi101, + float phi011, float phi111); +double volume_fraction(double phi000, double phi100, + double phi010, double phi110, + double phi001, double phi101, + double phi011, double phi111); + +#endif