main.cpp

//g++ main.cpp rotate.cpp -L/usr/X11R6/lib -lX11 -lXi -lXmu -lGL -lglut -lGLU ../../inmost.a -O5
// press space - explode mesh to see connection 
//#define OCTREECUTCELL_DEBUG
#define DISCR_DEBUG

#include "inmost.h"
#include "my_glut.h"
#include <iostream>
#include <sstream>
#include <algorithm>
#include "rotate.h"
#include <stdarg.h>


using namespace INMOST;
Mesh * mesh;
int drawgrid = 0;
int interactive = 0;
double zoom = 1;
int width = 800, height = 800;
double sleft = 1e20, sright = -1e20, sbottom = 1e20, stop = -1e20, sfar = -1e20, snear = 1e20;
Tag mat,colort; int maxmat = -1, maxcolor = 0;
Mesh::GeomParam table;
double reset_timer = Timer();
double shift[3] = {0,0,0};
int boundary = 2;
int color = 0;
int edges = 1;
bool perspective = false;
int text = 4;
int display_orphans = 1;
int badfaces = 1;
ElementSet problem_set, orphan_edges, orphan_faces;


#if defined(DISCR_DEBUG)
Tag avg_coord, tensor;

static Storage::real dot_prod(Storage::real x[3], Storage::real y[3])  { return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]; }
static void tensor_prod3(Storage::real * K, Storage::real v[3], Storage::real out[3])
{
	out[0] = v[0] * K[0] + v[1] * K[1] + v[2] * K[2];
	out[1] = v[0] * K[3] + v[1] * K[4] + v[2] * K[5];
	out[2] = v[0] * K[6] + v[1] * K[7] + v[2] * K[8];
}

void FindHarmonicPoint(Face & fKL, Cell & cK, Cell & cL, Tag tensor, Storage::real xK[3], Storage::real xL[3], Storage::real y[3], Storage::real & coef)
{
	Storage::real yK[3], yL[3], xKL[3], nKL[3], dK, dL, lK, lL, lKs[3], lLs[3], D, t;
	Storage::real coefK, coefL, coefQ, coefDiv;
	Storage::real_array KK = cK->RealArray(tensor), KL = cL->RealArray(tensor);
	fKL->OrientedUnitNormal(cK,nKL);
	fKL->Centroid(xKL);
	tensor_prod3(&KK[0],nKL,lKs);
	tensor_prod3(&KL[0],nKL,lLs);
	lK = dot_prod(nKL,lKs);
	lL = dot_prod(nKL,lLs);
	lKs[0] -= nKL[0]*lK;
	lKs[1] -= nKL[1]*lK;
	lKs[2] -= nKL[2]*lK;
	lLs[0] -= nKL[0]*lL;
	lLs[1] -= nKL[1]*lL;
	lLs[2] -= nKL[2]*lL;
	D = -dot_prod(xKL,nKL);
	dK = fabs(dot_prod(xK,nKL)+D);
	dL = fabs(dot_prod(xL,nKL)+D);
	yK[0] = xK[0] + dK*nKL[0];
	yK[1] = xK[1] + dK*nKL[1];
	yK[2] = xK[2] + dK*nKL[2];
	yL[0] = xL[0] - dL*nKL[0];
	yL[1] = xL[1] - dL*nKL[1];
	yL[2] = xL[2] - dL*nKL[2];
	coefDiv = (lL*dK+lK*dL);
	coefL = lL*dK/coefDiv;
	coefK = lK*dL/coefDiv;
	coefQ = dK*dL/coefDiv;
	y[0] = yK[0]*coefK + yL[0]*coefL + (lKs[0]-lLs[0])*coefQ;
	y[1] = yK[1]*coefK + yL[1]*coefL + (lKs[1]-lLs[1])*coefQ;
	y[2] = yK[2]*coefK + yL[2]*coefL + (lKs[2]-lLs[2])*coefQ;
	coef = coefL;
}


void FindBoundaryPoint(Face & fK, Cell & cK, Tag tensor, Storage::real xK[3], Storage::real y[3])
{
	Storage::real nK[3], lKs[3], lK, xfK[3], t;
	Storage::real_array KK = cK->RealArray(tensor);
	fK->Centroid(xfK);
	fK->OrientedUnitNormal(cK,nK);
	tensor_prod3(&KK[0],nK,lKs);
	lK = dot_prod(nK,lKs);
	t = (dot_prod(nK,xfK) - dot_prod(nK,xK))/lK;
	y[0] = xK[0] + t*lKs[0];
	y[1] = xK[1] + t*lKs[1];
	y[2] = xK[2] + t*lKs[2];
}
#endif


void printtext(const char * fmt, ... )
{

	unsigned int i;
	char stext[1024];
	va_list ap;
	if ( fmt == NULL ) return;
	va_start(ap,fmt);
	vsprintf(stext,fmt,ap);
	va_end(ap);
	for(i=0;i<strlen(stext);i++)
	{
		glutBitmapCharacter(GLUT_BITMAP_HELVETICA_10,
							stext[i]);
	}
}

void set_matrix3d()
{
	double aspect = (double)width/(double)height;
	double side = std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5;
	//double center[3] = { (sleft+sright)*0.5, (sbottom+stop)*0.5, (sfar+snear)*0.5};
	const double sc = 2;
	glMatrixMode (GL_PROJECTION);
	glLoadIdentity ();
	//~ glOrtho(center[0]-sc*side*zoom*aspect,center[0]+sc*zoom*side*aspect,
		//~ center[1]-sc*side*zoom,center[1]+sc*zoom*side,
		//~ center[2]-sc*side*100,center[2]+sc*side*100);
	if( !perspective )
	{
		glOrtho(-sc*side*zoom*aspect,sc*side*zoom*aspect,
				-sc*side*zoom,sc*side*zoom,
				-sc*side*100,sc*side*100);
	}
	else
		gluPerspective(60.0, aspect, 0.00001, 1000.0);
	glMatrixMode (GL_MODELVIEW);
}

