src/utils.hpp

Go to the documentation of this file.
00001 //
00002 // Ephi - simulation of magnetic fields and particles
00003 // Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
00004 // For more information please see http://www.mare.ee/indrek/ephi/
00005 //
00006 // This program is free software; you can redistribute it and/or modify
00007 // it under the terms of the GNU General Public License as published by
00008 // the Free Software Foundation; either version 2 of the License, or
00009 // (at your option) any later version.
00010 //
00011 // This program is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU General Public License along
00017 // with this program; if not, write to the Free Software Foundation, Inc.,
00018 // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
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&); // we are a dumb class, go away
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 

Generated on Thu Dec 6 20:31:14 2007 for Ephi by  doxygen 1.5.0