mech.cpp

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

#include <assert.h>
#include <float.h>
#include "mech.h"

#define MAX_VALS		16

mech::simple_grid::simple_grid(const grid &g)
:a(g.array),width(g.width),height(g.height),bound(g.bound),
 minvalue(g.minvalue),maxvalue(g.maxvalue),dx(g.dx),dy(g.dy)
{
}

mech::mech()
:grids(0),links(0),count(0)
{
	R  = 0.1;
	mu = 100;
	Cv = 0.1;
	lambda = 0.1;

	dx = 4e-2f;
	dy = 4e-2f;
	dt = 5e-7f;
	total_mass = 0;
}

mech::~mech()
{
	if(grids)delete[] grids;
	if(links)delete[] links;
}

int mech::construct_grids( list<interim_rect> &g,const list<interim_bound> &b,
						  const bound_cond &defbounds,const node &init ) 
{
	if(grids)delete[] grids;
	if(links)delete[] links;

	int idx, subidx, nodes_total = 0;
	list<interim_rect>::iterator i, j;

	count = g.size();
	grids = new grid[count];
	links = new int[count*count];
	memset(links,0,sizeof(int)*count*count);

	// calculate bounding rect for mech
	rect<real> phis_bounds = *g.begin();
	for(i = g.begin(); i != g.end(); i++)
	{
		if(phis_bounds.a.x > i->a.x)phis_bounds.a.x = i->a.x;
		if(phis_bounds.a.y > i->a.y)phis_bounds.a.y = i->a.y;
		if(phis_bounds.b.x < i->b.x)phis_bounds.b.x = i->b.x;
		if(phis_bounds.b.y < i->b.y)phis_bounds.b.y = i->b.y;
	}

	bounds.a.x = 0;
	bounds.a.y = 0;
	bounds.b.x = int((phis_bounds.b.x - phis_bounds.a.x)/dx);
	bounds.b.y = int((phis_bounds.b.y - phis_bounds.a.y)/dy);

	// setup grids
	for(i = g.begin(), idx = 0; i != g.end(); i++, idx++){
		rect<int> r;
		r.a.x = (i->a.x - phis_bounds.a.x) / dx;
		r.a.y = (i->a.y - phis_bounds.a.y) / dy;
		r.b.x = (i->b.x - phis_bounds.a.x) / dx;
		r.b.y = (i->b.y - phis_bounds.a.y) / dy;
		grids[idx].set_dims(r,dx,dy,init);
		grids[idx].set_cond(defbounds);
		nodes_total += grids[idx].width * grids[idx].height;
	}

	for(list<interim_bound>::const_iterator i = b.begin(); i != b.end(); i++){
		if(i->id >= 0 || i->id < count)switch(i->o){
		case left:
			grids[i->id].left_cond = *i;
			break;
		case right:
			grids[i->id].right_cond = *i;
			break;
		case top:
			grids[i->id].top_cond = *i;
			break;
		case bottom:
			grids[i->id].bottom_cond = *i;
			break;
		default:
			printf("mech::construct_grids warning: invalid bound cond orientation %d\n",i->o);
			break;
		}else printf("mech::construct_grids warning: invalid bound id %d\n",i->id);
	}

	// make table of links
	for(i = g.begin(), idx = 0; i != g.end(); i++, idx++){
		for(j = i, subidx = idx; j != g.end(); j++, subidx++)
		{
			if(i == j)continue;
			links[idx*count + subidx] =
			links[subidx*count + idx] =
				(i->inclusive_inside(j->a)
				|| i->inclusive_inside(j->b)
				|| i->inclusive_inside(point<real>(j->a.x,j->b.y))
				|| i->inclusive_inside(point<real>(j->b.x,j->a.y))
				|| j->inclusive_inside(i->a)
				|| j->inclusive_inside(i->b)
				|| j->inclusive_inside(point<real>(i->a.x,i->b.y))
				|| j->inclusive_inside(point<real>(i->b.x,i->a.y)));
		}
	}

	printf("mech::construct_grids completed size %d nodes\n",nodes_total);

	return 0;
}

list<const char *> parse(char *s)
{
	list<const char *> words;
	do{
		while(*s && iswspace(*s))s++;
		if(*s){
			words.push_back(s);
			while(*s && !iswspace(*s))s++;
			if(*s)*s++ = 0;
		}
	}while(*s);
	return words;
}

int parse_floats(list<const char *> &words,real *vals,int maxcount)
{
	int i;
	list<const char *>::iterator it;
	for(i = 0, it = words.begin();i < maxcount && it != words.end();i++, it++){
		sscanf(*it,SCANF_FORMAT,vals + i);
	}
	return i;
}

int parse_boundaries(list<const char *> &words,bound_cond *bc)
{
	int i;
	list<const char *>::iterator it;
	for(i = 0, it = words.begin();i < 4 && it != words.end();i++, it++){
		const char *word = *it;
		if(*word == 'd'){
			sscanf(word + 1,SCANF_FORMAT,&bc->d[i]);
			bc->t[i] = bound_cond::deriv;
		}else if(*word == 'D'){
			sscanf(word + 1,SCANF_FORMAT,&bc->d[i]);
			bc->t[i] = bound_cond::deriv2;
		}else{
			sscanf(word,SCANF_FORMAT,&bc->n[i]);
			bc->t[i] = bound_cond::value;
		}
	}
	return i;
}

