#ifndef NOTCL

#include <mystdlib.h>

#include <myadt.hpp>
#include <meshing.hpp>
#include <csg.hpp>
#include <stlgeom.hpp>

// #include <parallel.hpp>

#include <visual.hpp>

namespace netgen
{
  extern AutoPtr<Mesh> mesh;
  extern STLGeometry * stlgeometry;
  VisualSceneMesh vsmesh;

  VisualSceneMesh :: VisualSceneMesh ()
    : VisualScene()
  {
    filledlist = 0;
    linelist = 0;
    edgelist = 0;
    badellist = 0;
    tetlist = 0;
    prismlist = 0;
    hexlist = 0;
    pyramidlist = 0;
    identifiedlist = 0;
    pointnumberlist = 0;
    domainsurflist = 0;

    vstimestamp = GetTimeStamp();
    selecttimestamp = GetTimeStamp();
    filledtimestamp = GetTimeStamp();
    linetimestamp = GetTimeStamp();
    edgetimestamp = GetTimeStamp();
    pointnumbertimestamp = GetTimeStamp();

    tettimestamp = GetTimeStamp();
    prismtimestamp = GetTimeStamp();
    hextimestamp = GetTimeStamp();
    pyramidtimestamp = GetTimeStamp();

    badeltimestamp = GetTimeStamp();
    identifiedtimestamp = GetTimeStamp();
    domainsurftimestamp = GetTimeStamp();


    selface = -1;
    selelement = -1;
    locpi = 1;
    selpoint = -1;
    selpoint2 = -1;
    seledge = -1;
  }

  VisualSceneMesh :: ~VisualSceneMesh ()
  {
    ;
  }


