Extras/sph/common/geomx.cpp

   1 /*
   2   FLUIDS v.1 - SPH Fluid Simulator for CPU and GPU
   3   Copyright (C) 2008. Rama Hoetzlein, http://www.rchoetzlein.com
   4
   5   ZLib license
   6   This software is provided 'as-is', without any express or implied
   7   warranty.  In no event will the authors be held liable for any damages
   8   arising from the use of this software.
   9
  10   Permission is granted to anyone to use this software for any purpose,
  11   including commercial applications, and to alter it and redistribute it
  12   freely, subject to the following restrictions:
  13
  14   1. The origin of this software must not be misrepresented; you must not
  15      claim that you wrote the original software. If you use this software
  16      in a product, an acknowledgment in the product documentation would be
  17      appreciated but is not required.
  18   2. Altered source versions must be plainly marked as such, and must not be
  19      misrepresented as being the original software.
  20   3. This notice may not be removed or altered from any source distribution.
  21 */
  22
  23
  24 #include <vector.h>
  25 #include "geomx.h"
  26 #include "mdebug.h"
  27
  28 #include <assert.h>
  29
  30 GeomX::GeomX ()
  31 {
  32         mHeap = 0x0;
  33 }
  34
  35 void GeomX::FreeBuffers ()
  36 {
  37         for (int n=0; n < (int) mBuf.size(); n++) {
  38                 free ( mBuf[n].data );
  39                 mBuf[n].data = 0x0;
  40         }
  41         mBuf.clear ();
  42
  43         if ( mHeap != 0x0 ) free ( mHeap );
  44         mHeap = 0x0;
  45         mHeapNum = 0;
  46         mHeapMax = 0;
  47 }
  48
  49 void GeomX::ResetHeap ()
  50 {
  51         mHeapNum = 0;
  52         mHeapFree = -1;
  53 }
  54
  55 int GeomX::GetSize ()
  56 {
  57         int sum = mHeapNum * sizeof(hval) ;
  58         for (int n=0; n < (int) mBuf.size(); n++)
  59                 sum += mBuf[n].size;
  60         sum += (int) mAttribute.size() * sizeof(GeomAttr);
  61         return sum;
  62 }
  63
  64 void GeomX::ClearAttributes ()
  65 {
  66         mAttribute.clear ();
  67 }
  68
  69 int GeomX::CopyBuffer ( uchar bdest, uchar bsrc, GeomX& src )
  70 {
  71         if ( bsrc >= src.GetNumBuf() ) return -1;
  72
  73         if ( bdest >= GetNumBuf() ) {
  74                 for (int n=0; n <= bdest - GetNumBuf(); n++ )
  75                         AddBuffer ( 0, 0, 0 );
  76         }
  77         if ( mBuf[bdest].data != 0x0 ) {
  78                 free ( mBuf[bdest].data );
  79                 mBuf[bdest].data = 0x0;
  80         }
  81
  82         GeomBuf* buf = src.GetBuffer( bsrc );
  83         mBuf[bdest].dtype = buf->dtype;
  84         mBuf[bdest].max = buf->max;
  85         mBuf[bdest].num = buf->num;
  86         mBuf[bdest].stride = buf->stride;
  87         mBuf[bdest].size = buf->size;
  88         mBuf[bdest].data = (char*) malloc (  mBuf[bdest].max * mBuf[bdest].stride );
  89         memcpy ( mBuf[bdest].data, buf->data, mBuf[bdest].num * mBuf[bdest].stride );
  90
  91         return bdest;
  92 }
  93
  94
  95 void GeomX::CopyBuffers ( GeomX& src )
  96 {
  97         FreeBuffers ();
  98         for (int n = 0; n < src.GetNumBuf(); n++)
  99                 CopyBuffer ( n, n, src );
 100         CopyHeap ( src );
 101 }
 102
 103 void GeomX::CopyHeap ( GeomX& src )
 104 {
 105         hpos num, max, freepos;
 106         hval* src_data = src.GetHeap ( num, max, freepos );
 107
 108         if ( mHeap != 0x0 ) {
 109                 free ( mHeap );
 110                 mHeap = 0x0;
 111         }
 112
 113         mHeap = (hval*) malloc ( max * sizeof(hval) );
 114         mHeapMax = max;
 115         mHeapNum = num;
 116         mHeapFree = freepos;
 117
 118         memcpy ( mHeap, src_data, mHeapNum * sizeof(hval) );
 119 }
 120
 121 hval* GeomX::GetHeap ( hpos& num, hpos& max, hpos& free )
 122 {
 123         num = mHeapNum;
 124         max = mHeapMax;
 125         free = mHeapFree;
 126         return mHeap;
 127 }
 128
 129 void GeomX::CopyAttributes ( GeomX& src )
 130 {
 131         GeomAttr* attr;
 132         ClearAttributes ();
 133         int a;
 134         for (int n = 0; n < src.GetNumAttr(); n++) {
 135                 attr = src.GetAttribute ( n );
 136                 a = AddAttribute ( (uchar) attr->buf, attr->name, attr->stride, false );
 137                 mAttribute[a].offset = attr->offset;
 138         }
 139 }
 140
 141
 142 void GeomX::AddHeap ( int max )
 143 {
 144         mHeap = (hval*) malloc ( max * sizeof(hval ) );
 145         mHeapMax = max;
 146         mHeapNum = 0;
 147         mHeapFree = -1;
 148 }
 149
 150 int GeomX::AddAttribute ( uchar b, std::string name, ushort stride )
 151 {
 152         return AddAttribute ( b, name, stride, true );
 153 }
 154
 155 int GeomX::AddAttribute ( uchar b, std::string name, ushort stride, bool bExtend )
 156 {
 157         GeomBuf* buf = &mBuf[b];
 158         GeomAttr attr;
 159
 160         attr.buf = b;
 161         attr.name = name;
 162         attr.offset = buf->stride;
 163         attr.stride = stride;
 164
 165         if ( bExtend ) {
 166                 int new_stride = buf->stride + attr.stride;
 167                 char* new_data = (char*) malloc ( buf->max * new_stride );
 168                 char* old_data = buf->data;
 169                 for (int n=0; n < buf->num; n++) {
 170                         memcpy ( new_data, old_data, buf->stride );
 171                         old_data += buf->stride;
 172                         new_data += new_stride;
 173                 }
 174                 free ( buf->data );
 175                 buf->data = new_data;
 176                 buf->stride = new_stride;
 177         }
 178         mAttribute.push_back ( attr );
 179         return (int) mAttribute.size()-1;
 180 }
 181
 182 int GeomX::GetAttribute ( std::string name )
 183 {
 184         for (int n=0; n < (int) mAttribute.size(); n++) {
 185                 if ( mAttribute[n].name.compare ( name ) == 0 )
 186                         return n;
 187         }
 188         return -1;
 189 }
 190
 191 int GeomX::GetAttrOffset ( std::string name )
 192 {
 193         for (int n=0; n < (int) mAttribute.size(); n++) {
 194                 if ( mAttribute[n].name.compare ( name ) == 0 )
 195                         return mAttribute[n].offset;
 196         }
 197         return -1;
 198 }
 199
 200 int GeomX::AddBuffer ( uchar typ, ushort stride, int max )
 201 {
 202         GeomBuf buf;
 203         buf.dtype = typ;
 204         buf.stride = stride;
 205         buf.max = max;
 206         buf.num = 0;
 207         buf.size = 0;
 208         buf.data = (char*) malloc ( buf.max * buf.stride );
 209         mBuf.push_back ( buf );
 210         return (int) mBuf.size()-1;
 211 }
 212 void GeomX::ResetBuffer ( uchar b, int n )
 213 {
 214         mBuf[b].max = n;
 215
 216         if ( mBuf[b].data != 0x0 ) free ( mBuf[b].data );
 217
 218         char* new_data = (char*) malloc ( mBuf[b].max * mBuf[b].stride );
 219
 220         mBuf[b].data = new_data;
 221
 222         mBuf[b].num = 0;
 223         mBuf[b].size = mBuf[b].num*mBuf[b].stride;
 224 }
 225
 226 char* GeomX::AddElem ( uchar b, href& ndx )
 227 {
 228         if ( mBuf[b].num >= mBuf[b].max ) {
 229                 if ( long(mBuf[b].max) * 2 > ELEM_MAX ) {
 230                         error.PrintF ( "geom", "Maximum number of elements reached.\n" );
 231                         error.Exit ();
 232                 }
 233                 mBuf[b].max *= 2;
 234                 char* new_data = (char*) malloc ( mBuf[b].max * mBuf[b].stride );
 235                 memcpy ( new_data, mBuf[b].data, mBuf[b].num*mBuf[b].stride );
 236                 free ( mBuf[b].data );
 237                 mBuf[b].data = new_data;
 238         }
 239         mBuf[b].num++;
 240         mBuf[b].size += mBuf[b].stride;
 241         ndx = mBuf[b].num-1;
 242         return mBuf[b].data + ndx*mBuf[b].stride;
 243 }
 244
 245 char* GeomX::RandomElem ( uchar b, href& ndx )
 246 {
 247         ndx = mBuf[b].num * rand() / RAND_MAX;
 248         return mBuf[b].data + ndx*mBuf[b].stride;
 249 }
 250
 251 int GeomX::AddElem ( uchar b, char* data )
 252 {
 253         if ( mBuf[b].num >= mBuf[b].max ) {
 254                 mBuf[b].max *= 2;
 255                 char* new_data = (char*) malloc ( mBuf[b].max * mBuf[b].stride );
 256                 memcpy ( new_data, mBuf[b].data, mBuf[b].num*mBuf[b].stride );
 257                 free ( mBuf[b].data );
 258                 mBuf[b].data = new_data;
 259         }
 260         memcpy ( mBuf[b].data + mBuf[b].num*mBuf[b].stride, data, mBuf[b].stride );
 261         mBuf[b].num++;
 262         mBuf[b].size += mBuf[b].stride;
 263         return mBuf[b].num-1;
 264 }
 265
 266 bool GeomX::DelElem ( uchar b, int n )
 267 {
 268         if ( n >= 0 && n < mBuf[b].num ) {
 269                 memcpy ( mBuf[b].data + n*mBuf[b].stride, mBuf[b].data + (mBuf[b].num-1)*mBuf[b].stride, mBuf[b].stride );
 270                 mBuf[b].num--;
 271                 mBuf[b].size -= mBuf[b].stride;
 272                 return true;
 273         }
 274         return false;
 275 }
 276
 277 hpos GeomX::HeapExpand ( ushort size, ushort& ret  )
 278 {
 279         mHeapMax *= 2;
 280         if ( mHeapMax > HEAP_MAX ) {
 281                 error.PrintF ( "geom", "Geom heap size exceeds range of index.\n" );
 282                 error.Exit ();
 283         }
 284         hval* pNewHeap = (hval*) malloc ( mHeapMax * sizeof(hval) );
 285         if ( pNewHeap == 0x0 ) {
 286                 error.PrintF ( "geom", "Geom heap out of memory.\n" );
 287                 error.Exit ();
 288         }
 289         memcpy ( pNewHeap, mHeap, mHeapNum*sizeof(hval) );
 290         free ( mHeap );
 291         mHeap = pNewHeap;
 292         ret = size;
 293         assert ( mHeapNum >= 0 && mHeapNum < mHeapMax );
 294         return mHeapNum;
 295 }
 296
 297 #define LIST_INIT               4
 298
 299 hpos GeomX::HeapAlloc ( ushort size, ushort& ret  )
 300 {
 301         hval* pPrev = 0x0;
 302         hval* pCurr = mHeap + mHeapFree;
 303         hpos pos = -1;
 304
 305         if ( mHeapNum + size < mHeapMax ) {
 306                 // Heap not yet full.
 307                 pos = mHeapNum;
 308                 ret = size;
 309                 mHeapNum += size;
 310         } else {
 311                 // Heap full, search free space
 312                 if ( mHeapFree == -1 ) {
 313                         pos = HeapExpand ( size, ret );
 314                         mHeapNum += ret;
 315                 } else {
 316                         while ( *pCurr < size && pCurr != mHeap-1 ) {
 317                                 pPrev = pCurr;
 318                                 pCurr = mHeap + * (hpos*) (pCurr + FPOS);
 319                         }
 320                         if ( pCurr != mHeap-1 ) {
 321                                 // Found free space.
 322                                 if ( pPrev == 0x0 ) {
 323                                         mHeapFree = * (hpos*) (pCurr + FPOS);
 324                                         assert ( mHeapFree >= -1 );
 325                                 } else {
 326                                         * (hpos*) (pPrev+FPOS) = * (hpos*) (pCurr+FPOS);
 327                                 }
 328                                 pos = pCurr-mHeap;
 329                                 ret = *pCurr;
 330                         } else {
 331                                 // Heap full, no free space. Expand heap.
 332                                 pos = HeapExpand ( size, ret );
 333                                 mHeapNum += ret;
 334                         }
 335                 }
 336         }
 337
 338         assert ( pos >= 0 && pos <= mHeapNum );
 339         memset ( mHeap+pos, 0, size*sizeof(hval) );
 340         return pos;
 341 }
 342
 343 void GeomX::HeapAddFree ( hpos pos, int size )
 344 {
 345         // memset ( mHeap+pos, 0xFF, size*sizeof(hval) );
 346
 347         if ( pos < mHeapFree || mHeapFree == -1 ) {
 348                 // Lowest position. Insert at head of heap.
 349                 * (hpos*) (mHeap+pos) = size;
 350                 * (hpos*) (mHeap+pos + FPOS) = mHeapFree;
 351                 mHeapFree = pos;
 352                 if ( mHeapFree != -1 && *(hpos*) (mHeap+mHeapFree+FPOS) == 0x0 ) { error.PrintF ( "geomx", "ERROR: Heap pointer 0. pHeapFree, %d\n", pos ); error.Exit (); }
 353
 354                 assert ( mHeapFree >= -1 );
 355         } else {
 356                 hval* pCurr = mHeap + pos;
 357                 hval* pPrev = 0x0;
 358                 hval* pNext = mHeap + mHeapFree;
 359
 360                 if ( pCurr == pNext ) {         // Freeing the first block
 361                         mHeapFree = * (hpos*) (pCurr + FPOS);
 362                         * (hpos*) (mHeap+mHeapFree + FPOS) = 0xFFFFFFFF;
 363                         * (hpos*) (pCurr + FPOS) = 0xFFFFFFFF;
 364                         * (hpos*) pCurr = 0xFFFFFFFF;
 365                         return;
 366
 367                 }
 368
 369                 // Find first block greater than new free pos.
 370                 while ( pNext < pCurr && pNext != mHeap-1 ) {
 371                         pPrev = pNext;
 372                         pNext = mHeap + * (hpos*) (pNext + FPOS);
 373                 }
 374
 375                 int x = 0;
 376                 if ( pPrev + *(pPrev) == pCurr ) {                                                      // Prev block touches current one (expand previous)
 377                         * (hpos*) pPrev += size;                                                                // - prev.size = prev.size + curr.size
 378                         x=1;
 379
 380                 } else if ( pCurr + size == pNext && pNext != mHeap-1 ) {       // Curr block touches next one (expand current block)
 381                         * (hpos*) (pPrev + FPOS) = pos;                                                 // - prev.next = curr
 382                         * (hpos*) (pCurr + FPOS) = * (hpos*) (pNext + FPOS);    // - curr.next = next.next
 383                         * (hpos*) pCurr = size + * (hpos*) pNext;                               // - curr.size = size + next.size
 384                         * (hpos*) (pNext) = 0xFFFFFFFF;                                                 // - curr = null
 385                         * (hpos*) (pNext + FPOS) = 0xFFFFFFFF;                                  // - curr.next = null
 386                         x=2;
 387
 388                 } else {                                                                                                        // Otherwise (linked list insert)
 389                         * (hpos*) (pPrev + FPOS) = pos;                                                 // - prev.next = curr
 390                         if ( pNext != mHeap-1 )
 391                                 * (hpos*) (pCurr + FPOS) = pNext - mHeap;                       // - curr.next = next
 392                         else
 393                                 * (hpos*) (pCurr + FPOS) = 0xFFFFFFFF;                          // - curr.next = null (no next node)
 394                         * (hpos*) pCurr = size;                                                                 // - curr.size = size
 395                         x=3;
 396                 }
 397                 if ( pCurr !=mHeap-1 && *(hpos*) (pCurr+FPOS) == 0x0 ) { error.PrintF ( "geomx", "ERROR: Heap pointer 0. pCurr, %d, %d\n", x, pos ); error.Exit (); }
 398                 if ( pPrev !=mHeap-1 && *(hpos*) (pPrev+FPOS) == 0x0 ) { error.PrintF ( "geomx", "ERROR: Heap pointer 0. pPrev, %d, %d\n", x, pos ); error.Exit (); }
 399                 if ( pNext !=mHeap-1 && *(hpos*) (pNext+FPOS) == 0x0 ) { error.PrintF ( "geomx", "ERROR: Heap pointer 0. pNext, %d, %d\n", x, pos ); error.Exit (); }
 400                 if ( *(hpos*) (mHeap+mHeapFree+FPOS) == 0x0 ) { error.PrintF ( "geomx", "ERROR: Heap pointer 0. pHeapFree, %d, %d\n", x, pos ); error.Exit (); }
 401
 402                 // -- check for bugs (remove eventually)
 403                 pNext = mHeap + mHeapFree;
 404                 while ( pNext != mHeap-1 ) {
 405                         pPrev = pNext;
 406                         pNext = mHeap + * (hpos*) (pNext + FPOS);
 407                         if ( pNext < pPrev && pNext != mHeap-1 ) {
 408                                 error.PrintF ( "geomx", "ERROR: Heap free space out of order. %d, %d\n", x, pos );
 409                                 error.Exit ();
 410                         }
 411                 }
 412                 //---
 413         }
 414  }
 415
 416 void GeomX::ClearRefs ( hList& list )
 417 {
 418         list.cnt = 0;
 419         list.max = 0;
 420         list.pos = 0;
 421 }
 422
 423 hval GeomX::AddRef ( hval r, hList& list, hval delta )
 424 {
 425         if ( list.max == 0 ) {
 426                 list.cnt = 1;
 427                 list.pos = HeapAlloc ( LIST_INIT, list.max );
 428                 *(mHeap + list.pos) = r+delta;
 429         } else {
 430                 if ( list.cnt >= list.max ) {
 431                         int siz = list.max;
 432                         hpos new_pos = HeapAlloc ( siz+LIST_INIT, list.max );                           // Alloc new location
 433                         //printf ( "MOVE %d -> %d\n", list.pos, new_pos );
 434                         memcpy ( mHeap+new_pos, mHeap+list.pos, list.cnt*sizeof(hval) );        // Copy data to new location
 435                         HeapAddFree ( list.pos, siz );                                                                          // Free old location
 436                         list.pos = new_pos;
 437
 438                 }
 439                 *(mHeap + list.pos + list.cnt) = r+delta;
 440                 list.cnt++;
 441         }
 442         return list.pos;
 443 }