grid.h

this a Navier-Stokes equations solver. It support grids contains of multiple connected rectangles. S

#include <vector>
#include <list>
#include <assert.h>
#include <math.h>
#include "geometry.h"

#define PRESICION	4

#if PRESICION == 4
#define SCANF_FORMAT	"%f"
typedef float real;
#elif PRESICION == 8
#define SCANF_FORMAT	"%lf"
typedef double real;
#endif


using namespace std;

//enum rect_side{ left = 0, top, right, bottom };
enum Orientation{ left = 0, right = 1, top = 2, bottom = 3 };

struct node
{
	node():u(0),v(0),rho(0),T(0){}

	node(real _vx,real _vy,real den,real t):u(_vx),v(_vy),rho(den),T(t){}

	const real& operator[](int i) const {
		return *(reinterpret_cast<const real *>(this) + i);
	}
	real& operator[](int i){
		return *(reinterpret_cast<real *>(this) + i);
	}

	real velocity() const { return sqrt(u*u + v*v); }
	real rhouv() const { return rho*u*v; }
	real rhou2() const { return rho*u*u; }
	real rhov2() const { return rho*v*v; }
	real p(real R) const { return rho*T*R; }
	real e0(real Cv) const { return Cv*T + (u*u + v*v)/2; }

	real u;
	real v;
	real rho;
	real T;

	static const node dummy;
};

inline node operator - (const node &a,const node &b)
{
	return node(a.u - b.u,a.v - b.v, a.rho - b.rho,a.T - b.T);
}

struct bound_cond
{
	enum type {value = 0, deriv, deriv2};

	bound_cond(){
		t[0] = t[1] = t[2] = t[3] = value;
	}

	node get(const node &a,const node &b = node::dummy) const
	{
		node r;
		for(int i = 0; i < 4; i++)switch(t[i]){
		case value:
			r[i] = n[i];
			break;
		case deriv:
			r[i] = a[i] - d[i];
			break;
		case deriv2:
			assert(false);
			//r[i] = 2.f*b[i] - a[i] + d[i];
			break;
		}
		return r;
	}

	node n;
	node d;
	type t[4];
};

struct grid
{
	grid();
	~grid();

	void set_dims(const rect<int> &r,real _dx,real _dy,const node &init);
	void set_cond(const bound_cond &n);
	void update_minmax();

	const node* from_local_coords(int x,int y) const;
	const node* from_global_coords(int x,int y) const;
	const node* from_global_coords(const point<int> &pt) const
		{ return from_global_coords(pt.x,pt.y); }
	point<int> global_coords(int x,int y) const
		{ return point<int>(bound.a.x + x,bound.a.y + y); }

	node *array;
	node *prev_array;
	node *bound_arrays[4];
	node minvalue, maxvalue;
	real mass;
	int width, height;
	real dx, dy;
	rect<int> bound;
	bound_cond left_cond, right_cond, top_cond, bottom_cond;
};