  void VisualSceneMesh :: DrawScene ()
  {
    if (!mesh)
      {
	VisualScene::DrawScene();
	return;
      }

    lock = NULL;

    clock_t starttime, endtime;
    starttime = clock();

    BuildScene();

    glClearColor(backcolor, backcolor, backcolor, 1.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glEnable (GL_COLOR_MATERIAL);
    glColor3f (1.0f, 1.0f, 1.0f);
    glLineWidth (1.0f);

    SetLight();

    glPushMatrix();
    glMultMatrixf (transformationmat);

    GLdouble projmat[16];
    glGetDoublev (GL_PROJECTION_MATRIX, projmat);


#ifdef PARALLEL
   glEnable (GL_BLEND);
   glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
#endif


    glInitNames ();
    glPushName (0);

    //    glEnable (GL_LINE_SMOOTH);
//         glEnable (GL_BLEND);
//         glEnable (GL_POLYGON_SMOOTH);
//         glDisable (GL_DEPTH_TEST);
//         glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    //    glHint (GL_LINE_SMOOTH_HINT, GL_DONT_CARE);

    glDisable (GL_COLOR_MATERIAL);

    GLfloat matcol0[] = { 0, 0, 0, 1 };
    GLfloat matcol1[] = { 1, 1, 1, 1 };
    GLfloat matcolf[] = { 0, 1, 0, 1 };
    GLfloat matcolb[] = { 0.5, 0, 0, 1 };
    // GLfloat matcolblue[] = { 0, 0, 1, 1 };

    glMatrixMode (GL_MODELVIEW);

    glMaterialfv(GL_FRONT, GL_EMISSION, matcol0);
    glMaterialfv(GL_BACK, GL_EMISSION, matcol0);
    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, matcol1);
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolf);
    glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, matcolb);

    glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);



    // glPolygonOffset (1,10);
    glPolygonOffset (2,2);
    glEnable (GL_POLYGON_OFFSET_FILL);

    SetClippingPlane ();

    if (vispar.drawfilledtrigs)
      {
	if (filledtimestamp < mesh->GetTimeStamp () ||
	    filledtimestamp < selecttimestamp)
	  {
	    BuildFilledList ();
	  }


#ifdef PARALLELGL
	if (ntasks > 1 && vispar.drawtetsdomain > 0 && vispar.drawtetsdomain < ntasks)
	  glCallList (par_filledlists[vispar.drawtetsdomain]);
	else
#endif
	  glCallList (filledlist);
      }

    if (vispar.drawbadels)
      glCallList (badellist);

    if (vispar.drawprisms)
      {
	BuildPrismList ();
	glCallList (prismlist);
      }

    if (vispar.drawpyramids)
      {
	BuildPyramidList ();
	glCallList (pyramidlist);
      }

    if (vispar.drawhexes)
      {
	BuildHexList ();
	glCallList (hexlist);
      }

    if (vispar.drawtets)
      {
	BuildTetList ();
	glCallList (tetlist);
      }

    if (vispar.drawdomainsurf)
      {
	BuildDomainSurfList();
	glCallList (domainsurflist);
      }

    glDisable (GL_POLYGON_OFFSET_FILL);

    // draw lines

    glMatrixMode (GL_MODELVIEW);

    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcol0);
    glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, matcol0);
    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, matcol0);

    glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
    glLineWidth (1.0f);
    glColor3f (0.0f, 0.0f, 0.0f);
    glDisable (GL_LINE_SMOOTH);





    if (vispar.drawoutline)
      {
	glPolygonOffset (1, 1);
	glEnable (GL_POLYGON_OFFSET_LINE);

	if (linetimestamp < mesh->GetTimeStamp ())
	  BuildLineList ();

#ifdef PARALLELGL
	if (ntasks > 1 && vispar.drawtetsdomain > 0 && vispar.drawtetsdomain < ntasks)
	  // for (int dest = 1; dest < ntasks; dest++)
	  // if (vispar.drawtetsdomain == dest)
	  glCallList (par_linelists[vispar.drawtetsdomain]);
	else
#endif
	  glCallList (linelist);


	glDisable (GL_POLYGON_OFFSET_LINE);
      }



    if (vispar.drawidentified)
      {
	glPolygonOffset (1, -1);
	glEnable (GL_POLYGON_OFFSET_LINE);
	glCallList (identifiedlist);
	glDisable (GL_POLYGON_OFFSET_LINE);
      }

    if (vispar.drawpointnumbers ||
	vispar.drawedgenumbers ||
	vispar.drawfacenumbers ||
	vispar.drawelementnumbers)
      glCallList (pointnumberlist);


    glPopName();

    if (vispar.drawedges)
      {
	BuildEdgeList();
	glCallList (edgelist);
      }


    if (selpoint > 0 && selpoint <= mesh->GetNP())
      {
	/*
	glPointSize (3.0);
	glColor3d (0, 0, 1);
	glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolblue);
	glBegin (GL_POINTS);

	const Point3d p = mesh->Point(selpoint);
	glVertex3f (p.X(), p.Y(), p.Z());
	glEnd();
	*/

	glColor3d (0, 0, 1);

	static GLubyte cross[] =
	  {
	    0xc6, 0xee, 0x7c, 0x38, 0x7c, 0xee, 0xc6
	  };
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);

	glDisable (GL_COLOR_MATERIAL);
	glDisable (GL_LIGHTING);
	glDisable (GL_CLIP_PLANE0);

	const Point3d p = mesh->Point(selpoint);
	glRasterPos3d (p.X(), p.Y(), p.Z());
	glBitmap (7, 7, 3, 3, 0, 0, &cross[0]);
      }


    glDisable(GL_CLIP_PLANE0);

    glPopMatrix();

    if (vispar.colormeshsize)
      DrawColorBar (minh, maxh, 1);

    DrawCoordinateCross ();
    DrawNetgenLogo ();

    if (lock)
      {
	lock -> UnLock();
	delete lock;
	lock = NULL;
      }

    glFinish();


    endtime = clock();
    // cout << 1.0 / (double(endtime - starttime)/CLOCKS_PER_SEC) << " frames/sec" << endl;
  }


  void VisualSceneMesh :: BuildScene (int zoomall)
  {
    if (!mesh)
      {
	VisualScene::BuildScene (zoomall);
	return;
      }

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    int i, j;


    Point3d pmin, pmax;
    static double oldrad = 0;

    Array<Element2d> faces;

    int meshtimestamp = mesh->GetTimeStamp();
    if (meshtimestamp > vstimestamp || zoomall)
      {
        if (mesh->GetDimension() == 2)
          {
            // works in NGSolve, mesh view
            mesh->GetBox (pmin, pmax);
          }
        else
          {
            // otherwise strange zooms douring mesh generation
            mesh->GetBox (pmin, pmax, SURFACEPOINT);
          }

	if (vispar.use_center_coords && zoomall == 2)
	  {
	    center.X() = vispar.centerx; center.Y() = vispar.centery; center.Z() = vispar.centerz;
	  }
	else if (selpoint >= 1 && zoomall == 2)
	  center = mesh->Point (selpoint);
	else if (vispar.centerpoint >= 1 && zoomall == 2)
	  center = mesh->Point (vispar.centerpoint);
	else
	  center = Center (pmin, pmax);
	rad = 0.5 * Dist (pmin, pmax);
	if(rad == 0) rad = 1e-6;

	if (rad > 1.2 * oldrad ||
	    mesh->GetMajorTimeStamp() > vstimestamp ||
	    zoomall)
	  {
	    CalcTransformationMatrices();
	    oldrad = rad;
	  }
      }

    glEnable (GL_NORMALIZE);

    if (pointnumberlist)
      {
	glDeleteLists (pointnumberlist, 1);
	pointnumberlist = 0;
      }

    if (badellist)
      {
	glDeleteLists (badellist, 1);
	badellist = 0;
      }
    /*
      if (prismlist)
      {
      glDeleteLists (prismlist, 1);
      prismlist = 0;
      }

      if (pyramidlist)
      {
      glDeleteLists (pyramidlist, 1);
      pyramidlist = 0;
      }

      if (hexlist)
      {
      glDeleteLists (hexlist, 1);
      hexlist = 0;
      }
    */
    if (identifiedlist)
      {
	glDeleteLists (identifiedlist, 1);
	identifiedlist = 0;
      }

    pointnumberlist = glGenLists (1);
    glNewList (pointnumberlist, GL_COMPILE);

    if (vispar.drawpointnumbers ||
	vispar.drawedgenumbers ||
	vispar.drawfacenumbers ||
	vispar.drawelementnumbers)
      {
	//     	glEnable (GL_COLOR_MATERIAL);
	GLfloat textcol[3] = { 1 - backcolor,
			       1 - backcolor,
			       1 - backcolor };
	glColor3fv (textcol);
	glNormal3d (0, 0, 1);
	glPushAttrib (GL_LIST_BIT);
	// glListBase (fontbase);

	char buf[30];

	if (vispar.drawpointnumbers)
	  for (i = 1; i <= mesh->GetNP(); i++)
	    {
	      const Point3d & p = mesh->Point(i);
	      glRasterPos3d (p.X(), p.Y(), p.Z());

	      sprintf (buf, "%d", i);

	      // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
	      MyOpenGLText (buf);
	    }

	if (vispar.drawedgenumbers)
	  {
	    /*
	    for (SegmentIndex i = 0; i < mesh->GetNSeg(); i++)
	      {
		const Segment & seg = (*mesh)[i];

		const Point3d & p1 = mesh->Point(seg[0]);
		const Point3d & p2 = mesh->Point(seg[1]);
		const Point3d p = Center (p1, p2);
		glRasterPos3d (p.X(), p.Y(), p.Z());

		sprintf (buf, "%d", seg.edgenr);
		glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
	      }
	    */

	    const MeshTopology & top = mesh->GetTopology();
	    for (i = 1; i <= top.GetNEdges(); i++)
	      {
		int v1, v2;
		top.GetEdgeVertices (i, v1, v2);
		const Point3d & p1 = mesh->Point(v1);
		const Point3d & p2 = mesh->Point(v2);
		const Point3d p = Center (p1, p2);
		glRasterPos3d (p.X(), p.Y(), p.Z());

		sprintf (buf, "%d", i);
		// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
		MyOpenGLText (buf);

	      }

	  }


	if (vispar.drawfacenumbers)
	  {
	    const MeshTopology & top = mesh->GetTopology();
	    Array<int> v;
	    for (i = 1; i <= top.GetNFaces(); i++)
	      {
		top.GetFaceVertices (i, v);
		const Point3d & p1 = mesh->Point(v.Elem(1));
		const Point3d & p2 = mesh->Point(v.Elem(2));
		const Point3d & p3 = mesh->Point(v.Elem(3));
		Point3d p;
		if (v.Elem(4) == 0)
		  {
		    p = Center (p1, p2, p3);
		  }
		else
		  {
		    const Point3d & p4 = mesh->Point(v.Elem(4));
		    Point3d hp1 = Center (p1, p2);
		    Point3d hp2 = Center (p3, p4);
		    p = Center (hp1, hp2);
		  }

		glRasterPos3d (p.X(), p.Y(), p.Z());
		sprintf (buf, "%d", i);
		// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
		MyOpenGLText (buf);
	      }
	  }



	if (vispar.drawelementnumbers)
	  {
	    Array<int> v;
	    for (i = 1; i <= mesh->GetNE(); i++)
	      {
		// const ELEMENTTYPE & eltype = mesh->ElementType(i);
		Array<int> pnums;

		Point3d p;
		const Element & el = mesh->VolumeElement (i);

		if ( ! el.PNum(5)) //  eltype == TET )
		  {

		    pnums.SetSize(4);
		    for( int j = 0; j < pnums.Size(); j++)
		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);


		    const Point3d & p1 = mesh->Point(pnums[0]);
		    const Point3d & p2 = mesh->Point(pnums[1]);
		    const Point3d & p3 = mesh->Point(pnums[2]);
		    const Point3d & p4 = mesh->Point(pnums[3]);
		    p = Center (p1, p2, p3, p4);
		  }
		else if ( ! el.PNum(6)) // eltype == PYRAMID
		  {
		    pnums.SetSize(5);
		    for( int j = 0; j < pnums.Size(); j++)
		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);


		    const Point3d & p1 = mesh->Point(pnums[0]);
		    const Point3d & p2 = mesh->Point(pnums[1]);
		    const Point3d & p3 = mesh->Point(pnums[2]);
		    const Point3d & p4 = mesh->Point(pnums[3]);
		    const Point3d & p5 = mesh->Point(pnums[4]);

		    p.X()  = 0.3 * p5.X() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . X();
		    p.Y()  = 0.3 * p5.Y() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . Y();
		    p.Z()  = 0.3 * p5.Z() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . Z();

		  }
		else if ( ! el.PNum(7) ) // eltype == PRISM
		  {
		    pnums.SetSize(6);
		    for( int j = 0; j < pnums.Size(); j++)
		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);

		    const Point3d & p1 = mesh->Point(pnums[0]);
		    const Point3d & p2 = mesh->Point(pnums[1]);
		    const Point3d & p3 = mesh->Point(pnums[2]);
		    const Point3d & p11 = mesh->Point(pnums[3]);
		    const Point3d & p12 = mesh->Point(pnums[4]);
		    const Point3d & p13 = mesh->Point(pnums[5]);
		    p = Center (  Center (p1, p2, p3) , Center(p11, p12, p13) )  ;

		  }
		else if (! el.PNum(9) ) // eltype == HEX
		  {
		    pnums.SetSize(8);
		    for( int j = 0; j < pnums.Size(); j++)
		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);

      		    const Point3d & p1 = mesh->Point(pnums[0]);
		    const Point3d & p2 = mesh->Point(pnums[1]);
		    const Point3d & p3 = mesh->Point(pnums[2]);
		    const Point3d & p4 = mesh->Point(pnums[3]);
		    const Point3d & p5 = mesh->Point(pnums[4]);
		    const Point3d & p6 = mesh->Point(pnums[5]);
		    const Point3d & p7 = mesh->Point(pnums[6]);
		    const Point3d & p8 = mesh->Point(pnums[7]);

		    p = Center ( Center ( Center(p1, p3), Center(p2, p4) ) , Center( Center(p5, p7) , Center(p6, p8 ) ) );
		  }

		glRasterPos3d (p.X(), p.Y(), p.Z());
		sprintf (buf, "%d", i);
		// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
		MyOpenGLText (buf);

	      }
	  }


	glPopAttrib ();
	//      	glDisable (GL_COLOR_MATERIAL);
      }
    glEndList ();






    badellist = glGenLists (1);
    glNewList (badellist, GL_COMPILE);

    if (vispar.drawbadels)
      {
	//  SetClippingPlane ();

	static float badelcol[] = { 1.0f, 0.0f, 1.0f, 1.0f };
	glLineWidth (1.0f);

	for (i = 1; i <= mesh->GetNE(); i++)
	  {
	    if (mesh->VolumeElement(i).flags.badel ||
		mesh->VolumeElement(i).flags.illegal ||
		(i == vispar.drawelement))
	      {
		// copy to be thread-safe
		Element el = mesh->VolumeElement (i);
		el.GetSurfaceTriangles (faces);

		glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, badelcol);


		//	  if ( (el.GetNP() == 4) || (el.GetNP() == 10))
		if (el.PNum(1))
		  {
		    glBegin (GL_TRIANGLES);

		    for (j = 1; j <= faces.Size(); j++)
		      {
			Element2d & face = faces.Elem(j);
			const Point3d & lp1 = mesh->Point (el.PNum(face.PNum(1)));
			const Point3d & lp2 = mesh->Point (el.PNum(face.PNum(2)));
			const Point3d & lp3 = mesh->Point (el.PNum(face.PNum(3)));
			Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
			n /= (n.Length()+1e-12);
			glNormal3d (n.X(), n.Y(), n.Z());
			glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
			glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
			glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
		      }

		    glEnd();
		  }
	      }
	  }



	for (i = 1; i <= mesh->GetNE(); i++)
	  {
	    if (mesh->VolumeElement(i).flags.badel)
	      {
		// copy to be thread-safe
		Element el = mesh->VolumeElement (i);
		if ( (el.GetNP() == 4) || (el.GetNP() == 10))
		  {
		    glBegin (GL_LINES);
		    glVertex3d (0,0,0);
		    const Point3d & p = mesh->Point(el.PNum(1));
		    glVertex3d (p.X(), p.Y(), p.Z());
		    glEnd();
		  }
	      }
	  }


	for (i = 1; i <= mesh->GetNE(); i++)
	  {
	    Element el = mesh->VolumeElement (i);
	    int hascp = 0;
	    for (j = 1; j <= el.GetNP(); j++)
	      if (el.PNum(j) == vispar.centerpoint)
		hascp = 1;

	    if (hascp)
	      {
		(*testout) << "draw el " << i << " : ";
		for (j = 1; j <= el.GetNP(); j++)
		  (*testout) << el.PNum(j) << " ";
		(*testout) << endl;

		if (el.GetNP() == 4)
		  {
		    int et[6][2] =
		      { { 1, 2 },
			{ 1, 3 },
			{ 1, 4 },
			{ 2, 3 },
			{ 2, 4 },
			{ 3, 4 } } ;

		    for (j = 0; j < 6; j++)
		      {
			glBegin (GL_LINES);
			const Point3d & p1 = mesh->Point (el.PNum(et[j][0]));
			const Point3d & p2 = mesh->Point (el.PNum(et[j][1]));
			glVertex3d (p1.X(), p1.Y(), p1.Z());
			glVertex3d (p2.X(), p2.Y(), p2.Z());
			glEnd ();
		      }
		  }


		if (el.GetNP() == 10)
		  {
		    int et[12][2] =
		      { { 1, 5 },
			{ 2, 5 },
			{ 1, 6 },
			{ 3, 6 },
			{ 1, 7 },
			{ 4, 7 },
			{ 2, 8 },
			{ 3, 8 },
			{ 2, 9 },
			{ 4, 9 },
			{ 3, 10 },
			{ 4, 10 } };

		    for (j = 0; j < 12; j++)
		      {
			glBegin (GL_LINES);
			const Point3d & p1 = mesh->Point (el.PNum(et[j][0]));
			const Point3d & p2 = mesh->Point (el.PNum(et[j][1]));
			glVertex3d (p1.X(), p1.Y(), p1.Z());
			glVertex3d (p2.X(), p2.Y(), p2.Z());
			glEnd ();
		      }
		  }
	      }
	  }


	for (i = 1; i <= mesh->GetNSE(); i++)
	  {
	    Element2d el = mesh->SurfaceElement(i);
	    if (!el.BadElement())
	      continue;

	    int drawel = 1;
	    for (j = 1; j <= el.GetNP(); j++)
	      if (!el.PNum(j))
		drawel = 0;

	    if (!drawel)
	      continue;

	    cout << int (el.GetType()) << " " << flush;
	    switch (el.GetType())
	      {
	      case TRIG:
		{
		  glBegin (GL_TRIANGLES);

		  Point3d lp1 = mesh->Point (el.PNum(1));
		  Point3d lp2 = mesh->Point (el.PNum(2));
		  Point3d lp3 = mesh->Point (el.PNum(3));
		  Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
		  n /= (n.Length() + 1e-12);
		  glNormal3dv (&n.X());
		  glVertex3dv (&lp1.X());
		  glVertex3dv (&lp2.X());
		  glVertex3dv (&lp3.X());
		  glEnd();
		  break;
		}
	      case QUAD:
		{
		  glBegin (GL_QUADS);

		  const Point3d & lp1 = mesh->Point (el.PNum(1));
		  const Point3d & lp2 = mesh->Point (el.PNum(2));
		  const Point3d & lp3 = mesh->Point (el.PNum(4));
		  const Point3d & lp4 = mesh->Point (el.PNum(3));
		  Vec3d n = Cross (Vec3d (lp1, lp2),
				   Vec3d (lp1, Center (lp3, lp4)));
		  n /= (n.Length() + 1e-12);
		  glNormal3d (n.X(), n.Y(), n.Z());
		  glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		  glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		  glVertex3d (lp4.X(), lp4.Y(), lp4.Z());
		  glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
		  glEnd();
		  break;
		}
	      case TRIG6:
		{
		  int lines[6][2] = {
		    { 1, 6 }, { 2, 6 },
		    { 1, 5 }, { 3, 5 },
		    { 2, 4 }, { 3, 4 } };

		  glBegin (GL_LINES);
		  for (j = 0; j < 6; j++)
		    {
		      glVertex3dv ( mesh->Point (el.PNum(lines[j][0])) );
		      glVertex3dv ( mesh->Point (el.PNum(lines[j][0])) );
		    }
		  glEnd();
		  break;
		}

	      case QUAD6:
		{
		  int lines[6][2] = {
		    { 1, 5 }, { 2, 5 },
		    { 3, 6 }, { 4, 6 },
		    { 1, 4 }, { 2, 3 } };

		  glBegin (GL_LINES);

		  for (j = 0; j < 6; j++)
		    {
		      const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0]));
		      const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1]));

		      glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		      glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		    }
		  glEnd ();
		  break;
		}
	      default:
		PrintSysError ("Cannot draw surface element of type ",
			       int(el.GetType()));
	      }
	  }
	glLoadName (0);

      }
    glEndList ();


    if (1)
      {

	identifiedlist = glGenLists (1);
	glNewList (identifiedlist, GL_COMPILE);

	GLfloat identifiedcol[] = { 1, 0, 1, 1 };

	glLineWidth (3);

	//  for (i = 1; i <= mesh->GetNSeg(); i++)

        if (& mesh -> GetIdentifications() )
          {
            INDEX_2_HASHTABLE<int> & idpts =
              mesh->GetIdentifications().GetIdentifiedPoints();
            if (&idpts)
              {
                for (i = 1; i <= idpts.GetNBags(); i++)
                  for (j = 1; j <= idpts.GetBagSize(i); j++)
                    {
                      INDEX_2 pts;
                      int val;

                      idpts.GetData (i, j, pts, val);
                      const Point3d & p1 = mesh->Point(pts.I1());
                      const Point3d & p2 = mesh->Point(pts.I2());

                      glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
                                    identifiedcol);

                      glBegin (GL_LINES);
                      glVertex3f (p1.X(), p1.Y(), p1.Z());
                      glVertex3f (p2.X(), p2.Y(), p2.Z());
                      glEnd();
                    }
              }
          }

	glEndList ();
      }

    if (lock)
      {
	lock -> UnLock();
	delete lock;
	lock = NULL;
      }

    vstimestamp = meshtimestamp;
  }




  void VisualSceneMesh :: BuildFilledList()
  {
    static int timer = NgProfiler::CreateTimer ("Mesh::BuildFilledList");
    NgProfiler::RegionTimer reg (timer);

    // cout << "buildilled" << endl;

#ifdef PARALLELGL
    cout << "buildfillelist, id = " << id << ", nse = " << mesh -> GetNSE() << endl;

    if (id == 0 && ntasks > 1)
      {
	InitParallelGL();

	par_filledlists.SetSize (ntasks);

	for ( int dest = 1; dest < ntasks; dest++ )
	  {
	    MyMPI_Send ("redraw", dest);
	    MyMPI_Send ("filledlist", dest);
	  }
	for ( int dest = 1; dest < ntasks; dest++ )
          {
            MyMPI_Recv (par_filledlists[dest], dest);
            cout << "proc " << dest << " has drawn to list " << par_filledlists[dest] << endl;
          }

	if (filledlist)
	  glDeleteLists (filledlist, 1);

	filledlist = glGenLists (1);
	glNewList (filledlist, GL_COMPILE);

	for ( int dest = 1; dest < ntasks; dest++ )
	  glCallList (par_filledlists[dest]);

	glEndList();



	filledtimestamp = NextTimeStamp();
	return;
      }

#endif


    // clock_t starttime, endtime;
    // starttime = clock();

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    filledtimestamp = NextTimeStamp();

    if (filledlist)
      glDeleteLists (filledlist, 1);

    filledlist = glGenLists (1);
    glNewList (filledlist, GL_COMPILE);

    // cout << "I am p " << id << " and got filledlist " << filledlist << endl;

    bool checkvicinity = (stlgeometry != NULL) && stldoctor.showvicinity;

    glEnable (GL_NORMALIZE);

    glLineWidth (1.0f);

    Vector locms;

    if (vispar.colormeshsize)
      {
	glEnable (GL_COLOR_MATERIAL);
	locms.SetSize (mesh->GetNP());
	double maxh = -1;
	double minh = 1e99;
	for (int i = 1; i <= locms.Size(); i++)
	  {
	    Point3d p = mesh->Point(i);
	    locms(i-1) = mesh->GetH (p);
	    if (locms(i-1) > maxh) maxh = locms(i-1);
	    if (locms(i-1) < minh) minh = locms(i-1);
	  }
	if (!locms.Size())
	  { minh = 1; maxh = 10; }
      }
    else
      glDisable (GL_COLOR_MATERIAL);


    GLfloat matcol[] = { 0, 1, 0, 1 };
    GLfloat matcolsel[] = { 1, 0, 0, 1 };
#ifdef PARALLEL
    GLfloat mat_coll_transp[] = { 0, 1, 0, 0.3 };
    GLfloat mat_coll_transp_sel[] = { 1, 0, 0, 0.3 };
#endif

    CurvedElements & curv = mesh->GetCurvedElements();
    int hoplotn = 1 << vispar.subdivisions;

    for (int col = 1; col <= 2; col++)
      {
	if (col == 2)
	  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolsel);
	else
	  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol);

	for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++)
	  {
	    const Element2d & el = (*mesh)[sei];

#ifdef PARALLEL
	    if ( el.IsGhost() )
	      {
		if ( col == 2 )
		  glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_coll_transp_sel);
		else
		  glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_coll_transp);
	      }
	    else
	      {
		if (col == 2)
		  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolsel);
		else
		  glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol);
	      }