void set_matrix2d()
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-1,1,-1,1,-1,1);
	glMatrixMode(GL_MODELVIEW);
}

void reshape(int w, int h)
{
	width = w;
	height = h;
	set_matrix3d();
	glViewport(0, 0, w, h);
}

class Input
{
public:
	enum InputType {Double,Integer};
private:
	std::string str;
	void * input_link;
	InputType type;
	bool done;


public:
	Input(int * val) {input_link = val; type = Integer; done = false; str = "";}
	Input(double * val) {input_link = val; type = Double; done = false; str = "";}
	Input(void * link, InputType type) : input_link(link), type(type) {done = false; str = "";}
	Input(const Input & other):input_link(other.input_link), str(other.str), type(other.type), done(other.done) {}
	Input & operator =(Input const & other) {input_link = other.input_link; str = other.str; type = other.type; done = other.done; return *this;}
	~Input() {}
	void KeyPress(char c)
	{
		if( c == 13 )
		{

			done = true;
			if( type == Double ) *((double *)input_link) = atof(str.c_str());
			else if( type == Integer ) *((int *)input_link) = atoi(str.c_str());
			glutPostRedisplay();
		}
		else if( c == 8 )
		{
			if( !str.empty() ) str.erase(str.size()-1);
			glutPostRedisplay();
		}
		else if( (c >= '0' && c <= '9') || c=='+' || c=='-' || c=='.' || c=='e' || c == 'E' )
		{
			str += c;
			glutPostRedisplay();
		}
		else if( c == 27 )
		{

			done = true;
			glutPostRedisplay();
		}
	}
	bool Done() {return done;}
	void Draw()
	{
		float h = 24.0f/(float)height;

		glColor3f(1,1,1);
		glBegin(GL_QUADS);
		glVertex3f(-0.99,-0.99,1);
		glVertex3f(-0.99,-0.99+h,1);
		glVertex3f(0.99,-0.99+h,1);
		glVertex3f(0.99,-0.99,1);
		glEnd();
		glColor3f(0,0,0);
		glBegin(GL_LINE_LOOP);
		glVertex3f(-0.99,-0.99,1);
		glVertex3f(-0.99,-0.99+h,1);
		glVertex3f(0.99,-0.99+h,1);
		glVertex3f(0.99,-0.99,1);
		glEnd();

		glColor4f(0,0,0,1);
		glRasterPos2d(-0.985,-0.985);
		//printtext(str.c_str());
		printtext("input number (%s): %s",type == Integer ? "integer":"double", str.c_str());
	}
} * CommonInput = NULL;





void keyboard2(unsigned char key, int x, int y)
{
	if( key == '=' || key == '+' ||  key == '_' || key == '-' || key == 'w' || key == 's' || key == 'a' || key == 'd' || key == 'r' || key == 'p' || key == 'z')
	{
		interactive = false;
		glutPostRedisplay();
	}

}

int matfilter = 0;
int show_cell = -1, show_face = -1;
int parentfilter = -1;
std::vector<int> show_cells, show_faces;

