/*
 $Id: triang2.cc,v 1.9 1996/11/07 10:56:04 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "triang2.h"
#if NOKASKPLOT == 0
#include "feplot2.h"
#endif

#include "utils.h"
#include "numerics.h"

#include "mzibutil.h"
#include "alloc.h"
#include "triangtempl.h"

#include "cmdpars.h"
extern CmdPars Cmd;

Vector<PT*>    EDG2::p(2);
Vector<PT*>    TR2::p(3);
Vector<EDG*>   TR2::e(3);
Vector<PATCH*> TR2::f(3);
//-------------------------------------------------------------------------


MESH2:: MESH2 (const char* inFileName, Bool readMesh)
		 
    : MESH(inFileName), ptList(0,25), edgList(0,25), trList(0,25),

	ptAlloc(ElementsInBlock), trAlloc(2*ElementsInBlock),  
	edgAlloc(3*ElementsInBlock), oldGreenEdges(25), oldGreenTriangles(25)
{

    noOfPoints    = 0;
    noOfEdges     = 0;
    noOfTriangles = 0;

    // if (Cmd.isTrue("Circle")) circle = True;
    // else 		      circle = False;

    if (readMesh)
    {
	if (Cmd.isTrue("ReadOldFormat")) readTriangulationOldFormat(inFileName);
	else readTriangulation(inFileName);
    }
}
//-------------------------------------------------------------------------

MESH2:: ~MESH2()
{
    int i;

    if (Cmd.isTrue("ReadOldFormat"))
    {
	delete []  trList[0]->first;
	delete []  ptList[0]->first;
	delete [] edgList[0]->first;

	//RemoveGreen();
	//FORALL( trList,i)  trList[i]->deleteAll();
	//FORALL( ptList,i)  ptList[i]->deleteAll();
	//FORALL(edgList,i) edgList[i]->deleteAll();
    }

    FORALL( trList,i) delete  trList[i];
    FORALL( ptList,i) delete  ptList[i];
    FORALL(edgList,i) delete edgList[i];
}
//-------------------------------------------------------------------------

FEPlot* MESH2:: newFEPlot(int plotType, char* caption, float size)
{
#if NOKASKPLOT == 0
    return new FEPlotMESH2(this, plotType, caption, size);
#else
    return 0;
#endif
}
//-------------------------------------------------------------------------


void MESH2:: readTriangulation(const char* fileName0)
{
    int i, maxIndex;
    TR2 *t;

    DList<PT2>*   ptL = new DList<PT2>;    ptList.push(ptL);
    DList<TR2>*   trL = new DList<TR2>;    trList.push(trL);
    DList<EDG2>* edgL = new DList<EDG2>;  edgList.push(edgL);

    BufferedParser parser(fileName0, "end", CommentFlag);


    if (parser.findWord("Poin") != parser.ok) missingItem("points",fileName);
    if (parser.findWord("max")  != parser.ok) missingItem("maxIndex",fileName);

    parser.getValue(&maxIndex);
    checkMaxIndex(maxIndex);

    Vector<PT2*> ptAdr(0,maxIndex);		// for point addresses
    FORALL(ptAdr,i) ptAdr[i] = 0;

    readPoints(parser, ptAdr);


    if (parser.findWord("Elem") != parser.ok) missingItem("triangles",fileName);
    readElements(parser, ptAdr);


    // --		compute edges:

    Vector<Stack<EDG2*>*> edsOfPts(noOfPoints);
    FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG2*>(10);

    DListIter<TR2> iter(triangles());
    while ((t=iter.all()))  CompleteTriangle(t, edsOfPts);


    if (parser.findWord("Boun") != parser.ok)  missingItem("boundary",fileName);
    readBoundary(parser, ptAdr, edsOfPts);


    FORALL(edsOfPts,i) delete edsOfPts[i];   

    MakeGreen();

    iter.reset();
    while ((t=iter.all())) t->setBoundary();

    if (infoRefinement) Info();
    if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------

void MESH2:: readPoints(BufferedParser& parser, Vector<PT2*>& ptAdr)
{
    int  ptNo, count=0;
    PT2* p;

    while (parser.getValue(&ptNo) == parser.ok)  
    {
	if (ptNo > ptAdr.high()) pointIndexError(ptNo);
	
	p = ptAlloc.Get();			// new PT2;
	ptAdr[ptNo] = p;
	parser.getValue(&p->x);
	parser.getValue(&p->y);
	p->node = ++count;			// used in Find...

	ptList[0]->add(p);
	++noOfPoints;
    }
}
//-------------------------------------------------------------------------

void MESH2:: readElements(BufferedParser& parser, Vector<PT2*>& ptAdr)
{
    int pNo1, pNo2, pNo3, classA;
    TR2 *t;

    while (parser.getValue(&pNo1) == parser.ok)
    {
	parser.getValue(&pNo2);
	parser.getValue(&pNo3);

	t = getTR2();			// new TR2;
	t->p1 = ptAdr[pNo1];
	t->p2 = ptAdr[pNo2];
	t->p3 = ptAdr[pNo3];

	parser.getValue(&classA);
	t->classA = classA;			// int -> char

	trList[0]->add(t);
	++noOfTriangles;
    }
}
//-------------------------------------------------------------------------


void MESH2:: readBoundary(BufferedParser& parser, Vector<PT2*>& ptAdr, 
			  Vector<Stack<EDG2*>*>& edsOfPts)
{
    int  i, pNo1, pNo2, pNo;
    PT2  *pt1, *pt2, *p;
    EDG2 *e=0;


    while (parser.getValue(&pNo1) == parser.ok)	   // edges on boundaries
    {
	parser.getValue(&pNo2);

	pt1 = ptAdr[pNo1];
	pt2 = ptAdr[pNo2];


	Stack<EDG2*>& edges = *edsOfPts[pt1->node];

	FORALL(edges,i) 
	{
	    e = edges[i];
	    if (pointOnEdge(pt2,e)) break;
	}

	if (!e->onBoundary()) boundaryWarning(pNo1,pNo2);
	e->readBC(parser);
    }


    if (parser.findWord("Poin") == parser.ok)    // additional BCs on points
    {
	while (parser.getValue(&pNo) == parser.ok)	
	{
	    p = ptAdr[pNo];
	    p->readBC(parser);
	}
    }
}
//-------------------------------------------------------------------------


void MESH2:: CompleteTriangle(TR2* t, Vector<Stack<EDG2*>*>& edsOfPts)
{
    setRightHand(t);

    PT2 *p1 = t->p1, *p2 = t->p2, *p3 = t->p3;

    t->e1 = FindEdgByPts(p2,p3,edsOfPts);
    t->e2 = FindEdgByPts(p3,p1,edsOfPts);
    t->e3 = FindEdgByPts(p1,p2,edsOfPts);
    SetTE(t);
}
//-------------------------------------------------------------------------

EDG2* MESH2:: FindEdgByPts(PT2* p1, PT2* p2, Vector<Stack<EDG2*>*>& edsOfPts)
{
    int   i;
    EDG2* ed;

    Stack<EDG2*>& edges = *edsOfPts[p1->node];
    FORALL(edges,i) 
    {
	ed = edges[i];
	if (pointOnEdge(p2,ed)) return ed;
    }
    
    ed = edgAlloc.Get();			// new EDG2;
    ed->p1 = p1;
    ed->p2 = p2;

    edsOfPts[p1->node]->push(ed);
    edsOfPts[p2->node]->push(ed);
    
    edgList[0]->add(ed);
    ++noOfEdges;
    return ed;
}
//-------------------------------------------------------------------------

void MESH2:: setRightHand(TR2* t)
{
    Jacobian J;
    t->compJ(J,False);
    if (J.detJ < 0.0)  { PT2* p=t->p2;  t->p2=t->p3;  t->p3=p; }
}
//-------------------------------------------------------------------------

Bool MESH2:: pointOnEdge(PT2* p, EDG2* e)
{
    return (e->p1 == p || e->p2 == p);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

#define GETLINE() line=tp;lineNo++;while(((*tp)!='\n')&&((*tp)!=NUL))tp++;if ((*tp)!=NUL){*tp=NUL;tp++;}


Bool MESH2:: readTriangulationOldFormat(const char* inFileName)
{
    int ptdim, edgdim, trdim, bounddim;
    int countP, countE, countT;

    double x,y;
    int rc, i, k1, k2, k3,
    lineNo = 0, classA, partner;
    char ch, *line, *tp;
    char globBuf[256];
    PT2  *p;
    EDG2 *ed;
    TR2  *t;


    char* inFileText = 0;
    ZIBReadFile(inFileName, &inFileText);
    
    tp = inFileText;
        
    GETLINE();
    
    GETLINE();
    rc = sscanf(line,"Dimension:(%d,%d,%d,%d)", &ptdim, &edgdim,
                &trdim, &bounddim);
    if (rc<=0)
    {
	cout << "\n*** ReadTriangulation: Error at top of file\n";
	return False;
    }


    // allocate lists for depth 0, reset and push them onto the stacks:

    DList<PT2>* ptL = new DList<PT2>; 
    ptList.push(ptL);
    PT2*  PArray = new PT2[ptdim];
    for (i=0; i<ptdim; ++i) 
    { 
	PArray[i].reset();  
	ptList[0]->add(&PArray[i]); 
    }

    DList<TR2>* trL  = new DList<TR2>; 
    trList.push(trL);
    TR2*  TArray = new TR2[trdim];
    for (i=0; i<trdim; ++i) 
    { 
	TArray[i].reset();  
	trList[0]->add(&TArray[i]); 
    }

    DList<EDG2>* edgL = new DList<EDG2>;
    edgList.push(edgL);
    EDG2*  EArray = new EDG2[edgdim];
    for (i=0; i<edgdim; ++i) 
    { 
	EArray[i].reset();  
	edgList[0]->add(&EArray[i]); 
    }

    if (PArray==0 || EArray==0 || TArray==0) 
        cout << "\n*** Read Triangulation: error in new for Arrays\n";
    
    countP = 0; countE = 0; countT = 0;

    while (True)
    {
        GETLINE();
	for (i=0; i<256; ++i) globBuf[i] = '\0';
        rc = sscanf(line, "%d:(%le,%le),%c%d", &i, &x, &y, &ch, &classA);
        if (rc==0) break;
	
        p = &PArray[i];
        p->x = x;
        p->y = y;
        p->node = countP + 1;
        p->boundP = ((ch=='D')||(ch=='d'))?DIRICHLET:
			(((ch=='I')||(ch=='i'))?INTERIOR:
	 		(((ch=='N')||(ch=='n'))?NEUMANN:CAUCHY));
        if (rc==5) p->classA = classA;

        ++countP;
    }
    if (strncmp(line,"END",3) != 0)
    { 
        printf("\n*** ReadTri: syntax error reading points (line %d)\n", lineNo);
        return False;
    }
    noOfPoints = countP;

    
			/* edges */

    while (True)
    {
        GETLINE();
	for (i=0; i<256; ++i) globBuf[i] = '\0';
        rc = sscanf(line,"%d:(%d,%d),%c%s", &i,&k1,&k2,&ch,globBuf);
        if (rc==0) break;
        if ((k1>=ptdim)||(k2>=ptdim))
	{
            printf("\n *** ReadTri: bound check error reading edge (line %d)\n",
		   lineNo);
            return False;
	}
        countE++;
        ed = &EArray[i];
        ed->p1 = &PArray[k1];
        ed->p2 = &PArray[k2];
        ed->boundP = ((ch=='D')||(ch=='d'))?DIRICHLET:
		(((ch=='I')||(ch=='i'))?INTERIOR:
		 (((ch=='N')||(ch=='n'))?NEUMANN:CAUCHY));

	if ((rc==5)&&((*globBuf)!=NUL))
	{
            rc = sscanf(globBuf,"%d,%c%d", &classA, &ch, &partner);
            if (rc==0)
	    { 
		rc = sscanf(globBuf,",%c%d", &ch, &partner);
		if (rc==0)
		    printf("\n Not a number after boundary type (line %d)\n",
			  lineNo);
	    }
            else ed->classA = classA;
            if (rc>1)
	    {
                if ((ch=='S')||(ch=='s'))
		ed->firstSon = &EArray[partner];
                else if ((ch=='F')||(ch=='f'))
		{
                    EArray[i+1].father = ed->father = &EArray[partner];
                    EArray[partner].pm =  &PArray[k2];
		}
                else
		{
                    printf("\n *** Unknown Character '%c' (line %d)\n",
			    ch, lineNo);
		}
	    }
	}
    }
    if (strncmp(line,"END",3)!=0)
    { 
        printf("\n*** ReadTri: syntax error reading edges (line %d)\n", lineNo);
        return False;
    }
    noOfEdges=countE;


	    		/* triangles */
    while (True)
    {
        GETLINE();
	for (i=0; i<256; ++i) globBuf[i] = '\0';

        rc = sscanf(line,"%d:(%d,%d,%d)%s",&i,&k1,&k2,&k3,globBuf);
        if (rc==0) break;
        if ((k1>=edgdim)||(k2>=edgdim)||(k3>=edgdim))
	{
	    printf(
	      "\n *** ReadTri: bound check error reading triangle (line %d)\n", 
	       lineNo);
            cout.flush(); abort();
	}
        countT++;
        t = &TArray[i];
        t->e1 = &EArray[k1];
        t->e2 = &EArray[k2];
        t->e3 = &EArray[k3];
        // t->next = &TArray[i+1];
        // t->prev = &TArray[i-1];

        // CallProcs(NewTriangle,(ptr)t,sizeof(TR2),0);

        SetTP(t, 0);
        SetTE(t);

        if ((rc==5)&&((*globBuf)!=NUL))
	{
            rc = sscanf(globBuf,",%d,%d", &classA, &partner);
            if (rc==0)
	    {
                printf("\n *** Not a number after triangle type (line %d)\n", 
			lineNo);
	    }
            else t->classA = classA;

	    /*
	      if (rc>1)
              {
	      if ((partner>=0)&&(partner<trdim))
	      t->partner = &TArray[partner];
	      else
	      {
	      sprintf(globBuf, "ReadTri: bound check error reading 
		triangle (line %d)\n", lineNo);
	      ZIBStdOut(globBuf);
	      CleanUp(inFile);
	      return False;
	      }
              }
	    */
	}
    }
    if (strncmp(line,"END",3)!=0)
    { 
        printf("\n *** ReadTri: syntax error reading triangles (line %d)\n", 
	       lineNo);
        return False;
    }
    noOfTriangles = countT;
    


    		/* boundary */
    /*
    if (bounddim!=0)
    {
        for (k = 0;k<bounddim;k++)
	{
            GETLINE();
            rc = sscanf(line,"(%le,%le):%d\n",&x,&y,&k1);
            midPoint = NewArcVec();
            if (midPoint==0) break;
            midPoint[0] = x;
            midPoint[1] = y;
            if (rc==0) break;
            for (i=0;i<k1;i++)
	    {
                GETLINE();
                rc = sscanf(line,"%d",&k2);
                AE(&EArray[k2],arcData) = (ptr)midPoint;
	    }
	}
        GETLINE();
        if (strncmp(line,"END",3)!=0)
	{ 
            printf("\n *** ReadTri: syntax error reading bounds (line %d)\n", 
		   lineNo);
            return False;
	}
    }
    */

    MakeGreen();

    DListIter<TR2> iter(triangles());
    while ((t=iter.all())) t->setBoundary();

    if (infoRefinement)  Info();

    delete inFileText;
    return True;
}
//-------------------------------------------------------------------------

