1 
2 
3 // DO NOT EDIT !
4 // This file is generated using the MantaFlow preprocessor (prep generate).
5 
6 /******************************************************************************
7  *
8  * MantaFlow fluid solver framework
9  * Copyright 2011 Tobias Pfaff, Nils Thuerey
10  *
11  * This program is free software, distributed under the terms of the
12  * Apache License, Version 2.0
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Fast marching
16  *
17  ******************************************************************************/
18 
19 #ifndef _FASTMARCH_H
20 #define _FASTMARCH_H
21 
22 #include <queue>
23 #include "levelset.h"
24 
25 namespace Manta {
26 
27 //! Fast marching. Transport certain values
28 //  This class exists in two versions: for scalar, and for vector values - the only difference are
29 //  flag checks i transpTouch (for simplicity in separate classes)
30 
31 template<class GRID, class T>
fmInterpolateNeighbors(GRID * mpVal,int x,int y,int z,Real * weights)32 inline T fmInterpolateNeighbors(GRID *mpVal, int x, int y, int z, Real *weights)
33 {
34   T val(0.);
35   if (weights[0] > 0.0)
36     val += mpVal->get(x + 1, y + 0, z + 0) * weights[0];
37   if (weights[1] > 0.0)
38     val += mpVal->get(x - 1, y + 0, z + 0) * weights[1];
39   if (weights[2] > 0.0)
40     val += mpVal->get(x + 0, y + 1, z + 0) * weights[2];
41   if (weights[3] > 0.0)
42     val += mpVal->get(x + 0, y - 1, z + 0) * weights[3];
43   if (mpVal->is3D()) {
44     if (weights[4] > 0.0)
45       val += mpVal->get(x + 0, y + 0, z + 1) * weights[4];
46     if (weights[5] > 0.0)
47       val += mpVal->get(x + 0, y + 0, z - 1) * weights[5];
48   }
49   return val;
50 }
51 
52 template<class GRID, class T> class FmValueTransportScalar {
53  public:
FmValueTransportScalar()54   FmValueTransportScalar() : mpVal(0), mpFlags(0){};
~FmValueTransportScalar()55   ~FmValueTransportScalar(){};
initMarching(GRID * val,FlagGrid * flags)56   void initMarching(GRID *val, FlagGrid *flags)
57   {
58     mpVal = val;
59     mpFlags = flags;
60   }
isInitialized()61   inline bool isInitialized()
62   {
63     return mpVal != 0;
64   }
65 
66   //! cell is touched by marching from source cell
transpTouch(int x,int y,int z,Real * weights,Real time)67   inline void transpTouch(int x, int y, int z, Real *weights, Real time)
68   {
69     if (!mpVal || !mpFlags->isEmpty(x, y, z))
70       return;
71     T val = fmInterpolateNeighbors<GRID, T>(mpVal, x, y, z, weights);
72     (*mpVal)(x, y, z) = val;
73   };
74 
75  protected:
76   GRID *mpVal;
77   FlagGrid *mpFlags;
78 };
79 
80 template<class GRID, class T> class FmValueTransportVec3 {
81  public:
FmValueTransportVec3()82   FmValueTransportVec3() : mpVal(0), mpFlags(0){};
~FmValueTransportVec3()83   ~FmValueTransportVec3(){};
isInitialized()84   inline bool isInitialized()
85   {
86     return mpVal != 0;
87   }
initMarching(GRID * val,const FlagGrid * flags)88   void initMarching(GRID *val, const FlagGrid *flags)
89   {
90     mpVal = val;
91     mpFlags = flags;
92   }
93 
94   //! cell is touched by marching from source cell
transpTouch(int x,int y,int z,Real * weights,Real time)95   inline void transpTouch(int x, int y, int z, Real *weights, Real time)
96   {
97     if (!mpVal || !mpFlags->isEmpty(x, y, z))
98       return;
99 
100     T val = fmInterpolateNeighbors<GRID, T>(mpVal, x, y, z, weights);
101 
102     // set velocity components if adjacent is empty
103     if (mpFlags->isEmpty(x - 1, y, z))
104       (*mpVal)(x, y, z).x = val.x;
105     if (mpFlags->isEmpty(x, y - 1, z))
106       (*mpVal)(x, y, z).y = val.y;
107     if (mpVal->is3D()) {
108       if (mpFlags->isEmpty(x, y, z - 1))
109         (*mpVal)(x, y, z).z = val.z;
110     }
111   };
112 
113  protected:
114   GRID *mpVal;
115   const FlagGrid *mpFlags;
116 };
117 
118 class FmHeapEntryOut {
119  public:
120   Vec3i p;
121   // quick time access for sorting
122   Real time;
compare(const Real x,const Real y)123   static inline bool compare(const Real x, const Real y)
124   {
125     return x > y;
126   }
127 
128   inline bool operator<(const FmHeapEntryOut &o) const
129   {
130     const Real d = fabs((time) - ((o.time)));
131     if (d > 0.)
132       return (time) > ((o.time));
133     if (p.z != o.p.z)
134       return p.z > o.p.z;
135     if (p.y != o.p.y)
136       return p.y > o.p.y;
137     return p.x > o.p.x;
138   };
139 };
140 
141 class FmHeapEntryIn {
142  public:
143   Vec3i p;
144   // quick time access for sorting
145   Real time;
compare(const Real x,const Real y)146   static inline bool compare(const Real x, const Real y)
147   {
148     return x < y;
149   }
150 
151   inline bool operator<(const FmHeapEntryIn &o) const
152   {
153     const Real d = fabs((time) - ((o.time)));
154     if (d > 0.)
155       return (time) < ((o.time));
156     if (p.z != o.p.z)
157       return p.z < o.p.z;
158     if (p.y != o.p.y)
159       return p.y < o.p.y;
160     return p.x < o.p.x;
161   };
162 };
163 
164 //! fast marching algorithm wrapper class
165 template<class T, int TDIR> class FastMarch {
166 
167  public:
168   // MSVC doesn't allow static const variables in template classes
InvalidTime()169   static inline Real InvalidTime()
170   {
171     return -1000;
172   }
InvtOffset()173   static inline Real InvtOffset()
174   {
175     return 500;
176   }
177 
178   enum SpecialValues { FlagInited = 1, FlagIsOnHeap = 2 };
179 
180   FastMarch(const FlagGrid &flags,
181             Grid<int> &fmFlags,
182             Grid<Real> &levelset,
183             Real maxTime,
184             MACGrid *velTransport = NULL);
~FastMarch()185   ~FastMarch()
186   {
187   }
188 
189   //! advect level set function with given velocity */
190   void performMarching();
191 
192   //! test value for invalidity
isInvalid(Real v)193   inline bool isInvalid(Real v) const
194   {
195     return (v <= InvalidTime());
196   }
197 
198   void addToList(const Vec3i &p, const Vec3i &src);
199 
200   //! convert phi to time value
phi2time(Real phival)201   inline Real phi2time(Real phival)
202   {
203     return (phival - InvalidTime() + InvtOffset()) * -1.0;
204   }
205 
206   //! ... and back
time2phi(Real tval)207   inline Real time2phi(Real tval)
208   {
209     return (InvalidTime() - InvtOffset() - tval);
210   }
211 
_phi(int i,int j,int k)212   inline Real _phi(int i, int j, int k)
213   {
214     return mLevelset(i, j, k);
215   }
216 
217  protected:
218   Grid<Real> &mLevelset;
219   const FlagGrid &mFlags;
220   Grid<int> &mFmFlags;
221 
222   //! velocity extrpolation
223   FmValueTransportVec3<MACGrid, Vec3> mVelTransport;
224 
225   //! maximal time to march for
226   Real mMaxTime;
227 
228   //! fast marching list
229   std::priority_queue<T, std::vector<T>, std::less<T>> mHeap;
230   Real mReheapVal;
231 
232   //! weights for touching points
233   Real mWeights[6];
234 
235   template<int C> inline Real calcWeights(int &okCnt, int &invcnt, Real *v, const Vec3i &idx);
236 
237   inline Real calculateDistance(const Vec3i &pos);
238 };
239 
240 }  // namespace Manta
241 #endif
242