void keyboard(unsigned char key, int x, int y)
{
	if( CommonInput != NULL )
	{
		CommonInput->KeyPress(key);
		return;
	}
	if( key == 't' )
	{
		text = (text+1)%6;
		glutPostRedisplay();
	}
	else if( key == 'b' )
	{
		boundary = (boundary+1)%3;
		glutPostRedisplay();
	}
	else if( key == 'p' )
	{
		perspective = !perspective;
		reshape(width,height);
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 'w' )
	{
		double shiftmod[3] = {0,0,0};
		shiftmod[1] -= 0.03f*expf(zoom-1)*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5;
		rotatevector((double*)shiftmod);
		shift[0] += shiftmod[0];
		shift[1] += shiftmod[1];
		shift[2] += shiftmod[2];
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 's' )
	{
		double shiftmod[3] = {0,0,0};
		shiftmod[1] += 0.03f*expf(zoom-1)*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5;
		rotatevector((double*)shiftmod);
		shift[0] += shiftmod[0];
		shift[1] += shiftmod[1];
		shift[2] += shiftmod[2];
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 'a' )
	{
		double shiftmod[3] = {0,0,0};
		shiftmod[0] += 0.03f*expf(zoom-1)*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5;
		rotatevector((double*)shiftmod);
		shift[0] += shiftmod[0];
		shift[1] += shiftmod[1];
		shift[2] += shiftmod[2];
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 'd' )
	{
		double shiftmod[3] = {0,0,0};
		shiftmod[0] -= 0.03f*expf(zoom-1)*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5;
		rotatevector((double*)shiftmod);
		shift[0] += shiftmod[0];
		shift[1] += shiftmod[1];
		shift[2] += shiftmod[2];
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 'r' )
	{
		shift[0] = 0.0f;
		shift[1] = 0.0f;
		shift[2] = 0.0f;
		zoom = 1;
		quatinit();
		glutPostRedisplay();
		interactive = true;
	}
	else if( key == 'e' )
	{
		edges = !edges;
		glutPostRedisplay();
	}
	else if( key == 27 )
	{
		//delete mesh;
		exit(-1);
	}
	else if( key == '=' || key == '+')
	{
		zoom /= 1.1;
		reshape(width,height);
		glutPostRedisplay();
		interactive = true;
		//reset_timer = Timer();
	}
	else if( key == '_' || key == '-')
	{
		zoom *= 1.1;
		reshape(width,height);
		glutPostRedisplay();
		interactive = true;
		//reset_timer = Timer();
	}
	else if( key == ' ')
	{
		show_cells.clear();
		show_faces.clear();
		//mesh->Save("out.vtk");
		//drawgrid = (drawgrid+1)%2;
		//glutPostRedisplay();
	}
	else if( key == 'q' )
	{
		mesh->PrepareGeometricData(table);
	}
	else if( key == 'w' )
	{
		mesh->RemoveGeometricData(table);
	}
	else if( key-'0' >= 0 && key-'9' <= 0 )
	{
		matfilter = key-'0';
		glutPostRedisplay();
	}
	else if( key == 'm' )
	{
		if( parentfilter == -1 )
		{
			if( CommonInput == NULL ) CommonInput = new Input(&parentfilter);
		}
		else parentfilter = -1;
		glutPostRedisplay();
	}
	else if( key == 'c' )
	{
		if( CommonInput == NULL ) CommonInput = new Input(&show_cell);
		glutPostRedisplay();
	}
	else if( key == 'f' )
	{
		if( CommonInput == NULL ) CommonInput = new Input(&show_face);
		glutPostRedisplay();
	}
	else if( key == 'z' )
	{
		if( !colort.isValid() )
		{
			unsigned visited = 0;
			colort = mesh->CreateTag("COLOR",DATA_INTEGER,CELL,NONE,1);
			std::vector<Cell> queue;
			ElementArray<Element> around;
			queue.reserve(256);
			while(visited != mesh->NumberOfCells() )
			{
				maxcolor++;
				std::cout << "new part " << maxcolor << std::endl;
				for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
				{
					if( it->Integer(colort) == 0 )
					{
						queue.push_back(it->self());
						visited++;
						queue.back()->Integer(colort) = maxcolor;
						break;
					}
				}
				while(!queue.empty())
				{
					around = queue.back()->BridgeAdjacencies(FACE,CELL);

					queue.pop_back();
					for(ElementArray<Element>::iterator it = around.begin(); it != around.end(); ++it)
						if( it->Integer(colort) == 0 )
						{
							queue.push_back(it->getAsCell());
							visited++;
							queue.back()->Integer(colort) = maxcolor;
						}
				}
			}
		}
		color = !color;
	}
	else if( key == 'o' )
	{
		display_orphans = !display_orphans;
		glutPostRedisplay();
	}

	else if( key == 'l' )
	{
		badfaces = !badfaces;
		glutPostRedisplay();
	}
}




class face2gl
{
	double dist;
	double c[4];
	std::vector<double> verts;
public:
	face2gl():verts() {dist = 0; memset(c,0,sizeof(double)*4);}
	face2gl(const face2gl & other) :verts(other.verts) {dist = other.dist; memcpy(c,other.c,sizeof(double)*4);}
	face2gl & operator =(face2gl const & other) { verts = other.verts; dist = other.dist; memmove(c,other.c,sizeof(double)*4); return *this;}
	~face2gl() {}
	void draw() const
	{
		if( boundary == 2 ) glColor3dv(c);
		else glColor4dv(c);
		glBegin(GL_POLYGON);
		for(unsigned k = 0; k < verts.size(); k+=3)
			glVertex3dv(&verts[k]);
		glEnd();
	}
	bool operator <(const face2gl & other) const {return dist < other.dist;}
	void set_color(double r, double g, double b, double a) {c[0] = r; c[1] = g; c[2] = b; c[3] = a;}
	void add_vert(double x, double y, double z) {unsigned s = verts.size(); verts.resize(s+3); verts[s] = x; verts[s+1] = y; verts[s+2] = z;}
	void add_vert(double v[3]) {verts.insert(verts.end(),v,v+3);}
	void set_center(double cnt[3], double cam[3])
	{
		dist = sqrt((cnt[0]-cam[0])*(cnt[0]-cam[0])+(cnt[1]-cam[1])*(cnt[1]-cam[1])+(cnt[2]-cam[2])*(cnt[2]-cam[2]));
	}
};