void MESH2:: writeTriangulation(const char* /*outFileText*/, const Vector<Real>& u)
									   const
{    
    if (Cmd.isTrue("writeGrapeFormat")) writeGrapeFormat (fileName, u);
    if (Cmd.isTrue("writeZIBFormat"))   writeZIBFormat   (fileName, u);
}
//-------------------------------------------------------------------------


void MESH2:: writeZIBFormat(const char* outFileText, const Vector<Real>& u)
									const 
{
   int  i;
   int  node;
   TR2  *t;
   EDG2 *ed;
   PT2  *p;
   FILE *outFile;

   Vector<int> PNodes(noOfPoints);

   // saving the node numbers for later reuse
   i = 0;
   DListIter<PT2> PIter(points());
   while ((p=PIter.all()))  PNodes[++i] = p->node;


   char* outText = new char[strlen(outFileText)+10];
   strcpy(outText, outFileText);
   strcpy(strchr(outText,'.'),".zib.geo");

   outFile = fopen(outText, "w");
   if (outFile==0)
     {
       cout << "\n*** writeTriangulation: could not open file " << outText 
       << "\n";
       cout.flush(); abort();
     }

   fprintf(outFile, "#  %s\n\n", fileName); 
  
   fprintf(outFile, "Points: maxIndex %d\n", noOfPoints);

   node = 0;
   DListIter<PT2> iterPts(points());
   while ((p=iterPts.all()))  {
      p->node = node+1;
      fprintf(outFile, " %d %e %e\n", p->node, p->x, p->y);
      node++;
   }
   fprintf(outFile, "END\n");


   fprintf(outFile, "\nElements:\n");

   DListIter<TR2> iterTrs(triangles());
   while ((t=iterTrs.all()))  {
      fprintf(outFile, "%d %d %d %d\n",
                    (t->p1)->node, (t->p2)->node,
                    (t->p3)->node, t->classA );
   }
   fprintf(outFile, "END\n");


   fprintf(outFile, "\nBoundary: \n");

   DListIter<EDG2> iterEds(edges());
   while ((ed=iterEds.all()))  {
      if (!ed->onBoundary()) continue;
      if (ed->isDirichlet())
                  fprintf(outFile, "%d %d d %d\n",
                  (ed->p1)->node, (ed->p2)->node,ed->classA );
      else if (ed->isNeumann())
                  fprintf(outFile, "%d %d n %d\n",
                  (ed->p1)->node, (ed->p2)->node,ed->classA );
      else if (ed->isCauchy())
                  fprintf(outFile, "%d %d c %d\n",
                  (ed->p1)->node, (ed->p2)->node, ed->classA );
   }

   fprintf(outFile, "END\n");


   // output of values at the nodes located in the points 
   fprintf(outFile, "\nValues in nodes: \n");
   i = 0;
   iterPts.reset();
   while ((p=iterPts.all()))  fprintf(outFile, "%e\n", u[PNodes[++i]]);
   fprintf(outFile, "END\n");


   cout << "\n  WriteZIB-Data ok.; written to " << outText << "\n";
   delete outText;
   fclose(outFile);

   // restore of node numbers
   i = 0;
   PIter.reset();
   while ((p=PIter.all()))  p->node = PNodes[++i];
}
//-------------------------------------------------------------------------


