[80a6a52] | 1 | #ifndef __GEOWORLDVMPDG_H__
|
---|
| 2 | #define __GEOWORLDVMPDG_H__
|
---|
| 3 |
|
---|
| 4 | #include <mars_grids.h>
|
---|
| 5 | #include <mars_ptr.h>
|
---|
| 6 | #include <mars_util.h>
|
---|
| 7 |
|
---|
| 8 | #include <vector>
|
---|
| 9 | #include <set>
|
---|
| 10 |
|
---|
| 11 | namespace geoworld
|
---|
| 12 | {
|
---|
| 13 | template< typename T >
|
---|
| 14 | class MPDGVecType
|
---|
| 15 | {
|
---|
| 16 | public:
|
---|
| 17 | typedef T Precision;
|
---|
| 18 |
|
---|
| 19 | T altitude;
|
---|
| 20 | unsigned short mindex;
|
---|
| 21 |
|
---|
| 22 | MPDGVecType ()
|
---|
| 23 | : altitude(0), mindex(0) {}
|
---|
| 24 |
|
---|
| 25 | MPDGVecType (const MPDGVecType & copy)
|
---|
| 26 | : altitude(copy.altitude), mindex(copy.mindex) {}
|
---|
| 27 | MPDGVecType (const T altitude, const unsigned short nIdx)
|
---|
| 28 | : altitude(altitude), mindex(nIdx) {}
|
---|
| 29 |
|
---|
| 30 | inline MPDGVecType & operator += (const MPDGVecType & v)
|
---|
| 31 | { return operator += (v.altitude); }
|
---|
| 32 | inline MPDGVecType & operator -= (const MPDGVecType & v)
|
---|
| 33 | { return operator -= (v.altitude); }
|
---|
| 34 | inline MPDGVecType operator + (const MPDGVecType & v) const
|
---|
| 35 | { return operator + (v.altitude); }
|
---|
| 36 | inline MPDGVecType operator - (const MPDGVecType & v) const
|
---|
| 37 | { return operator - (v.altitude); }
|
---|
| 38 | inline MPDGVecType & operator += (const T & n)
|
---|
| 39 | {
|
---|
| 40 | altitude += n;
|
---|
| 41 | return *this;
|
---|
| 42 | }
|
---|
| 43 | inline MPDGVecType & operator -= (const T & n)
|
---|
| 44 | {
|
---|
| 45 | altitude -= n;
|
---|
| 46 | return *this;
|
---|
| 47 | }
|
---|
| 48 | inline MPDGVecType operator + (const T & n) const
|
---|
| 49 | {
|
---|
| 50 | return MPDGVecType(altitude + n, mindex);
|
---|
| 51 | }
|
---|
| 52 | inline MPDGVecType operator - (const T & n) const
|
---|
| 53 | {
|
---|
| 54 | return MPDGVecType(altitude - n, mindex);
|
---|
| 55 | }
|
---|
| 56 | inline MPDGVecType operator << (const unsigned short n) const
|
---|
| 57 | {
|
---|
| 58 | return MPDGVecType(altitude << n, mindex);
|
---|
| 59 | }
|
---|
| 60 | inline MPDGVecType operator >> (const unsigned short n) const
|
---|
| 61 | {
|
---|
| 62 | return MPDGVecType(altitude >> n, mindex);
|
---|
| 63 | }
|
---|
| 64 |
|
---|
| 65 | inline bool operator < (const MPDGVecType & b) const
|
---|
| 66 | { return mindex < b.mindex || (mindex == b.mindex && altitude < b.altitude); }
|
---|
| 67 | inline bool operator > (const MPDGVecType & b) const
|
---|
| 68 | { return mindex > b.mindex || (mindex == b.mindex && altitude > b.altitude); }
|
---|
| 69 |
|
---|
| 70 | inline bool operator <= (const MPDGVecType & b) const
|
---|
| 71 | { return mindex <= b.mindex || (mindex == b.mindex && altitude <= b.altitude); }
|
---|
| 72 | inline bool operator >= (const MPDGVecType & b) const
|
---|
| 73 | { return mindex >= b.mindex || (mindex == b.mindex && altitude >= b.altitude); }
|
---|
| 74 | };
|
---|
| 75 |
|
---|
| 76 | template< typename T >
|
---|
| 77 | struct MPDGVecTypeNumerics
|
---|
| 78 | {
|
---|
| 79 | typedef long double UP;
|
---|
| 80 | typedef double POS;
|
---|
| 81 |
|
---|
| 82 | static inline geoworld::MPDGVecType< T > zero () { return geoworld::MPDGVecType< T > (); }
|
---|
| 83 | static inline geoworld::MPDGVecType< T > step() { return geoworld::MPDGVecType< T > (mars::Numerics< T >::step(), 0); }
|
---|
| 84 | static inline geoworld::MPDGVecType< T > min () { return geoworld::MPDGVecType< T >(std::numeric_limits< T >::lowest(), std::numeric_limits< unsigned short >::lowest()); }
|
---|
| 85 | static inline geoworld::MPDGVecType< T > max () { return geoworld::MPDGVecType< T >(std::numeric_limits< T >::max(), std::numeric_limits< unsigned short >::max()); }
|
---|
| 86 | };
|
---|
| 87 |
|
---|
| 88 | template< typename E, typename T >
|
---|
| 89 | class ElementalFieldView : public mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >
|
---|
| 90 | {
|
---|
| 91 | public:
|
---|
| 92 | typedef std::vector< E > ElementList;
|
---|
| 93 |
|
---|
| 94 | using ScalarFieldView = typename mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
|
---|
| 95 | using Offset = typename ScalarFieldView::Offset;
|
---|
| 96 |
|
---|
| 97 | private:
|
---|
| 98 | ElementList * _pElements;
|
---|
| 99 |
|
---|
| 100 | protected:
|
---|
| 101 | ElementalFieldView (mars::MemoryBlock< MPDGVecType< T > > * pBlock, const mars::BBox< unsigned short > & bbox, ElementList * pElemList, const mars::GridType engt = mars::Normal)
|
---|
| 102 | : _pElements(pElemList), ScalarFieldView (pBlock, bbox, engt) {}
|
---|
| 103 | ElementalFieldView (mars::MemoryBlock< MPDGVecType< T > > * pBlock, const unsigned short nLeft, const unsigned short nTop, const unsigned short nWidth, const unsigned short nHeight, ElementList * pElemList, const mars::GridType engt = mars::Normal)
|
---|
| 104 | : _pElements(pElemList), ScalarFieldView (pBlock, nLeft, nTop, nWidth, nHeight, engt) {}
|
---|
| 105 |
|
---|
| 106 | // Reserved for use by the ElementalField and subclasses, used to duplicate the memory block and element list
|
---|
| 107 | ElementalFieldView (const ElementalFieldView & copy, mars::MemoryBlock< E > * const pBlock, ElementList * const pElements)
|
---|
| 108 | : ScalarFieldView(copy, pBlock), _pElements(pElements) {}
|
---|
| 109 |
|
---|
| 110 | // Copy constructor
|
---|
| 111 | ElementalFieldView (const ElementalFieldView & copy)
|
---|
| 112 | : ScalarFieldView(copy), _pElements(copy._pElements) {}
|
---|
| 113 |
|
---|
| 114 | ElementList * elements() { return _pElements; }
|
---|
| 115 |
|
---|
| 116 | static ElementalFieldView< E, T > createHosted (const unsigned short nWidth, const unsigned short nHeight, const mars::GridType engt)
|
---|
| 117 | { return ElementalFieldView< E, T > (new mars::MemoryBlock< MPDGVecType< T > > (nWidth, nHeight), 0, 0, nWidth, nHeight, new ElementList(), engt); }
|
---|
| 118 |
|
---|
| 119 | mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > rangeImpl (const int x0, const int y0, const int xN, const int yN) const
|
---|
| 120 | {
|
---|
| 121 | T
|
---|
| 122 | fAltMax = std::numeric_limits <T>::min(),
|
---|
| 123 | fAltMin = std::numeric_limits <T>::max();
|
---|
| 124 | unsigned short
|
---|
| 125 | nIdxMax = std::numeric_limits <unsigned short>::lowest(),
|
---|
| 126 | nIdxMin = std::numeric_limits <unsigned short>::max();
|
---|
| 127 |
|
---|
| 128 | int x, y;
|
---|
| 129 | int iti;
|
---|
| 130 |
|
---|
| 131 | #ifdef _OPENMP
|
---|
| 132 | #define RANGEIMPL_THREADCOUNT 4
|
---|
| 133 | #else
|
---|
| 134 | #define RANGEIMPL_THREADCOUNT 1
|
---|
| 135 | #endif
|
---|
| 136 | unsigned short vnIdxMin4[RANGEIMPL_THREADCOUNT], vnIdxMax4[RANGEIMPL_THREADCOUNT];
|
---|
| 137 | T vfAltMin4[RANGEIMPL_THREADCOUNT], vfAltMax4[RANGEIMPL_THREADCOUNT];
|
---|
| 138 |
|
---|
| 139 | for (iti = 0; iti < RANGEIMPL_THREADCOUNT; ++iti)
|
---|
| 140 | {
|
---|
| 141 | vfAltMin4[iti] = fAltMin;
|
---|
| 142 | vfAltMax4[iti] = fAltMax;
|
---|
| 143 | vnIdxMin4[iti] = nIdxMin;
|
---|
| 144 | vnIdxMax4[iti] = nIdxMax;
|
---|
| 145 | }
|
---|
| 146 |
|
---|
| 147 | #ifdef _OPENMP
|
---|
| 148 | #pragma omp parallel for private (x, y) num_threads(RANGEIMPL_THREADCOUNT)
|
---|
| 149 | #endif
|
---|
| 150 | for (x = x0; x <= xN; ++x)
|
---|
| 151 | {
|
---|
| 152 | const int ti =
|
---|
| 153 | #ifdef _OPENMP
|
---|
| 154 | omp_get_thread_num();
|
---|
| 155 | #else
|
---|
| 156 | 0;
|
---|
| 157 | #endif
|
---|
| 158 | for (y = y0; y <= yN; ++y)
|
---|
| 159 | {
|
---|
| 160 | const MPDGVecType< T > & elem = ScalarFieldView::get(x, y);
|
---|
| 161 |
|
---|
| 162 | if (elem.altitude > vfAltMax4[ti])
|
---|
| 163 | vfAltMax4[ti] = elem.altitude;
|
---|
| 164 | if (elem.altitude < vfAltMin4[ti])
|
---|
| 165 | vfAltMin4[ti] = elem.altitude;
|
---|
| 166 |
|
---|
| 167 | if (elem.mindex > vnIdxMax4[ti])
|
---|
| 168 | vnIdxMax4[ti] = elem.mindex;
|
---|
| 169 | if (elem.mindex < vnIdxMin4[ti])
|
---|
| 170 | vnIdxMin4[ti] = elem.mindex;
|
---|
| 171 | }
|
---|
| 172 | }
|
---|
| 173 |
|
---|
| 174 | for (iti = 0; iti < RANGEIMPL_THREADCOUNT; ++iti)
|
---|
| 175 | {
|
---|
| 176 | if (vfAltMax4[iti] > fAltMax)
|
---|
| 177 | fAltMax = vfAltMax4[iti];
|
---|
| 178 | if (vfAltMin4[iti] < fAltMin)
|
---|
| 179 | fAltMin = vfAltMin4[iti];
|
---|
| 180 |
|
---|
| 181 | if (vnIdxMax4[iti] > nIdxMax)
|
---|
| 182 | nIdxMax = vnIdxMax4[iti];
|
---|
| 183 | if (vnIdxMin4[iti] < nIdxMin)
|
---|
| 184 | nIdxMin = vnIdxMin4[iti];
|
---|
| 185 | }
|
---|
| 186 |
|
---|
| 187 | return mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics < T > > (MPDGVecType< T >(fAltMin, nIdxMin), MPDGVecType< T >(fAltMax, nIdxMax));
|
---|
| 188 | }
|
---|
| 189 |
|
---|
| 190 | friend class mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
|
---|
| 191 |
|
---|
| 192 | public:
|
---|
| 193 | typedef E Element;
|
---|
| 194 | typedef ElementalFieldView< E, T > View;
|
---|
| 195 | typedef T AltitudePrecision;
|
---|
| 196 |
|
---|
| 197 | ElementalFieldView(ElementalFieldView && move) : ScalarFieldView(move), _pElements(move._pElements) {}
|
---|
| 198 |
|
---|
| 199 | ElementalFieldView < E, T > * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom);
|
---|
| 200 | const ElementalFieldView < E, T > * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const;
|
---|
| 201 |
|
---|
| 202 | inline const E & getElement (const signed int x, const signed int y) const
|
---|
| 203 | {
|
---|
| 204 | return _pElements->at(
|
---|
| 205 | ScalarFieldView::get(
|
---|
| 206 | x + static_cast< const ScalarFieldView * const > (this) ->left,
|
---|
| 207 | y + static_cast< const ScalarFieldView * const > (this) ->top
|
---|
| 208 | ).mindex
|
---|
| 209 | );
|
---|
| 210 | }
|
---|
| 211 | inline size_t elementCount () const
|
---|
| 212 | { return _pElements->size(); }
|
---|
| 213 | inline const E & getElement (const size_t nIdx) const
|
---|
| 214 | { return _pElements->at(nIdx); }
|
---|
| 215 | inline E & getElement (const size_t nIdx)
|
---|
| 216 | { return _pElements->at(nIdx); }
|
---|
| 217 |
|
---|
| 218 | signed int findElementIndex (const E & elem) const
|
---|
| 219 | {
|
---|
| 220 | unsigned short nIdx = 0;
|
---|
| 221 | for (typename ElementList::const_iterator i = _pElements->begin(); i != _pElements->end(); ++i, ++nIdx)
|
---|
| 222 | if (*i == elem)
|
---|
| 223 | return nIdx;
|
---|
| 224 |
|
---|
| 225 | return -1;
|
---|
| 226 | }
|
---|
| 227 |
|
---|
| 228 | size_t countViewElements () const
|
---|
| 229 | {
|
---|
| 230 | std::set< unsigned short > setMIdx;
|
---|
| 231 |
|
---|
| 232 | for (size_t j = 0; j < ScalarFieldView::height; ++j)
|
---|
| 233 | for (size_t i = 0; i < ScalarFieldView::width; ++i)
|
---|
| 234 | setMIdx.insert(ScalarFieldView::get(i, j).mindex);
|
---|
| 235 |
|
---|
| 236 | return setMIdx.size();
|
---|
| 237 | }
|
---|
| 238 |
|
---|
| 239 | void fetchViewElements (std::set< E > & st) const
|
---|
| 240 | {
|
---|
| 241 | for (size_t j = 0; j < ScalarFieldView::height; ++j)
|
---|
| 242 | for (size_t i = 0; i < ScalarFieldView::width; ++i)
|
---|
| 243 | st.insert(getElement(i, j));
|
---|
| 244 | }
|
---|
| 245 | void fetchViewIndices (std::set< unsigned short > & st) const
|
---|
| 246 | {
|
---|
| 247 | for (size_t j = 0; j < ScalarFieldView::height; ++j)
|
---|
| 248 | for (size_t i = 0; i < ScalarFieldView::width; ++i)
|
---|
| 249 | st.insert(ScalarFieldView::operator () (i, j) .mindex);
|
---|
| 250 | }
|
---|
| 251 |
|
---|
| 252 | inline typename ElementalFieldView< E, T >::ElementList::const_iterator beginElements() const
|
---|
| 253 | { return _pElements->begin(); }
|
---|
| 254 | inline typename ElementalFieldView< E, T >::ElementList::const_iterator endElements() const
|
---|
| 255 | { return _pElements->end(); }
|
---|
| 256 |
|
---|
| 257 | template< typename L >
|
---|
| 258 | inline void reg (const typename L::const_iterator iBegin, const typename L::const_iterator iEnd, const L * = NULL)
|
---|
| 259 | { _pElements->insert(_pElements->end(), iBegin, iEnd); }
|
---|
| 260 |
|
---|
| 261 | inline void reg (const E & e)
|
---|
| 262 | { _pElements->push_back(e); }
|
---|
| 263 |
|
---|
| 264 | inline void clearRegs ()
|
---|
| 265 | { _pElements->clear(); }
|
---|
| 266 |
|
---|
| 267 | inline const T & altitude (const signed int x, const signed int y) const
|
---|
| 268 | { return ScalarFieldView::get(x, y).altitude; }
|
---|
| 269 | inline T & altitude (const signed int x, const signed int y)
|
---|
| 270 | { return ScalarFieldView::get(x, y).altitude; }
|
---|
| 271 |
|
---|
| 272 | inline void operator += (const T & n)
|
---|
| 273 | { ScalarFieldView::operator += (MPDGVecType< T >(n, 0)); }
|
---|
| 274 | inline void operator -= (const T & n)
|
---|
| 275 | { ScalarFieldView::operator -= (MPDGVecType< T >(n, 0)); }
|
---|
| 276 |
|
---|
| 277 | inline virtual mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > range () const
|
---|
| 278 | { return rangeImpl (0, 0, ScalarFieldView::right - ScalarFieldView::left, ScalarFieldView::bottom - ScalarFieldView::top); }
|
---|
| 279 |
|
---|
| 280 | inline void operator << (const ElementalFieldView< E, T > & copy)
|
---|
| 281 | {
|
---|
| 282 | _pElements->clear();
|
---|
| 283 | _pElements->insert(_pElements->end(), copy.beginElements(), copy.endElements());
|
---|
| 284 | ScalarFieldView::operator << (copy);
|
---|
| 285 | }
|
---|
| 286 | inline void operator >> (ElementalFieldView< E, T > & dest) const
|
---|
| 287 | { dest << *this; }
|
---|
| 288 |
|
---|
| 289 | inline void operator >> (Offset & dest) const
|
---|
| 290 | { ScalarFieldView::operator >> (dest); }
|
---|
| 291 | inline void operator << (const Offset & src)
|
---|
| 292 | { ScalarFieldView::operator << (src); }
|
---|
| 293 | };
|
---|
| 294 |
|
---|
| 295 | template< typename E, typename T >
|
---|
| 296 | ElementalFieldView < E, T > * ElementalFieldView< E, T >::createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
|
---|
| 297 | { return ScalarFieldView::template createTypedView< ElementalFieldView< E, T > > (nLeft, nTop, nRight, nBottom, _pElements); }
|
---|
| 298 |
|
---|
| 299 | template< typename E, typename T >
|
---|
| 300 | const ElementalFieldView < E, T > * ElementalFieldView< E, T >::createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
|
---|
| 301 | { return ScalarFieldView::template createTypedView< ElementalFieldView< E, T > > (nLeft, nTop, nRight, nBottom, _pElements); }
|
---|
| 302 |
|
---|
| 303 | template <typename E, typename T> class VarriedMPDG;
|
---|
| 304 |
|
---|
| 305 | template <typename E, typename T>
|
---|
| 306 | class ElementalField : public ElementalFieldView< E, T >
|
---|
| 307 | {
|
---|
| 308 | public:
|
---|
| 309 | using MemoryBlock = mars::MemoryBlock< MPDGVecType< T > >;
|
---|
| 310 | using ScalarFieldView = mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
|
---|
| 311 | using VarriedMPDGType = VarriedMPDG< E, T >;
|
---|
| 312 | using ElementList = typename ElementalFieldView< E, T >::ElementList;
|
---|
| 313 |
|
---|
| 314 | ElementalField (const unsigned short width, const unsigned short height, const mars::ScalarField< MPDGVecType< T > > & copy, const unsigned short x, const unsigned short y, const mars::GridType engt = mars::Normal)
|
---|
| 315 | : ElementalFieldView< E, T >(VarriedMPDGType::createHosted(width, height, engt))
|
---|
| 316 | { blitFrom (copy, x, y); }
|
---|
| 317 |
|
---|
| 318 | ElementalField (const unsigned short width, const unsigned short height, const mars::GridType engt = mars::Normal)
|
---|
| 319 | : ElementalFieldView< E, T >(ElementalFieldView< E, T >::createHosted(width, height, engt))
|
---|
| 320 | { memset(this->data(), 0, this->block() -> bytesize); }
|
---|
| 321 | ElementalField (const ElementalFieldView< E, T > & copy)
|
---|
| 322 | : ElementalFieldView< E, T >(copy, new MemoryBlock (copy.getMemoryBlock()), new ElementList (*copy._pElements)) {}
|
---|
| 323 |
|
---|
| 324 | virtual ~ElementalField()
|
---|
| 325 | { delete this->elements(); }
|
---|
| 326 |
|
---|
| 327 | inline std::istream & operator << (std::istream & ins)
|
---|
| 328 | { return *this->block() << ins; }
|
---|
| 329 | inline std::ostream & operator >> (std::ostream & outs) const
|
---|
| 330 | { return *this->block() >> outs; }
|
---|
| 331 |
|
---|
| 332 | inline void operator >> (typename ScalarFieldView::Offset & dest) const
|
---|
| 333 | { ElementalFieldView< E, T >::operator >> (dest); }
|
---|
| 334 | inline void operator << (const typename ScalarFieldView::Offset & src)
|
---|
| 335 | { ElementalFieldView< E, T >::operator << (src); }
|
---|
| 336 | };
|
---|
| 337 |
|
---|
| 338 | template <typename E, typename T>
|
---|
| 339 | class VarriedMPDG : public mars::MidpointDisplacementGrid< MPDGVecType< T >, MPDGVecTypeNumerics< T > >, public ElementalFieldView< E, T >
|
---|
| 340 | {
|
---|
| 341 | private:
|
---|
| 342 | const enum StorageMode
|
---|
| 343 | {
|
---|
| 344 | Concrete,
|
---|
| 345 | Virtual
|
---|
| 346 | } _ensmStorageMode;
|
---|
| 347 |
|
---|
| 348 | protected:
|
---|
| 349 | using MPDG = mars::MidpointDisplacementGrid< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
|
---|
| 350 | using Elems = ElementalFieldView< E, T >;
|
---|
| 351 |
|
---|
| 352 | class SeedFunctor : public MPDG::PinkFn
|
---|
| 353 | {
|
---|
| 354 | private:
|
---|
| 355 | const unsigned short _maxidx;
|
---|
| 356 | float _base;
|
---|
| 357 |
|
---|
| 358 | public:
|
---|
| 359 | inline SeedFunctor (unsigned short maxidx, const unsigned short nDim, const float fCoarseness)
|
---|
| 360 | : MPDG::PinkFn(nDim, fCoarseness), _maxidx(maxidx) {}
|
---|
| 361 |
|
---|
| 362 | inline void iteration (const unsigned short nStride, const unsigned short c)
|
---|
| 363 | {
|
---|
| 364 | MPDG::PinkFn::iteration(nStride, c);
|
---|
| 365 | _base = this->computeRandomAmount(0) / 2;
|
---|
| 366 | }
|
---|
| 367 |
|
---|
| 368 | inline MPDGVecType< T > operator () (int & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c)
|
---|
| 369 | {
|
---|
| 370 | MPDGVecType< T > result;
|
---|
| 371 |
|
---|
| 372 | result.mindex = mars::RAND(_maxidx);
|
---|
| 373 | result.altitude = static_cast< T > (this->applyFlux (_base));
|
---|
| 374 | return result;
|
---|
| 375 | }
|
---|
| 376 | inline MPDGVecType< T > operator () (int & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c, const MPDGVecType< T > & d)
|
---|
| 377 | {
|
---|
| 378 | return operator() (state, x, y, a, b, c);
|
---|
| 379 | }
|
---|
| 380 | };
|
---|
| 381 |
|
---|
| 382 | struct LocalFunctorState
|
---|
| 383 | {
|
---|
| 384 | T sumalt[3];
|
---|
| 385 | unsigned char ctalt[3];
|
---|
| 386 | MPDGVecType< T > result;
|
---|
| 387 | };
|
---|
| 388 |
|
---|
| 389 | class LocalFunctor : public MPDG::template AbstractFunctor<LocalFunctorState>
|
---|
| 390 | {
|
---|
| 391 | public:
|
---|
| 392 | inline LocalFunctor (const unsigned short nDim, const float fCoarseness)
|
---|
| 393 | : MPDG::template AbstractFunctor<LocalFunctorState>(nDim, fCoarseness) {}
|
---|
| 394 |
|
---|
| 395 | MPDGVecType< T > operator () (LocalFunctorState & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c)
|
---|
| 396 | {
|
---|
| 397 | // Count cells matching 'a'
|
---|
| 398 | state.ctalt[0] = 0;
|
---|
| 399 | state.sumalt[0] = a.altitude;
|
---|
| 400 | if (a.mindex == b.mindex)
|
---|
| 401 | {
|
---|
| 402 | state.sumalt[0] += b.altitude;
|
---|
| 403 | ++state.ctalt[0];
|
---|
| 404 | }
|
---|
| 405 | if (a.mindex == c.mindex)
|
---|
| 406 | {
|
---|
| 407 | state.sumalt[0] += c.altitude;
|
---|
| 408 | ++state.ctalt[0];
|
---|
| 409 | }
|
---|
| 410 |
|
---|
| 411 | // Count cells matching 'b'
|
---|
| 412 | state.ctalt[1] = 0;
|
---|
| 413 | state.sumalt[1] = b.altitude;
|
---|
| 414 | if (b.mindex == c.mindex)
|
---|
| 415 | {
|
---|
| 416 | state.sumalt[1] += c.altitude;
|
---|
| 417 | ++state.ctalt[1];
|
---|
| 418 | }
|
---|
| 419 |
|
---|
| 420 | // Flip a coin to pick between a dominant alpha or auxiliary beta cell
|
---|
| 421 | const unsigned int pick = mars::RAND(3);
|
---|
| 422 |
|
---|
| 423 | if (pick < 2)
|
---|
| 424 | state.result.altitude = this->applyFlux(state.sumalt[pick] / (state.ctalt[pick] + 1));
|
---|
| 425 | else
|
---|
| 426 | state.result.altitude = this->applyFlux(c.altitude);
|
---|
| 427 |
|
---|
| 428 | switch (pick)
|
---|
| 429 | {
|
---|
| 430 | case 0:
|
---|
| 431 | state.result.mindex = a.mindex;
|
---|
| 432 | break;
|
---|
| 433 | case 1:
|
---|
| 434 | state.result.mindex = b.mindex;
|
---|
| 435 | break;
|
---|
| 436 | case 2:
|
---|
| 437 | state.result.mindex = c.mindex;
|
---|
| 438 | break;
|
---|
| 439 | }
|
---|
| 440 |
|
---|
| 441 | return state.result;
|
---|
| 442 | }
|
---|
| 443 | MPDGVecType< T > operator () (LocalFunctorState & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c, const MPDGVecType< T > & d)
|
---|
| 444 | {
|
---|
| 445 | // Count cells matching 'a'
|
---|
| 446 | state.ctalt[0] = 0;
|
---|
| 447 | state.sumalt[0] = a.altitude;
|
---|
| 448 | if (a.mindex == b.mindex)
|
---|
| 449 | {
|
---|
| 450 | state.sumalt[0] += b.altitude;
|
---|
| 451 | ++state.ctalt[0];
|
---|
| 452 | }
|
---|
| 453 | if (a.mindex == c.mindex)
|
---|
| 454 | {
|
---|
| 455 | state.sumalt[0] += c.altitude;
|
---|
| 456 | ++state.ctalt[0];
|
---|
| 457 | }
|
---|
| 458 | if (a.mindex == d.mindex)
|
---|
| 459 | {
|
---|
| 460 | state.sumalt[0] += d.altitude;
|
---|
| 461 | ++state.ctalt[0];
|
---|
| 462 | }
|
---|
| 463 |
|
---|
| 464 | // Count cells matching 'b'
|
---|
| 465 | state.ctalt[1] = 0;
|
---|
| 466 | state.sumalt[1] = b.altitude;
|
---|
| 467 | if (b.mindex == c.mindex)
|
---|
| 468 | {
|
---|
| 469 | state.sumalt[1] += c.altitude;
|
---|
| 470 | ++state.ctalt[1];
|
---|
| 471 | }
|
---|
| 472 | if (b.mindex == d.mindex)
|
---|
| 473 | {
|
---|
| 474 | state.sumalt[1] += d.altitude;
|
---|
| 475 | ++state.ctalt[1];
|
---|
| 476 | }
|
---|
| 477 |
|
---|
| 478 | // Count cells matching 'c'
|
---|
| 479 | state.ctalt[2] = 0;
|
---|
| 480 | state.sumalt[2] = c.altitude;
|
---|
| 481 | if (c.mindex == d.mindex)
|
---|
| 482 | {
|
---|
| 483 | state.sumalt[2] += d.altitude;
|
---|
| 484 | ++state.ctalt[2];
|
---|
| 485 | }
|
---|
| 486 |
|
---|
| 487 | // Flip a coin to pick between a dominant alpha or auxiliary beta cell
|
---|
| 488 |
|
---|
| 489 | const unsigned int pick = mars::RAND(4);
|
---|
| 490 |
|
---|
| 491 | if (pick < 3)
|
---|
| 492 | state.result.altitude = this->applyFlux(state.sumalt[pick] / (state.ctalt[pick] + 1));
|
---|
| 493 | else
|
---|
| 494 | state.result.altitude = this->applyFlux(d.altitude);
|
---|
| 495 |
|
---|
| 496 | switch (pick)
|
---|
| 497 | {
|
---|
| 498 | case 0:
|
---|
| 499 | state.result.mindex = a.mindex;
|
---|
| 500 | break;
|
---|
| 501 | case 1:
|
---|
| 502 | state.result.mindex = b.mindex;
|
---|
| 503 | break;
|
---|
| 504 | case 2:
|
---|
| 505 | state.result.mindex = c.mindex;
|
---|
| 506 | break;
|
---|
| 507 | case 3:
|
---|
| 508 | state.result.mindex = d.mindex;
|
---|
| 509 | break;
|
---|
| 510 | }
|
---|
| 511 |
|
---|
| 512 | return state.result;
|
---|
| 513 | }
|
---|
| 514 | };
|
---|
| 515 |
|
---|
| 516 | VarriedMPDG (const short nIterCount, ElementalFieldView< E, T > v, const StorageMode ensm)
|
---|
| 517 | : MPDG(nIterCount, v), ElementalFieldView<E, T>(v), _ensmStorageMode(ensm) {}
|
---|
| 518 |
|
---|
| 519 | static VarriedMPDG * createHosted (const short nIterCount, const unsigned short nWidth, const unsigned short nHeight, const mars::GridType engt = mars::Normal)
|
---|
| 520 | { return new VarriedMPDG(nIterCount, ElementalFieldView< E, T >::createHosted(nWidth, nHeight, engt), Concrete); }
|
---|
| 521 |
|
---|
| 522 | public:
|
---|
| 523 | using ScalarFieldView = typename mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
|
---|
| 524 |
|
---|
| 525 | virtual ~VarriedMPDG()
|
---|
| 526 | {
|
---|
| 527 | switch (_ensmStorageMode)
|
---|
| 528 | {
|
---|
| 529 | case Concrete: delete this->elements(); break;
|
---|
| 530 | }
|
---|
| 531 | }
|
---|
| 532 |
|
---|
| 533 | inline static VarriedMPDG * createInstance (const unsigned short nMinWidth, const unsigned short nMinHeight, const mars::GridType engt)
|
---|
| 534 | {
|
---|
| 535 | assert(nMinHeight > 0 && nMinWidth > 0);
|
---|
| 536 | return createHosted (
|
---|
| 537 | MPDG::calcIterCount(std::min(nMinHeight, nMinWidth)),
|
---|
| 538 | MPDG::dimensionFor(nMinWidth, engt),
|
---|
| 539 | MPDG::dimensionFor(nMinHeight, engt),
|
---|
| 540 | engt
|
---|
| 541 | );
|
---|
| 542 | }
|
---|
| 543 | inline static VarriedMPDG * createInstance (Elems & view)
|
---|
| 544 | {
|
---|
| 545 | return new VarriedMPDG (
|
---|
| 546 | calcIterCount(std::min(view.width, view.height) - (view.gridtype == mars::Normal ? 1 : 0)),
|
---|
| 547 | view,
|
---|
| 548 | Virtual
|
---|
| 549 | );
|
---|
| 550 | }
|
---|
| 551 |
|
---|
| 552 | inline virtual mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > range () const
|
---|
| 553 | { return Elems::range(); }
|
---|
| 554 | inline virtual void clear()
|
---|
| 555 | { Elems::clear(); }
|
---|
| 556 | inline void stitch(const int nStitch)
|
---|
| 557 | { Elems::stitch(nStitch); }
|
---|
| 558 | void stitch(const int nStitch, const ScalarFieldView * pOther)
|
---|
| 559 | { Elems::stitch(nStitch, pOther); }
|
---|
| 560 | void stitchAll(const int nStitch)
|
---|
| 561 | { Elems::stitchAll(nStitch); }
|
---|
| 562 | inline void stitchAll(const int nStitch, const ScalarFieldView * pLeft, const ScalarFieldView * pTop, const ScalarFieldView * pRight, const ScalarFieldView * pBottom)
|
---|
| 563 | { Elems::stitchAll(nStitch, pLeft, pTop, pRight, pBottom); }
|
---|
| 564 |
|
---|
| 565 | // TODO: Type hierarchy seems broken (bogus diamond pattern)
|
---|
| 566 | inline void operator -= (const T & n)
|
---|
| 567 | {
|
---|
| 568 | Elems::operator -= (n);
|
---|
| 569 | }
|
---|
| 570 |
|
---|
| 571 | void seed (unsigned short nIterations, const float fCoarseNess, const signed int nFlux = -1)
|
---|
| 572 | {
|
---|
| 573 | SeedFunctor f (Elems::elementCount(), nFlux < 0 ? this->mindim : static_cast< unsigned short > (nFlux), fCoarseNess);
|
---|
| 574 |
|
---|
| 575 | this->iterations(nIterations, f);
|
---|
| 576 | }
|
---|
| 577 |
|
---|
| 578 | void run (const float fCoarseNess, const int nStitch = 0)
|
---|
| 579 | {
|
---|
| 580 | LocalFunctor f (this->mindim, fCoarseNess);
|
---|
| 581 |
|
---|
| 582 | this->iterations (this->itercount - this->index(), f, nStitch);
|
---|
| 583 | }
|
---|
| 584 | };
|
---|
| 585 | }
|
---|
| 586 |
|
---|
| 587 | #endif
|
---|