face2gl DrawFace(Element & f, int mmat, double campos[3])
{
	double cnt[3];
	face2gl ret;
	f->Centroid(cnt);
	ret.set_center(cnt,campos);
	ElementArray<Node> nodes = f->getNodes();

	if( f->nbAdjElements(CELL) == 0 )
		ret.set_color(1,0,0,0.25);
	else if( mmat == 0 )
	{
		ret.set_color(0.2,0,1,0.25);
	}
	else
	{
		double r = color ? mmat / (double)maxcolor : mmat / (double) maxmat;
		ret.set_color(1-r,r,1-r,0.25);
	}

	for( ElementArray<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
		ret.add_vert(&(kt->Coords()[0]));


	if( edges )
	{
		//glLineWidth(2.0);
		if( boundary == 1 ) glColor3f(0,0,1); else glColor3f(0,0,0);
		glBegin(GL_LINE_LOOP);
		for( ElementArray<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
			glVertex3dv(&(kt->Coords()[0]));
		glEnd();
		//glLineWidth(1.0);
	}
	return ret;
}

Tag problem_tag, mats_tag,parent_tag,cell_material_tag;

void fill_mat_str(Storage::integer id, Storage::integer_array mats, char str[2048])
{
	str[0] = '\0';
	sprintf(str,"%d",id);
	for(unsigned q = 0; q < mats.size(); q++)
		sprintf(str,"%s,%d",str,mats[q]+1);
}


void fill_str(Storage::integer id, char str[2048])
{
	str[0] = '\0';
	sprintf(str,"%d",id);
}

void whereami(double & cx, double & cy, double & cz)
{
   // Get the viewing matrix
   GLdouble modelview[16],projection[16];
   GLint viewport[4] = {0,0,1,1};
   glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
   glGetDoublev(GL_PROJECTION_MATRIX, projection);

   GLdouble outx, outy, outz;  // Var's to save the answer in

   gluUnProject(0.5, 0.5, 0.,
               modelview, projection, viewport,
               &outx, &outy, &outz);

   // Return the result.
   cx = outx;
   cy = outy;
   cz = outz;
}

void draw()
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	set_matrix3d();

	if( perspective )
		glTranslated(0,0,-zoom*2*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5);
	//glTranslated((sleft+sright)*0.5,(sbottom+stop)*0.5,(snear+sfar)*0.5);
	rotate();

	//glPointSize(5);
	//glColor3f(0,1,1);
	glLineWidth(3.0);
	glBegin(GL_LINES);
	glColor3f(1,0,0);
	glVertex3d(0,0,0);
	glVertex3d(zoom*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5*0.1,0,0);
	glColor3f(0,1,0);
	glVertex3d(0,0,0);
	glVertex3d(0,zoom*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5*0.1,0);
	glColor3f(0,0,1);
	glVertex3d(0,0,0);
	glVertex3d(0,0,zoom*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.5*0.1);
	glEnd();
	glLineWidth(1.0);
	//glPointSize(1);

	//glTranslated(-(sleft+sright)*0.5+shift[0],-(sbottom+stop)*0.5 + shift[1],-(snear+sfar)*0.5 + shift[2]);

	glTranslated(shift[0],shift[1],shift[2]);

	double campos[3] = {0.5,0.5,0};
	//~ double matrix[16], inverse[16];
	//~ glGetDoublev(GL_MODELVIEW_MATRIX,matrix);
	//~ inverse_v1(inverse,matrix);
	//~ campos[0] = inverse[12];
	//~ campos[1] = inverse[13];
	//~ campos[2] = inverse[14];

	//~ campos[0] = inverse[3];
	//~ campos[1] = inverse[7];
	//~ campos[2] = inverse[11];

	whereami(campos[0],campos[1],campos[2]);


	//glTranslated((l+r)*0.5,(b+t)*0.5,(near+far)*0.5);
	int pacef = std::max(1,mesh->FaceLastLocalID()/10000);

	std::vector<face2gl> polygons;

	if( badfaces )
	{
		for(INMOST_DATA_INTEGER_TYPE it = 0; it < mesh->FaceLastLocalID(); it += (interactive ? pacef : 1)) if( mesh->isValidFace(it) )
		{
			Face f = mesh->FaceByLocalID(it);
			if( !f->CheckNormalOrientation() )
				polygons.push_back(DrawFace(f,1,campos));
		}
	}


	char str[2048];

#if defined(DISCR_DEBUG)
	if(avg_coord.isValid())
	{
		glLineWidth(2.0);
		glColor3f(0,1,0);
		glBegin(GL_LINES);
		for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
		{
			Storage::real cnt[3];
			it->Centroid(cnt);
			ElementArray<Face> faces = it->getFaces();
			for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); ++jt) if(!jt->Boundary())
			{
				Storage::real_array cc = jt->RealArray(avg_coord);
				glVertex3dv(cnt);
				glVertex3dv(&cc[0]);
			}
		}
		glEnd();
		glLineWidth(1.0);
		glPointSize(4.0);
		glColor3f(1,0,0);
		glBegin(GL_POINTS);
		for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) if(!it->Boundary() )
		{
			Storage::real_array cc = it->RealArray(avg_coord);
			glVertex3dv(&cc[0]);
		}
		glEnd();
		glPointSize(1.0);
	}

	if(tensor.isValid())
	{
		glColor3f(0,0,0);
		glBegin(GL_LINES);
		Storage::real h = 0, nh = 0;
		for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) if(!it->Boundary())
		{
			h += (it->BackCell()->Volume() + it->FrontCell()->Volume());
			nh += it->Area();
		}
		h /= nh;
		for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) if(!it->Boundary() )
		{
			Storage::real nrm[3], scale, f1[3], f2[3], l, cnt[3], v[3];
			Cell c1 = it->BackCell(), c2 = it->FrontCell();
			if( c1.isValid() && c2.isValid() )
			{
				it->Centroid(cnt);
				it->OrientedUnitNormal(c1,nrm);
				scale = h;
				Storage::real_array K1 = c1->RealArray(tensor), K2 = c2->RealArray(tensor);
				f1[0] = K1[0]*nrm[0] + K1[1]*nrm[1] + K1[2]*nrm[2];
				f1[1] = K1[3]*nrm[0] + K1[4]*nrm[1] + K1[5]*nrm[2];
				f1[2] = K1[6]*nrm[0] + K1[7]*nrm[1] + K1[8]*nrm[2];

				l = 0.15*scale/sqrt(f1[0]*f1[0]+f1[1]*f1[1]+f1[2]*f1[2]);

				f1[0] *= l;
				f1[1] *= l;
				f1[2] *= l;

				v[0] = f1[0] + cnt[0];
				v[1] = f1[1] + cnt[1];
				v[2] = f1[2] + cnt[2];

				glVertex3dv(cnt);
				glVertex3dv(v);

				v[0] = -f1[0] + cnt[0];
				v[1] = -f1[1] + cnt[1];
				v[2] = -f1[2] + cnt[2];

				glVertex3dv(cnt);
				glVertex3dv(v);

				f2[0] = K2[0]*nrm[0] + K2[1]*nrm[1] + K2[2]*nrm[2];
				f2[1] = K2[3]*nrm[0] + K2[4]*nrm[1] + K2[5]*nrm[2];
				f2[2] = K2[6]*nrm[0] + K2[7]*nrm[1] + K2[8]*nrm[2];



				l = 0.15*scale/sqrt(f2[0]*f2[0]+f2[1]*f2[1]+f2[2]*f2[2]);

				f2[0] *= l;
				f2[1] *= l;
				f2[2] *= l;

				v[0] = f2[0] + cnt[0];
				v[1] = f2[1] + cnt[1];
				v[2] = f2[2] + cnt[2];

				glVertex3dv(cnt);
				glVertex3dv(v);

				v[0] = -f2[0] + cnt[0];
				v[1] = -f2[1] + cnt[1];
				v[2] = -f2[2] + cnt[2];

				glVertex3dv(cnt);
				glVertex3dv(v);
			}
		}

		glEnd();

	}
