/contrib/Netgen/libsrc/meshing/bisect.cpp

Large files files are truncated, but you can click here to view the full file

#include <mystdlib.h>
#include "meshing.hpp"

#define noDEBUG


namespace netgen
{
  //#include "../interface/writeuser.hpp"
  class MarkedTet;
  class MarkedPrism;
  class MarkedIdentification;
  class MarkedTri;
  class MarkedQuad;
  
  typedef MoveableArray<MarkedTet> T_MTETS;
  typedef MoveableArray<MarkedPrism> T_MPRISMS;
  typedef MoveableArray<MarkedIdentification> T_MIDS;
  typedef MoveableArray<MarkedTri> T_MTRIS;
  typedef MoveableArray<MarkedQuad> T_MQUADS;

  
  
  class MarkedTet
  {
  public:
    /// pnums of tet
    PointIndex pnums[4];
    /// material number
    int matindex;
    /// element marked for refinement
    /// marked = 1: marked by element marker, marked = 2 due to closure
    unsigned int marked:2;
    /// flag of Arnold-Mukherjee algorithm
    unsigned int flagged:1;
    /// tetedge (local coordinates 0..3)
    unsigned int tetedge1:3;
    unsigned int tetedge2:3;
    // marked edge of faces
    // face_j : face without node j,
    // mark_k : edge without node k
    
    char faceedges[4];
    // unsigned char faceedges[4];
    bool incorder;
    unsigned int order:6;

    MarkedTet()
    { 
      for (int i = 0; i < 4; i++) { faceedges[i] = 255; }
    }
  };

  ostream & operator<< (ostream & ost, const MarkedTet & mt)
  {
    for(int i=0; i<4; i++)
      ost << mt.pnums[i] << " ";

    ost << mt.matindex << " " << int(mt.marked) << " " << int(mt.flagged) << " " << int(mt.tetedge1) << " " << int(mt.tetedge2) << " ";
    
    ost << "faceedges = ";
    for(int i=0; i<4; i++)
      ost << int(mt.faceedges[i]) << " ";

    ost << " order = ";
    ost << mt.incorder << " " << int(mt.order) << "\n";
    return ost;
  }
  istream & operator>> (istream & ost, MarkedTet & mt)
  {
    for(int i=0; i<4; i++)
      ost >> mt.pnums[i];

    ost >> mt.matindex;

    int auxint;
    ost >> auxint;
    mt.marked = auxint;
    ost >> auxint;
    mt.flagged = auxint;
    ost >> auxint;
    mt.tetedge1 = auxint;
    ost >> auxint;
    mt.tetedge2 = auxint;
    
    char auxchar;

    for(int i=0; i<4; i++)
      {
	ost >> auxchar;
	mt.faceedges[i] = auxchar;
      }

    ost >> mt.incorder;
    ost >> auxint;
    mt.order = auxint;
    return ost;
  }

  class MarkedPrism
  {
  public:
    /// 6 point numbers
    PointIndex pnums[6];
    /// material number
    int matindex;
    /// marked for refinement
    int marked;
    /// edge without node k (0,1,2)
    int markededge;

    bool incorder;
    unsigned int order:6;
  };

  
  ostream & operator<< (ostream & ost, const MarkedPrism & mp)
  {
    for(int i=0; i<6; i++)
      ost << mp.pnums[i] << " ";

    ost << mp.matindex << " " << mp.marked << " " << mp.markededge << " " << mp.incorder << " " << int(mp.order) << "\n";
    return ost;
  }
  istream & operator>> (istream & ist, MarkedPrism & mp)
  {
    for(int i=0; i<6; i++)
      ist >> mp.pnums[i];

    ist >> mp.matindex >> mp.marked >> mp.markededge >> mp.incorder;
    int auxint;
    ist >> auxint;
    mp.order = auxint;
    return ist;
  }


  class MarkedIdentification
  {
  public:
    // number of points of one face (3 or 4)
    int np;
    /// 6 or 8 point numbers
    PointIndex pnums[8];
    /// marked for refinement
    int marked;
    /// edge starting with node k (0,1,2, or 3)
    int markededge;

    bool incorder;
    unsigned int order:6;
  };
    
  
  ostream & operator<< (ostream & ost, const MarkedIdentification & mi)
  {
    ost << mi.np << " ";
    for(int i=0; i<2*mi.np; i++)
      ost << mi.pnums[i] << " ";
    ost << mi.marked << " " << mi.markededge << " " << mi.incorder << " " << int(mi.order) << "\n";
    return ost;
  }
  istream & operator>> (istream & ist, MarkedIdentification & mi)
  {
    ist >> mi.np;
    for(int i=0; i<2*mi.np; i++)
      ist >> mi.pnums[i];
    ist >> mi.marked >> mi.markededge >> mi.incorder;
    int auxint;
    ist >> auxint;
    mi.order = auxint;
    return ist;
  }
  




  class MarkedTri
  {
  public:
    /// three point numbers
    PointIndex pnums[3];
    /// three geominfos
    PointGeomInfo pgeominfo[3];
    /// marked for refinement
    int marked;
    /// edge without node k
    int markededge;
    /// surface id
    int surfid;

    bool incorder;
    unsigned int order:6;
  };
  
  ostream & operator<< (ostream & ost, const MarkedTri & mt)
  {
    for(int i=0; i<3; i++)
      ost << mt.pnums[i] << " ";
    for(int i=0; i<3; i++)
      ost << mt.pgeominfo[i] << " ";
    ost << mt.marked << " " << mt.markededge << " " << mt.surfid << " " << mt.incorder << " " << int(mt.order) << "\n";
    return ost;
  } 
  istream & operator>> (istream & ist, MarkedTri & mt)
  {
    for(int i=0; i<3; i++)
      ist >> mt.pnums[i];
    for(int i=0; i<3; i++)
      ist >> mt.pgeominfo[i];
    ist >> mt.marked >> mt.markededge >> mt.surfid >> mt.incorder;
    int auxint;
    ist >> auxint;
    mt.order = auxint;
    return ist;
  }
    


  class MarkedQuad
  {
  public:
    /// point numbers
    PointIndex pnums[4];
    ///
    PointGeomInfo pgeominfo[4];
    /// marked for refinement
    int marked;
    /// marked edge: 0/2 = vertical, 1/3 = horizontal
    int markededge;
    /// surface id
    int surfid;

    bool incorder;
    unsigned int order:6;
  };

  ostream & operator<< (ostream & ost, const MarkedQuad & mt)
  {
    for(int i=0; i<4; i++)
      ost << mt.pnums[i] << " ";
    for(int i=0; i<4; i++)
      ost << mt.pgeominfo[i] << " ";
    ost << mt.marked << " " << mt.markededge << " " << mt.surfid << " " << mt.incorder << " " << int(mt.order) << "\n";
    return ost;
  } 
  istream & operator>> (istream & ist, MarkedQuad & mt)
  {
    for(int i=0; i<4; i++)
      ist >> mt.pnums[i];
    for(int i=0; i<4; i++)
      ist >> mt.pgeominfo[i];
    ist >> mt.marked >> mt.markededge >> mt.surfid >> mt.incorder;
    int auxint;
    ist >> auxint;
    mt.order = auxint;
    return ist;
  }




  void PrettyPrint(ostream & ost, const MarkedTet & mt)
  {
    int te1 = mt.tetedge1;
    int te2 = mt.tetedge2;
    int order = mt.order;

    ost << "MT: " << mt.pnums[0] << " - " << mt.pnums[1] << " - " 
	<< mt.pnums[2] << " - " << mt.pnums[3] << endl
	<< "marked edge: " << te1 << " - " << te2
	<< ", order = " << order << endl;
    //for (int k = 0; k < 4; k++)
    //  ost << int(mt.faceedges[k]) << "  ";
    for (int k = 0; k < 4; k++)
      {
	ost << "face";
	for (int j=0; j<4; j++)
	  if(j != k)
	    ost << " " << mt.pnums[j];
	for(int i=0; i<3; i++)
	  for(int j=i+1; j<4; j++)
	    if(i != k && j != k && int(mt.faceedges[k]) == 6-k-i-j)
	      ost << " marked edge " << mt.pnums[i] << " " << mt.pnums[j] << endl;
      }
    ost << endl;
  }




