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