#endif


	//glColor3f(0,0,0);

	if( !show_cells.empty() )
	{
		for(int qq = 0; qq < show_cells.size(); qq++)
		{
			if( show_cells[qq] < 0 || show_cells[qq] >= mesh->CellLastLocalID() ) continue;
			Element it = mesh->ElementByLocalID(CELL,show_cells[qq]);
			if( it.isValid() )
			{
				Storage::integer_array mats;
				ElementArray<Edge> edges = it->getEdges();
				Storage::real cnt[3];
				for(unsigned k = 0; k < edges.size(); k++)
				{
					edges[k].Centroid(cnt);
					if( mats_tag.isValid() )mats = edges[k].IntegerArray(mats_tag);

					ElementArray<Node> nodes = edges[k].getNodes();

					if( matfilter == 0 || (mats_tag.isValid() && std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{
						if( mats_tag.isValid() ) fill_mat_str(edges[k].LocalID(), mats,str);
						else fill_str(edges[k].LocalID(),str);
						glColor3f(0,0,1);
						glRasterPos3dv(cnt);
						if( text == 2 || text == 5 ) if(text) printtext(str);
						glColor3f(0,0,0);
						glBegin(GL_LINES);
						glVertex3dv(&nodes[0].Coords()[0]);
						glVertex3dv(&nodes[1].Coords()[0]);
						glEnd();
					}

					glColor3f(1,0,0);
					if(mats_tag.isValid())mats = nodes[0].IntegerArray(mats_tag);

					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{

						if( mats_tag.isValid() ) fill_mat_str(nodes[0].LocalID(), mats,str);
						else fill_str(nodes[0].LocalID(),str);
						glRasterPos3dv(&nodes[0].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}

					if(mats_tag.isValid()) mats = nodes[1].IntegerArray(mats_tag);

					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{
						if( mats_tag.isValid() ) fill_mat_str(nodes[1].LocalID(),mats,str);
						else fill_str(nodes[1].LocalID(),str);
						glRasterPos3dv(&nodes[1].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}
				}

				ElementArray<Face> faces = it->getFaces();

				for(unsigned k = 0; k < faces.size(); k++)
				{
					if(mats_tag.isValid())mats = faces[k].IntegerArray(mats_tag);
					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{
						if( faces[k].Boundary() )
							glColor3f(0.15,0.45,0);
						else glColor3f(0,1,0);
						faces[k].Centroid(cnt);
						if( mats_tag.isValid() ) fill_mat_str(faces[k].LocalID(),mats,str);
						else fill_str(faces[k].LocalID(), str);
						glRasterPos3dv(cnt);
						if( text == 3 || text == 5 ) if(text) printtext(str);

						//~ glColor3f(0,0,0);
						//~ Storage::real nrm[3];
						//~ faces[k].OrientedNormal(it->getAsCell(),nrm);
						//~ glLineWidth(2.0);
						//~ glBegin(GL_LINES);
						//~ glVertex3dv(cnt);
						//~ glVertex3d(cnt[0]+nrm[0]*20,cnt[1]+nrm[1]*20,cnt[2]+nrm[2]*20);
						//~ glEnd();
						//~ glLineWidth(1.0);
					}
					//polygons.push_back(DrawFace(&faces[k],0,campos));
				}


				if( matfilter == 0 || (mat.isValid() && it->Integer(mat) == matfilter-1) )
				{
					glColor3f(0,0,0);
					it->Centroid(cnt);
					sprintf(str,"%d",it->LocalID());
					if(parent_tag.isValid()) sprintf(str,"%s, %d",str,it->Integer(parent_tag));
					if(cell_material_tag.isValid()) sprintf(str,"%s, %d",str,it->Integer(cell_material_tag));
					if( problem_tag.isValid() ) sprintf(str,"%s, %d", str, it->Integer(problem_tag));
					glRasterPos3dv(cnt);
					if( text == 4 || text == 5 ) if(text) printtext(str);

					glPointSize(3);
					glColor3f(1,0,0);
					glBegin(GL_POINTS);
					glVertex3dv(cnt);
					glEnd();
					glPointSize(1);
				}
			}
		}
	}
	std::vector<face2gl> inside_face;
	if( !show_faces.empty() )
	{
		for(int qq = 0; qq < show_faces.size(); qq++)
		{
			Element it = InvalidElement();
			if( show_faces[qq] < 0 || show_faces[qq] >= mesh->LastLocalID(FACE) ) continue;

			it = mesh->ElementByLocalID(FACE,show_faces[qq]);
			if( it.isValid() )
			{
				Storage::integer_array mats;
				Storage::real cnt[3], cntf[3];
				ElementArray<Edge> edges = it->getEdges();
				it->Centroid(cntf);
				for(unsigned k = 0; k < edges.size(); k++)
				//if( edges[k].nbAdjElements(CELL) == 0 )
				{
					if(mats_tag.isValid() ) mats = edges[k].IntegerArray(mats_tag);

					ElementArray<Node> nodes = edges[k].getNodes();

					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{
						edges[k].Centroid(cnt);
						if( mats_tag.isValid() ) fill_mat_str(edges[k].LocalID(), mats,str);
						else fill_str(edges[k].LocalID(),str);
						glColor3f(0,0,1);
						glRasterPos3dv(cnt);
						if( text == 2 || text == 5 ) if(text) printtext(str);

						glColor3f(0,0,0);
						glBegin(GL_LINES);
						glVertex3dv(&nodes[0].Coords()[0]);
						glVertex3dv(&nodes[1].Coords()[0]);
						glEnd();

						face2gl f;
						f.set_color(0.6,0.6,0.6,0.1);
						f.add_vert(cntf);
						f.add_vert(&nodes[0].Coords()[0]);
						f.add_vert(&nodes[1].Coords()[0]);
						f.set_center(cntf,campos);
						inside_face.push_back(f);
					}

					glColor3f(1,0,0);
					 if( mats_tag.isValid() ) mats = nodes[0].IntegerArray(mats_tag);
					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{
						if( mats_tag.isValid() ) fill_mat_str(nodes[0].LocalID(), mats,str);
						else fill_str(nodes[0].LocalID(),str);
						glRasterPos3dv(&nodes[0].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}
					if( mats_tag.isValid() ) mats = nodes[1].IntegerArray(mats_tag);
					if( matfilter == 0 || (mats_tag.isValid() &&std::binary_search(mats.begin(),mats.end(),matfilter-1)) )
					{

						if( mats_tag.isValid() ) fill_mat_str(nodes[1].LocalID(),mats,str);
						else fill_str(nodes[1].LocalID(),str);
						glRasterPos3dv(&nodes[1].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}

				}

				if( matfilter == 0 || (mat.isValid() &&it->Integer(mat) == matfilter-1) )
				{
					glColor3f(0,1,0);
					if( mats_tag.isValid() ) fill_mat_str(it->LocalID(),it->IntegerArray(mats_tag),str);
					else fill_str(it->LocalID(),str);
					glRasterPos3dv(cntf);
					if( text == 3 || text == 5 ) if(text) printtext(str);
				}
#if 0
				ElementArray<Cell> cells = it->getCells();

				for(ElementArray<Cell>::iterator jt = cells.begin(); jt != cells.end(); ++jt)
				{
					glColor3f(0,0,0);
					jt->Centroid(cnt);
					fill_str(jt->LocalID(),str);
					glRasterPos3dv(cnt);
					if( text == 4 || text == 5 ) if(text) printtext(str);
				}

	#if defined(DISCR_DEBUG)
				Storage::real cnt0[3], cnt2[3], nrm[3], f1[3],f2[3],l;

				it->Centroid(cnt0);
				/*
				glColor3f(0,0,1);
				glBegin(GL_LINES);
				it->getAsFace()->UnitNormal(nrm);
				tensor_prod3(&it->getAsFace()->BackCell()->RealArrayDF(tensor)[0],
							 nrm,f1);
				l = dot_prod(f1,f1);
				glVertex3dv(cnt0);
				glVertex3d(cnt[0]-f1[0],cnt[1]-f1[1],cnt[2]-f1[2]);

				if(it->getAsFace()->FrontCell() != NULL )
				{
					tensor_prod3(&it->getAsFace()->BackCell()->RealArrayDF(tensor)[0],
								 nrm,f2);
					glVertex3dv(cnt0);
					glVertex3d(cnt[0]+f2[0],cnt[1]+f2[1],cnt[2]+f2[2]);
				}
				glEnd();
				*/
				for(ElementArray<Cell>::iterator jt = cells.begin(); jt != cells.end(); ++jt)
				{
					jt->Centroid(cnt);

					glColor3f(0,0,0);
					glBegin(GL_LINES);
					glVertex3dv(cnt0);
					glVertex3dv(cnt);
					glEnd();

					ElementArray<Face> faces = jt->getFaces();

					for(ElementArray<Face>::iterator qt = faces.begin(); qt != faces.end(); ++qt)
					{
						qt->Centroid(cnt);
						glColor3f(0,0,1);
						glBegin(GL_LINES);
						glVertex3dv(cnt0);
						glVertex3dv(cnt);
						glEnd();
						sprintf(str,"%d",qt->LocalID());
						glRasterPos3dv(cnt);
						printtext(str);

						if( !qt->Boundary() )
						{
							Storage::real stub;
							Cell * xL = qt->FrontCell();
							Cell * xK = qt->BackCell();
							Cell * xQ = xK == &*jt ? xL : xK;
							xQ->Centroid(cnt2);
							FindHarmonicPoint(&*qt,&*jt,xQ,tensor,cnt0,cnt2,cnt,stub);
							glColor3f(1,0,0);
							glBegin(GL_LINES);
							glVertex3dv(cnt0);
							glVertex3dv(cnt);
							glVertex3dv(cnt);
							glVertex3dv(cnt2);
							glEnd();

							sprintf(str,"%d",xQ->LocalID());
							glRasterPos3dv(cnt2);
							printtext(str);
						}
						else
						{
							Cell * xK = qt->BackCell();
							xK->Centroid(cnt2);
							FindBoundaryPoint(&*qt,xK,tensor,cnt2,cnt);

							glColor3f(0,1,0);
							glBegin(GL_LINES);
							glVertex3dv(cnt0);
							glVertex3dv(cnt);
							glEnd();

							sprintf(str,"%d",qt->LocalID());
							glRasterPos3dv(cnt);
							printtext(str);
						}
					}
				}

				glEnd();
	#endif
#endif
				//mats = it->IntegerArray(mats_tag);
				//if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
				//	polygons.push_back(DrawFace(&*it,0,campos));
			}
		}
	}
	std::sort(inside_face.rbegin(),inside_face.rend());

	glEnable(GL_BLEND);

	for(int q = 0; q < inside_face.size() ; q++) inside_face[q].draw();

	glDisable(GL_BLEND);

	//for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); it++)
#if defined(OCTREECUTCELL_DEBUG)
	if( show_cell == -1 && show_face == -1 )
	{
		int end = mesh->MaxLocalID(CELL);
		//for(int i = 0; i < end; ++i)
		if( problem_set != NULL ) for(ElementSet::iterator iit = problem_set->begin(); iit != problem_set->end(); ++iit)
		{
			Element * it = &*iit;//mesh->ElementByLocalID(CELL,i);
			//if( it == NULL ) continue;
			//if(it->HaveData(problem_tag) && it->Integer(problem_tag) != 0 )
			{
				if( parentfilter != -1 && parent_tag.isValid() && it->Integer(parent_tag) != parentfilter ) continue;
				Storage::integer_array mats;
				ElementArray<Edge> edges = it->getEdges();
				Storage::real cnt[3];
				for(unsigned k = 0; k < edges.size(); k++)
				{
					edges[k].Centroid(cnt);
					mats = edges[k].IntegerArray(mats_tag);

					ElementArray<Node> nodes = edges[k].getNodes();

					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					{
						fill_mat_str(edges[k].LocalID(), mats,str);
						glColor3f(0,0,1);
						glRasterPos3dv(cnt);
						if( text == 2 || text == 5 ) if(text) printtext(str);
						glColor3f(0,0,0);
						glBegin(GL_LINES);
						glVertex3dv(&nodes[0].Coords()[0]);
						glVertex3dv(&nodes[1].Coords()[0]);
						glEnd();
					}

					glColor3f(1,0,0);
					mats = nodes[0].IntegerArray(mats_tag);

					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					{

						fill_mat_str(nodes[0].LocalID(), mats,str);
						glRasterPos3dv(&nodes[0].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}

					mats = nodes[1].IntegerArray(mats_tag);

					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					{
						fill_mat_str(nodes[1].LocalID(),mats,str);
						glRasterPos3dv(&nodes[1].Coords()[0]);
						if( text == 1 || text == 5 ) if(text) printtext(str);
					}
				}

				ElementArray<Face> faces = it->getFaces();

				for(unsigned k = 0; k < faces.size(); k++)
				{
					mats = faces[k].IntegerArray(mats_tag);
					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					{
						if( faces[k].Boundary() )
							glColor3f(0.15,0.45,0);
						else glColor3f(0,1,0);
						faces[k].Centroid(cnt);
						fill_mat_str(faces[k].LocalID(),mats,str);
						glRasterPos3dv(cnt);
						if( text == 3 || text == 5 ) if(text) printtext(str);


						//~ glColor3f(0,0,0);
						//~ Storage::real nrm[3];
						//~ faces[k].OrientedNormal(it->getAsCell(),nrm);
						//~ glLineWidth(3.0);
						//~ glBegin(GL_LINES);
						//~ glVertex3dv(cnt);
						//~ glVertex3d(cnt[0]+nrm[0]*20,cnt[1]+nrm[1]*20,cnt[2]+nrm[2]*20);
						//~ glEnd();
						//~ glLineWidth(1.0);
					}
					//polygons.push_back(DrawFace(&faces[k],0,campos));
				}


				if( matfilter == 0 || it->Integer(mat) == matfilter-1 )
				{
					glColor3f(0,0,0);
					it->Centroid(cnt);
					sprintf(str,"%d, %d, %d",it->LocalID(),it->Integer(parent_tag),it->Integer(cell_material_tag)+1);
					if( problem_tag.isValid() ) sprintf(str,"%s, %d", str, it->Integer(problem_tag));
					glRasterPos3dv(cnt);
					if( text == 4 || text == 5 ) if(text) printtext(str);

					glPointSize(3);
					glColor3f(1,0,0);
					glBegin(GL_POINTS);
					glVertex3dv(cnt);
					glEnd();
					glPointSize(1);
				}






			}
		}
		if( display_orphans )
		{
			//end = mesh->MaxLocalID(EDGE);
			//for(Mesh::iteratorEdge it = mesh->BeginEdge(); it != mesh->EndEdge(); it++)
			//for(int i = 0; i < end; i++)
			if( orphan_edges != NULL ) for(ElementSet::iterator iit = orphan_edges.Begin(); iit != orphan_edges.End(); ++iit)
			{
				Element * it = &*iit;//mesh->ElementByLocalID(EDGE,i);
				//if( it == NULL ) continue;
				//if( it->nbAdjElements(CELL) == 0 )
				{
					Storage::integer_array mats = it->IntegerArray(mats_tag);
					Storage::real cnt[3];

					ElementArray<Node> nodes = it->getNodes();

					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					{
						it->Centroid(cnt);
						fill_mat_str(it->LocalID(), mats,str);
						glColor3f(1,0,1);
						glRasterPos3dv(cnt);
						if( text == 2 || text == 5 ) if(text) printtext(str);


						glBegin(GL_LINES);
						glVertex3dv(&nodes[0].Coords()[0]);
						glVertex3dv(&nodes[1].Coords()[0]);
						glEnd();


						glColor3f(1,0,1);
						mats = nodes[0].IntegerArray(mats_tag);


					}


					//~ glColor3f(1,0,0);
					//~ mats = nodes[0].IntegerArray(mats_tag);
					//~ fill_mat_str(nodes[0].LocalID(), mats,str);
					//~ glRasterPos3dv(&nodes[0].Coords()[0]);
					//~ printtext(str);
					//~
					//~ mats = nodes[1].IntegerArray(mats_tag);
					//~ fill_mat_str(nodes[1].LocalID(),mats,str);
					//~ glRasterPos3dv(&nodes[1].Coords()[0]);
					//~ printtext(str);

				}
			}
			//end = mesh->MaxLocalID(FACE);
			//for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); it++)
			//for(int i = 0; i < end; i++)
			if( orphan_faces != NULL ) for(ElementSet::iterator iit = orphan_faces->begin(); iit != orphan_faces->end(); ++iit)
			{
				Element * it = &*iit;//mesh->ElementByLocalID(FACE,i);
				//if( it == NULL ) continue;
				//if( it->nbAdjElements(CELL) == 0 )
				{

					Storage::integer_array mats;
					Storage::real cnt[3];
					ElementArray<Edge> edges = it->getEdges();

					for(unsigned k = 0; k < edges.size(); k++)
					//if( edges[k].nbAdjElements(CELL) == 0 )
					{
						mats = edges[k].IntegerArray(mats_tag);

						ElementArray<Node> nodes = edges[k].getNodes();

						if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
						{
							edges[k].Centroid(cnt);
							fill_mat_str(edges[k].LocalID(), mats,str);
							glColor3f(0,0,1);
							glRasterPos3dv(cnt);
							if( text == 2 || text == 5 ) if(text) printtext(str);

							glColor3f(0,0,0);
							glBegin(GL_LINES);
							glVertex3dv(&nodes[0].Coords()[0]);
							glVertex3dv(&nodes[1].Coords()[0]);
							glEnd();
						}

						glColor3f(1,0,0);
						mats = nodes[0].IntegerArray(mats_tag);
						if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
						{
							fill_mat_str(nodes[0].LocalID(), mats,str);
							glRasterPos3dv(&nodes[0].Coords()[0]);
							if( text == 1 || text == 5 ) if(text) printtext(str);
						}

						if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
						{
							mats = nodes[1].IntegerArray(mats_tag);
							fill_mat_str(nodes[1].LocalID(),mats,str);
							glRasterPos3dv(&nodes[1].Coords()[0]);
							if( text == 1 || text == 5 ) if(text) printtext(str);
						}
					}



					mats = it->IntegerArray(mats_tag);

					if( text == 3 || text == 5 )
					{
						it->Centroid(cnt);
						glColor3f(0,0.5,0.5);
						glRasterPos3dv(cnt);
						fill_mat_str(it->LocalID(), mats,str);
						if(text) printtext(str);
					}

			//~ EDGE 647252
					if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
						polygons.push_back(DrawFace(&*it,0,campos));
				}
			}
		}
	}
#endif
	if( boundary )
	{
		for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype) )
		{
			int end = mesh->LastLocalID(etype);
			//for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); it++)
			if( colort.isValid() )
			{
				for(int i = 0; i < end; i+= interactive ? pacef : 1)
				{
					Element it = mesh->ElementByLocalID(etype,i);
					if( !it.isValid() ) continue;
					if( it->Boundary() )
					{
						int mmat = color ? (!it->getAsFace()->BackCell().isValid() ? 0 : it->getAsFace()->BackCell()->Integer(colort)) : (!it->getAsFace()->BackCell().isValid() ? 0 : it->getAsFace()->BackCell()->Integer(mat));
						if( matfilter == 0 || matfilter-1 == mmat )
							polygons.push_back(DrawFace(it->self(),mmat,campos));
					}
				}
			}
			else
			{
				for(int i = 0; i < end; i+= interactive ? pacef : 1)
				{
					Element it = mesh->ElementByLocalID(etype,i);
					if( !it.isValid() ) continue;
					if( it->Boundary() )
						polygons.push_back(DrawFace(it->self(),0,campos));
				}
			}
		}
	}

	if( display_orphans )
	{
		if( orphan_edges.isValid()) for(ElementSet::iterator iit = orphan_edges.Begin(); iit != orphan_edges.End(); ++iit)
		{
			Element it = iit->self();//mesh->ElementByLocalID(EDGE,i);
			//if( it == NULL ) continue;
			//if( it->nbAdjElements(CELL) == 0 )
			{
				Storage::real cnt[3];
				
				ElementArray<Node> nodes = it->getNodes();

				{
					glBegin(GL_LINES);
					glVertex3dv(&nodes[0].Coords()[0]);
					glVertex3dv(&nodes[1].Coords()[0]);
					glEnd();

				}
			}
		}
		if( orphan_faces.isValid() ) for(ElementSet::iterator iit = orphan_faces.Begin(); iit != orphan_faces.End(); ++iit)
		{
			Element it = iit->self();
			{

				Storage::integer_array mats;
				Storage::real cnt[3];
				ElementArray<Edge> edges = it->getEdges();

				for(unsigned k = 0; k < edges.size(); k++)
				{
					ElementArray<Node> nodes = edges[k].getNodes();
					{
						glBegin(GL_LINES);
						glVertex3dv(&nodes[0].Coords()[0]);
						glVertex3dv(&nodes[1].Coords()[0]);
						glEnd();
					}
				}


				if( matfilter == 0 || std::binary_search(mats.begin(),mats.end(),matfilter-1) )
					polygons.push_back(DrawFace(it->self(),0,campos));
			}
		}
	}


	std::sort(polygons.begin(),polygons.end());


	//glDisable(GL_DEPTH_TEST);
	//if( boundary == 1 ) glEnable(GL_BLEND);
	//else glDisable(GL_BLEND);
	glEnable(GL_BLEND);

	for(int q = polygons.size(); q > 0 ; q--) polygons[q-1].draw();

	glDisable(GL_BLEND);
//~
	//~ for(int q = 0; q < polygons.size() ; q++)
		//~ polygons[q].draw();

	//glEnable(GL_DEPTH_TEST);