  int BTSortEdges (const Mesh & mesh,
		   const Array< Array<int,PointIndex::BASE>* > & idmaps,
		   INDEX_2_CLOSED_HASHTABLE<int> & edgenumber)
  {
    PrintMessage(4,"sorting ... ");

    //  if (mesh.PureTetMesh())
    if (1)
      {
	// new, fast version
      
	Array<INDEX_2> edges;
	Array<int> eclasses;
      
	int i, j, k;
	int cntedges = 0;
	int go_on;
	int ned(0);
      
	// enumerate edges:
	for (i = 1; i <= mesh.GetNE(); i++)
	  {
	    const Element & el = mesh.VolumeElement (i);
	    static int tetedges[6][2] =
	      { { 1, 2 },
		{ 1, 3 },
		{ 1, 4 },
		{ 2, 3 },
		{ 2, 4 },
		{ 3, 4 } } ;
	    static int prismedges[9][2] =
	      { { 1, 2 },
		{ 1, 3 },
		{ 2, 3 },
		{ 4, 5 },
		{ 4, 6 },
		{ 5, 6 },
		{ 1, 4 },
		{ 2, 5 },
		{ 3, 6 } };
	    int pyramidedges[6][2] =
	      { { 1, 2 },
		{ 3, 4 },
		{ 1, 5 },
		{ 2, 5 },
		{ 3, 5 },
		{ 4, 5 } };
	  
	    int (*tip)[2] = NULL;
	  
	    switch (el.GetType())
	      {
	      case TET:
	      case TET10:
		{
		  tip = tetedges;
		  ned = 6;
		  break;
		}
	      case PRISM:
	      case PRISM12:
		{
		  tip = prismedges;
		  ned = 6;
		  break;
		}
	      case PYRAMID:
		{
		  tip = pyramidedges;
		  ned = 6;
		  break;
		}
              default:
                throw NgException("Bisect, element type not handled in switch");
	      }
	      
	    for (j = 0; j < ned; j++)
	      {
		INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1]));
		i2.Sort();
		//(*testout) << "edge " << i2 << endl;
		if (!edgenumber.Used(i2))
		  {
		    cntedges++;
		    edges.Append (i2);
		    edgenumber.Set(i2, cntedges);
		  }
	      }
	  }
      
	// additional surface edges:
	for (i = 1; i <= mesh.GetNSE(); i++)
	  {
	    const Element2d & el = mesh.SurfaceElement (i);
	    static int trigedges[3][2] =
	      { { 1, 2 },
		{ 2, 3 },
		{ 3, 1 } };

	    static int quadedges[4][2] =
	      { { 1, 2 },
		{ 2, 3 },
		{ 3, 4 },
		{ 4, 1 } };


	    int (*tip)[2] = NULL;
	  
	    switch (el.GetType())
	      {
	      case TRIG:
	      case TRIG6:
		{
		  tip = trigedges;
		  ned = 3;
		  break;
		}
	      case QUAD:
	      case QUAD6:
		{
		  tip = quadedges;
		  ned = 4;
		  break;
		}
	      default:
		{
		  cerr << "Error: Sort for Bisect, SE has " << el.GetNP() << " points" << endl;
		  ned = 0;
		}
	      }
	      
	    for (j = 0; j < ned; j++)
	      {
		INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1]));
		i2.Sort();
		if (!edgenumber.Used(i2))
		  {
		    cntedges++;
		    edges.Append (i2);
		    edgenumber.Set(i2, cntedges);
		  }
	      }
	  }





	eclasses.SetSize (cntedges);
	for (i = 1; i <= cntedges; i++)
	  eclasses.Elem(i) = i;

	// identify edges in element stack
	do
	  {
	    go_on = 0;
	    for (i = 1; i <= mesh.GetNE(); i++)
	      {
		const Element & el = mesh.VolumeElement (i);	     
	      
		if (el.GetType() != PRISM &&
		    el.GetType() != PRISM12 &&
		    el.GetType() != PYRAMID)
		  continue;

		int prismpairs[3][4] =
		  { { 1, 2, 4, 5 },
		    { 2, 3, 5, 6 },
		    { 1, 3, 4, 6 } };
	      
		int pyramidpairs[3][4] =
		  { { 1, 2, 4, 3 },
		    { 1, 5, 4, 5 },
		    { 2, 5, 3, 5 } };
                
		int (*pairs)[4] = NULL;
		switch (el.GetType())
		  {
		  case PRISM:
		  case PRISM12:
		    {
		      pairs = prismpairs;
		      break;
		    }
		  case PYRAMID:
		    {
		      pairs = pyramidpairs;
		      break;
		    }
                  default:
                    throw NgException("Bisect, element type not handled in switch, 2");
		  }

		for (j = 0; j < 3; j++)
		  {
		    INDEX_2 e1 (el.PNum(pairs[j][0]), 
				el.PNum(pairs[j][1]));
		    INDEX_2 e2 (el.PNum(pairs[j][2]), 
				el.PNum(pairs[j][3]));
		    e1.Sort();
		    e2.Sort();
		      
		    int eclass1 = edgenumber.Get (e1);
		    int eclass2 = edgenumber.Get (e2);

		    //		  (*testout) << "identify edges " << eclass1 << "-" << eclass2 << endl;

		    if (eclasses.Get(eclass1) >
			eclasses.Get(eclass2))
		      {
			eclasses.Elem(eclass1) = 
			  eclasses.Get(eclass2);
			go_on = 1;
		      }
		    else if (eclasses.Get(eclass2) >
			     eclasses.Get(eclass1))
		      {
			eclasses.Elem(eclass2) = 
			  eclasses.Get(eclass1);
			go_on = 1;
		      }
		  }
	      }

	    for(SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
	      {
		const Element2d & el2d = mesh[sei];

		for(i = 0; i < el2d.GetNP(); i++)
		  {
		    INDEX_2 e1(el2d[i], el2d[(i+1) % el2d.GetNP()]);
		    e1.Sort();
		    INDEX_2 e2;
		    
		    for(k = 0; k < idmaps.Size(); k++)
		      {
			e2.I1() = (*idmaps[k])[e1.I1()];
			e2.I2() = (*idmaps[k])[e1.I2()];
			
			if(e2.I1() == 0 || e2.I2() == 0 ||
			   e1.I1() == e2.I1() || e1.I2() == e2.I2())
			  continue;
			
			e2.Sort();
			if(!edgenumber.Used(e2))
			  continue;
			

			int eclass1 = edgenumber.Get (e1);
			int eclass2 = edgenumber.Get (e2);
			
			if (eclasses.Get(eclass1) >
			    eclasses.Get(eclass2))
			  {
			    eclasses.Elem(eclass1) = 
			      eclasses.Get(eclass2);


			    go_on = 1;
			  }
			else if (eclasses.Get(eclass2) >
				 eclasses.Get(eclass1))
			  {
			    eclasses.Elem(eclass2) = 
			      eclasses.Get(eclass1);
			    go_on = 1;
			  }
		      }		      
		  }
		
	      }

	  }
	while (go_on);

// 	for (i = 1; i <= cntedges; i++)
// 	  {
// 	    (*testout) << "edge " << i << ": " 
// 		       << edges.Get(i).I1() << "-" << edges.Get(i).I2()
// 		       << ", class = " << eclasses.Get(i) << endl;
// 	  }
	
	// compute classlength:
	Array<double> edgelength(cntedges);

	/*
	for (i = 1; i <= cntedges; i++)
	  edgelength.Elem(i) = 1e20;
	*/

	for (i = 1; i <= cntedges; i++)
	  {
	    INDEX_2 edge = edges.Get(i);
	    double elen = Dist (mesh.Point(edge.I1()),
				mesh.Point(edge.I2()));
	    edgelength.Elem (i) = elen;
	  }

	/*
	  for (i = 1; i <= mesh.GetNE(); i++)
	  {
	  const Element & el = mesh.VolumeElement (i);
	  
	  if (el.GetType() == TET)
	  {
	  for (j = 1; j <= 3; j++)
	  for (k = j+1; k <= 4; k++)
	  {
	  INDEX_2 i2(el.PNum(j), el.PNum(k));
	  i2.Sort();
		    
	  int enr = edgenumber.Get(i2);
	  double elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2()));
	  if (elen < edgelength.Get(enr))
	  edgelength.Set (enr, elen);
	  }
	  }
	  else if (el.GetType() == PRISM)
	  {
	  for (j = 1; j <= 3; j++)
	  {
	  k = (j % 3) + 1;
		  
	  INDEX_2 i2(el.PNum(j), el.PNum(k));
	  i2.Sort();
		  
	  int enr = edgenumber.Get(i2);
	  double elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2()));
	  if (elen < edgelength.Get(enr))
	  edgelength.Set (enr, elen);
		  
	  i2 = INDEX_2(el.PNum(j+3), el.PNum(k+3));
	  i2.Sort();
		  
	  enr = edgenumber.Get(i2);
	  elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2()));
	  if (elen < edgelength.Get(enr))
	  edgelength.Set (enr, elen);
		  
	  if (!edgenumber.Used(i2))
	  {
	  cntedges++;
	  edgenumber.Set(i2, cntedges);
	  }
	  i2 = INDEX_2(el.PNum(j), el.PNum(j+3));
	  i2.Sort();
		  
	  enr = edgenumber.Get(i2);
	  elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2()));
	  if (elen < edgelength.Get(enr))
	  edgelength.Set (enr, elen);
	  }
	  }
	  }
	*/

      
	for (i = 1; i <= cntedges; i++)
	  {
	    if (eclasses.Get(i) != i)
	      {
		if (edgelength.Get(i) < edgelength.Get(eclasses.Get(i)))
		  edgelength.Elem(eclasses.Get(i)) = edgelength.Get(i);
		edgelength.Elem(i) = 1e20;
	      }
	  }


	TABLE<int> eclasstab(cntedges);
	for (i = 1; i <= cntedges; i++)
	  eclasstab.Add1 (eclasses.Get(i), i);


	// sort edges:
	Array<int> sorted(cntedges);
      
	QuickSort (edgelength, sorted);
      
	int cnt = 0;
	for (i = 1; i <= cntedges; i++)
	  {
	    int ii = sorted.Get(i);
	    for (j = 1; j <= eclasstab.EntrySize(ii); j++)
	      {
		cnt++;
		edgenumber.Set (edges.Get(eclasstab.Get(ii, j)), cnt); 
	      }
	  }
	return cnt;
      }

    else
    
      {
	// old version
      
	int i, j;
	int cnt = 0;
	int found;
	double len2, maxlen2;
	INDEX_2 ep;
      
	// sort edges by length, parallel edges (on prisms)
	// are added in blocks
      
	do
	  {
	    found = 0;
	    maxlen2 = 1e30;
	  
	    for (i = 1; i <= mesh.GetNE(); i++)
	      {
		const Element & el = mesh.VolumeElement (i);
		int ned;
		int tetedges[6][2] =
		  { { 1, 2 },
		    { 1, 3 },
		    { 1, 4 },
		    { 2, 3 },
		    { 2, 4 },
		    { 3, 4 } };
		int prismedges[6][2] =
		  { { 1, 2 },
		    { 1, 3 },
		    { 2, 4 },
		    { 4, 5 },
		    { 4, 6 },
		    { 5, 6 } };
		int pyramidedges[6][2] =
		  { { 1, 2 },
		    { 3, 4 },
		    { 1, 5 },
		    { 2, 5 },
		    { 3, 5 },
		    { 4, 5 } };

		int (*tip)[2];

		switch (el.GetType())
		  {
		  case TET:
		    {
		      tip = tetedges;
		      ned = 6;
		      break;
		    }
		  case PRISM:
		    {
		      tip = prismedges;
		      ned = 6;
		      break;
		    }
		  case PYRAMID:
		    {
		      tip = pyramidedges;
		      ned = 6;
		      break;
		    }
                  default:
                    throw NgException("Bisect, element type not handled in switch, 3");
		  }
	      
		for (j = 0; j < ned; j++)
		  {
		    INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1]));
		    i2.Sort();
		    if (!edgenumber.Used(i2))
		      {
			len2 = Dist (mesh.Point (i2.I1()),
				     mesh.Point (i2.I2()));
			if (len2 < maxlen2)
			  {
			    maxlen2 = len2;
			    ep = i2;
			    found = 1;
			  }
		      }
		  }
	      }
	    if (found)
	      {
		cnt++;
		edgenumber.Set (ep, cnt);
	      
	      
		// find connected edges:
		int go_on = 0;
		do
		  {
		    go_on = 0;
		    for (i = 1; i <= mesh.GetNE(); i++)
		      {
			const Element & el = mesh.VolumeElement (i);	      
			if (el.GetNP() != 6) continue;

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

			int pyramidpairs[3][4] =
			  { { 1, 2, 4, 3 },
			    { 1, 5, 4, 5 },
			    { 2, 5, 3, 5 } };
		      
			int (*pairs)[4];
			switch (el.GetType())
			  {
			  case PRISM:
			    {
			      pairs = prismpairs;
			      break;
			    }
			  case PYRAMID:
			    {
			      pairs = pyramidpairs;
			      break;
			    }
                          default:
                            throw NgException("Bisect, element type not handled in switch, 3a");
			  }

			for (j = 0; j < 3; j++)
			  {
			    INDEX_2 e1 (el.PNum(pairs[j][0]), 
					el.PNum(pairs[j][1]));
			    INDEX_2 e2 (el.PNum(pairs[j][2]), 
					el.PNum(pairs[j][3]));
			    e1.Sort();
			    e2.Sort();
			  
			    int used1 = edgenumber.Used (e1);
			    int used2 = edgenumber.Used (e2);
			  
			    if (used1 && !used2)
			      {
				cnt++;
				edgenumber.Set (e2, cnt);
				go_on = 1;
			      }
			    if (used2 && !used1)
			      {
				cnt++;
				edgenumber.Set (e1, cnt);
				go_on = 1;
			      }
			  }
		      }
		  }
		while (go_on);
	      }
	  }
	while (found);

	return cnt;
      }
  }




  void BTDefineMarkedTet (const Element & el,
			  INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
			  MarkedTet & mt)
  {
    int i, j, k;
    for (i = 0; i < 4; i++)
      mt.pnums[i] = el[i];

    mt.marked = 0;
    mt.flagged = 0;

    mt.incorder = 0;
    mt.order = 1;
  
    int val = 0;
    // find marked edge of tet:
    for (i = 0; i < 3; i++)
      for (j = i+1; j < 4; j++)
	{
	  INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
	  i2.Sort();
	  int hval = edgenumber.Get(i2);
	  if (hval > val)
	    {
	      val = hval;
	      mt.tetedge1 = i;
	      mt.tetedge2 = j;    
	    }
	}


    // find marked edges of faces:
    for (k = 0; k < 4; k++)
      {
	val = 0;
	for (i = 0; i < 3; i++)
	  for (j = i+1; j < 4; j++)
	    if (i != k && j != k)
	      {
		INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
		i2.Sort();
		int hval = edgenumber.Get(i2);
		if (hval > val)
		  {
		    val = hval;
                    int hi = 6 - k - i - j;
                    mt.faceedges[k] = char(hi);
		  }
	      }
      }
  }




  void BTDefineMarkedPrism (const Element & el,
			    INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
			    MarkedPrism & mp)
  {
    int i, j;

    if (el.GetType() == PRISM ||
	el.GetType() == PRISM12)
      {
	for (i = 0; i < 6; i++)
	  mp.pnums[i] = el[i];
      }
    else if (el.GetType() == PYRAMID)
      {
	static int map[6] = 
	  { 1, 2, 5, 4, 3, 5 };
	for (i = 0; i < 6; i++)
	  mp.pnums[i] = el.PNum(map[i]);
      }
    else if (el.GetType() == TET ||
	     el.GetType() == TET10)
      {
	static int map[6] = 
	  { 1, 4, 3, 2, 4, 3 };
	for (i = 0; i < 6; i++)
	  mp.pnums[i] = el.PNum(map[i]);
      
      }
    else
      {
	PrintSysError ("Define marked prism called for non-prism and non-pyramid");
      }
  


    mp.marked = 0;

    mp.incorder = 0;
    mp.order = 1;

    int val = 0;
    for (i = 0; i < 2; i++)
      for (j = i+1; j < 3; j++)
	{
	  INDEX_2 i2(mp.pnums[i], mp.pnums[j]);
	  i2.Sort();
	  int hval = edgenumber.Get(i2);
	  if (hval > val)
	    {
	      val = hval;
	      mp.markededge = 3 - i - j;
	    }
	}
  }



  bool BTDefineMarkedId(const Element2d & el, 
			INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 
			const Array<int,PointIndex::BASE> & idmap,
			MarkedIdentification & mi)
  {

    bool identified = true;
    mi.np = el.GetNP();
    int min1(0),min2(0);
    for(int j = 0; identified && j < mi.np; j++)
      {
	mi.pnums[j] = el[j];
	mi.pnums[j+mi.np] = idmap[el[j]];

	if(j == 0 || el[j] < min1)
	  min1 = el[j];
	if(j == 0 || mi.pnums[j+mi.np] < min2)
	  min2 = mi.pnums[j+mi.np];

	identified = (mi.pnums[j+mi.np] != 0 && mi.pnums[j+mi.np] != mi.pnums[j]);
      }

    identified = identified && (min1 < min2);

    if(identified)
      {
	mi.marked = 0;
	
	mi.incorder = 0;
	mi.order = 1;

	int val = 0;
	for (int i = 0; i < mi.np; i++)
	  {
	    INDEX_2 i2(mi.pnums[i], mi.pnums[(i+1)%mi.np]);
	    i2.Sort();
	    int hval = edgenumber.Get(i2);
	    if (hval > val)
	      {
		val = hval;
		mi.markededge = i;
	      }
	  }
      }

    return identified;
  }


  void BTDefineMarkedTri (const Element2d & el,
			  INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
			  MarkedTri & mt)
  {
    int i, j;
    for (i = 0; i < 3; i++)
      {
	mt.pnums[i] = el[i];
	mt.pgeominfo[i] = el.GeomInfoPi (i+1);
      }

    mt.marked = 0;
    mt.surfid = el.GetIndex();

    mt.incorder = 0;
    mt.order = 1;

    int val = 0;
    for (i = 0; i < 2; i++)
      for (j = i+1; j < 3; j++)
	{
	  INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
	  i2.Sort();
	  int hval = edgenumber.Get(i2);
	  if (hval > val)
	    {
	      val = hval;
	      mt.markededge = 3 - i - j;
	    }
	}
  }
  

  
  void PrettyPrint(ostream & ost, const MarkedTri & mt)
  {
    ost << "MarkedTrig: " << endl;
    ost << "  pnums = "; for (int i=0; i<3; i++) ost << mt.pnums[i] << " "; ost << endl; 
    ost << "  marked = " << mt.marked << ", markededge=" << mt.markededge << endl;
    for(int i=0; i<2; i++)
      for(int j=i+1; j<3; j++)
	if(mt.markededge == 3-i-j)
	  ost << "  marked edge pnums = " << mt.pnums[i] << " " << mt.pnums[j] << endl;
  }


  void PrettyPrint(ostream & ost, const MarkedQuad & mq)
  {
    ost << "MarkedQuad: " << endl;
    ost << "  pnums = "; for (int i=0; i<4; i++) ost << mq.pnums[i] << " "; ost << endl; 
    ost << "  marked = " << mq.marked << ", markededge=" << mq.markededge << endl;
  }





  void BTDefineMarkedQuad (const Element2d & el,
			   INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
			   MarkedQuad & mq)
  {
    int i;
    for (i = 0; i < 4; i++)
      mq.pnums[i] = el[i];
    Swap (mq.pnums[2], mq.pnums[3]);  

    mq.marked = 0;
    mq.markededge = 0;
    mq.surfid = el.GetIndex();
  }




  // mark elements due to local h
  int BTMarkTets (T_MTETS & mtets,
		  T_MPRISMS & mprisms,
		  const Mesh & mesh)
  {
    int marked = 0;

    int np = mesh.GetNP();
    Vector hv(np);
    for (int i = 0; i < np; i++)
      hv(i) = mesh.GetH (mesh.Point(i+1));

    double hfac = 1;
  
    for (int step = 1; step <= 2; step++)
      {
	for (int i = 1; i <= mtets.Size(); i++)
	  {
	    double h = 0;
	  
	    for (int j = 0; j < 3; j++)
	      for (int k = j+1; k < 4; k++)
		{
		  const Point<3> & p1 = mesh.Point (mtets.Get(i).pnums[j]);
		  const Point<3> & p2 = mesh.Point (mtets.Get(i).pnums[k]);
		  double hh = Dist2 (p1, p2);
		  if (hh > h) h = hh;
		}
	    h = sqrt (h);
	  
	    double hshould = 1e10;
	    for (int j = 0; j < 4; j++)
	      {
		double hi = hv (mtets.Get(i).pnums[j]-1);
		if (hi < hshould)
		  hshould = hi;
	      }
	  
	
	    if (step == 1)
	      {
		if (h / hshould > hfac)
		  hfac = h / hshould;
	      }
	    else
	      {
		if (h > hshould * hfac)
		  {
		    mtets.Elem(i).marked = 1;
		    marked = 1;
		  }
		else
		  mtets.Elem(i).marked = 0;
	      }
	  
	  }
	for (int i = 1; i <= mprisms.Size(); i++)
	  {
	    double h = 0;
	  
	    for (int j = 0; j < 2; j++)
	      for (int k = j+1; k < 3; k++)
		{
		  const Point<3> & p1 = mesh.Point (mprisms.Get(i).pnums[j]);
		  const Point<3> & p2 = mesh.Point (mprisms.Get(i).pnums[k]);
		  double hh = Dist2 (p1, p2);
		  if (hh > h) h = hh;
		}
	    h = sqrt (h);
	  
	    double hshould = 1e10;
	    for (int j = 0; j < 6; j++)
	      {
		double hi = hv (mprisms.Get(i).pnums[j]-1);
		if (hi < hshould)
		  hshould = hi;
	      }
	  
	
	    if (step == 1)
	      {
		if (h / hshould > hfac)
		  hfac = h / hshould;
	      }
	    else
	      {
		if (h > hshould * hfac)
		  {
		    mprisms.Elem(i).marked = 1;
		    marked = 1;
		  }
		else
		  mprisms.Elem(i).marked = 0;
	      }
	  
	  }



	if (step == 1)
	  {
	    if (hfac > 2)
	      hfac /= 2;
	    else
	      hfac = 1;
	  }

      }
    return marked;
  }














  void BTBisectTet (const MarkedTet & oldtet, int newp, 
		    MarkedTet & newtet1, MarkedTet & newtet2)
  {
#ifdef DEBUG
    *testout << "bisect tet " << oldtet << endl;
#endif    
    
    int i, j, k;
  
  
    // points vis a vis from tet-edge
    int vis1, vis2;
    vis1 = 0;
    while (vis1 == oldtet.tetedge1 || vis1 == oldtet.tetedge2)
      vis1++;
    vis2 = 6 - vis1 - oldtet.tetedge1 - oldtet.tetedge2;


    


    // is tet of type P ?
    int istypep = 0;
    for (i = 0; i < 4; i++)
      {
	int cnt = 0;
	for (j = 0; j < 4; j++)
	  if (oldtet.faceedges[j] == i)
	    cnt++;
	if (cnt == 3)
	  istypep = 1;
      }


  
    for (i = 0; i < 4; i++)
      {
	newtet1.pnums[i] = oldtet.pnums[i];
	newtet2.pnums[i] = oldtet.pnums[i];
      }
    newtet1.flagged = istypep && !oldtet.flagged;
    newtet2.flagged = istypep && !oldtet.flagged;

    int nm = oldtet.marked - 1;
    if (nm < 0) nm = 0;
    newtet1.marked = nm;
    newtet2.marked = nm;

#ifdef DEBUG
    *testout << "newtet1,before = " << newtet1 << endl;
    *testout << "newtet2,before = " << newtet2 << endl;
#endif

    for (i = 0; i < 4; i++)
      {
	if (i == oldtet.tetedge1)
	  {
	    newtet2.pnums[i] = newp;
	    newtet2.faceedges[i] = oldtet.faceedges[i];  // inherited face
	    newtet2.faceedges[vis1] = i;        // cut faces
	    newtet2.faceedges[vis2] = i;

	    j = 0;
	    while (j == i || j == oldtet.faceedges[i])
	      j++;
	    k = 6 - i - oldtet.faceedges[i] - j;
	    newtet2.tetedge1 = j;                        // tet-edge
	    newtet2.tetedge2 = k;         

	    // new face:
	    if (istypep && oldtet.flagged)
              {
                int hi = 6 - oldtet.tetedge1 - j - k;
                newtet2.faceedges[oldtet.tetedge2] = char(hi);
              }
	    else
	      newtet2.faceedges[oldtet.tetedge2] = oldtet.tetedge1;
            
#ifdef DEBUG
            *testout << "i = " << i << ", j = " << j << " k = " << k 
                     << " oldtet.tetedge1 = " << oldtet.tetedge1 
                     << " oldtet.tetedge2 = " << oldtet.tetedge2
                     << "   6-oldtet.tetedge1-j-k = " <<  6 - oldtet.tetedge1 - j - k 
                     << "   6-oldtet.tetedge1-j-k = " <<  short(6 - oldtet.tetedge1 - j - k)
                     << endl;
            *testout << "vis1 = " << vis1 << ", vis2 = " << vis2 << endl;
            for (int j = 0; j < 4; j++)
              if (newtet2.faceedges[j] > 3)
                {
                  *testout << "ERROR1" << endl;
                }
#endif
	  }

	if (i == oldtet.tetedge2)
	  {
	    newtet1.pnums[i] = newp;
	    newtet1.faceedges[i] = oldtet.faceedges[i];  // inherited face
	    newtet1.faceedges[vis1] = i;
	    newtet1.faceedges[vis2] = i;
	    j = 0;
	    while (j == i || j == oldtet.faceedges[i])
	      j++;
	    k = 6 - i - oldtet.faceedges[i] - j;
	    newtet1.tetedge1 = j;        
	    newtet1.tetedge2 = k;

	    // new face:
	    if (istypep && oldtet.flagged)
              {
                int hi = 6 - oldtet.tetedge2 - j - k;
                newtet1.faceedges[oldtet.tetedge1] = char(hi);
              }
	    else
	      newtet1.faceedges[oldtet.tetedge1] = oldtet.tetedge2;

#ifdef DEBUG
            for (int j = 0; j < 4; j++)
              if (newtet2.faceedges[j] > 3)
                {
                  *testout << "ERROR2" << endl;
                }
#endif
	  }
      }

    newtet1.matindex = oldtet.matindex;
    newtet2.matindex = oldtet.matindex;
    newtet1.incorder = 0;
    newtet1.order = oldtet.order;
    newtet2.incorder = 0;
    newtet2.order = oldtet.order;

    *testout << "newtet1 =  " << newtet1 << endl;
    *testout << "newtet2 =  " << newtet2 << endl;
  }


  

  void BTBisectPrism (const MarkedPrism & oldprism, int newp1, int newp2,
		      MarkedPrism & newprism1, MarkedPrism & newprism2)
  {
    int i;

    for (i = 0; i < 6; i++)
      {
	newprism1.pnums[i] = oldprism.pnums[i];
	newprism2.pnums[i] = oldprism.pnums[i];
      }  
    
    int pe1, pe2;
    pe1 = 0;
    if (pe1 == oldprism.markededge)
      pe1++;
    pe2 = 3 - oldprism.markededge - pe1;

    newprism1.pnums[pe2] = newp1;
    newprism1.pnums[pe2+3] = newp2;
    newprism1.markededge = pe2;
    newprism2.pnums[pe1] = newp1;
    newprism2.pnums[pe1+3] = newp2;
    newprism2.markededge = pe1;

    newprism1.matindex = oldprism.matindex;
    newprism2.matindex = oldprism.matindex;

    int nm = oldprism.marked - 1;
    if (nm < 0) nm = 0;
    newprism1.marked = nm;
    newprism2.marked = nm;

    newprism1.incorder = 0;
    newprism1.order = oldprism.order;
    newprism2.incorder = 0;
    newprism2.order = oldprism.order;
  }


  void BTBisectIdentification (const MarkedIdentification & oldid,
			       Array<int> & newp,
			       MarkedIdentification & newid1,
			       MarkedIdentification & newid2)
  {
    for(int i=0; i<2*oldid.np; i++)
      {
	newid1.pnums[i] = oldid.pnums[i];
	newid2.pnums[i] = oldid.pnums[i];
      }
    newid1.np = newid2.np = oldid.np;

    if(oldid.np == 3)
      {
	newid1.pnums[(oldid.markededge+1)%3] = newp[0];
	newid1.pnums[(oldid.markededge+1)%3+3] = newp[1];
	newid1.markededge = (oldid.markededge+2)%3;

	newid2.pnums[oldid.markededge] = newp[0];
	newid2.pnums[oldid.markededge+3] = newp[1];
	newid2.markededge = (oldid.markededge+1)%3;
      }
    else if(oldid.np == 4)
      {
	newid1.pnums[(oldid.markededge+1)%4] = newp[0];
	newid1.pnums[(oldid.markededge+2)%4] = newp[2];
	newid1.pnums[(oldid.markededge+1)%4+4] = newp[1];
	newid1.pnums[(oldid.markededge+2)%4+4] = newp[3];
	newid1.markededge = (oldid.markededge+3)%4;

	newid2.pnums[oldid.markededge] = newp[0];
	newid2.pnums[(oldid.markededge+3)%4] = newp[2];
	newid2.pnums[oldid.markededge+4] = newp[1];
	newid2.pnums[(oldid.markededge+3)%4+4] = newp[3];
	newid2.markededge = (oldid.markededge+1)%4;
      }

    
    int nm = oldid.marked - 1;
    if (nm < 0) nm = 0;
    newid1.marked = newid2.marked = nm;

    newid1.incorder = newid2.incorder = 0;
    newid1.order = newid2.order = oldid.order;
  }



  void BTBisectTri (const MarkedTri & oldtri, int newp, const PointGeomInfo & newpgi,
		    MarkedTri & newtri1, MarkedTri & newtri2)
  {
    int i;

    for (i = 0; i < 3; i++)
      {
	newtri1.pnums[i] = oldtri.pnums[i];
	newtri1.pgeominfo[i] = oldtri.pgeominfo[i];
	newtri2.pnums[i] = oldtri.pnums[i];
	newtri2.pgeominfo[i] = oldtri.pgeominfo[i];
      }  

    int pe1, pe2;
    pe1 = 0;
    if (pe1 == oldtri.markededge)
      pe1++;
    pe2 = 3 - oldtri.markededge - pe1;

    newtri1.pnums[pe2] = newp;
    newtri1.pgeominfo[pe2] = newpgi;
    newtri1.markededge = pe2;

    newtri2.pnums[pe1] = newp;
    newtri2.pgeominfo[pe1] = newpgi;
    newtri2.markededge = pe1;


    newtri1.surfid = oldtri.surfid;
    newtri2.surfid = oldtri.surfid;

    int nm = oldtri.marked - 1;
    if (nm < 0) nm = 0;
    newtri1.marked = nm;
    newtri2.marked = nm;

    newtri1.incorder = 0;
    newtri1.order = oldtri.order;
    newtri2.incorder = 0;
    newtri2.order = oldtri.order;
    
    
  }


  void BTBisectQuad (const MarkedQuad & oldquad, 
		     int newp1, const PointGeomInfo & npgi1, 
		     int newp2, const PointGeomInfo & npgi2, 
		     MarkedQuad & newquad1, MarkedQuad & newquad2)
  {
    int i;

    for (i = 0; i < 4; i++)
      {
	newquad1.pnums[i] = oldquad.pnums[i];
	newquad1.pgeominfo[i] = oldquad.pgeominfo[i];
	newquad2.pnums[i] = oldquad.pnums[i];
	newquad2.pgeominfo[i] = oldquad.pgeominfo[i];
      }  

/*    if (oldquad.marked==1) // he/sz: 2d quads or 3d prism
    {   
      newquad1.pnums[1] = newp1;
      newquad1.pgeominfo[1] = npgi1;
      newquad1.pnums[3] = newp2;
      newquad1.pgeominfo[3] = npgi2;

      newquad2.pnums[0] = newp1;
      newquad2.pgeominfo[0] = npgi1;
      newquad2.pnums[2] = newp2;
      newquad2.pgeominfo[2] = npgi2;
    }
      
    else if (oldquad.marked==2) // he/sz: 2d quads only
    {
      newquad1.pnums[0] = newp1;
      newquad1.pnums[1] = newp2;
      newquad1.pnums[3] = oldquad.pnums[2];  
      newquad1.pnums[2] = oldquad.pnums[0]; 
      newquad1.pgeominfo[0] = npgi1;
      newquad1.pgeominfo[1] = npgi2;
      newquad1.pgeominfo[3] = oldquad.pgeominfo[2]; 
      newquad1.pgeominfo[2] = oldquad.pgeominfo[0];

      newquad2.pnums[0] = newp2;
      newquad2.pnums[1] = newp1;
      newquad2.pnums[3] = oldquad.pnums[1];  
      newquad2.pnums[2] = oldquad.pnums[3]; 
      newquad2.pgeominfo[0] = npgi2;
      newquad2.pgeominfo[1] = npgi1;
      newquad2.pgeominfo[3] = oldquad.pgeominfo[1]; 
      newquad2.pgeominfo[2] = oldquad.pgeominfo[3];
    }
      
    */
      
    if (oldquad.markededge==0 || oldquad.markededge==2)
    {
      newquad1.pnums[1] = newp1;
      newquad1.pgeominfo[1] = npgi1;
      newquad1.pnums[3] = newp2;
      newquad1.pgeominfo[3] = npgi2;

      newquad2.pnums[0] = newp1;
      newquad2.pgeominfo[0] = npgi1;
      newquad2.pnums[2] = newp2;
      newquad2.pgeominfo[2] = npgi2;
    }
    else // 1 || 3 
    {
      newquad1.pnums[2] = newp1;
      newquad1.pgeominfo[2] = npgi1;
      newquad1.pnums[3] = newp2;
      newquad1.pgeominfo[3] = npgi2;

      newquad2.pnums[0] = newp1;
      newquad2.pgeominfo[0] = npgi1;
      newquad2.pnums[1] = newp2;
      newquad2.pgeominfo[1] = npgi2;
    }
    newquad1.surfid = oldquad.surfid;
    newquad2.surfid = oldquad.surfid;

    int nm = oldquad.marked - 1;
    if (nm < 0) nm = 0;

    newquad1.marked = nm;
    newquad2.marked = nm;
    
    if (nm==1)
    {
      newquad1.markededge=1;
      newquad2.markededge=1;
    }
    else
    {
      newquad1.markededge=0;
      newquad2.markededge=0;
    }
    
  }


  int MarkHangingIdentifications(T_MIDS & mids, 
				 const INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i, j;
    
    int hanging = 0;
    for (i = 1; i <= mids.Size(); i++)
      {
	if (mids.Elem(i).marked)
	  {
	    hanging = 1;
	    continue;
	  }

	const int np = mids.Get(i).np;

	for(j = 0; j < np; j++)
	  {
	    INDEX_2 edge1(mids.Get(i).pnums[j],
			  mids.Get(i).pnums[(j+1) % np]);
	    INDEX_2 edge2(mids.Get(i).pnums[j+np],
			  mids.Get(i).pnums[((j+1) % np) + np]);

	    edge1.Sort();
	    edge2.Sort();
	    if (cutedges.Used (edge1) ||
		cutedges.Used (edge2))
	      {
		mids.Elem(i).marked = 1;
		hanging = 1;
	      }
	  }
      }

    return hanging;
  }


  /*
  void IdentifyCutEdges(Mesh & mesh,
			INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i,j,k;

    Array< Array<int,PointIndex::BASE>* > idmaps;
    for(i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++)
      {
	idmaps.Append(new Array<int,PointIndex::BASE>);
	mesh.GetIdentifications().GetMap(i,*idmaps.Last());
      }


    
    for(SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
      {
	const Element2d & el2d = mesh[sei];
	
	for(i = 0; i < el2d.GetNP(); i++)
	  {
	    INDEX_2 e1(el2d[i], el2d[(i+1) % el2d.GetNP()]);
	    e1.Sort();

	    if(!cutedges.Used(e1))
	      continue;

	    
	    for(k = 0; k < idmaps.Size(); k++)
	      {
		INDEX_2 e2((*idmaps[k])[e1.I1()],
			   (*idmaps[k])[e1.I2()]);
		
		if(e2.I1() == 0 || e2.I2() == 0 ||
		   e1.I1() == e2.I1() || e1.I2() == e2.I2())
		  continue;
		
		e2.Sort();

		if(cutedges.Used(e2))
		  continue;

		Point3d np = Center(mesh.Point(e2.I1()),
				    mesh.Point(e2.I2()));
		int newp = mesh.AddPoint(np);
		cutedges.Set(e2,newp);
		(*testout) << "DAAA" << endl;
	      }
	  }
      }

    
    for(i=0; i<idmaps.Size(); i++)
      delete idmaps[i];
    idmaps.DeleteAll();
  }
  */


  int MarkHangingTets (T_MTETS & mtets, 
		       const INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i, j, k;

    int hanging = 0;
    for (i = 1; i <= mtets.Size(); i++)
      {
	MarkedTet & teti = mtets.Elem(i);

	if (teti.marked)
	  {
	    hanging = 1;
	    continue;
	  }

	for (j = 0; j < 3; j++)
	  for (k = j+1; k < 4; k++)
	    {
	      INDEX_2 edge(teti.pnums[j],
			   teti.pnums[k]);
	      edge.Sort();
	      if (cutedges.Used (edge))
		{
		  teti.marked = 1;
		  hanging = 1;
		}
	    }
      }

    return hanging;
  }



  int MarkHangingPrisms (T_MPRISMS & mprisms, 
			 const INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i, j, k;

    int hanging = 0;
    for (i = 1; i <= mprisms.Size(); i++)
      {
	if (mprisms.Elem(i).marked)
	  {
	    hanging = 1;
	    continue;
	  }

	for (j = 0; j < 2; j++)
	  for (k = j+1; k < 3; k++)
	    {
	      INDEX_2 edge1(mprisms.Get(i).pnums[j],
			    mprisms.Get(i).pnums[k]);
	      INDEX_2 edge2(mprisms.Get(i).pnums[j+3],
			    mprisms.Get(i).pnums[k+3]);
	      edge1.Sort();
	      edge2.Sort();
	      if (cutedges.Used (edge1) ||
		  cutedges.Used (edge2))
		{
		  mprisms.Elem(i).marked = 1;
		  hanging = 1;
		}
	    }
      }
    return hanging;
  }



  int MarkHangingTris (T_MTRIS & mtris, 
		       const INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i, j, k;

    int hanging = 0;
    for (i = 1; i <= mtris.Size(); i++)
      {
	if (mtris.Get(i).marked)
	  {
	    hanging = 1;
	    continue;
	  }
	for (j = 0; j < 2; j++)
	  for (k = j+1; k < 3; k++)
	    {
	      INDEX_2 edge(mtris.Get(i).pnums[j],
			   mtris.Get(i).pnums[k]);
	      edge.Sort();
	      if (cutedges.Used (edge))
		{
		  mtris.Elem(i).marked = 1;
		  hanging = 1;
                }
	    }
      }  
      return hanging;
  }



  int MarkHangingQuads (T_MQUADS & mquads, 
			const INDEX_2_CLOSED_HASHTABLE<int> & cutedges)
  {
    int i;

    int hanging = 0;
    for (i = 1; i <= mquads.Size(); i++)
      {
	if (mquads.Elem(i).marked)
	  {
	    hanging = 1;
	    continue;
	  }

	INDEX_2 edge1(mquads.Get(i).pnums[0],
		      mquads.Get(i).pnums[1]);
	INDEX_2 edge2(mquads.Get(i).pnums[2],
		      mquads.Get(i).pnums[3]);
	edge1.Sort();
	edge2.Sort();
	if (cutedges.Used (edge1) ||
	    cutedges.Used (edge2))
	  {
	    mquads.Elem(i).marked = 1;
            mquads.Elem(i).markededge = 0;
	    hanging = 1;
            continue;
	  }
          
        // he/sz: second case: split horizontally
        INDEX_2 edge3(mquads.Get(i).pnums[1],
                      mquads.Get(i).pnums[3]);
        INDEX_2 edge4(mquads.Get(i).pnums[2],
                      mquads.Get(i).pnums[0]);

        edge3.Sort();
        edge4.Sort();
        if (cutedges.Used (edge3) ||
            cutedges.Used (edge4))
        {
          mquads.Elem(i).marked = 1;
          mquads.Elem(i).markededge = 1;
          hanging = 1; 
          continue; 
        }
    
      }
    return hanging;
  }



  void ConnectToNodeRec (int node, int tonode, 
			 const TABLE<int> & conto, Array<int> & connecttonode)
  {
    int i, n2;
    //  (*testout) << "connect " << node << " to " << tonode << endl;
    for (i = 1; i <= conto.EntrySize(node); i++)
      {
	n2 = conto.Get(node, i);
	if (!connecttonode.Get(n2))
	  {
	    connecttonode.Elem(n2) = tonode;
	    ConnectToNodeRec (n2, tonode, conto, connecttonode);
	  }
      }
  }




  T_MTETS mtets;
  T_MPRISMS mprisms;
  T_MIDS mids;
  T_MTRIS mtris;
  T_MQUADS mquads;


  void WriteMarkedElements(ostream & ost)
  {
    ost << "Marked Elements\n";

    ost << mtets.Size() << "\n";
    for(int i=0; i<mtets.Size(); i++)
      ost << mtets[i];

    ost << mprisms.Size() << "\n";
    for(int i=0; i<mprisms.Size(); i++)
      ost << mprisms[i];

    ost << mids.Size() << "\n";
    for(int i=0; i<mids.Size(); i++)
      ost << mids[i];

    ost << mtris.Size() << "\n";
    for(int i=0; i<mtris.Size(); i++)
      ost << mtris[i];

    ost << mquads.Size() << "\n";
    for(int i=0; i<mquads.Size(); i++)
      ost << mquads[i];
    ost << endl;
  }

  bool ReadMarkedElements(istream & ist, const Mesh & mesh)
  {
    string auxstring("");
    if(ist)
      ist >> auxstring;

    if(auxstring != "Marked")
      return false;

    if(ist)
      ist >> auxstring;

    if(auxstring != "Elements")
      return false;

    int size;

    ist >> size;
    mtets.SetSize(size);
    for(int i=0; i<size; i++)
      {
        ist >> mtets[i];
        if(mtets[i].pnums[0] > mesh.GetNV() || 
           mtets[i].pnums[1] > mesh.GetNV() || 
           mtets[i].pnums[2] > mesh.GetNV() || 
           mtets[i].pnums[3] > mesh.GetNV())
          return false;
      }

    ist >> size;
    mprisms.SetSize(size);
    for(int i=0; i<size; i++)
      ist >> mprisms[i];

    ist >> size;
    mids.SetSize(size);
    for(int i=0; i<size; i++)
      ist >> mids[i];

    ist >> size;
    mtris.SetSize(size);
    for(int i=0; i<size; i++)
      ist >> mtris[i];

    ist >> size;
    mquads.SetSize(size);
    for(int i=0; i<size; i++)
      ist >> mquads[i];

    return true;
  }





  void BisectTetsCopyMesh (Mesh & mesh, const class CSGeometry *,
			   BisectionOptions & opt,
			   const Array< Array<int,PointIndex::BASE>* > & idmaps,
			   const string & refinfofile)
  {
    mtets.SetName ("bisection, tets");
    mprisms.SetName ("bisection, prisms");
    mtris.SetName ("bisection, trigs");
    mquads.SetName ("bisection, quads");
    mids.SetName ("bisection, identifications");

    //int np = mesh.GetNP();
    int ne = mesh.GetNE();
    int nse = mesh.GetNSE();
    int i, j, k, l, m;

    /*
      if (mtets.Size() + mprisms.Size() == mesh.GetNE())
      return;
    */

    bool readok = false;

    if(refinfofile != "")
      {
	PrintMessage(3,"Reading marked-element information from \"",refinfofile,"\"");
	ifstream ist(refinfofile.c_str());

	readok = ReadMarkedElements(ist,mesh);

	ist.close();
      }

    if(!readok)
      {
	PrintMessage(3,"resetting marked-element information");
	mtets.SetSize(0);
	mprisms.SetSize(0);
	mids.SetSize(0);
	mtris.SetSize(0);
	mquads.SetSize(0);
	
	
	INDEX_2_HASHTABLE<int> shortedges(100);
	for (i = 1; i <= ne; i++)
	  {
	    const Element & el = mesh.VolumeElement(i);
	    if (el.GetType() == PRISM ||
		el.GetType() == PRISM12)
	      {
		for (j = 1; j <= 3; j++)
		  {
		    INDEX_2 se(el.PNum(j), el.PNum(j+3));
		    se.Sort();
		    shortedges.Set (se, 1);
		  }
	      }
	  }
	
	
	
	// INDEX_2_HASHTABLE<int> edgenumber(np);
	INDEX_2_CLOSED_HASHTABLE<int> edgenumber(9*ne+4*nse);  
	
	BTSortEdges (mesh, idmaps, edgenumber);
	
	
	for (i = 1; i <= ne; i++)
	  {
	    const Element & el = mesh.VolumeElement(i);
	    
	    switch (el.GetType())
	      {
	      case TET:
	      case TET10:
		{
		  // if tet has short edge, it is handled as degenerated prism
		  
		  int foundse = 0;
		  for (j = 1; j <= 3; j++)
		    for (k = j+1; k <= 4; k++)
		      {
			INDEX_2 se(el.PNum(j), el.PNum(k));
			se.Sort();
			if (shortedges.Used (se))
			  {
			    //		      cout << "tet converted to prism" << endl;
			    
			    foundse = 1;
			    int p3 = 1;
			    while (p3 == j || p3 == k)
			      p3++;
			    int p4 = 10 - j - k - p3;
			    
			    // even permutation ?
			    int pi[4];
			    pi[0] = j;
			    pi[1] = k;
			    pi[2] = p3;
			    pi[3] = p4;
			    int cnt = 0;
			    for (l = 1; l <= 4; l++)
			      for (m = 0; m < 3; m++)
				if (pi[m] > pi[m+1])
				  {
				    Swap (pi[m], pi[m+1]);
				    cnt++;
				  }
			    if (cnt % 2)
			      Swap (p3, p4);
			    
			    Element hel = el;
			    hel.PNum(1) = el.PNum(j);
			    hel.PNum(2) = el.PNum(k);
			    hel.PNum(3) = el.PNum(p3);
			    hel.PNum(4) = el.PNum(p4);
			    
			    MarkedPrism mp;
			    BTDefineMarkedPrism (hel, edgenumber, mp);
			    mp.matindex = el.GetIndex();
			    mprisms.Append (mp);
			  }
		      }
		  if (!foundse)
		    {
		      MarkedTet mt;
		      BTDefineMarkedTet (el, edgenumber, mt);
		      mt.matindex = el.GetIndex();
		      mtets.Append (mt);
		    }
		  break;
		}
	      case PYRAMID:
		{
		  // eventually rotate
		  MarkedPrism mp;
		  
		  INDEX_2 se(el.PNum(1), el.PNum(2));
		  se.Sort();
		  if (shortedges.Used (se))
		    {
		      Element hel = el;
		      hel.PNum(1) = el.PNum(2);
		      hel.PNum(2) = el.PNum(3);
		      hel.PNum(3) = el.PNum(4);
		      hel.PNum(4) = el.PNum(1);
		      BTDefineMarkedPrism (hel, edgenumber, mp);
		    }
		  else
		    {
		      BTDefineMarkedPrism (el, edgenumber, mp);
		    }
		  
		  mp.matindex = el.GetIndex();
		  mprisms.Append (mp);
		  break;
		}
	      case PRISM:
	      case PRISM12:
		{
		  MarkedPrism mp;
		  BTDefineMarkedPrism (el, edgenumber, mp);
		  mp.matindex = el.GetIndex();
		  mprisms.Append (mp);
		  break;
		}
              default:
                throw NgException("Bisect, element type not handled in switch, 4");
	      }
	  }
	
	for (i = 1; i <= nse; i++)
	  {
	    const Element2d & el = mesh.SurfaceElement(i);
	    if (el.GetType() == TRIG ||
		el.GetType() == TRIG6)
	      {
		MarkedTri mt;
		BTDefineMarkedTri (el, edgenumber, mt);
		mtris.Append (mt);
	      }
	    else
	      {
		MarkedQuad mq;
		BTDefineMarkedQuad (el, edgenumber, mq);
		mquads.Append (mq);
	      }
	    
	    MarkedIdentification mi;
	    for(j=0; j<idmaps.Size(); j++)
	      if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi))
		mids.Append(mi);
	  }
      }
	



    mesh.mlparentelement.SetSize(ne);
    for (i = 1; i <= ne; i++)
      mesh.mlparentelement.Elem(i) = 0;
    mesh.mlparentsurfaceelement.SetSize(nse);
    for (i = 1; i <= nse; i++)
      mesh.mlparentsurfaceelement.Elem(i) = 0;
  
    if (printmessage_importance>0)
    {
      ostringstream str1,str2;
      str1 << "copied " << mtets.Size() << " tets, " << mprisms.Size() << " prisms";
      str2 << "       " << mtris.Size() << " trigs, " << mquads.Size() << " quads";

      PrintMessage(4,str1.str());
      PrintMessage(4,str2.str());
    }
  }


  /*
  void UpdateEdgeMarks2(Mesh & mesh,
			const Array< Array<int,PointIndex::BASE>* > & idmaps)
  {
    Array< Array<MarkedTet>*,PointIndex::BASE > mtets_old(mesh.GetNP());
    Array< Array<MarkedPrism>*,PointIndex::BASE > mprisms_old(mesh.GetNP());
    Array< Array<MarkedIdentification>*,PointIndex::BASE > mids_old(mesh.GetNP());
    Array< Array<MarkedTri>*,PointIndex::BASE > mtris_old(mesh.GetNP());
    Array< Array<MarkedQuad>*,PointIndex::BASE > mquads_old(mesh.GetNP());

    for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
      mtets_old[i] = new Array<MarkedTet>;
    for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
      mprisms_old[i] = new Array<MarkedPrism>;
    for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
      mids_old[i] = new Array<MarkedIdentification>;
    for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
      mtris_old[i] = new Array<MarkedTri>;
    for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
      mquads_old[i] = new Array<MarkedQuad>;

    for(int i=0; i<mtets.Size(); i++)
      mtets_old[mtets[i].pnums[0]]->Append(mtets[i]);
    for(int i=0; i<mprisms.Size(); i++)
      mprisms_old[mprisms[i].pnums[0]]->Append(mprisms[i]);
    for(int i=0; i<mids.Size(); i++)
      mids_old[mids[i].pnums[0]]->Append(mids[i]);
    for(int i=0; i<mtris.Size(); i++)
      {
	(*testout) << "i " << i << endl;
	(*testout) << "mtris[i] " << mtris[i].pnums[0] << " " << mtris[i].pnums[1] << " " << mtris[i].pnums[2] << endl; 
	mtris_old[mtris[i].pnums[0]]->Append(mtris[i]);
      }
    for(int i=0; i<mquads.Size(); i++)
      mquads_old[mquads[i].pnums[0]]->Append(mquads[i]);

   
    
    int np = mesh.GetNP();
    int ne = mesh.GetNE();
    int nse = mesh.GetNSE();
    int i, j, k, l, m;


//       if (mtets.Size() + mprisms.Size() == mesh.GetNE())
//       return;

    

    mtets.SetSize(0);
    mprisms.SetSize(0);
    mids.SetSize(0);
    mtris.SetSize(0);
    mquads.SetSize(0);


    INDEX_2_HASHTABLE<int> shortedges(100);
    for (i = 1; i <= ne; i++)
      {
	const Element & el = mesh.VolumeElement(i);
	if (el.GetType() == PRISM ||
	    el.GetType() == PRISM12)
	  {
	    for (j = 1; j <= 3; j++)
	      {
		INDEX_2 se(el.PNum(j), el.PNum(j+3));
		se.Sort();
		shortedges.Set (se, 1);
	      }
	  }
      }



    // INDEX_2_HASHTABLE<int> edgenumber(np);
    INDEX_2_CLOSED_HASHTABLE<int> edgenumber(9*ne+4*nse);  

    BTSortEdges (mesh, idmaps, edgenumber);


    for (i = 1; i <= ne; i++)
      {
	const Element & el = mesh.VolumeElement(i);
	  
	switch (el.GetType())
	  {
	  case TET:
	  case TET10:
	    {
	      // if tet has short edge, it is handled as degenerated prism

	      int foundse = 0;
	      for (j = 1; j <= 3; j++)
		for (k = j+1; k <= 4; k++)
		  {
		    INDEX_2 se(el.PNum(j), el.PNum(k));
		    se.Sort();
		    if (shortedges.Used (se))
		      {
//		      cout << "tet converted to prism" << endl;

			foundse = 1;
			int p3 = 1;
			while (p3 == j || p3 == k)
			  p3++;
			int p4 = 10 - j - k - p3;

			// even permutation ?
			int pi[4];
			pi[0] = j;
			pi[1] = k;
			pi[2] = p3;
			pi[3] = p4;
			int cnt = 0;
			for (l = 1; l <= 4; l++)
			  for (m = 0; m < 3; m++)
			    if (pi[m] > pi[m+1])
			      {
				Swap (pi[m], pi[m+1]);
				cnt++;
			      }
			if (cnt % 2)
			  Swap (p3, p4);

			Element hel = el;
			hel.PNum(1) = el.PNum(j);
			hel.PNum(2) = el.PNum(k);
			hel.PNum(3) = el.PNum(p3);
			hel.PNum(4) = el.PNum(p4);

			MarkedPrism mp;

			BTDefineMarkedPrism (hel, edgenumber, mp);
			mp.matindex = el.GetIndex();
			mprisms.Append (mp);
		      }
		  }
	      if (!foundse)
		{
		  MarkedTet mt;
		  
		  int oldind = -1;
		  for(l = 0; oldind < 0 && l<mtets_old[el[0]]->Size(); l++)
		    if(el[1] == (*mtets_old[el[0]])[l].pnums[1] &&
		       el[2] == (*mtets_old[el[0]])[l].pnums[2] &&
		       el[3] == (*mtets_old[el[0]])[l].pnums[3])
		      oldind = l;

		  if(oldind >= 0)
		    mtets.Append((*mtets_old[el[0]])[oldind]);
		  else
		    {
		      BTDefineMarkedTet (el, edgenumber, mt);
		      mt.matindex = el.GetIndex();
		      mtets.Append (mt);
		    }
		}
	      break;
	    }
	  case PYRAMID:
	    {
	      // eventually rotate
	      MarkedPrism mp;
	    
	      INDEX_2 se(el.PNum(1), el.PNum(2));
	      se.Sort();
	      if (shortedges.Used (se))
		{
		  Element hel = el;
		  hel.PNum(1) = el.PNum(2);
		  hel.PNum(2) = el.PNum(3);
		  hel.PNum(3) = el.PNum(4);
		  hel.PNum(4) = el.PNum(1);
		  BTDefineMa…