int mech::read_file(const char *filename)
{
	FILE *f = fopen(filename,"rt");
	if(!f){
		printf("mech::read_file() fail to open %s file\n",filename);
		return -1;
	}
	node init;
	bound_cond default_bounds;
	list<interim_rect> g;
	list<interim_bound> bounds;
	int node_i = 0, lineno = 0;

	while(!feof(f)){
		char s[0x100];
		real vals[MAX_VALS];

		fgets(s,sizeof(s),f);
		if(!strncmp(s,"mu",2)){
			__asm nop;
		}
		list<const char *> words = parse(s);
		if(words.empty()){
			lineno++;
			continue;
		}
		const char *name = words.front();words.pop_front();

		if(*name == ';'){
			lineno++;
			continue;
		}

		if(!strcmp(name,"rect")){
			parse_floats(words,vals,MAX_VALS);
			interim_rect x;
			x.id = (int)vals[0];
			x.a.x = vals[1];
			x.a.y = vals[2];
			x.b.x = vals[3];
			x.b.y = vals[4];
			if(x.a.x > x.b.x)swap(x.a.x,x.b.x);
			if(x.a.y > x.b.y)swap(x.a.y,x.b.y);
			g.push_back(x);
		}else if(!strcmp(name,"bound")){
			parse_boundaries(words,&default_bounds);
		}else if(!strcmp(name,"init")){
			parse_floats(words,vals,MAX_VALS);
			init.u = vals[0];
			init.v = vals[1];
			init.rho = vals[2];
			init.T = vals[3];
		}else if(!strcmp(name,"bound_left")){
			interim_bound b;
			b.id = atoi(words.front());
			b.o = left;
			words.pop_front();
			parse_boundaries(words,&b);
			bounds.push_back(b);
		}else if(!strcmp(name,"bound_right")){
			interim_bound b;
			b.id = atoi(words.front());
			b.o = right;
			words.pop_front();
			parse_boundaries(words,&b);
			bounds.push_back(b);
		}else if(!strcmp(name,"bound_top")){
			interim_bound b;
			b.id = atoi(words.front());
			b.o = top;
			words.pop_front();
			parse_boundaries(words,&b);
			bounds.push_back(b);
		}else if(!strcmp(name,"bound_bottom")){
			interim_bound b;
			b.id = atoi(words.front());
			b.o = bottom;
			words.pop_front();
			parse_boundaries(words,&b);
			bounds.push_back(b);
		}else if(!strcmp(name,"dx")){
			sscanf(words.front(),SCANF_FORMAT,&dx);
		}else if(!strcmp(name,"dy")){
			sscanf(words.front(),SCANF_FORMAT,&dy);
		}else if(!strcmp(name,"dt")){
			sscanf(words.front(),SCANF_FORMAT,&dt);
		}else if(!strcmp(name,"lambda")){
			sscanf(words.front(),SCANF_FORMAT,&lambda);
		}else if(!strcmp(name,"mu")){
			sscanf(words.front(),SCANF_FORMAT,&mu);
		}else if(!strcmp(name,"R")){
			sscanf(words.front(),SCANF_FORMAT,&R);
		}else if(!strcmp(name,"Pr")){
			
		}else{
			printf("mech::read_file() line %d warning: unknown parameter %s\n",lineno,name);
		}
		lineno++;
	}
	fclose(f);
	if(g.size() > 0){
		construct_grids(g,bounds,default_bounds,init);
	}
	return lineno;
}

int mech::save_to_file(const char *filename)
{
	FILE *f = fopen(filename,"wb");
	if(!f)return -1;

	for(int i = 0; i < count; i++){
		fprintf(f,"grid %d\n",i);
		fprintf(f,"bound %f %f %f %f\n",
			grids[i].bound.a.x, grids[i].bound.a.y,
			grids[i].bound.b.x, grids[i].bound.b.y);
		fprintf(f,"dim %d %d\n",grids[i].width,grids[i].height);
		for(int y = 0; y < grids[i].height; y++){
			for(int x = 0; x < grids[i].width; x++){
				node *n = grids[i].array + x + y*grids[i].width;
				fprintf(f,"%f %f %f %f ",n->u,n->v,n->rho,n->T);
			}
			fprintf(f,"\n");
		}
	}

	fclose(f);
	return 0;
}

list<mech::simple_grid> mech::lock_grids()
{
	list<simple_grid> res;

	for(int i = 0; i < count; i++)res.push_back(grids[i]);

	return res;
}

void mech::release_grids()
{

}

const node& mech::at( int x,int y ) 
{
	for(int i = 0; i < count; i++)
		if(grids[i].bound.inclusive_inside(point<int>(x,y)))
			return *grids[i].from_global_coords(x,y);
	return node::dummy;
}

const node* mech::get_bound(const grid *g,int x,int y)
{
	int *line = links + count*(g - grids);
	point<int> pt = g->global_coords(x,y);
	for(size_t i = 0;i < count; i++){
		if(line[i] && grids[i].bound.inside(pt)){
			//if(x == -1 && y == 25)
			//	__asm nop;
			const node *v = grids[i].from_global_coords(pt);
			//printf("point %d %d get from %d rect T=%f\n",x,y,i,v->T);
			return v;
		}
	}
	return g->from_local_coords(x,y);
}

real mech::solve()