#endif

       // Philippose - 06/07/2009
       // Modified the colour system to integrate the face colours into 
       // the mesh data structure, rather than limit it to the OCC geometry 
       // structure... allows other geometry types to use face colours too
       if(col == 1)
       {
          matcol[0] = mesh->GetFaceDescriptor(el.GetIndex()).SurfColour().X();
          matcol[1] = mesh->GetFaceDescriptor(el.GetIndex()).SurfColour().Y();
          matcol[2] = mesh->GetFaceDescriptor(el.GetIndex()).SurfColour().Z();

          matcol[3] = 1.0;
          glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol);
       }

	    bool drawel = !el.IsDeleted();

	    if (checkvicinity)
	      for (int j = 0; j < el.GetNP(); j++)
		if (!stlgeometry->Vicinity(el.GeomInfoPi(j+1).trignum))
		  drawel = 0;

	    if (!drawel)
	      continue;

	    if (vispar.colormeshsize && col == 2)
	      continue;
	    if (!vispar.colormeshsize &&
		(col == 2) != (el.GetIndex() == selface))
	      continue;

            glLoadName (sei+1);

	    switch (el.GetType())
	      {
	      case TRIG:
		{
		  if (curv.IsHighOrder()) //  && curv.IsSurfaceElementCurved(sei))
		    {
                      if (hoplotn > 128) hoplotn = 128;
                      Point<3> xa[129];
                      Vec<3> na[129];

		      for (int i = 0; i < hoplotn; i++)
                        {
                          glBegin (GL_TRIANGLE_STRIP);

                          for (int j = 0; j <= hoplotn-i; j++)
                            for (int k = 0; k < 2; k++)
                              {
                                if (j == hoplotn-i && k == 1) continue;

                                if (i > 0 && k == 0)
                                  {
                                    glNormal3dv (na[j]);
                                    glVertex3dv (xa[j]);
                                    continue;
                                  }

                                Point<2> xref (double(j) / hoplotn, double(i+k) / hoplotn);
                                Point<3> xglob;
                                Mat<3,2> dxdxi;
                                Vec<3> dx, dy, n;

                                curv.CalcSurfaceTransformation (xref, sei, xglob, dxdxi);
                                for (int i = 0; i < 3; i++)
                                  {
                                    dx(i) = dxdxi(i,0);
                                    dy(i) = dxdxi(i,1);
                                  }
                                n = Cross (dx, dy);
                                glNormal3dv (n);
                                glVertex3dv (xglob);

                                if (k == 1)
                                  {
                                    na[j] = n;
                                    xa[j] = xglob;
                                  }
                              }
                          glEnd();
                        }
		    }
		  else // not high order
		    {
                      glBegin (GL_TRIANGLES);

		      const Point<3> & lp0 = (*mesh) [el[0]];
		      const Point<3> & lp1 = (*mesh) [el[1]];
		      const Point<3> & lp2 = (*mesh) [el[2]];

		      Vec<3> n = Cross (lp1-lp0, lp2-lp0).Normalize();
		      glNormal3dv (n);

		      if (vispar.colormeshsize)
			{
			  SetOpenGlColor  (locms(el[0]-1), minh, maxh, 1);
			  glVertex3dv (lp0);
			  SetOpenGlColor  (locms(el[1]-1), minh, maxh, 1);
			  glVertex3dv (lp1);
			  SetOpenGlColor  (locms(el[2]-1), minh, maxh, 1);
			  glVertex3dv (lp2);
			}
		      else
			{
			  glVertex3dv (lp0);
			  glVertex3dv (lp1);
			  glVertex3dv (lp2);
			}

		      glEnd();
		    }

		  break;
		}
	      case QUAD:
		{
		  // cout << "BuildFilledList: QUAD" << endl;
		  // CurvedElements & curv = mesh->GetCurvedElements();
		  if (curv.IsHighOrder()) //  && curv.IsSurfaceElementCurved(sei))
		    {
		      Point<2> xr[4];
		      Point<3> xg;
		      Vec<3> dx, dy, n;

		      glBegin (GL_QUADS);

		      for (int i = 0; i < hoplotn; i++)
			for (int j = 0; j < hoplotn; j++)
			  {
			    xr[0](0) = (double)    i/hoplotn; xr[0](1) = (double)    j/hoplotn;
			    xr[1](0) = (double)(i+1)/hoplotn; xr[1](1) = (double)    j/hoplotn;
			    xr[2](0) = (double)(i+1)/hoplotn; xr[2](1) = (double)(j+1)/hoplotn;
			    xr[3](0) = (double)    i/hoplotn; xr[3](1) = (double)(j+1)/hoplotn;

			    for (int l=0; l<4; l++)
			      {
				Mat<3,2> dxdxi;

				curv.CalcSurfaceTransformation (xr[l], sei, xg, dxdxi);
				for (int i = 0; i < 3; i++)
				  {
				    dx(i) = dxdxi(i,0);
				    dy(i) = dxdxi(i,1);
				  }

				n = Cross (dx, dy);
				n.Normalize();
				glNormal3d (n(0), n(1), n(2));
				glVertex3d (xg(0), xg(1), xg(2));
			      }

			  }

		      glEnd();
		    }

		  else // not high order

		    {
		      glBegin (GL_QUADS);

		      const Point<3> & lp1 = mesh->Point (el.PNum(1));
		      const Point<3> & lp2 = mesh->Point (el.PNum(2));
		      const Point<3> & lp3 = mesh->Point (el.PNum(4));
		      const Point<3> & lp4 = mesh->Point (el.PNum(3));

		      Vec<3> n = Cross (lp2-lp1,  Center (lp3, lp4)-lp1);
		      n.Normalize();
		      glNormal3dv (n);

		      glVertex3dv (lp1);
		      glVertex3dv (lp2);
		      glVertex3dv (lp4);
		      glVertex3dv (lp3);

		      glEnd ();
		    }
		  break;
		}

	      case TRIG6:
		{
		  glBegin (GL_TRIANGLES);

		  static int trigs[4][3] = {
		    { 1, 6, 5 },
		    { 2, 4, 6 },
		    { 3, 5, 4 },
		    { 4, 5, 6 } };

		  for (int j = 0; j < 4; j++)
		    {
                      const Point<3> & lp1 = mesh->Point (el.PNum(trigs[j][0]));
                      const Point<3> & lp2 = mesh->Point (el.PNum(trigs[j][1]));
                      const Point<3> & lp3 = mesh->Point (el.PNum(trigs[j][2]));
		      // Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
                      Vec<3> n = Cross (lp2-lp1, lp3-lp1);
		      glNormal3dv (n);

		      glVertex3dv (lp1);
		      glVertex3dv (lp2);
		      glVertex3dv (lp3);
		    }
		  glEnd();
		  break;
		}

	      case QUAD6:
		{
		  glBegin (GL_QUADS);
		  static int quads[2][4] = {
		    { 1, 5, 6, 4 },
		    { 5, 2, 3, 6 } };

		  for (int j = 0; j < 2; j++)
		    {
                      Point3d lp1 = mesh->Point (el.PNum(quads[j][0]));
                      Point3d lp2 = mesh->Point (el.PNum(quads[j][1]));
                      Point3d lp3 = mesh->Point (el.PNum(quads[j][2]));
		      Point3d lp4 = mesh->Point (el.PNum(quads[j][3]));
		      Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
		      n /= (n.Length() + 1e-12);
		      glNormal3dv (&n.X());
		      glVertex3dv (&lp1.X());
		      glVertex3dv (&lp2.X());
		      glVertex3dv (&lp3.X());
		      glVertex3dv (&lp4.X());
		    }
		  glEnd();
		  break;
		}

	      case QUAD8:
		{
		  glBegin (GL_TRIANGLES);
		  static int boundary[] =
		    { 1, 5, 2, 8, 3, 6, 4, 7, 1 };

		  Point3d c(0,0,0);
		  for (int j = 0; j < 4; j++)
		    {
		      const Point3d & hp = mesh->Point (el[j]);
		      c.X() -= 0.25 * hp.X();
		      c.Y() -= 0.25 * hp.Y();
		      c.Z() -= 0.25 * hp.Z();
		    }
		  for (int j = 4; j < 8; j++)
		    {
		      const Point3d & hp = mesh->Point (el[j]);
		      c.X() += 0.5 * hp.X();
		      c.Y() += 0.5 * hp.Y();
		      c.Z() += 0.5 * hp.Z();
		    }

		  for (int j = 0; j < 8; j++)
		    {
		      Point3d lp1 = mesh->Point (el.PNum(boundary[j]));
		      Point3d lp2 = mesh->Point (el.PNum(boundary[j+1]));

		      Vec3d n = Cross (Vec3d (c, lp1), Vec3d (c, lp2));
		      n /= (n.Length() + 1e-12);
		      glNormal3dv (&n.X());
		      glVertex3dv (&lp1.X());
		      glVertex3dv (&lp2.X());
		      glVertex3dv (&c.X());
		    }
		  glEnd();
		  break;
		}


	      default:
		PrintSysError ("Cannot draw (2) surface element of type ",
			       int(el.GetType()));
	      }



	  }
      }
    glLoadName (0);
    glEndList ();