void MESH2:: writeGrapeFormat(const char* outFileText, const Vector<Real>& u) 
									const 
{

   FILE *outFile;
   int  i;
   int  node;
   int  k[3];

   TR2 *t;
   PT2  *p;

   Vector<int> PNodes(noOfPoints), TNodes(noOfTriangles);

   i = 0;
   DListIter<PT2> PIter(points());
   while ((p=PIter.all()))  PNodes[++i] = p->node;

   i = 0;
   DListIter<TR2> TIter(triangles());
   while ((t=TIter.all()))  TNodes[++i] = t->node;

   char* outText = new char[strlen(outFileText)+10];
   strcpy(outText, outFileText);
   strcpy(strchr(outText,'.'),".grape");

   outFile = fopen(outText, "w");
   if (outFile==0)
     {
       cout << "\n*** grapeTriangulation: could not open file " << outText 
       << "\n";
       cout.flush(); abort();
     }

   fprintf(outFile, "%d %d\n",  noOfPoints, noOfTriangles);


   node = 0;
   DListIter<PT2> iterPts(points());
   while ((p=iterPts.all()))  {
      p->node = node;
      fprintf(outFile, "%e %e\n", p->x, p->y);
      node++;
   }

   node = 0;
   DListIter<TR2> iterT(triangles());
   while ((t=iterT.all()))  {
      t->node = node;

      fprintf(outFile, "%d %d %d \n",
                    (t->p1)->node, (t->p2)->node,
                    (t->p3)->node);

      node++;
   }

   iterT.reset();
   while ((t=iterT.all()))  
   {
       EDG2* ed;

       ed = t->e1;
       if ( (ed->t1 == nil) || (ed->t2 == nil) )
       k[0] = -1;
       else
       if ( (ed->t1 == t) )  
       k[0] = ed->t2->node;
       else k[0] = ed->t1->node;

       ed = t->e2;
       if ( (ed->t1 == nil) || (ed->t2 == nil) )
       k[1] = -1;
       else
       if ( (ed->t1 == t) )  k[1] = ed->t2->node;
         else k[1] = ed->t1->node;

       ed = t->e3;
       if ( (ed->t1 == nil) || (ed->t2 == nil) )
       k[2] = -1;
       else
       if ( (ed->t1 == t) )  k[2] = ed->t2->node;
       else k[2] = ed->t1->node;

       fprintf(outFile, "%d %d %d\n",
            k[0], k[1], k[2]);
   }

   // fprintf(outFile, "1\n");

   i = 0;
   iterPts.reset();
   while ((p=iterPts.all()))  fprintf(outFile, "%e\n", u[PNodes[++i]]);

   cout << "\n  GRAPE-Data ok.; written to " << outText << "\n";

   delete outText;
   fclose(outFile);

   i = 0;
   PIter.reset();
   while ((p=PIter.all()))  p->node = PNodes[++i];

   i = 0;
   TIter.reset();
   while ((t=TIter.all()))  t->node = TNodes[++i];
   
   return;
}
//-------------------------------------------------------------------------


