00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef __utils_hpp__
00022 #define __utils_hpp__
00023
00024 #include <vector>
00025 #include <list>
00026
00027 template <class T>
00028 class array16
00029 {
00030 public:
00031 array16 () : mem_ptr(0), mem_begin(0), mem_end(0), len(0) { }
00032
00033 void resize (size_t s)
00034 {
00035 unsigned char *old_mem_ptr = mem_ptr;
00036 T *old_mem_begin = mem_begin;
00037 size_t old_len = len;
00038
00039 len = s;
00040 if ( old_len < len )
00041 {
00042 mem_ptr = (unsigned char *)malloc (16 + (len + 1) * sizeof (T));
00043 size_t n = ((size_t)mem_ptr)%16;
00044 mem_begin = (T *)(mem_ptr + n);
00045 mem_end = mem_begin + len;
00046 if ( old_mem_ptr )
00047 {
00048 memcpy (mem_begin, old_mem_begin, (len < old_len ? len : old_len) * sizeof (T));
00049 free (old_mem_ptr);
00050 }
00051 }
00052 }
00053
00054 void erase (size_t index)
00055 {
00056 memcpy (mem_begin + index, mem_begin + index + 1, (len - index - 1) * sizeof (T));
00057 resize (len - 1);
00058 }
00059
00060 void release ()
00061 {
00062 if ( mem_ptr )
00063 free (mem_ptr);
00064 mem_ptr = 0;
00065 mem_begin = mem_end = 0;
00066 len = 0;
00067 }
00068
00069 T& operator[] (size_t i) { return mem_begin[i]; }
00070 const T& operator[] (size_t i) const { return mem_begin[i]; }
00071
00072 T* begin() { return mem_begin; }
00073 T* end() { return mem_end; }
00074 const T* begin() const { return mem_begin; }
00075 const T* end() const { return mem_end; }
00076
00077 size_t size () const { return len; }
00078
00079 ~array16 ()
00080 {
00081 release ();
00082 }
00083
00084 private:
00085 unsigned char *mem_ptr;
00086 T *mem_begin;
00087 T *mem_end;
00088 size_t len;
00089
00090 array16 (const array16&);
00091 array16& operator= (const array16&);
00092 };
00093
00094 template <class T, size_t CI = 32>
00095 class allocator
00096 {
00097 public:
00098 allocator() : it(0), end(0), memsize(0) { }
00099
00100 ~allocator ()
00101 {
00102 clear();
00103 }
00104
00105 void clear ()
00106 {
00107 while ( !cache.empty() )
00108 {
00109 free (cache.front());
00110 cache.pop_front();
00111 }
00112 it = end = 0;
00113 memsize = 0;
00114 }
00115
00116 void reserve (size_t n)
00117 {
00118 if ( (size_t)(end - it) < n )
00119 {
00120 size_t bc = sizeof (T) * n;
00121 memsize += bc;
00122 it = (T *)malloc (bc);
00123 end = it + n;
00124 cache.push_back (it);
00125 }
00126 }
00127
00128 T * get()
00129 {
00130 if ( it == end )
00131 reserve (CI);
00132 return it++;
00133 }
00134
00135 T * get(size_t n)
00136 {
00137 reserve (n);
00138 T *ret = it;
00139 it += n;
00140 return ret;
00141 }
00142
00143 T * getz()
00144 {
00145 T *ret = get();
00146 memset (ret, 0, sizeof(T));
00147 return ret;
00148 }
00149
00150 T * getz(size_t n)
00151 {
00152 reserve (n);
00153 T *ret = it;
00154 it += n;
00155 memset (ret, 0, sizeof(T) * n);
00156 return ret;
00157 }
00158
00159 size_t size() const { return memsize; }
00160
00161 private:
00162 typedef std::list<T *> cache_t;
00163 cache_t cache;
00164 T *it;
00165 T *end;
00166 size_t memsize;
00167 };
00168
00169 class FPWriter
00170 {
00171 public:
00172 FPWriter (FILE *fp) : fp(fp) { }
00173
00174 void write (const void *ptr, size_t len)
00175 {
00176 fwrite (ptr, len, 1, fp);
00177 }
00178
00179 void write_char (char c)
00180 {
00181 fwrite (&c, 1, 1, fp);
00182 }
00183
00184 void write (size_t v)
00185 {
00186 unsigned int n = v;
00187 fwrite (&n, sizeof (n), 1, fp);
00188 }
00189
00190 void write (float v)
00191 {
00192 fwrite (&v, sizeof (v), 1, fp);
00193 }
00194
00195 void write (double v)
00196 {
00197 fwrite (&v, sizeof (v), 1, fp);
00198 }
00199
00200 void write (const vect3d& v)
00201 {
00202 write (prec2double (v.x));
00203 write (prec2double (v.y));
00204 write (prec2double (v.z));
00205 }
00206
00207 ~FPWriter () { if ( fp ) fclose (fp); }
00208
00209 private:
00210 FILE *fp;
00211 };
00212
00213 class FPReader
00214 {
00215 public:
00216 FPReader (FILE *fp) : fp(fp) { }
00217
00218 void read (void *buf, size_t len)
00219 {
00220 fread (buf, len, 1, fp);
00221 }
00222
00223 void read (double& v)
00224 {
00225 fread (&v, sizeof(v), 1, fp);
00226 }
00227
00228 void read (float& v)
00229 {
00230 fread (&v, sizeof(v), 1, fp);
00231 }
00232
00233 void read (vect3d& v)
00234 {
00235 double x, y, z;
00236 read (x);
00237 read (y);
00238 read (z);
00239 v.x = x;
00240 v.y = y;
00241 v.z = z;
00242 }
00243
00244 void read (size_t& v)
00245 {
00246 unsigned int n;
00247 fread (&n, sizeof(n), 1, fp);
00248 v = n;
00249 }
00250
00251 void read_char (char& c)
00252 {
00253 fread (&c, 1, 1, fp);
00254 }
00255
00256 ~FPReader () { if ( fp ) fclose (fp); }
00257
00258 private:
00259 FILE *fp;
00260 };
00261
00262 template <int N>
00263 class differator
00264 {
00265 public:
00266 differator (const vect3d& resolution) : res(resolution) { }
00267 virtual ~differator () { }
00268
00269 virtual void getValues (const vect3d& pos, prec_t *vals) = 0;
00270
00271 void getDifference (const vect3d& pos, prec_t* dfdx, prec_t *dfdy, prec_t* dfdz,
00272 prec_t* d2fdxdy, prec_t* d2fdxdz, prec_t* d2fdydz,
00273 prec_t* d3fdxdydz)
00274 {
00275 vect3d p1, p2, p3, p4, p5, p6, p7, p8;
00276 prec_t h, h1, h2, h3;
00277 prec_t fp1[N], fp2[N], fp3[N], fp4[N], fp5[N], fp6[N], fp7[N], fp8[N];
00278
00279 p1 = pos - vect3d(res.x, 0, 0);
00280 p2 = pos + vect3d(res.x, 0, 0);
00281 getValues (p1, fp1);
00282 getValues (p2, fp2);
00283 h = p2.x - p1.x;
00284 for ( size_t i = 0; i < N; i++ )
00285 dfdx[i] = (fp2[i] - fp1[i]) / h;
00286
00287 p1 = pos - vect3d(0, res.y, 0);
00288 p2 = pos + vect3d(0, res.y, 0);
00289 getValues (p1, fp1);
00290 getValues (p2, fp2);
00291 h = p2.y - p1.y;
00292 for ( size_t i = 0; i < N; i++ )
00293 dfdy[i] = (fp2[i] - fp1[i]) / h;
00294
00295 p1 = pos - vect3d(0, 0, res.z);
00296 p2 = pos + vect3d(0, 0, res.z);
00297 getValues (p1, fp1);
00298 getValues (p2, fp2);
00299 h = p2.z - p1.z;
00300 for ( size_t i = 0; i < N; i++ )
00301 dfdz[i] = (fp2[i] - fp1[i]) / h;
00302
00303 p1 = pos + vect3d(-res.x, -res.y, 0);
00304 p2 = pos + vect3d(-res.x, res.y, 0);
00305 p3 = pos + vect3d(res.x, -res.y, 0);
00306 p4 = pos + vect3d(res.x, res.y, 0);
00307 getValues (p1, fp1);
00308 getValues (p2, fp2);
00309 getValues (p3, fp3);
00310 getValues (p4, fp4);
00311 h1 = p3.x - p1.x;
00312 h2 = p2.y - p1.y;
00313 h = h1*h2;
00314 for ( size_t i = 0; i < N; i++ )
00315 d2fdxdy[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;
00316
00317 p1 = pos + vect3d(-res.x, 0, -res.z);
00318 p2 = pos + vect3d(-res.x, 0, res.z);
00319 p3 = pos + vect3d(res.x, 0, -res.z);
00320 p4 = pos + vect3d(res.x, 0, res.z);
00321 getValues (p1, fp1);
00322 getValues (p2, fp2);
00323 getValues (p3, fp3);
00324 getValues (p4, fp4);
00325 h1 = p3.x - p1.x;
00326 h2 = p2.z - p1.z;
00327 h = h1*h2;
00328 for ( size_t i = 0; i < N; i++ )
00329 d2fdxdz[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;
00330
00331 p1 = pos + vect3d(0, -res.y, -res.z);
00332 p2 = pos + vect3d(0, -res.y, res.z);
00333 p3 = pos + vect3d(0, res.y, -res.z);
00334 p4 = pos + vect3d(0, res.y, res.z);
00335 getValues (p1, fp1);
00336 getValues (p2, fp2);
00337 getValues (p3, fp3);
00338 getValues (p4, fp4);
00339 h1 = p3.y - p1.y;
00340 h2 = p2.z - p1.z;
00341 h = h1*h2;
00342 for ( size_t i = 0; i < N; i++ )
00343 d2fdydz[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;
00344
00345 p1 = pos + vect3d(-res.x, -res.y, -res.z);
00346 p2 = pos + vect3d(-res.x, -res.y, res.z);
00347 p3 = pos + vect3d(-res.x, res.y, -res.z);
00348 p4 = pos + vect3d(-res.x, res.y, res.z);
00349 p5 = pos + vect3d(res.x, -res.y, -res.z);
00350 p6 = pos + vect3d(res.x, -res.y, res.z);
00351 p7 = pos + vect3d(res.x, res.y, -res.z);
00352 p8 = pos + vect3d(res.x, res.y, res.z);
00353 getValues (p1, fp1);
00354 getValues (p2, fp2);
00355 getValues (p3, fp3);
00356 getValues (p4, fp4);
00357 getValues (p5, fp5);
00358 getValues (p6, fp6);
00359 getValues (p7, fp7);
00360 getValues (p8, fp8);
00361 h1 = p5.x - p1.x;
00362 h2 = p3.y - p1.y;
00363 h3 = p2.z - p1.z;
00364 h = h1*h2*h3;
00365 for ( size_t i = 0; i < N; i++ )
00366 d3fdxdydz[i] = (((fp8[i] - fp7[i]) - (fp6[i] - fp5[i])) - ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i])))/h;
00367 }
00368
00369 private:
00370 vect3d res;
00371 };
00372
00373 #endif // !__utils_hpp__
00374