#ifdef PARALLELGL
    glFinish();
    if (id > 0)
      MyMPI_Send (filledlist, 0);
#endif

    // endtime = clock();
    // cout << "BuildFillList time = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl;
  }


  void VisualSceneMesh :: BuildLineList()
  {
    static int timer = NgProfiler::CreateTimer ("Mesh::BuildLineList");
    NgProfiler::RegionTimer reg (timer);

#ifdef PARALLELGL

    if (id == 0 && ntasks > 1)
      {
	InitParallelGL();

	par_linelists.SetSize (ntasks);

	for ( int dest = 1; dest < ntasks; dest++ )
	  {
	    MyMPI_Send ("redraw", dest);
	    MyMPI_Send ("linelist", dest);
	  }
	for ( int dest = 1; dest < ntasks; dest++ )
	  MyMPI_Recv (par_linelists[dest], dest);

	if (linelist)
	  glDeleteLists (linelist, 1);

	linelist = glGenLists (1);
	glNewList (linelist, GL_COMPILE);

	for ( int dest = 1; dest < ntasks; dest++ )
	  glCallList (par_linelists[dest]);

	glEndList();


	linetimestamp = NextTimeStamp();
	return;
      }

#endif

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    linetimestamp = NextTimeStamp();

    bool checkvicinity = (stlgeometry != NULL) && stldoctor.showvicinity;

    if (linelist)
      glDeleteLists (linelist, 1);

    linelist = glGenLists (1);
    glNewList (linelist, GL_COMPILE);

    // cout << "linelist = " << linelist << endl;

    glLineWidth (1.0f);


    int hoplotn = 1 << vispar.subdivisions;

    // PrintMessage (3, "nse = ", mesh->GetNSE());
    for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++)
      {
	const Element2d & el = (*mesh)[sei];

	bool drawel = !el.IsDeleted();
	if (checkvicinity)
	  for (int j = 0; j < el.GetNP(); j++)
	    if (!stlgeometry->Vicinity(el.GeomInfoPi(j+1).trignum))
	      drawel = 0;

	if (!drawel)
	  continue;

	switch (el.GetType())
	  {
	  case TRIG:
	    {
	      CurvedElements & curv = mesh->GetCurvedElements();
	      if (curv.IsHighOrder()) //  && curv.IsSurfaceElementCurved(sei))
		{
		  Point<3> xg;
		  glBegin (GL_LINE_LOOP);
		  for (int i = 0; i < hoplotn; i++)
		    {
		      Point<2> xr (double(i) / hoplotn, 0);
		      curv.CalcSurfaceTransformation (xr, sei, xg);
		      glVertex3dv (xg);
		    }
		  for (int i = 0; i < hoplotn; i++)
		    {
		      Point<2> xr (double(hoplotn-i) / hoplotn, double(i)/hoplotn);
		      curv.CalcSurfaceTransformation (xr, sei, xg);
		      glVertex3dv (xg);
		    }
		  for (int i = 0; i < hoplotn; i++)
		    {
		      Point<2> xr (0, double(hoplotn-i) / hoplotn);
		      curv.CalcSurfaceTransformation (xr, sei, xg);
		      glVertex3dv (xg);
		    }

		  glEnd();
		}
	      else
		{
		  glBegin (GL_TRIANGLES);

		  const Point<3> & lp0 = (*mesh) [el[0]];
		  const Point<3> & lp1 = (*mesh) [el[1]];
		  const Point<3> & lp2 = (*mesh) [el[2]];

		  glVertex3dv (lp0);
		  glVertex3dv (lp1);
		  glVertex3dv (lp2);

		  glEnd();
		}

	      break;

	    }

	  case QUAD:
	    {
	      CurvedElements & curv = mesh->GetCurvedElements();
	      if (curv.IsHighOrder()) //  && curv.IsSurfaceElementCurved(sei))
		{
		  Point<2> xr;
		  Point<3> xg;

		  glBegin (GL_LINE_STRIP);

		  for (int side = 0; side < 4; side++)
		    {
		      for (int i = 0; i <= hoplotn; i++)
			{
			  switch (side)
			    {
			    case 0:
			      xr(0) = (double) i/hoplotn;
			      xr(1) = 0.;
			      break;
			    case 1:
			      xr(0) = 1.;
			      xr(1) = (double) i/hoplotn;
			      break;
			    case 2:
			      xr(0) = (double) (hoplotn-i)/hoplotn;
			      xr(1) = 1.;
			      break;
			    case 3:
			      xr(0) = 0.;
			      xr(1) = (double) (hoplotn-i)/hoplotn;
			      break;
			    }

			  curv.CalcSurfaceTransformation (xr, sei, xg);
			  glVertex3d (xg(0), xg(1), xg(2));

			}

		    }
		  glEnd();

		} else {

		  glBegin (GL_QUADS);

		  const Point3d & lp1 = mesh->Point (el.PNum(1));
		  const Point3d & lp2 = mesh->Point (el.PNum(2));
		  const Point3d & lp3 = mesh->Point (el.PNum(4));
		  const Point3d & lp4 = mesh->Point (el.PNum(3));
		  Vec3d n = Cross (Vec3d (lp1, lp2),
				   Vec3d (lp1, Center (lp3, lp4)));
		  glNormal3d (n.X(), n.Y(), n.Z());
		  glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		  glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		  glVertex3d (lp4.X(), lp4.Y(), lp4.Z());
		  glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
		  glEnd();

		}

	      break;

	    }

	  case TRIG6:
	    {
	      int lines[6][2] = {
		{ 1, 6 }, { 2, 6 },
		{ 1, 5 }, { 3, 5 },
		{ 2, 4 }, { 3, 4 } };

	      glBegin (GL_LINES);
	      for (int j = 0; j < 6; j++)
		{
		  const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0]));
		  const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1]));

		  glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		  glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		}

	      glEnd();
	      break;
	    }

	  case QUAD6:
	    {
	      int lines[6][2] = {
		{ 1, 5 }, { 2, 5 },
		{ 3, 6 }, { 4, 6 },
		{ 1, 4 }, { 2, 3 } };

	      glBegin (GL_LINES);

	      for (int j = 0; j < 6; j++)
		{
		  const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0]));
		  const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1]));

		  glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		  glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		}
	      glEnd ();
	      break;
	    }

	  case QUAD8:
	    {
	      int lines[8][2] = {
		{ 1, 5 }, { 2, 5 }, { 3, 6 }, { 4, 6 },
		{ 1, 7 }, { 4, 7 }, { 2, 8 }, { 3, 8 }
	      };

	      glBegin (GL_LINES);

	      for (int j = 0; j < 8; j++)
		{
		  const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0]));
		  const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1]));

		  glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		  glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		}
	      glEnd ();
	      break;
	    }



	  default:
	    PrintSysError ("Cannot draw (4) surface element of type ",
			   int(el.GetType()));
	  }
      }

    glEndList ();


#ifdef PARALLELGL
    glFinish();
    if (id > 0)
      MyMPI_Send (linelist, 0);
