versor  3.0
C++11 library for Geometric algebra
vsr_field.h
1 //
2 // vsr_field.h
3 // Versor
4 //
5 //
6 // Created by Pablo Colapinto on 11/27/12.
7 // Copyright (c) 2012 __MyCompanyName__. All rights reserved.
8 
9 // TO DO: define on any lattice, allow for read/write from/into values
10 //
11 
12 #ifndef Versor_vsr_field_h
13 #define Versor_vsr_field_h
14 
15 #include "vsr_cubicLattice.h"
16 #include "util/vsr_math.h"
17 #include "gfx/gfx_patch.h"
18 
19 namespace vsr{
20 
21  #define ITN \
22  for (int i = 0; i < this->mNum; ++i) {
23 
24  #define END }
25 
26  #define ITER \
27  for (int i = 0; i < this->mWidth; ++i){\
28  for (int j = 0; j < this->mHeight; ++j){\
29  for (int k = 0; k < this->mDepth; ++k){\
30  int tidx = this->idx(i,j,k);
31 
32  #define ITEND }}}
33 
34  #define ITERED(x) \
35  for (int i = 0; i < x.w(); ++i){\
36  for (int j = 0; j < x.h(); ++j){\
37  for (int k = 0; k < x.d(); ++k){\
38  int tidx = idx(i,j,k);
39 
40  #define ITERV \
41  for (int i = 0; i < this->mWidth; ++i){\
42  double u = 1.0 * i/this->mWidth;\
43  for (int j = 0; j < this->mHeight; ++j){\
44  double v = 1.0 * j/this->mHeight;\
45  for (int k = 0; k < this->mDepth; ++k){\
46  double w = 1.0 * k/this->mDepth; \
47  int tidx = this->idx(i,j,k);
48 
49 
50  #define BOUNDITER\
51  for (int i = 1; i < this->mWidth-1; ++i){ \
52  for (int j = 1; j < this->mHeight-1; ++j){ \
53  for (int k = 1; k < this->mDepth-1; ++k){ \
54 
55  #define BOUNDED(x)\
56  for (int i = 1; i < x.w()-1; ++i){ \
57  for (int j = 1; j < x.h()-1; ++j){ \
58  for (int k = 1; k < x.d()-1; ++k){ \
59 
60  #define BOUNDEND \
61  }}}
62 
65 
66  template < class T >
67  class Field : public CubicLattice < typename T::space::point > {
68 
69  protected:
70 
71  T * mData;
72 
73  using GridType = typename T::space::point;
74 
75  public:
76 
77  //Vector Derivative in Euclidean Metric
78  //typedef typename ProductN<VEC,T::idx>::GP VecDeriv;
79 
80  T * dataPtr() const { return mData; }
81  void dataPtr( T* d ) { mData = d; }
82 
84  void zero() { ITER mData[tidx] = T(0); ITEND }
85 
86  Field( int w=1, int h=1, int d=1, double spacingX = 1.0, double spacingY = 1.0, double spacingZ = 1.0) :
87  CubicLattice<GridType>(w,h,d,spacingX,spacingY,spacingZ), mData( NULL )
88 
89  {
90  alloc();
91  init();
92  }
93 
94  Field& resize( int w, int h, int d, double spacingX = 1.0, double spacingY = 1.0, double spacingZ=1.0){
95 
96  CubicLattice<GridType>::resize(w,h,d,spacingX,spacingY,spacingZ);
97 
98  alloc(); init();
99  return *this;
100  }
101 
102  // Field& respace( double sx, ){
103  // this->spacing(s);
104  // init();
105  // return *this;
106  // }
107 
108  void onDestroy(){
109  if (mData) delete[] mData;
110  }
111 
113  void alloc(){
114  if (mData) delete[] mData;
115  mData = new T[this->mNum];
116  }
117 
119  T& operator [] (int i) { return ( mData[i] ); }
121  T operator [] (int i) const { return mData[i]; }
123  T& at(int w = 0, int h = 0, int d = 0) { return mData[ this->idx(w, h, d) ]; }
125  T at(int w = 0, int h = 0, int d = 0) const { return mData[ this->idx(w, h, d) ]; }
126 
127  //INITIALIZE
128  void basicInit(){
129  ITER
130  mData[ tidx ] = this->mPoint[tidx].template copy<T>();//T( this->mPoint[ tidx ] );
131  ITEND
132  }
133  //SPECIALIZE HERE
134  void init(){
135  basicInit();
136  }
137 
138  void reset() { init(); }
139 
141 // Vec rand(){
142 // return Vec( px(Rand::Int(mWidth)) , py(Rand::Int(mHeight)), pz(Rand::Int(mDepth)) );
143 // }
144 //
145 
146  //data at bottom front corner of vxl at p
147  template<class B>
148  T& dataAt( const B& p ){
149  int idx = vxlAt(p).a;
150  return mData[idx];
151  }
152 
153  //2d Euler (bounded)
154  template<class V>
155  T euler2d( const V& v){
156  V t = this->range2D(v); //bind to range of field
157  return surf(t[0], t[1]);
158  }
159 
160  //3d Euler (bounded) v an arbitrary point in space
161  template<class V>
162  T euler3d( const V& v) const{
163  V t= this->range(v);
164  return vol(t[0], t[1], t[2]);
165  }
166 
167 
169  T surf(VSR_PRECISION u, VSR_PRECISION v){
170  Patch p = this->surfIdx(u,v);
171 
172  T a = mData[ p.a ];//at ( iw, ih, 0 ) ;
173  T b = mData[ p.b ];//at ( iw + 1, ih, 0 );
174  T c = mData[ p.c ];//at ( iw + 1, ih + 1, 0 );
175  T d = mData[ p.d ];//at ( iw, ih + 1, 0 );
176 
177  return Interp::surface<T> (a,b,c,d, p.rw, p.rh);
178  }
179 
180  T vol(VSR_PRECISION u, VSR_PRECISION v, VSR_PRECISION w) const {
181  VPatch p = this->vidx(u,v,w);
182 
183  T a = mData[ p.a ]; T b = mData[ p.b ]; T c = mData[ p.c ]; T d = mData[ p.d ];
184  T e = mData[ p.e ]; T f = mData[ p.f ]; T g = mData[ p.g ]; T h = mData[ p.h ];
185  return Interp::volume<T> (a,b,c,d,e,f,g,h, p.rw, p.rh, p.rd );
186  }
187 
188  template<class A>
189  T vol(const A& a){ return vol(a[0],a[1],a[2]); }
190 
191  //Contour Integral
192  template<class V>
193  vector<V> contour(const V& v, int num, double force){
194  vector<V> vp;
195  V tv = v;
196  for (int i = 0; i < num; ++i){
197  vp.push_back(tv);
198  tv += euler3d(tv) * force;
199  }
200  return vp;
201  }
202 
203 
205 // T quadSurf(double u, double v){
206 // Patch p = surfIdx(u,v);
207 //
208 // }
209 
210 
211  //a 6-faced Kernel (3D 'plus" sign)
212 
213  T sumNbrs (int idx) const{
214  T tdx;
215  for (int i = 1; i < 7; ++i){
216  if (this->nbr(idx)[i] != -1 ) tdx += mData[ this->nbr(idx)[i] ];
217  }
218  return tdx;
219  }
220 
221  T diffNbrs (int idx) const {
222  T tdx;
223  tdx += mData[ this->nbr(idx).xr ] - mData [ this->nbr(idx).xl ]; //xr - xl
224  tdx += mData[ this->nbr(idx).yt ] - mData [ this->nbr(idx).yb ]; //yt - yb
225  tdx += mData[ this->nbr(idx).zb ] - mData [ this->nbr(idx).zf ]; //zb - zf
226  return tdx;
227  }
228 
229 
230  T diffXNbrs (int ix) const{
231  return mData[ this->nbr(ix).xr ] - mData [ this->nbr(ix).xl ];;
232  }
233 
234  T diffYNbrs ( int ix) const{
235  return mData[ this->nbr(ix).yt ] - mData [ this-> nbr(ix).yb ];
236  }
237 
238  T diffZNbrs (int ix) const{
239  return mData[ this->nbr(ix).zb ] - mData [ this->nbr(ix).zf ]; //xr - xl
240  }
241 
242  T dx(int ix) const{ return diffXNbrs(ix); } //or . . .
243  T dy(int ix) const{ return diffYNbrs(ix); }
244  T dz(int ix) const{ return diffZNbrs(ix); }
245 
246  //Create Vector Derivative field of another Field f. . . e1 * dx + e2 * dy + e3 * dz . . . in other dimensions use reciprocal frame
247  // template<class B>
248  // Field& der(const Field<B>& f){
249  // BOUNDITER
250  // int ix = idx(i,j,k);
251  // typename decltype(Vec() * T()) tdx;
252  // tdx += Vec::x * f.diffXNbrs(ix) * (mWidth-2);
253  // tdx += Vec::y * f.diffYNbrs(ix) * (mHeight-2);
254  // tdx += Vec::z * f.diffZNbrs(ix) * (mDepth-2);
255  // mData[ix] = tdx;
256  // BOUNDEND
257  // return *this;
258  // }
259 
260 
261  //Scalar Tensor
262  double tensNbrs (int idx) const {
263  double tdx;
264  tdx += mData[ this->nbr(idx).xr ][0] - mData [ this->nbr(idx).xl ][0]; //xr - xl
265  tdx += mData[ this->nbr(idx).yt ][1] - mData [ this->nbr(idx).yb ][1]; //yt - yb
266  tdx += mData[ this->nbr(idx).zb ][2] - mData [ this->nbr(idx).zf ][2]; //zb - zf
267  return tdx;
268  }
269 
270  //Scalar Tensor Weighted By 1.0/dim (? spacing * ?) vec for now . . .
271  double tensNbrsWt(int idx) const {
272  double tdx; double ww = 1.0 / this->w(); double wh = 1.0/this->h(); double wd = 1.0/this->d();
273  tdx += ( mData[ this->nbr(idx).xr ][0] - mData [ this->nbr(idx).xl ][0] ) * ww; //xr - xl
274  tdx += ( mData[ this->nbr(idx).yt ][1] - mData [ this->nbr(idx).yb ][1] ) * wh; //yt - yb
275  tdx += ( mData[ this->nbr(idx).zb ][2] - mData [ this->nbr(idx).zf ][2] ) * wd; //zb - zf
276  return tdx;
277  }
278 
279 
280  //Scalar Tensor "reverse weighted" By dim (? spacing * ?)
281  double tensNbrsRwt(int idx) const {
282  double tdx;
283  tdx += ( mData[ this->nbr(idx).xr ][0] - mData [ this->nbr(idx).xl ][0] ) * this-> w(); //xr - xl
284  tdx += ( mData[ this->nbr(idx).yt ][1] - mData [ this->nbr(idx).yb ][1] ) * this-> h(); //yt - yb
285  tdx += ( mData[ this->nbr(idx).zb ][2] - mData [ this->nbr(idx).zf ][2] ) * this-> d(); //zb - zf
286  return tdx;
287  }
288 
290  void gsSolver(const Field& prev){
291  //Iterative Pressure Solver substracts pressure tensor out
292  int it = 20;
293  for (int m = 0; m < it; ++m){
294  BOUNDITER
295  int ix = this->idx(i,j,k);
296  T td = sumNbrs(ix);
297  mData[ix] = (prev[ix] + td ) / 6.0;
298  BOUNDEND
299 
300  boundaryConditions(0);
301  }
302  }
303 
305  void diffuse(const Field& prev, double diffRate, bool bounded, bool ref){
306 
307  static int it = 20;
308 
309  double rate = diffRate * .001 * this->mNum;
310  if (bounded) {
311  //field is padded ? bounded
312  //iterate
313  for (int k = 0; k < it; ++k){
314  //for each data member
315  for (int i = 0; i < this->num(); ++i){
316 
317  //sum up neighboring values
318  T td = sumNbrs(i);
319  //multiply by (rate / 1 + 6*rate ) (or multiply afterwards?
320  td *= rate ;
321  //add to old value
322  mData[i] = ( prev[i] + td ) / (1 + 6*rate);
323  }
324  }
325  //break;
326  } else {
327  //unbounded so set bounds in loop
328  //iterate
329  for (int n = 0; n < it; ++n){
330  //cout << "BOUNDED:";
331  BOUNDITER
332  int tdx = this->idx(i,j,k);
333  T td = sumNbrs(tdx);
334  //multiply by (rate )
335  td *= rate;
336  //add to old value and divide new result by (1 + 6 * rate)
337  mData[tdx] = ( prev[tdx] + td ) / (1 + 6*rate);
338  BOUNDEND
339 
340  boundaryConditions(ref);
341  }
342  //break;
343  }
344  }
345 
347  template<class B>
348  void advect(const Field& prev, const Field<B>& f, double dt, bool ref){
349 
350  double dt0 = dt;// * mWidth;
351  BOUNDITER
352  int tidx = this->idx(i,j,k);
353  auto p = this->mPoint[ tidx ] + f.euler3d( this->mPoint[tidx] ) * -dt0;// .trs( f.euler3d( mPoint[tidx] ) * -dt0 );//f[tidx] * -dt0 ); //Lattice Point
354  mData[tidx] = prev.euler3d( p );
355  BOUNDEND
356 
357  boundaryConditions(ref);
358  }
359 
360  //backwards Twist Advection?
361 
362  // Sets values of a scalar field to divergence of another Field b
363  template <class B>
364  Field& div(const Field<B>& f){
365  //Sum Differences of each Vxl Face (= DIVERGENCE TENSOR)
366  BOUNDITER
367  int ix = this->idx(i,j,k);
368  mData[ix] = f.tensNbrsWt(ix) * ( -.5 ); //sca
369  BOUNDEND
370  boundaryConditions(0);
371  return *this;
372  }
373 
374  //swap ptrs to data with another field f of same type
375  Field& swap( Field& f ){
376  T * tmp = dataPtr(); dataPtr( f.dataPtr() ); f.dataPtr( tmp );
377  return *this;
378  }
379 
380  //Operators
381  template<class N>
382  Field& operator * ( N val ) { ITN mData[i] *= val; END return *this; }
383 
384  //Operators Pairwise substraction
385  /* template< Bits::Type ... XS > */
386  /* Field& operator -= ( const Field< MV<XS...> >& f ) { ITN mData[i] -= f[i]; END return *this; } */
387 
388 
389  //here written fro a scalar fiedl: specialize below for vectors etcs
390  void boundaryConditions(bool ref){
391 
392  for (int i = 0; i < this->mFace.size(); ++i){
393  int ix = this->mFace[i];
394  static Nbr n;
395  n = this->mNbr[ ix ];
396  int type = n.type;
397 
398  mData[ix] = T(0);
399 
400  //negation for now treat as vectors
401  if (type & LEFT) mData[ix] += ref ? -mData[n.xr] : mData[n.xr];
402  if (type & RIGHT) mData[ix] += ref ? -mData[n.xl] : mData[n.xl];
403  if (type & TOP) mData[ix] += ref ? -mData[n.yb] : mData[n.yb];
404  if (type & BOTTOM) mData[ix] += ref ? -mData[n.yt] : mData[n.yt];
405  if (type & BACK) mData[ix] += ref ? -mData[n.zf] : mData[n.zf];
406  if (type & FRONT) mData[ix] += ref ? -mData[n.zb] : mData[n.zb];
407 
408  mData[ix] /= Math::bitcount(type);
409  }
410 
411 
412  }
413 
414 
415 
416  };
417 
418 // //SPECIALIZATIONS
419 
420  // template<> void Field<Vec>::boundaryConditions(bool ref){
421  //
422  // for (int i = 0; i < mFace.size(); ++i){
423  // int ix = mFace[i];
424  // static Nbr n;
425  // n = mNbr[ ix ];
426  // int type = n.type;
427  //
428  // //negation for now treat as vectors
429  // if (type & LEFT) mData[ix][0] = ref ? -mData[n.xr][0] : mData[n.xr][0];
430  // if (type & RIGHT) mData[ix][0] = ref ? -mData[n.xl][0] : mData[n.xl][0];
431  // if (type & TOP) mData[ix][1] = ref ? -mData[n.yb][1] : mData[n.yb][1];
432  // if (type & BOBits::TypeOM) mData[ix][1] = ref ? -mData[n.yt][1] : mData[n.yt][1];
433  // if (type & BACK) mData[ix][2] = ref ? -mData[n.zf][2] : mData[n.zf][2];
434  // if (type & FRONT) mData[ix][2] = ref ? -mData[n.zb][2] : mData[n.zb][2];
435  // }
436  // }
437  //
438  //
439  // template<> void Field< Vec > :: init() {
440  // ITER mData[tidx] = Vec(mPoint[tidx]).unit(); ITEND
441  // }
442  //
443  // template<> void Field< Sca > :: init() {
444  // ITN mData[i] = Sca(0); END
445  // }
446  //
447  // template<> void Field< Dll > :: init(){
448  // ITN mData[i] = Frame( mPoint[i] ).dll(); END
449  // }
450  //
451  // template<> void Field< Frame > :: init(){
452  // //printf("FRAME FIELD INIT\n");
453  // ITN mData[i].pos() = mPoint[i]; END
454  // }
455 
456 
457 // template<> void Field< Vec > :: draw(float r, float g, float b, float a) {
458 // //cout <<"vecdraw"<<endl;
459 // drawPush(r,g,b,a);
460 // }
461 //
462 //
463 // template<> void Field< Sca > :: draw(float r, float g, float b, float a) {
464 // ITN GL::push(); glColor4f(r,g,b,mData[i][0] * a); GL::translate(mPoint[i].w()); GL::Glyph::Cube( mSpacing ); GL::pop(); END
465 // }
466 // template<> Frame Field< Frame > :: surf(double u, double v){
467 //
468 // Patch p = surfIdx(u,v);
469 //
470 // Dll a = mData[p.a].dll();
471 // Dll b = mData[p.b].dll();
472 // Dll c = mData[p.c].dll();
473 // Dll d = mData[p.d].dll();
474 //
475 // return Gen::mot( Interp::surface<Dll> (a,b,c,d, p.rw, p.rh) );
476 // }
477 
478 
479 // template <>
480 // double Field < Frame > :: tensNbrs (int idx) {
481 // double tdx;
482 // tdx += mData[ nbr(idx).xr ].pos()[0] - mData [ nbr(idx).xl ].pos()[0]; //xr - xl
483 // tdx += mData[ nbr(idx).yt ].pos()[1] - mData [ nbr(idx).yb ].pos()[1]; //yt - yb
484 // tdx += mData[ nbr(idx).zb ].pos()[2] - mData [ nbr(idx).zf ].pos()[2]; //zb - zf
485 // return tdx;
486 // }
487 
488 
489 
490 
491 
492 // template<>
493 // void Field<Dll> :: init(){
494 // ITER
495 // mGrid[ tidx ] = Ro::null( px(i), py(j), pz(k) );
496 // mData[ tidx ] = T( mGrid[ tidx ] );
497 // ITEND
498 // }
499 
500 // template< class T >
501 // void Field< T > ::
502 
503 
504 
505  //still working on this
506  /*
507  template <class T>
508  class Fluid {
509 
510  Field<T> mField; // Current Vector Veloctiy field
511  Field<T> mPrev; // Previously evaluated Vector field
512 
513  Field<Sca> mDensity; // Current Density Field
514  Field<Sca> mDensityPrev; // Previously Evaluated Density Field
515 
516  Field<Sca> mPressure;
517  Field<Sca> mDivergence;
518 
519  //Vector Derivative of Scalar Pressure Field is a Vector Field
520  Field< typename Sca::VecDeriv > mDerivative;
521 
522  public:
523 
524  Fluid(int _w, int _h, int _d)
525  : mField(_w,_h,_d), mPrev(_w,_h,_d),
526  mPressure(_w,_h,_d), mDivergence(_w,_h,_d),
527  mDensity(_w*2,_h*2,_d*2,.5), mDensityPrev(_w*2,_h*2,_d*2,.5),
528  mDerivative(_w,_h,_d)
529  {
530 
531  }
532 
533  void reset() {
534 
535  mField.reset(); mPrev.reset();
536  mDensity.reset(); mDensityPrev.reset();
537  mPressure.reset();
538  mDivergence.reset();
539  }
540 
541 
542  Field<T>& velocity() { return mField; }
543  Field<Sca>& density() { return mDensity; }
544  Field<Sca>& pressure() { return mPressure; }
545 
546  //incompressibility of fluids
547  void project(){
548 
549  //Reset Pressure field
550  mPressure.zero();
551  mDerivative.zero();
552  mDivergence.zero();
553 
554  //Calculate Divergence tensor
555  mDivergence.div( mField );
556 
557  //Iterative Gauss-Seidel Solver using Divergence Field
558  mPressure.gsSolver( mDivergence );
559 
560  //Substract out Vector Derivative Gradient of Scalar Pressure
561  mField -= mDerivative.der( mPressure ) * .5;
562 
563  //Satisfy Boundary Conditions
564  mField.boundaryConditions(1);
565 
566  }
567 
568  void step(double visc, double diff, double advectRate){
569 
570  //VECTOR VELOCITY FIELD
571  mField.swap(mPrev);
572 
573  mField.diffuse( mPrev, visc, 0, 1);
574 
575  project();
576 
577  mField.swap(mPrev);
578 
579  //self-advection
580  mField.advect( mPrev, mPrev, advectRate, 1);
581 
582  project();
583 
584  //DENSITY FIELD
585  mDensity.swap(mDensityPrev);
586  mDensity.diffuse(mDensityPrev, diff, 0, 0);
587  mDensity.swap(mDensityPrev);
588 
589  mDensity.advect(mDensityPrev, mField, advectRate, 0);
590  }
591 
592  };
593 
594  */
595  #undef ITER
596  #undef ITEND
597  #undef BOUNDITER
598  #undef BOUNDEND
599  #undef ITN
600  #undef END
601 
602 } //vsr::
603 
604 #endif
605 
606 
void gsSolver(const Field &prev)
Guass Siedel Relaxation Solver Using a Previous Field State.
Definition: vsr_field.h:290
VPatch vidx(double u, double v, double w) const
Volume Index at u,v,w [0-1].
Definition: vsr_cubicLattice.h:412
void diffuse(const Field &prev, double diffRate, bool bounded, bool ref)
Backwards Diffusion Using a Previous Field State.
Definition: vsr_field.h:305
T & operator[](int i)
Set Data by Index.
Definition: vsr_field.h:119
Volume Patch Info Container for Euler integration of a 3d Field.
Definition: vsr_patch.h:54
A Basic 3D Field (slowly porting this over from the now defunct vsr_lattice class) Use to Evaluate Ne...
Definition: vsr_field.h:67
double td() const
Total Depth.
Definition: vsr_cubicLattice.h:144
Info Container for Euler integration of a 2d Field.
Definition: vsr_patch.h:44
T & at(int w=0, int h=0, int d=0)
Set Data by Coordinate.
Definition: vsr_field.h:123
void alloc()
Allocate Memory.
Definition: vsr_field.h:113
void advect(const Field &prev, const Field< B > &f, double dt, bool ref)
Backwards Advection Using a Previous Field State prev and Based on a Velocity Frame f...
Definition: vsr_field.h:348
T at(int w=0, int h=0, int d=0) const
Get Data by Coordinate.
Definition: vsr_field.h:125
Discretized Volume Indexing (Isometric Cubic Lattice w/o data)
Definition: vsr_cubicLattice.h:61
int idx(int i, int j, int k=0) const
Get Index In Array.
Definition: vsr_cubicLattice.h:113
T & dataAt(const B &p)
Random Vector in Field.
Definition: vsr_field.h:148
the versor library namespace
Definition: vsr_algebra.h:29
Patch surfIdx(double u, double v)
Indicex of surface at u, v [0, 1].
Definition: vsr_cubicLattice.h:389
Vxl vxlAt(const V &tv) const
Voxel of Vector v.
Definition: vsr_cubicLattice.h:219
T sumNbrs(int idx) const
Get QUADRIC Interpolated Data at eval u,v [0-1].
Definition: vsr_field.h:213
void zero()
Zero Out All Data.
Definition: vsr_field.h:84
T surf(VSR_PRECISION u, VSR_PRECISION v)
Get BILINEAR Interpolated Data at eval u,v [0-1.0].
Definition: vsr_field.h:169