void MESH2:: Refine()
{
    TR2* t;
    int k, noPt[3], noEdg[3], noTr[3];
    Timer timer, accTimer;

    noPt[0]  = noOfPoints;
    noTr[0]  = noOfTriangles;
    noEdg[0] = noOfEdges;
    
    RemoveGreen();				// delete all green triangless

    noPt[1]  = noOfPoints;
    noTr[1]  = noOfTriangles;
    noEdg[1] = noOfEdges;

    DList<TR2>*  trL = new DList<TR2>;    
    trList.push(trL);

    FORALL(oldGreenTriangles,k)  RefineTriangle(oldGreenTriangles[k]);

    DListIter<TR2> iter(triangles());
    while ((t=iter.all()))  RefineTriangle(t);

    noPt[2]  = noOfPoints;
    noTr[2]  = noOfTriangles;
    noEdg[2] = noOfEdges;

    MakeGreen();

    ++refStep;

    if (infoRefinement) refinementInfo (noPt, noEdg, noTr);
    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------


void MESH2:: Flatten()
{
    int noReleasedEdges=0, noReleasedTriangles=0;
    int noPoints=0, noEdges=0, noTriangles=0;
    int k, greenTriangles, greenEdges;
    const int impossibleDepth=255;

    TR2* t;
    EDG2* ed;
    PT2* p;
    DList<TR2>*  trL = new DList<TR2>;    
    DList<EDG2>* edgL = new DList<EDG2>;
    DList<PT2>*  ptL = new DList<PT2>;

    Stack<TR2*>  trToRetain(25), trToRelease(25);    
    Stack<EDG2*> edgToRetain(25), edgToRelease(25);
    Stack<PT2*>  ptToRetain(25);

    DListIter<TR2>  iterTr(triangles());
    DListIter<EDG2> iterEdg(edges());
    DListIter<PT2>  iterPt(points());

    greenEdges     = noOfEdges;
    greenTriangles = noOfTriangles;

    RemoveGreen();				// delete all green triangless

    greenEdges -= noOfEdges;
    greenTriangles -= noOfTriangles;

    while ((ed=iterEdg.allOfHistory())) ed->mark=False;
    iterEdg.reset();

    for (k=oldGreenEdges.l; k<=oldGreenEdges.h; k++)
    { 
	ed = oldGreenEdges[k];
	ed->mark=False; 
    }

    while (!oldGreenTriangles.empty())
    {
	t = oldGreenTriangles.pop();
	if (t->refined())
	{
	    trToRelease.push(t);
	    noReleasedTriangles++;
	}
	else
	{
	    trToRetain.push(t);
	    noTriangles++;
	    (t->e1)->mark = True;
	    (t->e2)->mark = True;
	    (t->e3)->mark = True;
	}
    }

    while ((t=iterTr.allOfHistory()))
    {
	if (t->refined())
	{
	    trToRelease.push(t);
	    noReleasedTriangles++;
	}
	else
	{
	    trToRetain.push(t);
	    noTriangles++;
	    (t->e1)->mark = True;
	    (t->e2)->mark = True;
	    (t->e3)->mark = True;
	}
    }

    for (k=oldGreenEdges.l; k<=oldGreenEdges.h; k++)
    {
	ed = oldGreenEdges[k];
	if (ed->mark)
	{
	    edgToRetain.push(ed);
	    noEdges++;
	}
	else
	{
	    edgToRelease.push(ed);
	    noReleasedEdges++;
	}
    }
    while ((ed=iterEdg.allOfHistory()))
    {
	if (ed->mark)
	{
	    edgToRetain.push(ed);
	    noEdges++;
	}
	else
	{
	    edgToRelease.push(ed);
	    noReleasedEdges++;
	}
    }

    while (!oldGreenEdges.empty()) oldGreenEdges.pop();
    
    while ((p=iterPt.allOfHistory()))
    {
	ptToRetain.push(p);
	noPoints++;
    }

    while (!trList.empty())  { delete trList.pop(); }
    trList.push(trL);

    while (!edgList.empty()) { delete edgList.pop(); }
    edgList.push(edgL);

    while (!ptList.empty())  { delete ptList.pop(); }
    ptList.push(ptL);

    for (k=ptToRetain.l; k<=ptToRetain.h; k++)
    {
	p = ptToRetain[k];
	ptL->add(p);
	p->depth =0;
    }

    for (k=edgToRetain.l; k<=edgToRetain.h; k++)
    {
	ed = edgToRetain[k];
	if (ed->firstSon)
	{
	    ed->depth = impossibleDepth;
	    oldGreenEdges.push(ed);
	    noEdges--;
	    ed->father = nil;
	    continue;
	}
	if ((ed->father) && (!((ed->father)->mark))) ed->father = nil;
	ed->depth = 0;
	edgL->add(ed);
    }

    for (k=trToRetain.l; k<=trToRetain.h; k++)
    {
	t = trToRetain[k];
	t->father = nil;

	if (((t->e1)->firstSon)||((t->e2)->firstSon)||((t->e3)->firstSon))
	{ 
	    t->depth = impossibleDepth; 
	    t->mark = False;
	    oldGreenTriangles.push(t);
	}
	else  
	{ 
	    t->depth = 0; 
	    trL->add(t);
	}
    }
    
    while (!trToRelease.empty())
    {
	t = trToRelease.pop();
	DelTE(t,3);
	returnTR2(t);
    }
    while (!edgToRelease.empty())
    {
	ed = edgToRelease.pop();
	edgAlloc.Return(ed);
    }

    refStep = 0;
    maxDepth = 0;
    ptIterLevel = 0;
    iterLevel = 0;
    noOfPoints = noPoints;
    noOfEdges  = noEdges;
    noOfTriangles = noTriangles;

    greenEdges = noOfEdges;
    greenTriangles = noOfTriangles;

    MakeGreen();

    greenEdges -= noOfEdges;
    greenTriangles -= noOfTriangles;

    if (infoRefinement)
    {
	cout << "\n Flatten: " << noReleasedEdges << " edges, " <<
	            noReleasedTriangles << " triangles released\n";
	cout <<   "          " << noPoints << " points, " <<
		        noEdges+greenEdges << " edges, " <<
	            noTriangles+greenTriangles << " triangles retained\n";
    }
    return;
}
//-------------------------------------------------------------------------

#define NEIGH_TR2(ed,t) (((ed->t1)==t)?(ed->t2):(ed->t1))


/* ***************************************************************
   RefineTriangle (DoRefTr) internal routine to refine a triangle, loTrues for
   neighbors do be refined too.
   --------
   TR2 *t       Adress of Triangle
   return int  True or False
   *************************************************************** */


int MESH2:: RefineTriangle(TR2 *t)
{
    TR2 *neighbor, *bigNeighbor;
    EDG2 *ed, *edFather;
    EDG2* eds[1+3];
    int k;
    

    if (t->mark) t->mark = False;         	// return if not marked 
    else  return True;

    if (t->firstSon != 0) return True;  	// return if already refined 
    
    if (RefineTriangleRed(t) != True) 		// refine it 
    {
	cout << "\n*** RefineTriangle: triangle not refined \n"; cout.flush(); abort();
    }

    /*
      Algorithm: For all three edges:
      (1)  look if neighbor triangle has 2 refined neighbors
      and refine it this case,
      (2)  if no neighbor triangle exists, "father neighbor"
      triangle needs to be refined (green refined at
      the last level),
      This routine uses recursion.
      */

    t->getEdges(eds);

    for (k=1; k<=3; k++)
    {
        ed=eds[k];
        if ((ed->boundP)==INTERIOR)
	{
            neighbor = NEIGH_TR2(ed,t);
            if (neighbor!=0)
	    {
                if (Mark2(neighbor)==True)                      // (1) 
		{
                    neighbor->mark=True;
                    if (!RefineTriangle(neighbor)) return False;
		}
	    }
            else                                                // (2) 
	    {
                edFather = ed->father;
                if (edFather==0) continue;
                if ((edFather->t1)==0) continue;

                bigNeighbor = (edFather->t1->firstSon == 0) ?
					edFather->t1 : edFather->t2;
                if (bigNeighbor == 0) continue;  // return True; ???
                bigNeighbor->mark=True;
                if (!RefineTriangle(bigNeighbor))
			{ printf("\nHuch\n"); return False; }
	    }
        }
    }
    return True;
}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   Mark2(t)     tests if triangle has two refined neighbors
   ----------
   TR2 *t        Adress of the triangle
   return int   True or False
   ****************************************************************** */


int MESH2:: Mark2(TR2 *t)
{
    int k, anz=0;
    EDG2* eds[1+3];
    
    t->getEdges(eds);

    if (t->mark) return True;
    for (k=1; k<=3; k++) if ((eds[k]->firstSon)!=0)  anz++;
    return (anz > 1) ? True : False;
}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   RefineTriangleRed (RefineTr)  does the actual refinement
   ----------
   TR2 *t        Adress of the triangle
   return int   True or False
   ****************************************************************** */


int MESH2:: RefineTriangleRed(TR2 *t)
{
    TR2  *t1, *t2, *t3, *t4;
    EDG2 *e1 = t->e1,  *e2 = t->e2,  *e3 = t->e3, *ereg,
    *ie1,*ie2,*ie3, *e11, *e12, *e21, *e22, *e31, *e32;

    // Get 4 new triangles, 3 new edges:

    t1 = getTR2();  
    t2 = getTR2();  
    t3 = getTR2(); 
    t4 = getTR2();
    noOfTriangles += 3;

    ie1 = edgAlloc.Get(); 
    ie2 = edgAlloc.Get(); 
    ie3 = edgAlloc.Get();
    noOfEdges += 3;

    // Set the "son" edges eik, i=1,2,3; k=1,2;
    //  if neeccessary refine the edges

    if ((e1->firstSon)==0)         /* Edge not refined */
    {
        e11 = RefineEdge(e1);
        if (e11==0) return False; /* No space for Refinement */
    }
    else e11 = e1->firstSon;         /* Take the existing "son" edge */
    e12 = e11->next;
    if ((e11->p1)==(t->p3))
    { ereg = e11; e11 = e12; e12 = ereg; }
    if ((e2->firstSon)==0)
    {
        e21 = RefineEdge(e2);
        if (e21==0) return False;
    }
    else e21 = e2->firstSon;
    e22 = e21->next;
    if ((e21->p1)==(t->p1))
    { ereg = e21; e21 = e22; e22 = ereg; }
    if ((e3->firstSon)==0)
    {
        e31 = RefineEdge(e3);
        if (e31==0) return False;
    }
    else e31 = e3->firstSon;
    e32 = e31->next;
    if ((e31->p1)==(t->p2))
    { ereg = e31; e31 = e32; e32 = ereg; }


	    /* Build the new triangles */

    InnerEdges(t,e1,e2,e3,ie1,ie2,ie3);
    t1->e1 = e32; t1->e2 = e11; t1->e3 = ie1;
    t2->e1 = e12; t2->e2 = e21; t2->e3 = ie2;
    t3->e1 = e22; t3->e2 = e31; t3->e3 = ie3;
    t4->e3 = ie1; t4->e2 = ie3; t4->e1 = ie2;
    SetTP(t1,t); SetTE(t1);
    SetTP(t2,t); SetTE(t2);
    SetTP(t3,t); SetTE(t3);
    SetTP(t4,t); SetTE(t4);
    t->firstSon = t1;
    t->type = T_RED;

    ie1->depth = ie2->depth = ie3->depth = t1->depth;

    trList[trList.h]->add(t1);
    trList[trList.h]->add(t2);
    trList[trList.h]->add(t3);
    trList[trList.h]->add(t4);

    
    if (ie1->depth > edgList.h)
    {
	DList<EDG2>* edgL = new DList<EDG2>;
	edgList.push(edgL);
    }
    edgList[ie1->depth]->add(ie1);
    edgList[ie1->depth]->add(ie2);
    edgList[ie1->depth]->add(ie3);

    return True;
}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   RefineEdge(ed)  does the actual refinement of an edge
   ----------
   EDG2 *ed        Adress of the edge
   return int     True or False
   ****************************************************************** */


EDG2* MESH2:: RefineEdge(EDG2 *ed)
{
    EDG2 *e1, *e2;
    PT2 *pm;

    e1 = edgAlloc.Get();
    e2 = edgAlloc.Get();
    ++noOfEdges;

    pm = ptAlloc.Get();
    ++noOfPoints;

    pm->classA = ed->classA; 
    pm->boundP = ed->boundP;
    pm->depth  = ed->depth + 1;


    if ( (ed->isoP == CIRCULAR) && (ed->boundP != INTERIOR) ) 
       setXYpm(pm, ed);
    else
    {
	pm->x = 0.5*(ed->p1->x + ed->p2->x);
	pm->y = 0.5*(ed->p1->y + ed->p2->y);
    }


    e1->p1 = ed->p1; e1->p2 = pm;
    e1->boundP = ed->boundP; e1->father = ed; e1->type = ed->type;
    e2->p1 = pm; e2->p2 = ed->p2;
    e2->boundP = ed->boundP; e2->father = ed; e2->type = ed->type;
    e1->classA = e2->classA = ed->classA;
    e1->isoP = e2->isoP = ed->isoP;
    e1->isBlueDiagonal = e2->isBlueDiagonal = ed->isBlueDiagonal;
    ed->firstSon = e1; ed->pm =pm;
    e1->depth = e2->depth = ed->depth+1;
    

    if (e1->depth > edgList.h)  // create new list at depth if necessary
    {
	DList<EDG2>* edgL = new DList<EDG2>;
	edgList.push(edgL);
    }
    edgList[e1->depth]->add(e1); 	
    edgList[e1->depth]->add(e2); 	

    if (pm->depth > ptList.h)  	
    {
	DList<PT2>* ptL = new DList<PT2>;
	ptList.push(ptL);
    }
    ptList[pm->depth]->add(pm);

    return e1;
}
/*-------------------------------------------------------------------------*/
void MESH2:: setXYpm(PT2* pm, const EDG2* e)
{
    Real x1, x2, y1, y2, phi1, phi2, phi, radius;
    Real xx, yy;
    int quadrant;
    Real realPi=3.1415926535897932384,  realPi2=1.57079632679489661923;

    // dieser Teil ist von einem frueheren Programm uebernommen.
    // fuer Kreisboegen geht das sicher einfacher.

    x1 = e->p1->x; y1 = e->p1->y;
    x2 = e->p2->x; y2 = e->p2->y;

    radius = sqrt(x1*x1 + y1*y1);

    if (x1>=0.0)
       quadrant = (y1>0)?0:3;
    else
      quadrant = (y1>=0)?1:2;
    xx = fabs((quadrant==0)||(quadrant==2)?x1:y1);
    yy = fabs((quadrant==0)||(quadrant==2)?y1:x1);
    phi1 = atan2(yy,xx)+quadrant*realPi2;

    if (x2>=0.0)
       quadrant = (y2>0)?0:3;
    else
      quadrant = (y2>=0)?1:2;
    xx = fabs((quadrant==0)||(quadrant==2)?x2:y2);
    yy = fabs((quadrant==0)||(quadrant==2)?y2:x2);
    phi2 = atan2(yy,xx)+quadrant*realPi2;
    

    phi = 0.5*(phi1 + phi2);

    if (fabs(phi1 - phi2) > realPi) phi = phi - realPi;
    pm->x = radius * cos(phi);
    pm->y = radius * sin(phi);

}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   InnerEdges(..)      sets the points of the inner edges
   ----------
   EDG2 *e1,*e2,*e3     Edges of the triangle to be refined
   EDG2 *ie1,*ie2,*ie3  Inner edges
   ****************************************************************** */


void MESH2:: InnerEdges (TR2 *t,EDG2 *e1,EDG2 *e2,EDG2 *e3,EDG2 *ie1,
			 EDG2 *ie2,EDG2 *ie3)
{
    ie1->p1 = e1->pm; ie1->p2 = e3->pm; 
    ie1->type = T_RED; ie1->classA = t->classA;
    ie2->p1 = e2->pm; ie2->p2 = e1->pm; 
    ie2->type = T_RED; ie2->classA = t->classA;
    ie3->p1 = e3->pm; ie3->p2 = e2->pm; 
    ie3->type = T_RED; ie3->classA = t->classA;
    ie1->isBlueDiagonal = e2->isBlueDiagonal;
    ie2->isBlueDiagonal = e3->isBlueDiagonal;
    ie3->isBlueDiagonal = e1->isBlueDiagonal;
    return;
}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   MakeTriangleGreen(t)   generates two green triangles
   ----------
   TR2 *t        Adress of the triangle
   return int   True or False
   ****************************************************************** */


int MESH2:: MakeGreen()
{
    TR2* t;
    int k;

    DListIter<TR2> iter(triangles());
    while ((t=iter.all()))  MakeTriangleGreen(t);

    FORALL(oldGreenTriangles,k) 
	 if (!(oldGreenTriangles[k]->firstSon))
	    MakeTriangleGreen(oldGreenTriangles[k]);

    return True;
}
//-------------------------------------------------------------------------


int MESH2:: MakeTriangleGreen(TR2 *t)
{
    PT2 *p1, *p2;
    EDG2 *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew;
    TR2 *t1New,*t2New;
    
    // Testing if one neighbor triangle is refined and setting pi, ei
    //   corresponding to this neighbor:
  
    if ((e1->firstSon)!=0)
    {
        p1 = t->p1; p2 = t->p2;
        e1 = t->e1; e2 = t->e2; e3 = t->e3;
    }
    else
    {
        if ((e2->firstSon)!=0)
	{
            p1 = t->p2; p2 = t->p3;
            e1 = t->e2; e2 = t->e3; e3 = t->e1;
	}
        else
	{
            if ((e3->firstSon)!=0)
	    {
                p1 = t->p3; p2 = t->p1;
                e1 = t->e3; e2 = t->e1; e3 = t->e2;
	    }
            else return True;                // nothing to refine 
	}
    }

    // e1->pm->mark = T_GREEN;	// "green point"


    // Get one new edge and two new triangles

    eNew = edgAlloc.Get();
    noOfEdges++;

    t1New = getTR2();
    t2New = getTR2();
    noOfTriangles += 1;

    t1New->next = t2New;
    t2New->prev = t1New;

    //  Set the fields for the new objects

    eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;
    eNew->type = T_GREEN; 
    eNew->classA = t->classA;
    eNew->depth  = t->depth+1;

    t1New->e1 = e3; t1New->e3 = eNew; t1New->type = T_GREEN;  // !!! T_WHITE;
    t2New->e1 = e2; t2New->e2 = eNew; t2New->type = T_GREEN;  // !!! T_WHITE;

    if ((e1->p1)==p2)                 // test for orientation 
    {
        t1New->e2 = e1->firstSon; t2New->e3 = (e1->firstSon)->next;
    }
    else
    {
        t1New->e2 = (e1->firstSon)->next; t2New->e3 = e1->firstSon;
    }
    DelTE(t,3);
    SetTP(t1New,t); SetTP(t2New,t);
    SetTE(t1New);   SetTE(t2New);
    t->firstSon = t1New;

    // if (t->type == T_GREEN) cout << "*** double green triangle! \n"; // !!!

    t->type = T_GREEN;


    //  Append two new ones to their lists

    trList[trList.h]->add(t1New);
    trList[trList.h]->add(t2New);

    
    if (eNew->depth > edgList.h)
    {
	DList<EDG2>* edgL = new DList<EDG2>;
	edgList.push(edgL);
    }
    edgList[eNew->depth]->add(eNew);

    return True;
}
//-------------------------------------------------------------------------

/* *******************************************************************
   RemoveGreen()   remove all green triangles
   ----------
   ****************************************************************** */


void MESH2:: RemoveGreen()
{
    TR2  *t1old,*t2old,*tFather;
    EDG2 *eMid = 0;
    int k;
    

    // Loop for all triangles of last refinement step (green triangles only here)

    TR2* t = trList[refStep]->first;
	
    while (t != 0)
    {              
	k = 0;
	tFather = t->father;
	if (tFather == 0) 	      { t = t->next; continue; }
	if (tFather->type != T_GREEN) { t = t->next; continue; }

	// Remove two green triangles, print internal error message
	// if the green edge is not found

	t1old = t; t2old = t->next;

	if (t->e1->type==T_GREEN) { eMid = t->e1; k++;}
	if (t->e2->type==T_GREEN) { eMid = t->e2; k++;}
	if (t->e3->type==T_GREEN) { eMid = t->e3; k++;}

	if (k != 1)
	  cout << "\n*** RemoveGreen: Triangle with impossible state k="
	       << k << "\n";

	if 	(t1old->mark) tFather->mark = True;     	// refine father
	else if (t2old->mark) tFather->mark = True;     	// refine father
	else tFather->mark = False; 

	t = t2old->next;
	    
	//  -- Remove the two green triangles; return the space       
	    
	trList[refStep]->remove(t1old);
	trList[refStep]->remove(t2old);

	DelTE(t1old,3);   
	DelTE(t2old,3);
	returnTR2(t1old);		
	returnTR2(t2old);		
	
	tFather->firstSon = 0;
	tFather->type = T_WHITE;
	SetTE(tFather);

	edgList[eMid->depth]->remove(eMid);
	edgAlloc.Return(eMid);

	--noOfTriangles;
	--noOfEdges;
    }

    return;
}
/*-------------------------------------------------------------------------*/

/* *******************************************************************
   SetTP(t,f)   set the points of a triangle, computed from the
   edges
   ----------
   TR2 *t        Adress of the triangle
   TR2 *f        Father of t
   ****************************************************************** */


void MESH2:: SetTP(TR2 *t, TR2 *f)
{
    PT2 *P11 = (t->e1)->p1, *P12 = (t->e1)->p2,
    	*P21 = (t->e2)->p1, *P22 = (t->e2)->p2,
    	*P31 = (t->e3)->p1, *P32 = (t->e3)->p2;

    t->p1 = ((P31==P11)||(P31==P12)) ? P32:P31;
    t->p2 = ((P11==P21)||(P11==P22)) ? P12:P11;
    t->p3 = ((P21==P31)||(P21==P32)) ? P22:P21;
    t->father = f;
    if (f != 0)
    { 
	t->depth = (f->depth)+1; 
	t->classA = f->classA; 
    }
    if (t->depth > maxDepth && t->type != T_GREEN)  maxDepth = t->depth;
    return;
}
//-------------------------------------------------------------------------

/* *******************************************************************
   SetTE(t)     set the triangles of an edge, computed from a new
   triangle, makes checks for programming errors
   ----------
   TR2 *t        Adress of the triangle
   ****************************************************************** */

void MESH2:: SetTE(TR2 *t)
{
    EDG2 *ed;
    int k = 0;
    
    ed = t->e1;
    if ((ed->t1)==0)      { ed->t1 = t; k++; }
    else if ((ed->t2)==0) { ed->t2 = t; k++; }
    ed = t->e2;
    if ((ed->t1)==0)      { ed->t1 = t; k++; }
    else if ((ed->t2)==0) { ed->t2 = t; k++; }
    ed = t->e3;
    if ((ed->t1)==0)      { ed->t1 = t; k++; }
    else if ((ed->t2)==0) { ed->t2 = t; k++; }

    if (k!=3) 
    {
	cout << "\n *** Tri: impossible state (SetTE), k = " << k;
	cout << "  (type: " << (int) t->type << ")\n";

        printf("\n adress t:  %p (%d,%d,%d)",t,
		       t->p1->node,t->p2->node,t->p3->node);
        printf("\n adress e1:  %p (%d-%d)",t->e1,
		       t->e1->p1->node,t->e1->p2->node);
        printf("\n adress e2:  %p (%d-%d)",t->e2,
		       t->e2->p1->node,t->e2->p2->node);
        printf("\n adress e3:  %p (%d-%d)",t->e3,
		       t->e3->p1->node,t->e3->p2->node);
        printf("\n t->e1->t1: %p",(t->e1)->t1);
        printf("\n t->e1->t2: %p",(t->e1)->t2);
        printf("\n t->e2->t1: %p",(t->e2)->t1);
        printf("\n t->e2->t2: %p",(t->e2)->t2);
        printf("\n t->e3->t1: %p",(t->e3)->t1);
        printf("\n t->e3->t2: %p",(t->e3)->t2);

#if NOKASKPLOT == 0
	FEPlotMESH2* fePlot = new FEPlotMESH2(this, SCREEN);
	fePlot->plotElements();
	fePlot->plotTriangleNodes();
	
	fePlot->plot->set(PENSIZE, BIG);
	fePlot->plotEdge(t->e1);
	fePlot->plotEdge(t->e2);
	fePlot->plotEdge(t->e3);
	fePlot->plot->flush();
	Pause();
#endif
    }

    if ( t->e1->onBoundary() || t->e2->onBoundary() || t->e3->onBoundary() )
         t->boundP = BOUNDARY;

    return;
}
/* *******************************************************************
   DelTE(t, check)  delete a triangle on its edges
   ----------
   TR2 *t        Adress of the triangle
   ****************************************************************** */

void MESH2:: DelTE(TR2 *t, int check)
{
    EDG2 *ed;
    int k = 0;
    /*
      if triangle on edge, remove it, if not generate internal
      error message
      */
    ed = t->e1;
    if ((ed->t1)==t)      { ed->t1 = 0; k++; }
    else if ((ed->t2)==t) { ed->t2 = 0; k++; }
    ed = t->e2;
    if ((ed->t1)==t)      { ed->t1 = 0; k++; }
    else if ((ed->t2)==t) { ed->t2 = 0; k++; }
    ed = t->e3;
    if ((ed->t1)==t)      { ed->t1 = 0; k++; }
    else if ((ed->t2)==t) { ed->t2 = 0; k++; }
    
    if (k!=check) 
    {
	cout << "\n *** Tri: impossible state (DelTE), k = " << k;
	cout << "  (type: " << (int) t->type << ")\n";

        printf("\n t (%d,%d,%d):  %p",t->p1->node,t->p2->node,t->p3->node,t);
        printf("\n e1 (%d-%d):  %p",t->e1->p1->node,t->e1->p2->node,t->e1);
        printf("\n e2 (%d-%d):  %p",t->e2->p1->node,t->e2->p2->node,t->e2);
        printf("\n e3 (%d-%d):  %p",t->e3->p1->node,t->e3->p2->node,t->e3);
        printf("\n t->e1->t1: %p",(t->e1)->t1);
        printf("\n t->e1->t2: %p",(t->e1)->t2);
        printf("\n t->e2->t1: %p",(t->e2)->t1);
        printf("\n t->e2->t2: %p",(t->e2)->t2);
        printf("\n t->e3->t1: %p",(t->e3)->t1);
        printf("\n t->e3->t2: %p",(t->e3)->t2);

#if NOKASKPLOT == 0
	FEPlotMESH2* fePlot = new FEPlotMESH2(this, SCREEN);
	fePlot->plotElements();
	fePlot->plotTriangleNodes();

	fePlot->plot->set(PENSIZE, BIG);
	fePlot->plotEdge(t->e1);
	fePlot->plotEdge(t->e2);
	fePlot->plotEdge(t->e3);
	fePlot->plot->flush();
	Pause();
#endif
    }
    return;
}
//-------------------------------------------------------------------------

ostream& operator << (ostream& os, MESH2& t)
{
    PT2*  p;
    EDG2* e;
    TR2*  tr;

    os << "\nTriangulation";
    if (t.fileName) os << " "<<t.fileName;
    os <<":\n";

    DListIter<PT2>  iterPT (t.points());
    DListIter<TR2>  iterTR (t.triangles());
    DListIter<EDG2> iterEDG(t.edges());

    while ((p =iterPT.all()))  p->print(cout);
    while ((e =iterEDG.all())) e->print(cout);
    while ((tr=iterTR.all()))  tr->print(cout);

    os << "\nRefinement step " << t.refStep;
    os << "  --   points: " << t.noOfPoints;
    os << "  edges: " << t.noOfEdges;
    os << "  triangles: " << t.noOfTriangles; 

    os << "\n";
    return os;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MESH2:: Info()
{
    cout << "\n Triangulation:    points     edges triangles "
    << "  Depth=" << maxDepth << "  Steps=" << refStep
    << Form("\n%15s%10d%10d%10d", " ", noOfPoints, noOfEdges, noOfTriangles)
    << "\n";
}
//-------------------------------------------------------------------------

int MESH2:: MemSpace()
{
    int space = varAllocator.MemSpace();

    space += edgAlloc.MemSpace();
    space += ptAlloc.MemSpace();
    space += trAlloc.MemSpace();

    return space;
}
//-------------------------------------------------------------------------


void MESH2:: refinementInfo(int* noPt, int* noEdg, int* noTr)
{
    int space  = MemSpace();
    int blocks = varAllocator.NoOfBlocks();

    blocks += edgAlloc.NoOfBlocks();
    blocks += ptAlloc.NoOfBlocks();
    blocks += trAlloc.NoOfBlocks();

    cout << "\n    Refine:        points     edges triangles "
    << "  Depth=" << maxDepth << "  Steps=" << refStep
    << Form("\n%15s%10d%10d%10d", "previous", noPt[0], noEdg[0], noTr[0]);

    if (infoRefinement > 1)
    {
	cout << Form("\n%15s%10d%10d%10d", "removed", 
		     noPt[1]-noPt[0], noEdg[1]-noEdg[0], noTr[1]-noTr[0])
	<< Form("\n%15s%10d%10d%10d", "generated", 
		noPt[2]-noPt[1], noEdg[2]-noEdg[1], noTr[2]-noTr[1])
	<< Form("\n%15s%10d%10d%10d", "added",
		noOfPoints-noPt[2], noOfEdges-noEdg[2], noOfTriangles-noTr[2])
	<<  "   for green closure";
    }
    cout << Form("\n%15s%10d%10d%10d   Space : %1.6f MB","total",
		 noOfPoints, noOfEdges, noOfTriangles, float(space)/1e6);

    cout << Form("\n%15s%10.2f%10.2f%10.2f   Blocks: %1d\n","ratio",
		 float(noOfPoints)/noPt[0], float(noOfEdges)/noEdg[0],
		 float(noOfTriangles)/noTr[0], blocks);
	

    if (infoRefinement > 2) 
    {
	int npDepth, neDepth, step, np=0, ne=0, nt=0;
	
	cout << "\n    Grid Statistics:"
	<< Form("\n%15s%10s%10s%10s%10s", "depth/step", "points", "edges", 
		  "triangles","(total)");

	FORALL (trList, step)
	{
	    npDepth = (step<= ptList.h) ?  ptList[step]->noOfElements : 0;
	    neDepth = (step<=edgList.h) ? edgList[step]->noOfElements : 0;

	    cout << Form("\n%15d%10d%10d%10d", step,
			 npDepth, neDepth, trList[step]->noOfElements);
	    np += npDepth;
	    ne += neDepth;
	    nt += trList[step]->noOfElements;
	}
	
	cout << Form("\n%15s%10d%10d%10d", "total sum", np, ne, nt);
	cout << "\n";
    }
}
//-------------------------------------------------------------------------


void MESH2:: consistencyCheck()
{
/*
  TR2*  t;
  EDG2* e;
  PT2*  p;
  
  DListIter<PT2>  iterPT (points());
  DListIter<TR2>  iterTR (triangles());
  DListIter<EDG2> iterEDG(edges());
  
  while (p=iterPT.all()) 
  {
  if (p->depth != depth) 
  cout << "\n*** Check Pt: Point with wrong depth\n";
  }
  
  while (e=iterEDG.all())
  {
  if (e->depth != depth) 
  cout << "\n*** Check Edg: Edge with wrong depth\n";
  }
  
  while (t=iterTR.all())
  {
  if (t->depth != depth) 
  cout << "\n*** Check Tri: Tri with wrong depth\n";
  
  if (t->type != T_GREEN)
  {
  if (t->e1->depth != depth) 
  {
  cout << "\n*** Check Tri: edge 1 with wrong depth " 
  << t->e1->depth << " " << depth << "\n";
  FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
  fePlot->plotElements();
  fePlot->plot->set(PENSIZE, BIG);
  fePlot->plotEdge(t->e1);
  fePlot->plot->flush();
  Pause();
  }
  if (t->e2->depth != depth) 
  {
  cout << "\n*** Check Tri: edge 2 with wrong depth " 
  << t->e2->depth << " " << depth << "\n";
  FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
  fePlot->plotElements();
  fePlot->plot->set(PENSIZE, BIG);
  fePlot->plotEdge(t->e2);
  fePlot->plot->flush();
  Pause();
  }
  if (t->e3->depth != depth) 
  {
  cout << "\n*** Check Tri: edge 3 with wrong depth " 
  << t->e3->depth << " " << depth << "\n";
  FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
  fePlot->plotElements();
  fePlot->plot->set(PENSIZE, BIG);
  fePlot->plotEdge(t->e3);
  fePlot->plot->flush();
  Pause();
  }
  }
  
	  int k = 0;
	  
	  ed = t->e1;
	  if ((ed->t1)==t)  ++k;
	  if ((ed->t2)==t)  ++k; 
	  ed = t->e2;
	  if ((ed->t1)==t)  ++k;
	  if ((ed->t2)==t)  ++k;
	  ed = t->e3;
	  if ((ed->t1)==t)  ++k;
	  if ((ed->t2)==t)  ++k;
	  
	  if (k != 3) 
	  {
	  cout << "\n *** Check Tri: impossible state, k = " << k << "\n";
	  cout << "type: " << t->type << "\n";
	  printf("\n adress t:  %p",t);
	  printf("\n t->e1->t1: %p",(t->e1)->t1);
	  printf("\n t->e1->t2: %p",(t->e1)->t2);
	  printf("\n t->e2->t1: %p",(t->e2)->t1);
	  printf("\n t->e2->t2: %p",(t->e2)->t2);
	  printf("\n t->e3->t1: %p",(t->e3)->t1);
	  printf("\n t->e3->t2: %p",(t->e3)->t2);
	  
	  FEPlotMESH2* fePlot = new FEPlotMESH2(*this, X11);
	  fePlot->plotElements();
	  
	  fePlot->plot->set(PENSIZE, BIG);
	  fePlot->plotEdge(t->e1);
	  fePlot->plotEdge(t->e2);
	  fePlot->plotEdge(t->e3);
	  fePlot->plot->flush();
	    Pause();
    }
    */
}
//-------------------------------------------------------------------------


void MESH2:: Resolve(Real reqhMin, int nSteps)
{
    Real hMin, h;
    int   depth;
    PT2*   pt;
    TR2*   tr;
    EDG2*  e;
    Timer timer;


    if (nSteps==0)			// -> reqhMin will be used
    {
	// --		get minimal h: --

	DListIter<EDG2> iterEDG(edges());

	hMin = machMax(Real(0.0));

	while ((e=iterEDG.all()))
	{
	    h = e->lengthSqr();
	    if (h < hMin) hMin = h;
	}
	hMin = sqrt(hMin);
	nSteps = Round(hMin/reqhMin);
    }
    if (nSteps > 10) 
    {
	cout << "\n*** MESH2::Resolve: nSteps = " << nSteps << " too large\n";
	cout.flush(); abort();
    }

			// resolve mesh:

    RemoveGreen();	

    for (depth=0; depth < nSteps; ++depth)
    {
	DListIter<TR2>  iter(triangles(), depth);
	while ((tr=iter.all())) tr->mark = 1;
	iter.reset();
	while ((tr=iter.all())) RefineTriangle(tr);
    }

/*
		   // delete lower levels:

    for (depth=0; depth < nSteps; ++depth)
    {
	edgList[depth]->deleteAll();
	trList [depth]->deleteAll();
    }
    edgList[0] = edgList[nSteps]; 
    trList [0] = trList [nSteps]; 

    edgList.h = 0;
    trList.h  = 0;

    for (TrListProc tl(trList[0], all,&tr); tr; tr=tl.All()) 
    {
	tr->father = 0;
	tr->depth  = 0;
    }
    for (EdgListProc el(edgList[0],all,&e); e; e=el.All())
    {
	e->father = 0;
	e->depth  = 0;
    }

    for (depth=0; depth < nSteps; ++depth) ptList[0]->concat(*ptList[depth+1]);
    ptList.h = 0;
     
*/

			// set node numbers:

    DListIter<PT2> iterPT(points());
    int node = 0;
    while ((pt=iterPT.all()))  { pt->node  = ++node; } // pt->depth = 0; }

/* 
    noOfPoints     = node;
    noOfEdges      = edgList[0]->noOfElements;
    noOfTriangles  = trList [0]->noOfElements;
  
    maxDepth = 0;
*/

    MakeGreen();

    if (Cmd.isTrue("timeResolve")) { cout << "\tResolve: "; timer.cpu(); }
    if (infoRefinement) Info();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void PT2:: print(ostream& os) const
{
    os << "Point " << node << "  BC " << (int) boundP << "  class "
    << (int) classA << "  mark " << (int) mark 
    << " x: " << x << " y: " << y << "\n";
}

void PT2:: reset()
{
    boundP = INTERIOR;
    mark = classA = depth = NUL;
    node   = 0;
}
void PT2:: readBC(BufferedParser& parser)
{
    readBCValues(&boundP, &classA, parser);
}
//-------------------------------------------------------------------------

void PT2:: getNodeCoord(Vector<Real>& v) const  { v[1]=x;  v[2]=y; }

//-------------------------------------------------------------------------

void PT2:: getCoordinates(Vector<Real>& v) const { v[1]=x;  v[2]=y; }

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void EDG2:: print(ostream& os) const
{
    os << "Edge with points " << p1->node << " " << p2->node
    << "  type " << (int) type << "  BC " << (int) boundP  
    << "  class " << (int) classA << "  depth " << (int) depth << "\n";
}

void EDG2:: reset()
{
    boundP = INTERIOR;
    mark = classA = type = depth = isoP = NUL;
    isBlueDiagonal = NUL;
    p1 = p2 = pm = 0;
    t1 = t2 = 0;
    node = 0; 
    father = firstSon = 0;
}
//-------------------------------------------------------------------------

void EDG2:: compNormalVector(Vector<Real>& n, Real* length) const
{
    Vector<Real> u(3),v(3),t(3);

    p1->getCoordinates(u);
    p2->getCoordinates(v);

    t[1] = v[1] - u[1];
    t[2] = v[2] - u[2];

    n[1] =  t[2];			// rotate 90 deg. clockwise
    n[2] = -t[1];

    *length = sqrt(n[1]*n[1]+n[2]*n[2]);
    n[1] /= (*length);
    n[2] /= (*length);
}
//-------------------------------------------------------------------------

// --  		point at unit coordinate -> point at x

void EDG2:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
    Real L1 = 1.0 - u[1];
    Real L2 = u[1];
    x[1] = L1*p1->x + L2*p2->x; 
    x[2] = L1*p1->y + L2*p2->y; 
}
//-------------------------------------------------------------------------

void EDG2:: getNodeCoord(Vector<Real>& v) const  
{ 
    v[1] = 0.5*(p1->x + p2->x); 
    v[2] = 0.5*(p1->y + p2->y); 
}
//-------------------------------------------------------------------------

Real EDG2:: hMax() const
{
    return sqrt(sqr(p2->x - p1->x) + sqr(p2->y - p1->y));
}
//-------------------------------------------------------------------------

void EDG2:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2) const
{
    p1->getCoordinates(cP1);
    p2->getCoordinates(cP2);
}
//-------------------------------------------------------------------------



void EDG2:: readBC(BufferedParser& parser)
{
    readBCValues(&boundP, &classA, &isoP, parser);

    if (boundP == DIRICHLET)
    {
	p1->boundP = p2->boundP = boundP;
	p1->classA = p2->classA = classA;
    }
    else if (boundP == CAUCHY || boundP == NEUMANN || boundP == INTERIOR)
    {
	if (p1->boundP != DIRICHLET)
	{
	    p1->boundP = boundP;
	    p1->classA = classA;
	}
	if (p2->boundP != DIRICHLET) 
	{
	    p2->boundP = boundP;
	    p2->classA = classA;
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void TR2:: compJ(Jacobian& j, Bool checkDetJ) const
{
    j.J(1,1) = p2->x - p1->x;		// Jacobian j.J(i,k) = dx(k)/dr(i)
    j.J(2,1) = p3->x - p1->x;
    j.J(1,2) = p2->y - p1->y;
    j.J(2,2) = p3->y - p1->y;

    j.detJ   = j.J(1,1)*j.J(2,2) - j.J(1,2)*j.J(2,1);

    if (checkDetJ)
      if (j.detJ <= 0.0)
	{
	  cout << "\n*** Jacobi-Det.<= 0 " << j.detJ << " encountered\n";
	  cout.flush(); abort();
	}

}

void TR2:: compJinv(Jacobian& j) const
{
    int i,k;

    compJ(j);
    j.Jinv(1,1) =  j.J(2,2);   j.Jinv(1,2) = -j.J(1,2); 
    j.Jinv(2,1) = -j.J(2,1);   j.Jinv(2,2) =  j.J(1,1);

    for (i=1; i<=2; ++i)  
    for (k=1; k<=2; ++k) j.Jinv(i,k) /= j.detJ;
}
//-------------------------------------------------------------------------

void TR2:: print(ostream& os) const
{
    TR::print(os);
    os << "  type " << (int) type << "  class " << (int)classA
    << "  mark "<< (int) mark << "  depth " << (int) depth << "\n";
}

void TR2:: reset()
{
    p1 = p2 = p3 = 0;
    e1 = e2 = e3 = 0;
    mark = type = classA = depth = NUL;
    boundP = INTERIOR;
    father = firstSon = 0;
    node = 0;
    error = 0.0;
}

void TR2:: setBoundary()
{
    if ( e1->onBoundary() || e2->onBoundary() || e3->onBoundary() )
         boundP = BOUNDARY;
}
//-------------------------------------------------------------------------

int TR2:: NoOfSons() const
{
    if      (type == MESH::T_RED)   return 4;
    else if (type == MESH::T_GREEN) return 2;
    else
    {
	cout << "\n*** TR2::NoOfSons(): unknown type " << type << "\n";
	cout.flush(); abort();
	return 0;
    }
}
//-------------------------------------------------------------------------
void TR2:: getEdges(EDG2* eds[]) const 
{ 
    eds[1] = e1;  eds[2] = e2;  eds[3] = e3; 
}

//-------------------------------------------------------------------------

void TR2:: getEdgeOrientation(Vector<Bool>& orient)  const
{
    int i;
    FORALL(orient,i) orient[i] = False;

    if (p1->node == 0 || p2->node == 0 || p3->node == 0) // !!!
    {
	cout <<"\n*** TR2:: getEdgeOrientation: node == 0 !\n"; 
	cout.flush(); abort();
    }

    if (p2->node < p3->node) orient[1] = True;
    if (p3->node < p1->node) orient[2] = True;
    if (p1->node < p2->node) orient[3] = True;
}
//-------------------------------------------------------------------------

void TR2:: getAllFaces(Vector<PATCH*>& faces) const
{
    faces[1] = e1; faces[2] = e2; faces[3] = e3;
}
//-------------------------------------------------------------------------

void TR2:: getNeighbours(Vector<PATCH*>& neighbours) const
{
    if (e1->t1 != this) neighbours[1] = e1->t1;
    else		neighbours[1] = e1->t2;

    if (e2->t1 != this) neighbours[2] = e2->t1;
    else		neighbours[2] = e2->t2;

    if (e3->t1 != this) neighbours[3] = e3->t1;
    else		neighbours[3] = e3->t2;
}
//-------------------------------------------------------------------------

// --  		point at unit coordinate -> point at x

void TR2:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
    Real L1 = 1.0 - u[1] - u[2];
    Real L2 = u[1];
    Real L3 = u[2];

    x[1] = L1*p1->x + L2*p2->x + L3*p3->x; 
    x[2] = L1*p1->y + L2*p2->y + L3*p3->y; 
}
//-------------------------------------------------------------------------

// --  		 point at x -> unit coordinates

void TR2:: unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const
{
    int i, k;
    Jacobian J;
    Vector<Real> xn(2);

    compJinv(J);

    xn[1] = x[1] - p1->x;
    xn[2] = x[2] - p1->y;

    for (i=1; i<=2; ++i) 		// jacobian is transposed!
    {
	u[i] = 0.0;
	for (k=1; k<=2; ++k)  u[i] += J.Jinv(k,i) * xn[k];
    }
}
//-------------------------------------------------------------------------

void TR2:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2, 
			  Vector<Real>& cP3) const
{
    p1->getCoordinates(cP1);
    p2->getCoordinates(cP2);
    p3->getCoordinates(cP3);
}
//-------------------------------------------------------------------------


Real TR2:: volume() const
{
    Matrix<Real> J(2,2);
    
    J(1,1) = p2->x - p1->x;		// Jacobian J(i,k) = dx(k)/dr(i)
    J(2,1) = p3->x - p1->x;
    J(1,2) = p2->y - p1->y;
    J(2,2) = p3->y - p1->y;

    return 0.5 * Abs(J(1,1)*J(2,2) - J(1,2)*J(2,1));
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MESH2:: resetElemIter(Bool lastStep)
{
    if (elemIter != 0)
    {
	cout << "\n*** resetElemIter: elemIter != 0"
	     << " \t (previous operation not completed?)\n";
	cout.flush(); abort();
    }

    if (lastStep) iterLevel = trList.high();
    else	  iterLevel = 0;
}
//-------------------------------------------------------------------------


PATCH* MESH2:: elemIterAll()
{
    if (elemIter != 0)				// iterate rest of list
    {
	for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())  
	    if (!elemIter->refined()) return elemIter;	

	if (++iterLevel > trList.high()) return 0;
    }

    while (1)					  // iterate new lists on stack
    {
	for (elemIter=trList[iterLevel]->first; elemIter; 
	     elemIter=elemIter->Next()) 
	  if (!elemIter->refined()) return elemIter;

	if (++iterLevel > trList.high()) return 0;
    }
}
//-------------------------------------------------------------------------

PATCH* MESH2:: elemIterAllOfHistory() 
{
    if (elemIter != 0)				// iterate rest of list
    {
	for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next()) 
	  return elemIter;

	if (++iterLevel > trList.high()) return 0;	
    }

    while (1)					  // iterate new lists on stack
    {
	for (elemIter=trList[iterLevel]->first; elemIter; 
	     elemIter=elemIter->Next())  return elemIter;

	if (++iterLevel > trList.high()) return 0;		
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


PATCH* MESH2:: edgIterAll()
{
    if (edgIter != 0)				// iterate rest of list
    {
	for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next())  
	    if (!edgIter->refined()) return edgIter;	

	if (++edgIterLevel > edgList.high()) return 0;
    }

    while (1)					  // iterate new lists on stack
    {
	for (edgIter=edgList[edgIterLevel]->first; edgIter; 
	     edgIter=edgIter->Next()) 
	  if (!edgIter->refined()) return edgIter;

	if (++edgIterLevel > edgList.high()) return 0;
    }
}
//-------------------------------------------------------------------------

PATCH* MESH2:: edgIterAllOfHistory() 
{
    if (edgIter != 0)				// iterate rest of list
    {
	for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next()) 
	  return edgIter;

	if (++edgIterLevel > edgList.high()) return 0;	
    }

    while (1)					  // iterate new lists on stack
    {
	for (edgIter=edgList[edgIterLevel]->first; edgIter; 
	     edgIter=edgIter->Next())  return edgIter;

	if (++edgIterLevel > edgList.high()) return 0;		
    }
}
//-------------------------------------------------------------------------


PATCH* MESH2:: ptIterAll()
{
    if (ptIter != 0)				// iterate rest of list
    {
	for (ptIter=ptIter->Next(); ptIter; ptIter=ptIter->Next())  
	    if (!ptIter->refined()) return ptIter;	

	if (++ptIterLevel > ptList.high()) return 0;
    }

    while (1)					  // iterate new lists on stack
    {
	for (ptIter=ptList[ptIterLevel]->first; ptIter; 
	     ptIter=ptIter->Next()) 
	  if (!ptIter->refined()) return ptIter;

	if (++ptIterLevel > ptList.high()) return 0;
    }
}

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MESH2BISECT:: MESH2BISECT (const char* inFileName, Bool /*readMesh*/)
		 
    : MESH2(inFileName), newTriangles(25)
{ }
//-------------------------------------------------------------------------


void MESH2BISECT:: Refine()
{
    TR2   *t, *tNext;
    int   noPt[3], noEdg[3], noTr[3];
    Timer timer, accTimer;

    while (!newTriangles.empty()) newTriangles.pop();

    noPt[0]  = noOfPoints;
    noTr[0]  = noOfTriangles;
    noEdg[0] = noOfEdges;

    noPt[1]  = noOfPoints;
    noTr[1]  = noOfTriangles;
    noEdg[1] = noOfEdges;

    DListIter<TR2> iter(triangles());
    while ((t=iter.all())) EdgBisection(t);

    int depth;
    FORALL(trList,depth)	
    {
    	for (t=trList[depth]->first; t; t=tNext)
    	{
	   tNext = t->next;			// a tri may be removed!
	   RefineTriangle(t);
    	}
    }

    int moreTrianglesToRefine;

    while (1)
    {
	moreTrianglesToRefine = 0;

	iter.reset();
	while ((t=iter.all())) moreTrianglesToRefine += TestBisection(t);
	if (moreTrianglesToRefine == 0) break;

	if (infoRefinement>2) 
	  cout << "\n moreTrianglesToRefine = " << moreTrianglesToRefine << "\n";
       
	FORALL(trList,depth)
	{
	    for (t=trList[depth]->first; t; t=tNext)
	    {
	   	tNext = t->next;		// a tri may be removed!
		RefineTriangle(t);
	    }
	}
    }

    noPt[2]  = noOfPoints;
    noTr[2]  = noOfTriangles;
    noEdg[2] = noOfEdges;

    ++refStep;

    if (infoRefinement) refinementInfo (noPt, noEdg, noTr);
    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------

int MESH2BISECT:: RefineTriangle(TR2 *t)
{
    if (t->mark) t->mark = False;         	// return if not marked 
    else  return True;
    
    if (t->firstSon) return True;  		// return if already refined 
    
    if (!RefineTriangleRed(t))	 		// refine it 
    {
	cout << "\n*** RefineTriangle: triangle not bisected \n"; cout.flush(); abort();
    }
    return True;
}
//-------------------------------------------------------------------------

int MESH2BISECT:: TestBisection(TR2 *t)
{
    EDG2* eds[1+3];
    int k;

    t->getEdges(eds);
    for (k=1; k<=3; k++) 
      if (eds[k]->firstSon != 0)  { t->mark=True; return 1; }

    return 0;
}
//-------------------------------------------------------------------------

void MESH2BISECT:: EdgBisection(TR2 *t)
{
    EDG2* eds[1+3];
    int k;

    if (t->mark)
    {
	t->getEdges(eds);
	for (k=1; k<=3; k++) if (eds[k]->firstSon == 0) RefineEdge(eds[k]);
    }
}
//-------------------------------------------------------------------------

/* *******************************************************************
   MESH2BISECT::RefineTriangleRed   does the actual bisection of a triangle
   ----------
   TR2 *t       Adress of the triangle
   return int   True or False
   ****************************************************************** */


int MESH2BISECT:: RefineTriangleRed(TR2 *t)
{
    PT2 *p1=0, *p2=0;
    EDG2 *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew, *e11;
    TR2 *t1New, *t2New;
    TR2* tFather = t->father;
    int  longestEdge;

    if (!t->isVolatile())		// compute longest edge
    {
	Real lengthEdg1 = e1->hMax();
	Real lengthEdg2 = e2->hMax();
	Real lengthEdg3 = e3->hMax();
	
	longestEdge = 1;
	Real maxLength = lengthEdg1;
	if (lengthEdg2 > maxLength) { longestEdge=2; maxLength=lengthEdg2; }
	if (lengthEdg3 > maxLength) { longestEdge=3; }
    }
    else  	// a volatile son: the edge with lower depth MUST be refined
    {
	if      (e1->depth < t->depth) longestEdge = 1;
	else if (e2->depth < t->depth) longestEdge = 2;
	else 			       longestEdge = 3;
    }

    if (longestEdge==1)
    {
        p1 = t->p1; p2 = t->p2;
        e1 = t->e1; e2 = t->e2; e3 = t->e3;
    }
    else if (longestEdge==2)
    {
	p1 = t->p2; p2 = t->p3;
	e1 = t->e2; e2 = t->e3; e3 = t->e1;
    }
    else
    {
	p1 = t->p3; p2 = t->p1;
	e1 = t->e3; e2 = t->e1; e3 = t->e2;
    }

    // Get one new edge and two new triangles

    eNew = edgAlloc.Get();
    ++noOfEdges;

    t1New = getTR2();
    t2New = getTR2();
    ++noOfTriangles;

    //  Set the fields for the new objects

    if (e1->firstSon == 0) e11 = RefineEdge(e1);
    else 		   e11 = e1->firstSon;

    eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;
    eNew->type = T_RED; 
    eNew->classA = t->classA;
    eNew->depth  = e11->depth;

    t1New->e1 = e3; t1New->e3 = eNew; 
    t2New->e1 = e2; t2New->e2 = eNew; 

    if ((e1->p1)==p2)                 // test for orientation 
    {
        t1New->e2 = e1->firstSon; 
	t2New->e3 = (e1->firstSon)->next;
    }
    else
    {
        t1New->e2 = (e1->firstSon)->next; 
	t2New->e3 = e1->firstSon;
    }

    DelTE(t,3);

    if (!t->isVolatile())     // a non-volatile son creates two volatile sons
    {
	t1New->type = t2New->type = T_VOLATILE;
    	SetTP(t1New,t); 
	SetTP(t2New,t);

    	t->firstSon = t1New;
    	t->type = T_RED1;

    	if (t1New->depth > trList.h) trList.push(new DList<TR2>);
    	trList[t1New->depth]->add(t1New);
    	trList[t2New->depth]->add(t2New);
    }
    else 	// a volatile son creates two non-volatile sons and is removed
    {	
	t1New->type = t2New->type = T_WHITE;

       	SetTP(t1New,tFather); 
	SetTP(t2New,tFather);

	if      (tFather->type == T_RED1) tFather->type = T_RED2;
	else if (tFather->type == T_RED2) tFather->type = T_RED3;
    	else { cout << "\n*** error in tFather->type\n"; cout.flush(); abort(); }

	if (tFather->firstSon == t) tFather->firstSon = t1New;	
	trList[t->depth]->substitute(t,t1New,t2New);

	returnTR2(t);
    }

    SetTE(t1New);   
    SetTE(t2New);

    newTriangles.push(t1New); 
    newTriangles.push(t2New);

    if (eNew->depth > edgList.h) edgList.push(new DList<EDG2>);
    edgList[eNew->depth]->add(eNew);

    return True;
}