#endif
  }



  void VisualSceneMesh :: BuildEdgeList()
  {
    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    if (edgetimestamp > mesh->GetTimeStamp () && vispar.shrink == 1)
      return;

    edgetimestamp = NextTimeStamp();

    if (edgelist)
      glDeleteLists (edgelist, 1);

    edgelist = glGenLists (1);
    glNewList (edgelist, GL_COMPILE);


    GLfloat matcoledge[] = { 0, 0, 1, 1 };
    GLfloat matcolsingedge[] = { 1, 0, 1, 1 };

    glEnable (GL_POLYGON_OFFSET_LINE);
    glPolygonOffset (1, -1);

    glEnable (GL_COLOR_MATERIAL);
    glDisable (GL_LIGHTING);

    for (int i = 1; i <= mesh->GetNSeg(); i++)
      {
	const Segment & seg = mesh->LineSegment(i);
	const Point3d & p1 = (*mesh)[seg[0]];
	const Point3d & p2 = (*mesh)[seg[1]];

	if (seg.singedge_left || seg.singedge_right)
	  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
			matcolsingedge);
	else
	  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
			matcoledge);

	if (seg.singedge_left || seg.singedge_right)
	  glColor3fv (matcolsingedge);
	else
	  glColor3fv (matcoledge);

	if (seg.edgenr == seledge)
	  glLineWidth(5);
	else
	  glLineWidth(2);

	if (mesh->GetCurvedElements().IsHighOrder())
	  {
	    int hoplotn = 1 << vispar.subdivisions;
	    // mesh->GetCurvedElements().GetNVisualSubsecs();

	    Point<3> x;
	    glBegin (GL_LINE_STRIP);

	    for (int j = 0; j <= hoplotn; j++)
	      {
		mesh->GetCurvedElements().CalcSegmentTransformation ((double) j/hoplotn, i-1, x);
		glVertex3d (x(0), x(1), x(2));
                /*
                cout << "x = " << x(0) << ", " << x(1) << ", " << x(2)
                     << ", norm = 1+" << sqrt(x(0)*x(0)+x(1)*x(1))-1
                     << ", phi = " << atan2(x(1), x(0))/M_PI << endl;
                */
	      }

	    glEnd();

	  }
	else
	  {
	    glBegin (GL_LINES);
	    Point<3> hp1 = p1;
	    Point<3> hp2 = p2;
	    Point<3> c = Center(p1, p2);
	    if (vispar.shrink < 1)
	      {
		hp1 = c + vispar.shrink * (hp1 - c);
		hp2 = c + vispar.shrink * (hp2 - c);
	      }
	    glVertex3dv (hp1);
	    glVertex3dv (hp2); // p2.X(), p2.Y(), p2.Z());
	    glEnd();
	  }
      }

    glLineWidth (2);
    glDisable (GL_POLYGON_OFFSET_LINE);

    glDisable (GL_COLOR_MATERIAL);
    glEnable (GL_LIGHTING);

    glEndList();
  }




  void VisualSceneMesh :: BuildPointNumberList()
  {
    ;
  }



  // Bernstein Pol B_{n,i}(x) = n! / i! / (n-i)! (1-x)^{n-i} x^i
  static inline double Bernstein (int n, int i, double x)
  {
    double val = 1;
    for (int j = 1; j <= i; j++)
      val *= x;
    for (int j = 1; j <= n-i; j++)
      val *= (1-x) * (j+i) / j;
    return val;
  }

  void ToBernstein (int order, Point<3> * pts, int stride)
  {
    static DenseMatrix mat, inv;
    static Vector vec1, vec2;

    if (mat.Height () != order+1)
      {
	mat.SetSize (order+1);
	inv.SetSize (order+1);
	vec1.SetSize (order+1);
	vec2.SetSize (order+1);
	for (int i = 0; i <= order; i++)
	  {
	    double x = double(i) / order;
	    for (int j = 0; j <= order; j++)
              mat(i,j) = Bernstein (order, j, x);
	  }

	CalcInverse (mat, inv);
      }

    for (int i = 0; i < 3; i++)
      {
	for (int j = 0; j <= order; j++)
	  vec1(j) = pts[j*stride](i);

	inv.Mult (vec1, vec2);

	for (int j = 0; j <= order; j++)
	  pts[j*stride](i) = vec2(j);
      }
  }














  void VisualSceneMesh :: BuildTetList()
  {

    if (tettimestamp > mesh->GetTimeStamp () &&
	tettimestamp > vispar.clipplanetimestamp )
      return;

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    tettimestamp = NextTimeStamp();

    if (tetlist)
      glDeleteLists (tetlist, 1);


    tetlist = glGenLists (1);
    glNewList (tetlist, GL_COMPILE);



    Array<Element2d> faces;

    BitArray shownode(mesh->GetNP());
    if (vispar.clipenable)
      {
	shownode.Clear();
	for (int i = 1; i <= shownode.Size(); i++)
	  {
	    Point<3> p = mesh->Point(i);

	    double val =
	      p[0] * clipplane[0] +
	      p[1] * clipplane[1] +
	      p[2] * clipplane[2] +
	      clipplane[3];

	    if (val > 0) shownode.Set (i);
	  }
      }
    else
      shownode.Set();


#ifdef PARALLEL

    static float tetcols[][8] =
      {
	{ 1.0f, 1.0f, 0.0f, 1.0f },
	{ 1.0f, 0.0f, 0.0f, 1.0f },
	{ 0.0f, 1.0f, 0.0f, 1.0f },
	{ 0.0f, 0.0f, 1.0f, 1.0f },
	{ 1.0f, 1.0f, 0.0f, 0.3f },
	{ 1.0f, 0.0f, 0.0f, 0.3f },
	{ 0.0f, 1.0f, 0.0f, 0.3f },
	{ 0.0f, 0.0f, 1.0f, 0.3f }
      };

#else
    static float tetcols[][4] =
      {
	{ 1.0f, 1.0f, 0.0f, 1.0f },
	{ 1.0f, 0.0f, 0.0f, 1.0f },
	{ 0.0f, 1.0f, 0.0f, 1.0f },
	{ 0.0f, 0.0f, 1.0f, 1.0f }
      };
#endif


    CurvedElements & curv = mesh->GetCurvedElements();


    if (!curv.IsHighOrder())
      glShadeModel (GL_FLAT);
    else
      glShadeModel (GL_SMOOTH);

    int hoplotn = max (2, 1 << vispar.subdivisions);

    for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
      {
	if (vispar.drawtetsdomain > 0 &&
	    vispar.drawtetsdomain != (*mesh)[ei].GetIndex())
	  continue;

	const Element & el = (*mesh)[ei];

	if ((el.GetType() == TET || el.GetType() == TET10) && !el.IsDeleted())
	  {

	    bool drawtet = 1;
	    for (int j = 0; j < 4; j++)
	      if (!shownode.Test(el[j]))
		drawtet = 0;
	    if (!drawtet) continue;


	    int ind = el.GetIndex() % 4;


#ifdef PARALLEL
	    if (vispar.drawmetispartition && (el.GetPartition()!=-1))
	      ind = el.GetPartition() % 4;
	    // (*testout) << "ind (" << i << ") = " << ind << endl;


            if ( el.IsGhost() )
              {
                glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, tetcols[ind+4]);
              }
            else
              {
                glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, tetcols[ind]);
              }
#else
	    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, tetcols[ind]);
#endif




	    if (curv.IsHighOrder()) //  && curv.IsElementCurved(ei))
	      {
		const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (TET);
		const Point3d * vertices = MeshTopology :: GetVertices (TET);

		/*
		Point<3> grid[11][11];
		Point<3> fpts[3];
		int order = vispar.subdivisions+1;

		for (int trig = 0; trig < 4; trig++)
		  {
		    for (int j = 0; j < 3; j++)
		      fpts[j] = vertices[faces[trig][j]-1];

		    static Point<3> c(0.25, 0.25, 0.25);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 3; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[3] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      double(iy)/order };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l);

			  curv.CalcElementTransformation (xl, i-1, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 0.999, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }
		*/



		int order = curv.GetOrder();

		Array<Point<3> > ploc ( (order+1)*(order+1) );
		Array<Point<3> > pglob ( (order+1)*(order+1) );
		Point<3> fpts[3];

		for (int trig = 0; trig < 4; trig++)
		  {
		    for (int j = 0; j < 3; j++)
		      fpts[j] = vertices[faces[trig][j]-1];

		    static Point<3> c(0.25, 0.25, 0.25);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 3; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0, ii = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++, ii++)
			{
			  double lami[3] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      double(iy)/order };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l);

			  ploc[ii] = xl;
			}

		    curv.CalcMultiPointElementTransformation (&ploc, ei, &pglob, 0);

		    Point<3> grid[11][11];
		    for (int ix = 0, ii = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++, ii++)
			grid[ix][iy] = pglob[ii];

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(hoplotn, 0.0, 0.9999f, hoplotn, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, hoplotn, 0, hoplotn);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }
	      }

	    else

	      {
                Point<3> pts[4];
                for (int j = 0; j < 4; j++)
                  pts[j] = (*mesh)[el[j]];

                if (vispar.shrink < 1)
                  {
                    Point<3> c = Center (pts[0], pts[1], pts[2], pts[3]);
                    for (int j = 0; j < 4; j++)
                      pts[j] = c + vispar.shrink * (pts[j]-c);
                  }


                Vec<3> n;

                glBegin (GL_TRIANGLE_STRIP);

                n = Cross (pts[1]-pts[0], pts[2]-pts[0]);
                glNormal3dv (n);
                glVertex3dv (pts[0]);
                glVertex3dv (pts[1]);
                glVertex3dv (pts[2]);

                n = Cross (pts[3]-pts[1], pts[2]-pts[1]);
                glNormal3dv (n);
                glVertex3dv (pts[3]);

                n = Cross (pts[3]-pts[2], pts[0]-pts[2]);
                glNormal3dv (n);
                glVertex3dv (pts[0]);

                n = Cross (pts[1]-pts[3], pts[0]-pts[3]);
                glNormal3dv (n);
                glVertex3dv (pts[1]);

                glEnd();
	      }
	  }
      }

    glEndList ();
  }




  void VisualSceneMesh :: BuildPrismList()
  {
    if (prismtimestamp > mesh->GetTimeStamp () &&
	prismtimestamp > vispar.clipplanetimestamp )
      return;

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    prismtimestamp = NextTimeStamp();



    if (prismlist)
      glDeleteLists (prismlist, 1);

    prismlist = glGenLists (1);
    glNewList (prismlist, GL_COMPILE);

    static float prismcol[] = { 0.0f, 1.0f, 1.0f, 1.0f };
    glLineWidth (1.0f);

    Array<Element2d> faces;


    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, prismcol);

    for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
      {
	const Element & el = (*mesh)[ei];
	if (el.GetType() == PRISM && !el.IsDeleted())
	  {
	    int j;
	    int i = ei + 1;

	    CurvedElements & curv = mesh->GetCurvedElements();
	    if (curv.IsHighOrder()) //  && curv.IsElementCurved(ei))
	      {
		const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PRISM);
		const Point3d * vertices = MeshTopology :: GetVertices (PRISM);

		Point<3> grid[11][11];
		Point<3> fpts[4];
		int order = vispar.subdivisions+1;

		for (int trig = 0; trig < 2; trig++)
		  {
		    for (int j = 0; j < 3; j++)
		      fpts[j] = vertices[faces[trig][j]-1];

		    static Point<3> c(1.0/3.0, 1.0/3.0, 0.5);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 3; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[3] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      double(iy)/order };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l);

			  curv.CalcElementTransformation (xl, i-1, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 0.999f, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }

		for (int quad = 2; quad < 5; quad++)
		  {
		    for (int j = 0; j < 4; j++)
		      fpts[j] = vertices[faces[quad][j]-1];

		    static Point<3> c(1.0/3.0, 1.0/3.0, 0.5);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 4; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[4] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (  double(iy)/order),
			      (1-double(ix)/order) * (  double(iy)/order) };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) =
			      lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l) + lami[3] * fpts[3](l);

			  curv.CalcElementTransformation (xl, ei, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }





		/*
		int hoplotn = 1 << vispar.subdivisions;
		// int hoplotn = curv.GetNVisualSubsecs();

		const Point3d * facepoint = MeshTopology :: GetVertices (TRIG);
		const ELEMENT_FACE * elface = MeshTopology :: GetFaces(TRIG);

		glBegin (GL_TRIANGLES);

		for (int trig = 0; trig<2; trig++)
		  {

		    Vec<3> x0,x1,d0,d1;
		    x0 = facepoint[1] - facepoint[2];
		    x1 = facepoint[0] - facepoint[2];
		    x0.Normalize();
		    x1.Normalize();
		    if (trig == 1) swap (x0,x1);

		    Point<3> xr[3];
		    Point<3> xg;
		    Vec<3> dx, dy, dz, n;

		    for (int i1 = 0; i1 < hoplotn; i1++)
		      for (int j1 = 0; j1 < hoplotn-i1; j1++)
			for (int k = 0; k < 2; k++)
			  {
			    if (k == 0)
			      {
				xr[0](0) = (double)    i1/hoplotn; xr[0](1) = (double)    j1/hoplotn;
				xr[1](0) = (double)(i1+1)/hoplotn; xr[1](1) = (double)    j1/hoplotn;
				xr[2](0) = (double)    i1/hoplotn; xr[2](1) = (double)(j1+1)/hoplotn;
			      } else
				{
				  if (j1 == hoplotn-i1-1) continue;
				  xr[0](0) = (double)(i1+1)/hoplotn; xr[0](1) = (double)    j1/hoplotn;
				  xr[1](0) = (double)(i1+1)/hoplotn; xr[1](1) = (double)(j1+1)/hoplotn;
				  xr[2](0) = (double)    i1/hoplotn; xr[2](1) = (double)(j1+1)/hoplotn;
				};

			    for (int l=0; l<3; l++)
			      {
				Mat<3,3> dxdxi;
				xr[l](2) = (double) trig;
				curv.CalcElementTransformation (xr[l], i-1, xg, dxdxi);
				for (int i = 0; i < 3; i++)
				  {
				    dx(i) = dxdxi(i,0);
				    dy(i) = dxdxi(i,1);
				    dz(i) = dxdxi(i,2);
				  }

				Vec<3> d0 = x0(0)*dx + x0(1)*dy + x0(2)*dz;
				Vec<3> d1 = x1(0)*dx + x1(1)*dy + x1(2)*dz;
				n = Cross (d1, d0);
				glNormal3d (n(0), n(1), n(2));
				glVertex3d (xg(0), xg(1), xg(2));
			      }
			  }

		  }

		glEnd ();

		glBegin (GL_QUADS);

		for (int quad = 0; quad<3; quad++)
		  {
		    const Point3d * facepoint = MeshTopology :: GetVertices (PRISM);

		    Vec<3> x0,x1;
		    int xyz;

		    switch (quad)
		      {
		      case 0:
			x0 = facepoint[5] - facepoint[2];
			x1 = facepoint[0] - facepoint[2];
			xyz = 0;
			break;
		      case 1:
			x0 = facepoint[4] - facepoint[0];
			x1 = facepoint[1] - facepoint[0];
			xyz = 0;
			break;
		      case 2:
			x0 = facepoint[1] - facepoint[2];
			x1 = facepoint[5] - facepoint[2];
			xyz = 1;
			break;
		      }

		    x0.Normalize();
		    x1.Normalize();

		    swap (x0,x1);

		    Point<3> xr[4];
		    Point<3> xg;
		    Vec<3> dx, dy, dz, n;

		    for (int i1 = 0; i1 < hoplotn; i1++)
		      for (int j1 = 0; j1 < hoplotn; j1++)
			{
			  xr[0](xyz) = (double)    i1/hoplotn; xr[0](2) = (double)    j1/hoplotn;
			  xr[1](xyz) = (double)(i1+1)/hoplotn; xr[1](2) = (double)    j1/hoplotn;
			  xr[2](xyz) = (double)(i1+1)/hoplotn; xr[2](2) = (double)(j1+1)/hoplotn;
			  xr[3](xyz) = (double)    i1/hoplotn; xr[3](2) = (double)(j1+1)/hoplotn;

			  for (int l=0; l<4; l++)
			    {
			      switch (quad)
				{
				case 0: xr[l](1) = 0; break;
				case 1: xr[l](1) = 1-xr[l](0); break;
				case 2: xr[l](0) = 0; break;
				}

			      Mat<3,3> dxdxi;
			      curv.CalcElementTransformation (xr[l], i-1, xg, dxdxi);
			      for (int i = 0; i < 3; i++)
				{
				  dx(i) = dxdxi(i,0);
				  dy(i) = dxdxi(i,1);
				  dz(i) = dxdxi(i,2);
				}

			      Vec<3> d0 = x0(0)*dx + x0(1)*dy + x0(2)*dz;
			      Vec<3> d1 = x1(0)*dx + x1(1)*dy + x1(2)*dz;
			      n = Cross (d1, d0);
			      glNormal3d (n(0), n(1), n(2));
			      glVertex3d (xg(0), xg(1), xg(2));
			    }
			}
		  }
		glEnd ();
		*/
	      }
	    else
	      {
		Point3d c(0,0,0);
		if (vispar.shrink < 1)
		  {
		    for (j = 1; j <= 6; j++)
		      {
			Point3d p = mesh->Point(el.PNum(j));
			c.X() += p.X() / 6;
			c.Y() += p.Y() / 6;
			c.Z() += p.Z() / 6;
		      }
		  }

		el.GetSurfaceTriangles (faces);
		glBegin (GL_TRIANGLES);
		for (j = 1; j <= faces.Size(); j++)
		  {
		    Element2d & face = faces.Elem(j);
		    Point3d lp1 = mesh->Point (el.PNum(face.PNum(1)));
		    Point3d lp2 = mesh->Point (el.PNum(face.PNum(2)));
		    Point3d lp3 = mesh->Point (el.PNum(face.PNum(3)));
		    Vec3d n = Cross (Vec3d (lp1, lp3), Vec3d (lp1, lp2));
		    n /= (n.Length()+1e-12);
		    glNormal3d (n.X(), n.Y(), n.Z());
		    if (vispar.shrink < 1)
		      {
			lp1 = c + vispar.shrink * (lp1 - c);
			lp2 = c + vispar.shrink * (lp2 - c);
			lp3 = c + vispar.shrink * (lp3 - c);
		      }
		    glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		    glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		    glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
		  }

		glEnd();
	      }
	  }
      }
    glEndList ();
  }




  void VisualSceneMesh :: BuildHexList()
  {
    if (hextimestamp > mesh->GetTimeStamp () &&
	hextimestamp > vispar.clipplanetimestamp )
      return;

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    hextimestamp = NextTimeStamp();

    if (hexlist) glDeleteLists (hexlist, 1);

    hexlist = glGenLists (1);
    glNewList (hexlist, GL_COMPILE);


    static float hexcol[] = { 1.0f, 1.0f, 0.0f, 1.0f };
    glLineWidth (1.0f);
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, hexcol);

    Array<Element2d> faces;
    // int hoplotn = 1 << vispar.subdivisions;

    for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
      {
	const Element & el = (*mesh)[ei];
	if (el.GetType() == HEX && !el.IsDeleted())
	  {
	    CurvedElements & curv = mesh->GetCurvedElements();
	    if (curv.IsHighOrder()) //  && curv.IsElementCurved(ei))
	      {
		/* // classical
		   glBegin (GL_QUADS);

		   const ELEMENT_FACE * faces = MeshTopology :: GetFaces (HEX);
		   const Point3d * vertices = MeshTopology :: GetVertices (HEX);

		   Point<3> grid[33][33];
		   Vec<3> gridn[33][33];
		   Point<3> fpts[4];
		   for (int quad = 0; quad<6; quad++)
		   {
		   for (int j = 0; j < 4; j++)
		   fpts[j] = vertices[faces[quad][j]-1];

		   static Point<3> c(0.5, 0.5, 0.5);
		   if (vispar.shrink < 1)
		   for (int j = 0; j < 4; j++)
		   fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		   Vec<3> taux = fpts[1]-fpts[0];
		   Vec<3> tauy = fpts[3]-fpts[0];

		   for (int ix = 0; ix <= hoplotn; ix++)
		   for (int iy = 0; iy <= hoplotn; iy++)
		   {
		   Point<3> xl;
		   Mat<3,3> dxdxi;
		   double lami[4] =
		   { (1-double(ix)/hoplotn) * (1-double(iy)/hoplotn),
		   (  double(ix)/hoplotn) * (1-double(iy)/hoplotn),
		   (  double(ix)/hoplotn) * (  double(iy)/hoplotn),
		   (1-double(ix)/hoplotn) * (  double(iy)/hoplotn) };
		   for (int l = 0; l < 3; l++)
		   xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
		   lami[2] * fpts[2](l) + lami[3] * fpts[3](l);

		   curv.CalcElementTransformation (xl, ei, grid[ix][iy], dxdxi);

		   Vec<3> gtaux = dxdxi * taux;
		   Vec<3> gtauy = dxdxi * tauy;
		   gridn[ix][iy] = Cross (gtauy, gtaux).Normalize();
		   }

		   for (int ix = 0; ix < hoplotn; ix++)
		   for (int iy = 0; iy < hoplotn; iy++)
		   {
		   glNormal3dv (gridn[ix][iy]);
		   glVertex3dv (grid[ix][iy]);

		   glNormal3dv (gridn[ix+1][iy]);
		   glVertex3dv (grid[ix+1][iy]);

		   glNormal3dv (gridn[ix+1][iy+1]);
		   glVertex3dv (grid[ix+1][iy+1]);

		   glNormal3dv (gridn[ix][iy+1]);
		   glVertex3dv (grid[ix][iy+1]);
		   }
		   }

		   glEnd ();
		*/

		const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (HEX);
		const Point3d * vertices = MeshTopology :: GetVertices (HEX);

		Point<3> grid[11][11];
		Point<3> fpts[4];
		int order = vispar.subdivisions+1;

		for (int quad = 0; quad<6; quad++)
		  {
		    for (int j = 0; j < 4; j++)
		      fpts[j] = vertices[faces[quad][j]-1];

		    static Point<3> c(0.5, 0.5, 0.5);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 4; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[4] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (  double(iy)/order),
			      (1-double(ix)/order) * (  double(iy)/order) };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l) + lami[3] * fpts[3](l);

			  curv.CalcElementTransformation (xl, ei, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }
	      }
	    else
	      {
		Point3d c(0,0,0);
		if (vispar.shrink < 1)
		  {
		    for (int j = 1; j <= 8; j++)
		      {
			Point3d p = mesh->Point(el.PNum(j));
			c.X() += p.X();
			c.Y() += p.Y();
			c.Z() += p.Z();
		      }
		    c.X() /= 8;
		    c.Y() /= 8;
		    c.Z() /= 8;
		  }

		glBegin (GL_TRIANGLES);

		el.GetSurfaceTriangles (faces);
		for (int j = 1; j <= faces.Size(); j++)
		  {
		    Element2d & face = faces.Elem(j);
		    Point<3> lp1 = mesh->Point (el.PNum(face.PNum(1)));
		    Point<3> lp2 = mesh->Point (el.PNum(face.PNum(2)));
		    Point<3> lp3 = mesh->Point (el.PNum(face.PNum(3)));
		    Vec<3> n = Cross (lp3-lp1, lp2-lp1);
		    n.Normalize();
		    glNormal3dv (n);

		    if (vispar.shrink < 1)
		      {
			lp1 = c + vispar.shrink * (lp1 - c);
			lp2 = c + vispar.shrink * (lp2 - c);
			lp3 = c + vispar.shrink * (lp3 - c);
		      }

		    glVertex3dv (lp1);
		    glVertex3dv (lp2);
		    glVertex3dv (lp3);
		  }

		glEnd();
	      }
	  }
      }
    glEndList ();
  }









  void VisualSceneMesh :: BuildPyramidList()
  {
    if (pyramidtimestamp > mesh->GetTimeStamp () &&
	pyramidtimestamp > vispar.clipplanetimestamp )
      return;

    if (!lock)
      {
	lock = new NgLock (mesh->Mutex());
	lock -> Lock();
      }

    pyramidtimestamp = NextTimeStamp();


    if (pyramidlist)
      glDeleteLists (pyramidlist, 1);


    pyramidlist = glGenLists (1);
    glNewList (pyramidlist, GL_COMPILE);

    static float pyramidcol[] = { 1.0f, 0.0f, 1.0f, 1.0f };
    glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, pyramidcol);

    glLineWidth (1.0f);
    Array<Element2d> faces;

    for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
      {
	const Element & el = (*mesh)[ei];
	if (el.GetType() == PYRAMID && !el.IsDeleted())
	  {
	    int i = ei + 1;
            
	    CurvedElements & curv = mesh->GetCurvedElements();
	    if (curv.IsHighOrder()) //  && curv.IsElementCurved(ei))
	      {

		const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PYRAMID);
		const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID);

		Point<3> grid[11][11];
		Point<3> fpts[4];
		int order = vispar.subdivisions+1;

		for (int trig = 0; trig < 4; trig++)
		  {
		    for (int j = 0; j < 3; j++)
		      fpts[j] = vertices[faces[trig][j]-1];

		    static Point<3> c(0.375, 0.375, 0.25);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 3; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[3] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      double(iy)/order };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l);

			  curv.CalcElementTransformation (xl, i-1, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 0.999f, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }

		for (int quad = 4; quad < 5; quad++)
		  {
		    for (int j = 0; j < 4; j++)
		      fpts[j] = vertices[faces[quad][j]-1];

		    static Point<3> c(0.375, 0.375, 0.25);
		    if (vispar.shrink < 1)
		      for (int j = 0; j < 4; j++)
			fpts[j] += (1-vispar.shrink) * (c-fpts[j]);

		    for (int ix = 0; ix <= order; ix++)
		      for (int iy = 0; iy <= order; iy++)
			{
			  double lami[4] =
			    { (1-double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (1-double(iy)/order),
			      (  double(ix)/order) * (  double(iy)/order),
			      (1-double(ix)/order) * (  double(iy)/order) };

			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) =
			      lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
			      lami[2] * fpts[2](l) + lami[3] * fpts[3](l);

			  curv.CalcElementTransformation (xl, ei, grid[ix][iy]);
			}

		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
		    for (int j = 0; j <= order; j++)
		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);

		    glMap2d(GL_MAP2_VERTEX_3,
			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
			    &grid[0][0](0));
		    glEnable(GL_MAP2_VERTEX_3);
		    glEnable(GL_AUTO_NORMAL);

		    glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0);
		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);

		    glDisable (GL_AUTO_NORMAL);
		    glDisable (GL_MAP2_VERTEX_3);
		  }






		/*
		int hoplotn = 1 << vispar.subdivisions;

		const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PYRAMID);
		const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID);

		Point<3> grid[33][33];
		Vec<3> gridn[33][33];


		glBegin (GL_TRIANGLES);

		for (int trig = 0; trig < 4; trig++)
		  {
		    Point<3> p0 = vertices[faces[trig][0]-1];
		    Point<3> p1 = vertices[faces[trig][1]-1];
		    Point<3> p2 = vertices[faces[trig][2]-1];

		    if (vispar.shrink < 1)
		      {
			static Point<3> c(0.375, 0.375, 0.25);
			p0 = c + vispar.shrink * (p0 - c);
			p1 = c + vispar.shrink * (p1 - c);
			p2 = c + vispar.shrink * (p2 - c);
		      }


		    Vec<3> taux = p0-p2;
		    Vec<3> tauy = p1-p2;
		    Vec<3> gtaux, gtauy;

		    Point<3> xl;
		    Mat<3,3> dxdxi;

		    for (int ix = 0; ix <= hoplotn; ix++)
		      for (int iy = 0; iy <= hoplotn-ix; iy++)
			{
			  for (int l = 0; l < 3; l++)
			    xl(l) =
			      (1-double(ix+iy)/hoplotn) * p2(l) +
			      (double(ix)/hoplotn) * p0(l) +
			      (double(iy)/hoplotn) * p1(l);

			  curv.CalcElementTransformation (xl, i-1, grid[ix][iy], dxdxi);

			  gtaux = dxdxi * taux;
			  gtauy = dxdxi * tauy;
			  gridn[ix][iy] = Cross (gtauy, gtaux).Normalize();
			}

		    for (int ix = 0; ix < hoplotn; ix++)
		      for (int iy = 0; iy < hoplotn-ix; iy++)
			{
			  glNormal3dv (gridn[ix][iy]);
			  glVertex3dv (grid[ix][iy]);

			  glNormal3dv (gridn[ix+1][iy]);
			  glVertex3dv (grid[ix+1][iy]);

			  glNormal3dv (gridn[ix][iy+1]);
			  glVertex3dv (grid[ix][iy+1]);

			  if (iy < hoplotn-ix-1)
			    {
			      glNormal3dv (gridn[ix][iy+1]);
			      glVertex3dv (grid[ix][iy+1]);

			      glNormal3dv (gridn[ix+1][iy]);
			      glVertex3dv (grid[ix+1][iy]);

			      glNormal3dv (gridn[ix+1][iy+1]);
			      glVertex3dv (grid[ix+1][iy+1]);
			    }
			}
		  }

		glEnd ();




		glBegin (GL_QUADS);

		for (int quad = 4; quad < 5; quad++)
		  {
		    Point<3> p0 = vertices[faces[quad][0]-1];
		    Point<3> p1 = vertices[faces[quad][1]-1];
		    Point<3> p2 = vertices[faces[quad][2]-1];
		    Point<3> p3 = vertices[faces[quad][3]-1];

		    if (vispar.shrink < 1)
		      {
			static Point<3> c(0.375, 0.375, 0.25);
			p0 = c + vispar.shrink * (p0 - c);
			p1 = c + vispar.shrink * (p1 - c);
			p2 = c + vispar.shrink * (p2 - c);
			p3 = c + vispar.shrink * (p3 - c);
		      }

		    Vec<3> taux = p1-p0;
		    Vec<3> tauy = p3-p0;
		    Vec<3> gtaux, gtauy;

		    Point<3> xl, xg;
		    Mat<3,3> dxdxi;

		    for (int ix = 0; ix <= hoplotn; ix++)
		      for (int iy = 0; iy <= hoplotn; iy++)
			{
			  Point<3> xl;
			  for (int l = 0; l < 3; l++)
			    xl(l) =
			      (1-double(ix)/hoplotn)*(1-double(iy)/hoplotn) * p0(l) +
			      (  double(ix)/hoplotn)*(1-double(iy)/hoplotn) * p1(l) +
			      (  double(ix)/hoplotn)*(  double(iy)/hoplotn) * p2(l) +
			      (1-double(ix)/hoplotn)*(  double(iy)/hoplotn) * p3(l);

			  curv.CalcElementTransformation (xl, i-1, grid[ix][iy], dxdxi);

			  gtaux = dxdxi * taux;
			  gtauy = dxdxi * tauy;
			  gridn[ix][iy] = Cross (gtauy, gtaux).Normalize();
			}

		    for (int ix = 0; ix < hoplotn; ix++)
		      for (int iy = 0; iy < hoplotn; iy++)
			{
			  glNormal3dv (gridn[ix][iy]);
			  glVertex3dv (grid[ix][iy]);

			  glNormal3dv (gridn[ix+1][iy]);
			  glVertex3dv (grid[ix+1][iy]);

			  glNormal3dv (gridn[ix+1][iy+1]);
			  glVertex3dv (grid[ix+1][iy+1]);

			  glNormal3dv (gridn[ix][iy+1]);
			  glVertex3dv (grid[ix][iy+1]);
			}
		  }

		glEnd ();
		*/


	      }
	    else
	      {



		Point3d c(0,0,0);
		if (vispar.shrink < 1)
		  {
		    for (int j = 1; j <= 5; j++)
		      {
			Point3d p = mesh->Point(el.PNum(j));
			c.X() += p.X() / 5;
			c.Y() += p.Y() / 5;
			c.Z() += p.Z() / 5;
		      }
		  }


		el.GetSurfaceTriangles (faces);

		if (el.PNum(1))
		  {
		    glBegin (GL_TRIANGLES);

		    for (int j = 1; j <= faces.Size(); j++)
		      {
			Element2d & face = faces.Elem(j);
			Point3d lp1 = mesh->Point (el.PNum(face.PNum(1)));
			Point3d lp2 = mesh->Point (el.PNum(face.PNum(2)));
			Point3d lp3 = mesh->Point (el.PNum(face.PNum(3)));
			Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
			n /= (n.Length()+1e-12);
			n *= -1;
			glNormal3d (n.X(), n.Y(), n.Z());

			if (vispar.shrink < 1)
			  {
			    lp1 = c + vispar.shrink * (lp1 - c);
			    lp2 = c + vispar.shrink * (lp2 - c);
			    lp3 = c + vispar.shrink * (lp3 - c);
			  }

			glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
			glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
			glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
		      }

		    glEnd();
		  }
	      }
	  }
      }
    glEndList ();
  }

  void VisualSceneMesh :: BuildBadelList()
  {
    ;
  }

  void VisualSceneMesh :: BuildIdentifiedList()
  {
    ;
  }

  void VisualSceneMesh :: BuildDomainSurfList()
  {
    if (domainsurflist)
      glDeleteLists (domainsurflist, 1);

    domainsurflist = glGenLists (1);
    glNewList (domainsurflist, GL_COMPILE);

    int i, j;
    glLineWidth (1.0f);

    glDisable (GL_COLOR_MATERIAL);

    for (i = 1; i <= mesh->GetNSE(); i++)
      {
	Element2d el = mesh->SurfaceElement (i);

	int drawel = 1;
	for (j = 1; j <= el.GetNP(); j++)
	  {
	    if (!el.PNum(j))
	      drawel = 0;
	  }

	if (!drawel)
	  continue;

	if (el.GetIndex() < 1 || el.GetIndex() > mesh->GetNFD())
	  continue;
	int domin = mesh->GetFaceDescriptor(el.GetIndex()).DomainIn();
	int domout = mesh->GetFaceDescriptor(el.GetIndex()).DomainOut();

	int fac;
	if (domin == vispar.drawdomainsurf)
	  fac = 1;
	else if (domout == vispar.drawdomainsurf)
	  fac = -1;
	else
	  continue;


	GLfloat matcol[] = { 1, 0, 0, 1 };
	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol);


	if (el.GetNP() == 3)
	  {
	    glBegin (GL_TRIANGLES);

	    const Point3d & lp1 = mesh->Point (el.PNum(1));
	    const Point3d & lp2 = mesh->Point (el.PNum(2));
	    const Point3d & lp3 = mesh->Point (el.PNum(3));
	    Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
	    n /= ( fac * (n.Length()+1e-12));
	    glNormal3d (n.X(), n.Y(), n.Z());

	    if (!vispar.colormeshsize)
	      {
		glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
	      }
	    glEnd();
	  }
	else if (el.GetNP() == 4)
	  {
	    glBegin (GL_QUADS);

	    const Point3d & lp1 = mesh->Point (el.PNum(1));
	    const Point3d & lp2 = mesh->Point (el.PNum(2));
	    const Point3d & lp3 = mesh->Point (el.PNum(4));
	    const Point3d & lp4 = mesh->Point (el.PNum(3));
	    Vec3d n = Cross (Vec3d (lp1, lp2),
			     Vec3d (lp1, Center (lp3, lp4)));
	    n /= (fac * (n.Length()+1e-12));
	    glNormal3d (n.X(), n.Y(), n.Z());
	    glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
	    glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
	    glVertex3d (lp4.X(), lp4.Y(), lp4.Z());
	    glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
	    glEnd();
	  }
	else if (el.GetNP() == 6)
	  {
	    glBegin (GL_TRIANGLES);
	    static int trigs[4][3] = {
	      { 1, 6, 5 },
	      { 2, 4, 6 },
	      { 3, 5, 4 },
	      { 4, 5, 6 } };

	    for (j = 0; j < 4; j++)
	      {
		const Point3d & lp1 = mesh->Point (el.PNum(trigs[j][0]));
		const Point3d & lp2 = mesh->Point (el.PNum(trigs[j][1]));
		const Point3d & lp3 = mesh->Point (el.PNum(trigs[j][2]));
		Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3));
		n /= (fac * (n.Length() + 1e-12));
		glNormal3d (n.X(), n.Y(), n.Z());
		glVertex3d (lp1.X(), lp1.Y(), lp1.Z());
		glVertex3d (lp2.X(), lp2.Y(), lp2.Z());
		glVertex3d (lp3.X(), lp3.Y(), lp3.Z());
	      }
	    glEnd();
	  }
      }
    glEndList ();
  }








  void VisualSceneMesh :: MouseDblClick (int px, int py)
  {
     if(lock)
      {
	lock->UnLock();
	delete lock;
	lock = NULL;
      }



    MouseDblClickSelect(px,py,clipplane,backcolor,transformationmat,center,rad,
			filledlist,selelement,selface,seledge,selpoint,selpoint2,locpi);

    selecttimestamp = NextTimeStamp();

    /*
    int i, hits;

    // select surface triangle by mouse click

    GLuint selbuf[10000];
    glSelectBuffer (10000, selbuf);


    glRenderMode (GL_SELECT);

    GLint viewport[4];
    glGetIntegerv (GL_VIEWPORT, viewport);


    glMatrixMode (GL_PROJECTION);
    glPushMatrix();

    GLdouble projmat[16];
    glGetDoublev (GL_PROJECTION_MATRIX, projmat);

    glLoadIdentity();
    gluPickMatrix (px, viewport[3] - py, 1, 1, viewport);
    glMultMatrixd (projmat);



    glClearColor(backcolor, backcolor, backcolor, 1.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode (GL_MODELVIEW);

    glPushMatrix();
    glMultMatrixf (transformationmat);


    //  SetClippingPlane();

    glInitNames();
    glPushName (1);

    glPolygonOffset (1, 1);
    glEnable (GL_POLYGON_OFFSET_FILL);

    glDisable(GL_CLIP_PLANE0);

    if (vispar.clipenable)
      {
	Vec<3> n(clipplane[0], clipplane[1], clipplane[2]);
	double len = Abs(n);
	double mu = -clipplane[3] / (len*len);
	Point<3> p (mu * n);
	n /= len;
	Vec<3> t1 = n.GetNormal ();
	Vec<3> t2 = Cross (n, t1);

	double xi1mid = (center - p) * t1;
	double xi2mid = (center - p) * t2;

	glLoadName (0);
	glBegin (GL_QUADS);
	glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid-rad) * t2);
	glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid-rad) * t2);
	glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid+rad) * t2);
	glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid+rad) * t2);
	glEnd ();
      }

    //  SetClippingPlane();

    glCallList (filledlist);

    glDisable (GL_POLYGON_OFFSET_FILL);

    glPopName();

    glMatrixMode (GL_PROJECTION);
    glPopMatrix();

    glMatrixMode (GL_MODELVIEW);
    glPopMatrix();

    glFlush();


    hits = glRenderMode (GL_RENDER);

    //  cout << "hits = " << hits << endl;

    int minname = 0;
    GLuint mindepth = 0;

    // find clippingplane
    GLuint clipdepth = 0; // GLuint(-1);

    for (i = 0; i < hits; i++)
      {
	int curname = selbuf[4*i+3];
	if (!curname) clipdepth = selbuf[4*i+1];
      }

    for (i = 0; i < hits; i++)
      {
	int curname = selbuf[4*i+3];
	GLuint curdepth = selbuf[4*i+1];

	if (curname && (curdepth > clipdepth) &&
	    (curdepth < mindepth || !minname))
	  {
	    mindepth = curdepth;
	    minname = curname;
	  }
      }

    seledge = -1;
    if (minname)
      {
	const Element2d & sel = mesh->SurfaceElement(minname);


	cout << "select element " << minname
	     << " on face " << sel.GetIndex() << endl;
	cout << "Nodes: ";
	for (i = 1; i <= sel.GetNP(); i++)
	  cout << sel.PNum(i) << " ";
	cout << endl;

	selelement = minname;
	selface = mesh->SurfaceElement(minname).GetIndex();

	locpi = (locpi % sel.GetNP()) + 1;
	selpoint2 = selpoint;
	selpoint = sel.PNum(locpi);
	cout << "selected point " << selpoint
	     << ", pos = " << mesh->Point (selpoint)
	     << endl;

	for (i = 1; i <= mesh->GetNSeg(); i++)
	  {
	    const Segment & seg = mesh->LineSegment(i);
	    if (seg[0] == selpoint && seg[1] == selpoint2 ||
		seg[1] == selpoint && seg[0] == selpoint2)
	      {
		seledge = seg.edgenr;
		cout << "seledge = " << seledge << endl;
	      }
	  }

      }
    else
      {
	selface = -1;
	selelement = -1;
	selpoint = -1;
	selpoint2 = -1;
      }

    glDisable(GL_CLIP_PLANE0);

    selecttimestamp = NextTimeStamp();
    */

  }





  void MouseDblClickSelect (const int px, const int py,
			    const GLdouble * clipplane, const GLdouble backcolor,
			    const float * transformationmat,
			    const Point3d & center,
			    const double rad,
			    const int displaylist,
			    int & selelement, int & selface, int & seledge, int & selpoint,
			    int & selpoint2, int & locpi)
  {
    int i, hits;

    // select surface triangle by mouse click

    GLuint selbuf[10000];
    glSelectBuffer (10000, selbuf);


    glRenderMode (GL_SELECT);

    GLint viewport[4];
    glGetIntegerv (GL_VIEWPORT, viewport);


    glMatrixMode (GL_PROJECTION);
    glPushMatrix();

    GLdouble projmat[16];
    glGetDoublev (GL_PROJECTION_MATRIX, projmat);

    glLoadIdentity();
    gluPickMatrix (px, viewport[3] - py, 1, 1, viewport);
    glMultMatrixd (projmat);



    glClearColor(backcolor, backcolor, backcolor, 1.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode (GL_MODELVIEW);

    glPushMatrix();
    glMultMatrixf (transformationmat);


    //  SetClippingPlane();

    glInitNames();
    glPushName (1);

    glPolygonOffset (1, 1);
    glEnable (GL_POLYGON_OFFSET_FILL);

    glDisable(GL_CLIP_PLANE0);

    if (vispar.clipenable)
      {
	Vec<3> n(clipplane[0], clipplane[1], clipplane[2]);
	double len = Abs(n);
	double mu = -clipplane[3] / (len*len);
	Point<3> p (mu * n);
	n /= len;
	Vec<3> t1 = n.GetNormal ();
	Vec<3> t2 = Cross (n, t1);

	double xi1mid = (center - p) * t1;
	double xi2mid = (center - p) * t2;

	glLoadName (0);
	glBegin (GL_QUADS);
	glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid-rad) * t2);
	glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid-rad) * t2);
	glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid+rad) * t2);
	glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid+rad) * t2);
	glEnd ();
      }

    //  SetClippingPlane();

    glCallList (displaylist);

    glDisable (GL_POLYGON_OFFSET_FILL);

    glPopName();

    glMatrixMode (GL_PROJECTION);
    glPopMatrix();

    glMatrixMode (GL_MODELVIEW);
    glPopMatrix();

    glFlush();


    hits = glRenderMode (GL_RENDER);

    //cout << "hits = " << hits << endl;

    int minname = 0;
    GLuint mindepth = 0;

    // find clippingplane
    GLuint clipdepth = 0; // GLuint(-1);

    for (i = 0; i < hits; i++)
      {
	int curname = selbuf[4*i+3];
	if (!curname) clipdepth = selbuf[4*i+1];
      }

    for (i = 0; i < hits; i++)
      {
	int curname = selbuf[4*i+3];
	GLuint curdepth = selbuf[4*i+1];
	/*
	  cout << selbuf[4*i] << " " << selbuf[4*i+1] << " "
	  << selbuf[4*i+2] << " " << selbuf[4*i+3] << endl;
	*/
	if (curname && (curdepth > clipdepth) &&
	    (curdepth < mindepth || !minname))
	  {
	    mindepth = curdepth;
	    minname = curname;
	  }
      }

    seledge = -1;
    if (minname)
      {
	const Element2d & sel = mesh->SurfaceElement(minname);


	cout << "select element " << minname
	     << " on face " << sel.GetIndex() << endl;
	cout << "Nodes: ";
	for (i = 1; i <= sel.GetNP(); i++)
	  cout << sel.PNum(i) << " ";
	cout << endl;

	selelement = minname;
	selface = mesh->SurfaceElement(minname).GetIndex();

	locpi = (locpi % sel.GetNP()) + 1;
	selpoint2 = selpoint;
	selpoint = sel.PNum(locpi);
	cout << "selected point " << selpoint
	     << ", pos = " << mesh->Point (selpoint)
	     << endl;

	for (i = 1; i <= mesh->GetNSeg(); i++)
	  {
	    const Segment & seg = mesh->LineSegment(i);
	    if ( (seg[0] == selpoint && seg[1] == selpoint2) ||
		 (seg[1] == selpoint && seg[0] == selpoint2) )
	      {
		seledge = seg.edgenr;
		cout << "seledge = " << seledge << endl;
	      }
	  }

      }
    else
      {
	selface = -1;
	selelement = -1;
	selpoint = -1;
	selpoint2 = -1;
      }

    glDisable(GL_CLIP_PLANE0);



#ifdef PARALLELGL
    vsmesh.Broadcast ();
#endif
  }


  void VisualSceneMesh :: SetSelectedFace (int asf)
  {
    selface = asf;
    selecttimestamp = GetTimeStamp();
  }




}



#endif // NOTCL