/*
 $Id: triang3.cc,v 1.5 1996/10/11 16:13:22 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "triang3.h"
#include "utils.h"
#include "numerics.h"
#include "mzibutil.h"
#include "triangtempl.h"

#include "cmdpars.h"
extern CmdPars Cmd;

Vector<PT*>    EDG3::p(2);
Vector<PT*>    TR3::p(3);
Vector<EDG*>   TR3::e(3);
Vector<PATCH*> TR3::f(3);
Vector<PT*>    TET3::p(4);
Vector<EDG*>   TET3::e(6);
Vector<PATCH*> TET3::f(4);
//-------------------------------------------------------------------------

static const Real oneSixth = 1./6.;

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

StaticAllocator<TR3> TR3::alloc(ElementsInBlock);

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


MESH3:: MESH3(const char* inFileName, Bool readMesh)

    : MESH(inFileName), ptList(0,25), trList(0,25), edgList(0,25), tetList(0,25),

	ptAlloc(ElementsInBlock), edgAlloc(6*ElementsInBlock), 
	tetAlloc(5*ElementsInBlock) 
{

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

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

    // if (Cmd.isTrue("Cylinder")) cylinder = True;
    // else 		        cylinder = False;

    if (Cmd.isSet("refTetrahedron","BeyStabilDiag")) 
    {
	rightHand = False;
	refinementType = BeyStabilDiag;
    }
    else 
    {
	rightHand = True;			
	refinementType = ShortDiag;
    }

    if (readMesh) readTriangulation(inFileName);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


EFMESH3:: EFMESH3(const char* inFileName, Bool /*readMesh*/)

       : MESH3(inFileName, False), EFTetAlloc(5*ElementsInBlock) 
{
    readTriangulation(fileName);
}
//-------------------------------------------------------------------------

MESH3:: ~MESH3()
{
    int i;
    FORALL( trList,i) delete  trList[i];
    FORALL( ptList,i) delete  ptList[i];
    FORALL(edgList,i) delete edgList[i];
    FORALL(tetList,i) delete tetList[i];
}
//-------------------------------------------------------------------------


void MESH3:: readTriangulation(const char* fileName0)
{    
    if (Cmd.isTrue("readHyperFormat")) 	readHyperFormat(fileName0);
    else
    if (Cmd.isTrue("readGrapeFormat")) 	readGrapeFormat(fileName0);
    else
    if (Cmd.isTrue("readZIBFormat"))   	readZIBFormat(fileName0);
    else
    if (Cmd.isTrue("readOldFormat"))   	readOldZIBFormat(fileName0);
    else			 	readZIBFormat(fileName0);
}
//-------------------------------------------------------------------------


void MESH3:: readOldZIBFormat(const char* /*fileName0*/)
{
    int ptdim, tetdim, bounddim;

    int    rc, i, k;
    int    classA;
    int    dirTrDim, trNo;
    char   ch, line[256] ;

    PT3    *p;
    TR3    *t=0;
    TET3   *td;
    
    int    maxIndex;
    
    DList<PT3>* ptL = new DList<PT3>;  ptList.push(ptL);
    DList<EDG3>* edgL  = new DList<EDG3>;  edgList.push(edgL);
    DList<TR3>* trL  = new DList<TR3>;  trList.push(trL);
    DList<TET3>* tetL = new DList<TET3>; tetList.push(tetL);

//    BufferedParser parser(fileName, "end", CommentFlag);
Parser parser(fileName, "end", CommentFlag);
    if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
    
    if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
    rc = sscanf(line,"Dimension:(%d,%d,%d,%d)",
                      &ptdim,&tetdim,&dirTrDim,&bounddim);
    maxIndex = ptdim - 1;
    checkMaxIndex(maxIndex);
    noOfTetrahedra = tetdim;
   
    
    Vector<PT3*> ptAdr(0,maxIndex);		// for point addresses
    FORALL(ptAdr,i) ptAdr[i] = 0;
    
    int  ptNo, count=0;
    double x, y, z;
    
    for (i=1; i<= ptdim; i++)
    {
	if (parser.getLine(line) != parser.ok) missingItem("line",fileName);

        ch = ' ';
        rc = sscanf(line, "%d:(%le,%le,%le),%c%d", &ptNo, &x, &y, &z, &ch,
		    &classA);
        if (rc==0) break;
	
	if (ptNo > ptAdr.high()) pointIndexError(ptNo);
	
	p = ptAlloc.Get();			// new PT3;
	ptAdr[ptNo] = p;

	p->x = x;
	p->y = y;
	p->z = z;
	
        p->classA = 0;
        if (rc==6) p->classA = classA;

        if ((ch=='D')||(ch=='d')) p->boundP = DIRICHLET;
        else if ((ch=='I')||(ch=='i')) p->boundP = INTERIOR;
          else if ((ch=='N')||(ch=='n')) p->boundP = NEUMANN;
            else if ((ch=='C')||(ch=='c')) p->boundP = CAUCHY;
              else  p->boundP = INTERIOR;

	p->node = ++count;			// used in Find...

	ptList[0]->add(p);
    }
    
    noOfPoints=count;
    
    if (parser.findWord("end") != parser.ok) missingItem("end",fileName);
    parser.skipRestOfLine();
    
    int  pNo1, pNo2 ,pNo3, pNo4, tetNo;
    
    for (i=1; i<= tetdim; i++)
    {
	classA = 0;
	
	if (parser.getLine(line) != parser.ok) missingItem("line",fileName);

	rc = sscanf(line,"%d:(%d,%d,%d,%d),%d", &tetNo,&pNo1,&pNo2,&pNo3,&pNo4,
		    &classA);
	if (rc==0) break;
	if ((pNo1>maxIndex)||(pNo2>maxIndex)||(pNo3>maxIndex)||(pNo4>maxIndex))
	{
	    cout << "\n *** ReadTri: bound check error reading tetrahedron\n";
	    return;
	}
	
	td = getTET3();				// new TET3;
	td->p1 = ptAdr[pNo1];
	td->p2 = ptAdr[pNo2];
	td->p3 = ptAdr[pNo3];
	td->p4 = ptAdr[pNo4];
	
	td->classA = classA;
	SetTPts(td,0); 
	
	tetList[0]->add(td);
    }
    
    // --		compute edges and faces:
    
    Vector<Stack<EDG3*>*> edsOfPts(noOfPoints);
    FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG3*>(15);
    
    Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
    FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);
    
    DListIter<TET3> iterTet(tetras());
    while ((td=iterTet.all()))  completeTetrahedron(td, edsOfPts, trsOfPts);
    
    
    if (parser.findWord("end") != parser.ok) missingItem("end",fileName);
    parser.skipRestOfLine();
    
    
    PT3 *pt1, *pt2, *pt3;
    
    for (i=1; i<= dirTrDim; i++)
    {
	if (parser.getLine(line) != parser.ok) missingItem("line",fileName);

	ch = ' ';
	classA = 0;
	rc = sscanf(line,"%d:(%d,%d,%d),%c%d,(%le,%le,%le)", 
		       &trNo,&pNo1,&pNo2,&pNo3,&ch,&classA,&x,&y,&z);
	if (rc==0) break;
	   
	pt1 = ptAdr[pNo1];
	pt2 = ptAdr[pNo2];
	pt3 = ptAdr[pNo3];

	Stack<TR3*>& trs = *trsOfPts[pt1->node];

	FORALL(trs,k) 
	{
	    t = trs[k];
	    if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
	}

	if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);

	t->classA = classA;
	   
	if ((ch == 'd') || (ch == 'D')) 
	   {
	       t->boundP = DIRICHLET;
	       (t->e1)->boundP = DIRICHLET;
	       (t->e2)->boundP = DIRICHLET;
	       (t->e3)->boundP = DIRICHLET;
	       (t->e1)->classA = classA;
	       (t->e2)->classA = classA;
	       (t->e3)->classA = classA;
	   }
	else if ( (ch == 'c') || (ch == 'C'))
	{
	    t->boundP = CAUCHY;
	    
	    if ((t->e1)->boundP != DIRICHLET)
	    {
		(t->e1)->boundP = CAUCHY;
		(t->e1)->classA = classA;
	    }
	    if ((t->e2)->boundP != DIRICHLET) 
	    {
		(t->e2)->boundP = CAUCHY;
		(t->e2)->classA = classA;
	    }
	    if ((t->e3)->boundP != DIRICHLET) 
	    {
		(t->e3)->boundP = CAUCHY;
		(t->e3)->classA = classA;
	    }
	}
	else if ( (ch == 'n') || (ch == 'N'))
	{
	    t->boundP = NEUMANN;
	    if ((t->e1)->boundP != DIRICHLET) 
	    {
		(t->e1)->boundP = NEUMANN;
		(t->e1)->classA = classA;
	    }
	    if ((t->e2)->boundP != DIRICHLET) 
	    {
		(t->e2)->boundP = NEUMANN;
		(t->e2)->classA = classA;
	    }
	    if ((t->e3)->boundP != DIRICHLET) 
	    {
		(t->e3)->boundP = NEUMANN;
		(t->e3)->classA = classA;
	    }
	}
   }

    DListIter<TR3> iterTr(triangles());
    while ((t=iterTr.all()))  setNeighbsOfTetra(t);
    
    iterTr.reset();  
    while ((t=iterTr.all()))  
    if ( !(t->onBoundary()) ) 
    {
	trList[0]->remove(t);
	delete t;
    }
    
    iterTet.reset();
    while ((td=iterTet.all())) td->setBoundary();

    FORALL(edsOfPts,i) delete edsOfPts[i];
    FORALL(trsOfPts,i) delete trsOfPts[i];
    
    if (infoRefinement) Info();
    if (Cmd.isTrue("printMesh")) cout << *this;
    
    return;
}
//-------------------------------------------------------------------------


void MESH3:: readGrapeFormat(const char* /*fileName0*/)
{
    int   i, j, maxIndex, noOfTds;
    int   count=0;
    int   classA, pNo1, pNo2 ,pNo3, pNo4;
    PT3   *p, *pt1=0, *pt2=0, *pt3=0;
    TR3*  t=0;
    TET3* td;


    DList<PT3>*   ptL = new DList<PT3>;    ptList.push(ptL);
    DList<TET3>* tetL = new DList<TET3>;  tetList.push(tetL);
    DList<TR3>*   trL = new DList<TR3>;    trList.push(trL);
    DList<EDG3>* edgL = new DList<EDG3>;  edgList.push(edgL);

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


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

    checkMaxIndex(maxIndex);

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

    for (i=1; i <= maxIndex; i++)  
    {
	p = ptAlloc.Get();			// new PT3;
	ptAdr[i] = p;
	parser.getValue(&p->x);
	parser.getValue(&p->y);
	parser.getValue(&p->z);

	p->node = ++count;			// used in Find...
	ptList[0]->add(p);
	++noOfPoints;
    }


    rightHand = False;
    for (i=1; i <= noOfTds; i++)  
    {
	parser.getValue(&pNo1); 
	parser.getValue(&pNo2); 
	parser.getValue(&pNo3); 
	parser.getValue(&pNo4); 
	parser.getValue(&classA);   

	td = getTET3();				// new TET3;
	td->p1 = ptAdr[pNo1];
	td->p2 = ptAdr[pNo2];
	td->p3 = ptAdr[pNo3];
	td->p4 = ptAdr[pNo4];

	td->classA = 1;
	SetTPts(td,0); 

	tetList[0]->add(td);
	++noOfTetrahedra;
    }


    // --		compute edges and faces:

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

    Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
    FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);

    DListIter<TET3> iterTet(tetras());
    while ((td=iterTet.all()))  completeTetrahedron(td, edsOfPts, trsOfPts);

    int no, pNo[5];

    iterTet.reset();
    while ((td=iterTet.all())) 
    {
	parser.getValue(&no); pNo[1]=no;
	parser.getValue(&no); pNo[2]=no;
	parser.getValue(&no); pNo[3]=no;
	parser.getValue(&no); pNo[4]=no;

	for (j=1; j<=4; j++)
	{
	    if (pNo[j] != -1) continue;
	    if ( j == 1)
	       { pt1 = td->p2; pt2 = td->p3; pt3 = td->p4; }
	    else if ( j == 2)
	       { pt1 = td->p1; pt2 = td->p3; pt3 = td->p4; }
	    else if ( j == 3)
	       { pt1 = td->p1; pt2 = td->p2; pt3 = td->p4; }
	    else if ( j == 4)
	       { pt1 = td->p1; pt2 = td->p2; pt3 = td->p3; }

	    Stack<TR3*>& trs = *trsOfPts[pt1->node];

	    FORALL(trs,i) 
	    {
		t = trs[i];
		if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
	    }
	    if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);

	    t->boundP =  DIRICHLET; t->classA = 1;
    
	    t->e1->boundP = t->e2->boundP = t->e3->boundP = DIRICHLET;
	    t->e1->classA = t->e2->classA = t->e3->classA = 1;

	    t->p1->boundP = t->p2->boundP = t->p3->boundP = DIRICHLET;
	    t->p1->classA = t->p2->classA = t->p3->classA = 1;
	}
     
    }


    DListIter<TR3> iterTr(triangles());
    while ((t=iterTr.all()))  setNeighbsOfTetra(t);

    iterTr.reset();  
    while ((t=iterTr.all()))  
	if ( !(t->onBoundary()) ) 
	{
	   trList[0]->remove(t);
	   delete t;
	}

    iterTet.reset();
    while ((td=iterTet.all())) td->setBoundary();

    FORALL(edsOfPts,i) delete edsOfPts[i];
    FORALL(trsOfPts,i) delete trsOfPts[i];
   
    if (infoRefinement) Info();
    if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------


void MESH3:: readZIBFormat(const char* /*fileName0*/)
{
    int   i, maxIndex;
    TR3*  t;
    TET3* td;


    DList<PT3>*   ptL = new DList<PT3>;    ptList.push(ptL);
    DList<TET3>* tetL = new DList<TET3>;  tetList.push(tetL);
    DList<TR3>*   trL = new DList<TR3>;    trList.push(trL);
    DList<EDG3>* edgL = new DList<EDG3>;  edgList.push(edgL);

    BufferedParser parser(fileName, "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<PT3*> ptAdr(0,maxIndex);		// for point addresses
    FORALL(ptAdr,i) ptAdr[i] = 0;


    readPoints(parser, ptAdr);


    if (parser.findWord("Elem") != parser.ok) missingItem("tetrahedra",fileName);

    readElements(parser, ptAdr);


    // --		compute edges and faces:

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

    Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
    FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);

    DListIter<TET3> iterTet(tetras());
    while ((td=iterTet.all()))  completeTetrahedron(td, edsOfPts, trsOfPts);


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

    readBoundary(parser, ptAdr, trsOfPts);

    DListIter<TR3> iterTr(triangles());
    while ((t=iterTr.all()))  setNeighbsOfTetra(t);

    iterTr.reset();  
    while ((t=iterTr.all()))  
	if ( !(t->onBoundary()) ) 
	{
	   trList[0]->remove(t);
	   delete t;
	}

    iterTet.reset();
    while ((td=iterTet.all())) td->setBoundary();

    FORALL(edsOfPts,i) delete edsOfPts[i];
    FORALL(trsOfPts,i) delete trsOfPts[i];
   
    if (infoRefinement) Info();
    if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------

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

    while (parser.getValue(&ptNo) == parser.ok)  
    {
	if (ptNo > ptAdr.high()) pointIndexError(ptNo);
	
	p = ptAlloc.Get();			// new PT3;
	ptAdr[ptNo] = p;
	parser.getValue(&p->x);
	parser.getValue(&p->y);
	parser.getValue(&p->z);

	p->node = ++count;			// used in Find...

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

void MESH3:: readElements(BufferedParser& parser, Vector<PT3*>& ptAdr)
{
    int  pNo1, pNo2 ,pNo3, pNo4, classA;
    TET3 *td;

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

	td = getTET3();				// new TET3;
	td->p1 = ptAdr[pNo1];
	td->p2 = ptAdr[pNo2];
	td->p3 = ptAdr[pNo3];
	td->p4 = ptAdr[pNo4];

	parser.getValue(&classA);		// int -> char 
	td->classA = classA;
	SetTPts(td,0); 

	tetList[0]->add(td);
	++noOfTetrahedra;
    }
}
//-------------------------------------------------------------------------

void MESH3:: readBoundary(BufferedParser& parser, Vector<PT3*>& ptAdr, 
			  Vector<Stack<TR3*>*>& trsOfPts)
{
    int i, pNo1, pNo2 ,pNo3, pNo;
    PT3 *pt1, *pt2, *pt3, *p;
    TR3 *t = 0;

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

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

	Stack<TR3*>& trs = *trsOfPts[pt1->node];

	FORALL(trs,i) 
	{
	    t = trs[i];
	    if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
	}

	if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);
	t->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 MESH3:: readHyperFormat(const char* /*fileName0*/)
{
    int   i, nPoints, nEdges, nTriangles, nTetrahedrons;
    int   pointsSection, trianglesSection, tetrahedronsSection;
    int   pointDataSection, triangleDataSection, tetrahedronDataSection;
    TR3*  t;
    TET3* td;


    DList<PT3>*   ptL = new DList<PT3>;    ptList.push(ptL);
    DList<TET3>* tetL = new DList<TET3>;  tetList.push(tetL);
    DList<TR3>*   trL = new DList<TR3>;    trList.push(trL);
    DList<EDG3>* edgL = new DList<EDG3>;  edgList.push(edgL);

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

    if (parser.findWord("nPoints") != parser.ok) missingItem("nPoints",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    parser.getValue(&nPoints);  

    checkMaxIndex(nPoints);
    noOfPoints = nPoints;

    if (parser.findWord("nEdges") != parser.ok) missingItem("nEdges",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    parser.getValue(&nEdges);  

    checkMaxIndex(nEdges);

    if (parser.findWord("nTriangles") != parser.ok) 
      missingItem("nTriangles",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    parser.getValue(&nTriangles);  

    checkMaxIndex(nTriangles);

    if (parser.findWord("nTetrahedrons") != parser.ok)
                                   missingItem("nTetrahedrons",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    parser.getValue(&nTetrahedrons);  
    //cout << "\n nTetrahedrons = " << nTetrahedrons << "\n";
    checkMaxIndex(nTetrahedrons);
    noOfTetrahedra = nTetrahedrons;

    if (parser.findWord("Points") != parser.ok) missingItem("Points",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&pointsSection);  
    //cout << "\n pointsSection = " << pointsSection << "\n";
    checkMaxIndex(pointsSection);

    if (parser.findWord("Triangles") != parser.ok) 
    				missingItem("Triangles",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&trianglesSection);  
    //cout << "\n trianglesSection = " << trianglesSection << "\n";
    checkMaxIndex(trianglesSection);

    if (parser.findWord("Tetrahedrons") != parser.ok)
                                            missingItem("Tetrahedrons",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&tetrahedronsSection);  
    //cout << "\n tetrahedronsSection = " << tetrahedronsSection << "\n";
    checkMaxIndex(tetrahedronsSection);


    if (parser.findWord("PointData") != parser.ok) 
    				missingItem("PointData",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&pointDataSection);  
    //cout << "\n pointDataSection = " << pointDataSection << "\n";
    checkMaxIndex(pointDataSection);

    if (parser.findWord("TriangleData") != parser.ok) missingItem("TriangleData",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&triangleDataSection);  
    //cout << "\n triangleDataSection = " << triangleDataSection << "\n";
    checkMaxIndex(triangleDataSection);

    if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    parser.getValue(&tetrahedronDataSection);  
    //cout << "\n tetrahedronDataSection = " << tetrahedronDataSection << "\n";
    checkMaxIndex(tetrahedronDataSection);


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

    readHyperPoints(parser, pointsSection, ptAdr);

    readHyperElements(parser, tetrahedronsSection, ptAdr);

    readHyperClass(parser, tetrahedronDataSection);


    // --		compute edges and faces:

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

    Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
    FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);

    DListIter<TET3> iterTet(tetras());
    while ((td=iterTet.all()))  completeTetrahedron(td, edsOfPts, trsOfPts);

    //cout << "\n readHyperBoundary\n";
    readHyperBoundary(parser, trianglesSection, triangleDataSection, ptAdr, 
		      trsOfPts);

    DListIter<TR3> iterTr(triangles());
    while ((t=iterTr.all()))  setNeighbsOfTetra(t);

    iterTr.reset();  
    while ((t=iterTr.all()))  
	if ( !(t->onBoundary()) ) 
	{
	   trList[0]->remove(t);
	   delete t;
	}

    iterTet.reset();
    while ((td=iterTet.all())) td->setBoundary();

    FORALL(edsOfPts,i) delete edsOfPts[i];
    FORALL(trsOfPts,i) delete trsOfPts[i];
   
    if (infoRefinement) Info();
    if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------

void MESH3:: readHyperPoints(BufferedParser& parser, int pointsSection, Vector<PT3*>& ptAdr)
{
    int  count=0;
    PT3* p;

    int section = 0;
    while (section != pointsSection)
    {
 	if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    	parser.getValue(&section); 
    } 

    for (int i=1; i<=noOfPoints; i++)
    {
	
	p = ptAlloc.Get();			// new PT3;
	ptAdr[i] = p;
	parser.getValue(&p->x);
	parser.getValue(&p->y);
	parser.getValue(&p->z);
	p->node = ++count;			// used in Find...

	// cout << "\n p = " << p->x << " , " << p->y << " , " << p->z << "\n";
	ptList[0]->add(p);
    }
}
//-------------------------------------------------------------------------

void MESH3:: readHyperElements(BufferedParser& parser, int tetrahedronsSection, Vector<PT3*>& ptAdr)
{
    int  pNo1, pNo2 ,pNo3, pNo4;
    TET3 *td;

    int section = 0;
    while (section != tetrahedronsSection)
    {
 	if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    	parser.getValue(&section); 
    } 

    for (int i=1; i<=noOfTetrahedra; i++)
    {
	parser.getValue(&pNo1);
	parser.getValue(&pNo2);
	parser.getValue(&pNo3);
	parser.getValue(&pNo4);

	td = getTET3();				// new TET3;
	td->p1 = ptAdr[pNo1];
	td->p2 = ptAdr[pNo2];
	td->p3 = ptAdr[pNo3];
	td->p4 = ptAdr[pNo4];

	SetTPts(td,0); 

	// cout << "\n tet = " << pNo1 << ", " << pNo2 << ", " << pNo3 << ", " << pNo4 << "\n";

	tetList[0]->add(td);
    }
}
//-------------------------------------------------------------------------

void MESH3:: readHyperClass(BufferedParser& parser, int tetrahedronDataSection)
{
    int  classA;
    TET3 *td;

    int section = 0;
    while (section != tetrahedronDataSection)
    {
 	if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    	parser.getValue(&section); 
    } 


    DListIter<TET3> iterTet(tetras());
    while ((td=iterTet.all()))  
    {
	parser.getValue(&classA);		// int -> char 
	td->classA = classA;

	//cout << "\n tet->class = " << classA  << "\n";
    }
}

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

void MESH3:: readHyperBoundary(BufferedParser& parser, int trianglesSection, 
			       int triangleDataSection,
			       Vector<PT3*>& ptAdr, Vector<Stack<TR3*>*>& trsOfPts)
{
    int i, pNo1, pNo2 ,pNo3, pNo4, pNo5;
    PT3 *pt1, *pt2, *pt3;
    TR3 *t = 0;
    Vector<int> trClass(noOfTriangles), trBound(noOfTriangles);

    parser.rewind();

    if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);
    
    int section = 0;
    while (section != triangleDataSection)
    {
 	if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    	parser.getValue(&section); 
    } 

    for (i=1; i<=noOfTriangles; i++)		   // triangles on boundaries
    {
	parser.getValue(&pNo1);
	parser.getValue(&pNo2);
        trBound[i] = pNo1;
        trClass[i] = pNo2;
        // cout << "\n class, boundP = " << pNo1 << ", " << pNo2 << "\n";
    }
   
  
    parser.rewind();  

    if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
    if (parser.findWord("=")  != parser.ok) missingItem("=",fileName);

    section = 0;
    while (section != trianglesSection)
    {
 	if (parser.findChar('@')  != parser.ok) missingItem("@",fileName);
    	parser.getValue(&section); 
    } 


    for (int j=1; j<=noOfTriangles; j++)		   // triangles on boundaries
    {
	parser.getValue(&pNo1);
	parser.getValue(&pNo2);
	parser.getValue(&pNo3);
	parser.getValue(&pNo4);
	parser.getValue(&pNo5);

        // cout << "\ntriangles = " << pNo1 << pNo2 << pNo3 << pNo4 << pNo5 << "\n";
	pt1 = ptAdr[pNo1];
	pt2 = ptAdr[pNo2];
	pt3 = ptAdr[pNo3];

	Stack<TR3*>& trs = *trsOfPts[pt1->node];

	FORALL(trs,i) 
	{
	    t = trs[i];
	    if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
	}


	if ( t->onBoundary() ) 
        {
           t->classA = trClass[j];
           t->boundP = trBound[j];

	   t->readHyperBC(parser);
        }
    }

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

Bool MESH3:: pointOnEdge(PT3* p, EDG3* ed)
{
    return (ed->p1 == p || ed->p2 == p);
}
//-------------------------------------------------------------------------

Bool MESH3:: edgeOnTri(EDG3* ed, TR3* t)
{
    return (t->e1 == ed || t->e2 == ed || t->e3 == ed);
}
//-------------------------------------------------------------------------


EDG3* MESH3:: FindEdgByPts(PT3 *p1, PT3 *p2, Vector<Stack<EDG3*>*>& edsOfPts)
{ 
    int  i;
    EDG3 *ed;

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

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


TR3* MESH3:: FindTrByEds(PT3 *p1, PT3 *p2, PT3 *p3, EDG3 *e1, EDG3 *e2, EDG3 *e3,
			 Vector<Stack<TR3*>*>& trsOfPts)
{ 
    int i;
    TR3 *t;

    Stack<TR3*>& trs = *trsOfPts[p1->node];
    FORALL(trs,i) 
    {
	t = trs[i];
	if (edgeOnTri(e1,t) && edgeOnTri(e2,t) && edgeOnTri(e3,t)) return t;
    }

    t = new TR3;
    noOfTriangles++;

    t->e1 = e1; 
    t->e2 = e2; 
    t->e3 = e3;

    SetTP(t,0); 

    trsOfPts[p1->node]->push(t);
    trsOfPts[p2->node]->push(t);
    trsOfPts[p3->node]->push(t);

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


int MESH3:: completeTetrahedron(TET3 *td, Vector<Stack<EDG3*>*>& edsOfPts,
			                  Vector<Stack<TR3*>*>&  trsOfPts)
{ 
   PT3  *p1 = td->p1, *p2 = td->p2, *p3 = td->p3, *p4 = td->p4;
   EDG3 *e1, *e2, *e3, *e4, *e5, *e6;
   TR3  *t1, *t2, *t3, *t4;

   td->e1 = FindEdgByPts(p1,p3,edsOfPts);
   td->e2 = FindEdgByPts(p2,p3,edsOfPts);
   td->e3 = FindEdgByPts(p1,p2,edsOfPts);
   td->e4 = FindEdgByPts(p3,p4,edsOfPts);
   td->e5 = FindEdgByPts(p2,p4,edsOfPts);
   td->e6 = FindEdgByPts(p1,p4,edsOfPts);
 
   e1 = td->e1; e2 = td->e2; e3 = td->e3;
   e4 = td->e4; e5 = td->e5; e6 = td->e6;

   t1 = FindTrByEds(p2,p3,p4,e2,e4,e5,trsOfPts);
   t2 = FindTrByEds(p1,p3,p4,e1,e4,e6,trsOfPts);
   t3 = FindTrByEds(p1,p2,p4,e3,e5,e6,trsOfPts);
   t4 = FindTrByEds(p1,p2,p3,e1,e2,e3,trsOfPts);

   // neighborhood of this tetrahedron
   int k = 0;

   if ((t1->t31)==0)      { t1->t31 = td; k++; }
   else if ((t1->t32)==0) { t1->t32 = td; k++; }

   if ((t2->t31)==0)      { t2->t31 = td; k++; }
   else if ((t2->t32)==0) { t2->t32 = td; k++; }

   if ((t3->t31)==0)      { t3->t31 = td; k++; }
   else if ((t3->t32)==0) { t3->t32 = td; k++; }

   if ((t4->t31)==0)      { t4->t31 = td; k++; }
   else if ((t4->t32)==0) { t4->t32 = td; k++; }

   if (k!=4)
    {
        cout << "\n*** Tetra: impossible state (completeTetrahedron), k= " 
	     << k << "\n";
	cout.flush(); abort();
    }

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


void MESH3:: setNeighbsOfTetra(TR3* t)
{
   TET3 *tet1 = t->t31, *tet2 = t->t32;

   if (tet2 != nil) 
   {
     if (!PinTr(tet2->p1,t)) 
        {
          if (tet1) tet2->n1 = tet1; 
	  else { tet2->n1 = t; SetTP(tet2,t,2,3,4); }
        }
     else 
     if (!PinTr(tet2->p2,t)) 
        {
          if (tet1) tet2->n2 = tet1; 
	  else { tet2->n2 = t; SetTP(tet2,t,1,4,3); }
        }
     else 
     if (!PinTr(tet2->p3,t)) 
        {
          if (tet1) tet2->n3 = tet1;
	  else { tet2->n3 = t; SetTP(tet2,t,1,2,4); }
        }
     else 
     if (!PinTr(tet2->p4,t)) 
        {
          if (tet1) tet2->n4 = tet1; 
	  else { tet2->n4 = t; SetTP(tet2,t,1,3,2); }
        }
   }
   
   if (tet1 != nil) 
   {
     if (!PinTr(tet1->p1,t)) 
        {
          if (tet2) tet1->n1 = tet2; 
	  else { tet1->n1 = t; SetTP(tet1,t,2,3,4); }
        }
     else 
     if (!PinTr(tet1->p2,t)) 
        {
          if (tet2) tet1->n2 = tet2; 
	  else { tet1->n2 = t; SetTP(tet1,t,1,4,3); }
        }
     else 
     if (!PinTr(tet1->p3,t)) 
        {
          if (tet2) tet1->n3 = tet2; 
	  else { tet1->n3 = t; SetTP(tet1,t,1,2,4); }
        }
     else 
     if (!PinTr(tet1->p4,t)) 
        {
          if (tet2) tet1->n4 = tet2; 
	  else { tet1->n4 = t; SetTP(tet1,t,1,3,2); }
        }
   }
 
   return;
}
//-------------------------------------------------------------------------


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


void MESH3:: writeZIBFormat(const char* outFileText, const Vector<Real>& u)const 
{
   int  i;
   int  node;

   TET3 *tet;
   TR3  *t;
   PT3  *p;
   FILE *outFile;

   Vector<int> PNodes(noOfPoints);

   // saving the node numbers for later reuse
   i = 0;
   DListIter<PT3> 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<PT3> iterPts(points());
   while ((p=iterPts.all()))  {
      p->node = node+1;
      fprintf(outFile, " %d %e %e %e\n", p->node, p->x, p->y, p->z);
      node++;
   }
   fprintf(outFile, "END\n");

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

   DListIter<TET3> iterTet(tetras());
   while ((tet=iterTet.all()))  {
      fprintf(outFile, "%d %d %d %d %d\n",
                    (tet->p1)->node, (tet->p2)->node,
                    (tet->p3)->node, (tet->p4)->node,
                    tet->classA );
   }
   fprintf(outFile, "END\n");


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

   DListIter<TR3> iterTrs(triangles());
   while ((t=iterTrs.all()))  {
      if (!t->onBoundary()) continue;
      if (t->isDirichlet())
      {
         if (t->isoP == CIRCULAR)
           fprintf(outFile, "%d %d %d D %d S\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
         else if (t->isoP == CYLINDRICAL)
	      fprintf(outFile, "%d %d %d D %d C\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
	      else fprintf(outFile, "%d %d %d D %d P\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);

      }
      else if (t->isNeumann())
      {
         if (t->isoP == CIRCULAR)
           fprintf(outFile, "%d %d %d N %d S\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
         else if (t->isoP == CYLINDRICAL)
	  fprintf(outFile, "%d %d %d N %d C\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
	 else fprintf(outFile, "%d %d %d N %d P\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
      }
      else if (t->isCauchy())
      {
         if (t->isoP == CIRCULAR)
           fprintf(outFile, "%d %d %d C %d S\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
         else if (t->isoP == CYLINDRICAL)
	  fprintf(outFile, "%d %d %d C %d C\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
	 else fprintf(outFile, "%d %d %d C %d P\n",
                  (t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
      }
      else 
         fprintf(outFile, "%d %d %d I \n",
         (t->p1)->node, (t->p2)->node,(t->p3)->node );
   }

   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  ZIB-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 MESH3:: writeHyperFormat(const char* outFileText, const Vector<Real>& u)
									const 
{
   FILE *outFile;
   int  i;

   TET3 *tet;
   TR3  *t;
   EDG3 *ed;
   PT3  *p, *p1, *p2, *p3, *p4;
   PATCH *n1, *n2, *n3, *n4;

   Vector<int>  PNodes(noOfPoints),     ENodes(noOfEdges),
                TrNodes(noOfTriangles), TetNodes(noOfTetrahedra);
   Vector<char> TetMark(noOfTetrahedra);

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

   i = 0;
   DListIter<EDG3> EIter(edges());
   while ((ed=EIter.all()))  {
     ENodes[++i] = ed->node;
     ed->node = i;
   }
 
   i = 0;
   DListIter<TR3> TrIter(triangles());
   while ((t=TrIter.all()))  {
     TrNodes[++i] = t->node;
     t->node = i;
   }

   i = 0;
   DListIter<TET3> TetIter(tetras());
   while ((tet=TetIter.all()))  {
     TetNodes[++i] = tet->node;
     tet->node = i;
     TetMark[i] = tet->mark;
     tet->mark = 0;
   }


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

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

   fprintf(outFile, "# HyperMesh ASCII 0.2 \n"); 
   fprintf(outFile, "#                     \n"); 
   fprintf(outFile, "#                     \n"); 

   fprintf(outFile, "nPoints =  %d \n", noOfPoints);
   fprintf(outFile, "nEdges =  %d \n", noOfEdges);
   fprintf(outFile, "nTriangles =  %d \n", noOfTriangles);
   fprintf(outFile, "nTetrahedrons =  %d \n", noOfTetrahedra);

   fprintf(outFile, "#                     \n"); 

   fprintf(outFile, "Points { double[3] }  = @1  \n"); 
   fprintf(outFile, "Edges  { int[2] } = @2  \n"); 
   fprintf(outFile, "Triangles  { int[5] } = @3  \n"); 
   fprintf(outFile, "Tetrahedrons  { int[4] } = @4  \n"); 

   fprintf(outFile, "TriangleEncoding  = points tetrahedrons  \n"); 
   fprintf(outFile, "TetrahedronEncoding  = points  \n"); 

   fprintf(outFile, "PointData  { double } = @5  \n"); 
   fprintf(outFile, "TriangleData  { byte[2] } = @6  \n"); 
   fprintf(outFile, "TetrahedronData  { byte[2] } = @7  \n"); 
   fprintf(outFile, "#                     \n"); 
   fprintf(outFile, "#End of Header. Data follows.\n"); 
   fprintf(outFile, "#                     \n"); 

   fprintf(outFile, "\n@1  \n"); 
   DListIter<PT3> iterPts(points());		// coordinates of the points
   while ((p=iterPts.all()))  {
      fprintf(outFile, "%e %e %e\n", p->x, p->y, p->z);
   }


   fprintf(outFile, "\n@2  \n"); 
   DListIter<EDG3> iterEds(edges());
   while ((ed=iterEds.all()))  {
         fprintf(outFile, "%d %d\n", (ed->p1)->node, (ed->p2)->node);
   }



   fprintf(outFile, "\n@3  \n"); 
   int trsMark[1+4];
   Vector<char>  trBoundP(noOfTriangles), trClassA(noOfTriangles);

   DListIter<TET3> iterTets(tetras());
   TET3 *nbTet = 0;
   int trsM;
   i = 0;

   while ((tet=iterTets.all()))  {
         tet->getMarkedTriangles(trsMark);
	 p1=tet->p1; p2=tet->p2; p3=tet->p3; p4=tet->p4;
	 n1=tet->n1; n2=tet->n2; n3=tet->n3; n4=tet->n4;
	 if ( !(tet->orient()) ) {
	    p1 = tet->p2; p2 = tet->p1;
	    n1 = tet->n2; n2 = tet->n1;
	    trsM =trsMark[1]; trsMark[1] = trsMark[2]; trsMark[2] = trsM;
	 }

	 if (!trsMark[1]) 
	 {
           fprintf(outFile, "%d %d %d ", p3->node, p2->node, p4->node);
           nbTet = n1->castToTET3();
           if ( nbTet ) 
	   {
	     fprintf(outFile, " %d %d", tet->node, nbTet->node );
             nbTet->markTriangles(tet);
             trBoundP[++i] = INTERIOR; trClassA[i] = 1;
	   }
	   else        
	   {
	     fprintf(outFile, " %d 0", tet->node);
	     t = n1->castToTR3();
             trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
	   }
           fprintf(outFile, "\n");
	 }

	 if (!trsMark[2]) 
	 {
           fprintf(outFile, "%d %d %d ", p1->node, p3->node, p4->node);
           nbTet = n2->castToTET3();
           if ( nbTet ) 
	   {
	     fprintf(outFile, " %d %d", tet->node, nbTet->node);
             nbTet->markTriangles(tet);
             trBoundP[++i] = INTERIOR; trClassA[i] = 1;
	   }
	   else        
	   {
	     fprintf(outFile, " %d 0", tet->node);
	     t = n2->castToTR3();
             trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
	   }
           fprintf(outFile, "\n");
	 }

	 if (!trsMark[3]) 
	 {
           fprintf(outFile, "%d %d %d ",
                  p2->node, p1->node, p4->node);
           nbTet = n3->castToTET3();
           if ( nbTet ) 
	   {
	     fprintf(outFile, " %d %d", tet->node, nbTet->node);
             nbTet->markTriangles(tet);
             trBoundP[++i] = INTERIOR; trClassA[i] = 1;
	   }
	   else
	   {
             fprintf(outFile, " %d 0", tet->node);
	     t = n3->castToTR3();
             trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
	   }
           fprintf(outFile, "\n");
	 }

	 if (!trsMark[4]) 
	 {
           fprintf(outFile, "%d %d %d ",
                  p1->node, p2->node, p3->node);
           nbTet = n4->castToTET3();
           if ( nbTet ) 
	   {
	     fprintf(outFile, " %d %d", tet->node, nbTet->node);
             nbTet->markTriangles(tet);
             trBoundP[++i] = INTERIOR; trClassA[i] = 1;
	   }
	   else
	   {
             fprintf(outFile, " %d 0", tet->node);
	     t = n4->castToTR3();
             trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
	   }
           fprintf(outFile, "\n");
	 }

   }


   fprintf(outFile, "\n@4  \n"); 
   iterTets.reset();
   while ((tet=iterTets.all()))  
	{
	  p1=tet->p1; p2=tet->p2; p3=tet->p3; p4=tet->p4;
	  if ( !(tet->orient()) ) { p1 = tet->p2; p2 = tet->p1; }
      	  fprintf(outFile, "%d %d %d %d\n", p1->node, p2->node, p3->node, p4->node);
	}


   fprintf(outFile, "\n@5  \n"); 			// data at points
   iterPts.reset();
   while ((p=iterPts.all()))  fprintf(outFile, "%e\n", u[PNodes[p->node]]);


   fprintf(outFile, "\n@6  \n"); 			// data at triangles
   for (i=1; i<= noOfTriangles; i++)
   {
       fprintf(outFile, "%d %d \n", trBoundP[i], trClassA[i]);
   }

   fprintf(outFile, "\n@7  \n"); 			// data at tetrahedra
   iterTets.reset();
   while ((tet=iterTets.all()))  
   fprintf(outFile, "%d\n", tet->classA);

   fprintf(outFile, "\n#End of file.\n"); 			// data at tetrahedra

   cout << "\n  Hyper-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];

   i = 0;
   EIter.reset();
   while ((ed=EIter.all()))  ed->node = ENodes[++i];

   i = 0;
   TrIter.reset();
   while ((t=TrIter.all()))  t->node = TrNodes[++i];

   i = 0;
   TetIter.reset();
   while ((tet=TetIter.all()))  
	{
	  tet->node = TetNodes[++i];
          tet->mark = TetMark[i];
	}

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


void MESH3:: writeGrapeFormat(const char* outFileText, const Vector<Real>& u) 
									const 
{
   FILE *outFile;
   int  i;
   int  node;
   int  k[4];

   TET3 *tet;
   PT3  *p;

   Vector<int> PNodes(noOfPoints), TetNodes(noOfTetrahedra);

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

   i = 0;
   DListIter<TET3> TIter(tetras());
   while ((tet=TIter.all()))  TetNodes[++i] = tet->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\n", noOfTetrahedra);
   fprintf(outFile, "%d\n", noOfPoints);


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

   node = 0;
   DListIter<TET3> iterTet(tetras());
   while ((tet=iterTet.all()))  {
      tet->node = node+1;

      if (tet->orient())
         fprintf(outFile, "%d %d %d %d %d\n",
                    (tet->p1)->node, (tet->p2)->node,
                    (tet->p4)->node, (tet->p3)->node,
                    tet->classA );
      else
         fprintf(outFile, "%d %d %d %d %d\n",
                    (tet->p1)->node, (tet->p2)->node,
                    (tet->p3)->node, (tet->p4)->node,
                    tet->classA );

      node++;
   }

   iterTet.reset();
   while ((tet=iterTet.all()))  
   {
      TET3* n;
      n = (tet->n1)->castToTET3();
      if ( n == nil ) k[0] = -1;
      else k[0] = n->node;

      n = (tet->n2)->castToTET3();
      if ( n == nil ) k[1] = -1;
      else k[1] = n->node;

      if (tet->orient())
      {
        n = (tet->n3)->castToTET3();
        if ( n == nil ) k[3] = -1;
        else k[3] = n->node;

        n = (tet->n4)->castToTET3();
        if ( n == nil ) k[2] = -1;
        else k[2] = n->node;
      }
      else 
      {
        n = (tet->n3)->castToTET3();
        if ( n == nil ) k[2] = -1;
        else k[2] = n->node;

        n = (tet->n4)->castToTET3();
        if ( n == nil ) k[3] = -1;
        else k[3] = n->node;
      }
      fprintf(outFile, "%d %d %d %d\n", k[0], k[1], k[2], k[3]);
   }
    
   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 ((tet=TIter.all()))  tet->node = TetNodes[++i];
   
   return;
}
//-------------------------------------------------------------------------


void MESH3:: writeAVSFormat(const char* outFileText, const Vector<Real>& u)
									const 
{
   int  i;
   int  node, noOfPValues, noOfTetValues, noOfModValues;
   int  noOfComponents, noOfElems;

   TET3 *tet;
   PT3  *p;
   FILE *outFile;

   Vector<int> PNodes(noOfPoints);

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


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

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

   noOfPValues   = 1;
   noOfTetValues = 1;
   noOfModValues = 0;

   fprintf(outFile, "%d %d %d %d %d\n", 
           noOfPoints,noOfTetrahedra,noOfPValues,noOfTetValues, noOfModValues);


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

   node = 0;
   DListIter<TET3> iterTet(tetras());
   while ((tet=iterTet.all()))  {
      fprintf(outFile, "%d %d tet %d %d %d %d\n",
                    node, tet->classA,
                    (tet->p1)->node, (tet->p2)->node,
                    (tet->p3)->node, (tet->p4)->node
                   );
      node++;
   }


   noOfComponents = 1;
   noOfElems      = 1;

   fprintf(outFile, "%d %d\n", noOfComponents,noOfElems);
   fprintf(outFile, "Temperature, Degree\n");

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

   fprintf(outFile, "%d %d\n", noOfComponents,noOfElems);
   fprintf(outFile, "Material, ClassA\n");
   i = 0;
   iterTet.reset();
   while ((tet=iterTet.all()))  {
     fprintf(outFile, "%d %d\n", i, int(tet->classA));
     i++;
   }

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

   delete outText;
   fclose(outFile);

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

}
//-------------------------------------------------------------------------


void MESH3:: Refine()
{
    TET3* tet;
    int noPt[3], noEdg[3], noTr[3], noTet[3];
    Timer timer, accTimer;

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

    RemoveGreen();              // delete all green tetras

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

    DList<TET3>* tetL = new DList<TET3>;  
    tetList.push(tetL);

    moreRedTd = 0;
    DListIter<TET3> iter(tetras());
    while ((tet=iter.all()))  RefineTetra(tet);

    
    	// 	Algorithm: For all tetrahedra:
        // 	(1)  look if there are 3 refined edges, 
	//	two  of these on different faces,
        //    	and refine it red this case,
        //  	(2)  if no neighbor tetrahedron exists, "father neighbor"
        //     	tetrahedron needs to be red refined (green refined at
    	//	the last level).
	//	This routine does not use recursion.
    

    if (infoRefinement >= 2) cout << "\n    Refine: moreRedTd =";
    moreRedTd = 1;
    while (moreRedTd)
    {
	moreRedTd = 0;
	iter.reset();
	while ((tet=iter.all()))  CloseRedTet(tet);
	if (infoRefinement >= 2) cout << " " << moreRedTd;
	
	if (moreRedTd == 0) break;
	iter.reset();
	while ((tet=iter.all()))  RefineTetra(tet);
    }
    
    noPt[2]  = noOfPoints;
    noTr[2]  = noOfTriangles;
    noEdg[2] = noOfEdges;
    noTet[2] = noOfTetrahedra;


    MakeGreen();

    // now: tetrahedra with a boundary face are marked when generated
    // iter.reset();
    // while (tet=iter.all()) tet->setBoundary();

    ++refStep;

    if (Cmd.isTrue("triCheck"))
	consistencyCheck();


    if (infoRefinement) refinementInfo (noPt, noEdg, noTr, noTet);
    timeInfo(timer, accTimer);
    if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------

int MESH3:: CloseRedTet(TET3 *tet)
{
   TR3  *t;
   TR3  *trs[1+4];
   EDG3 *ed, *edSon1, *edSon2;
   EDG3 *eds[1+6];
   int i, noOfRefEds;
  
   //   request to refine a further tetrahedron tet

   if ( (tet->mark) == True)    
     {
       moreRedTd++;
       return True;
     }

   noOfRefEds = TestEdsInTd(tet);

   if ( noOfRefEds < 3) return True;
   if ( noOfRefEds > 3)
     { 
       tet->mark = True; 
       moreRedTd++;
       return True;
     }

   tet->mark = True;

   tet->getTriangles(trs);
   for (i=1; i<=4; i++)
      if (TestFace(trs[i]) == 3) {tet->mark = False; break;}          
   if ( tet->mark  == True)    
      {
        moreRedTd++;
	tet->returnTriangles(trs);
        return True;
      }

   tet->getEdges(eds);
   for (i=1; i<=6; i++)
       {
         ed = eds[i];
         edSon1 = ed->firstSon;
         if ( (edSon1) != 0)
           {
             edSon2 = edSon1->next;
             if ( (edSon1->firstSon != 0) ||
                  (edSon2->firstSon != 0)   )
               {
                tet->mark = True;
                moreRedTd++;
		tet->returnTriangles(trs);
                return True;
               }
           }
       }

   TET3 *nbTet;
   for (i=1; i<=4; i++)
   {
	t = trs[i]; nbTet = t->t32;
	if (nbTet == 0) continue;
	if (nbTet->refined()) 
        {
            if (testSons(t,nbTet) > 0)
            {
		tet->mark = True;
		moreRedTd++;
		tet->returnTriangles(trs);
		return True;
            }
        }
   }

   tet->returnTriangles(trs);

   return True;
}


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

void MESH3:: RemoveGreen()
{
   RelGreen3Td(); 
   RelGreen2ATd();
   RelGreen2BTd();
   RelGreen1Td(); 

   RelGreen();  
   RelGreenIrr(); 

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


int MESH3:: MakeGreen()
{
   TET3* tet;

   moreGreen1Td  = 0;
   moreGreen2ATd = 0;
   moreGreen2BTd = 0;
   moreGreen3Td  = 0;

   DListIter<TET3> iter(tetras());
   while ((tet=iter.all()))  MarkGreenTd(tet);

   if (infoRefinement >= 2)     
     cout << "\n    GREEN-I  -Closure: moreGreen-I-Td   = " 
          << moreGreen1Td<< "\n";

   iter.reset();
   while ((tet=iter.all()))  {
       if (!MkGreen1Td(tet)) {
          cout << "*** Tri: failed making GREEN-I tetrahedra (CloseRef)\n";
          cout.flush(); abort(); //return False;
       }
   }

   if (infoRefinement >= 2)    
     cout << "    GREEN-II -Closure: moreGreen-IIA-Td = " << moreGreen2ATd;

   iter.reset();
   while ((tet=iter.all()))  {
       if (!MkGreen2ATd(tet)) {
          cout << "Tri: failed making GREEN-II tetrahedra (CloseRef)\n";
          cout.flush(); abort(); // return False;
       }
   }

   if (infoRefinement >= 2) 
     cout << "  moreGreen-IIB-Td = " << moreGreen2BTd << "\n";
   iter.reset();
   while ((tet=iter.all()))  {
       if (!MkGreen2BTd(tet)) {
          cout << "Tri: failed making GREEN-II tetrahedra (CloseRef)\n";
          cout.flush(); abort(); // return False;
       }
   }

   if (infoRefinement >= 2)
     cout << "    GREEN-III-Closure: moreGreen-III-Td = " << moreGreen3Td <<"\n";

   iter.reset();
   while ((tet=iter.all()))  {
       if (!MkGreen3Td(tet)) {
          cout << "Tri: failed making GREEN-III tetrahedra (CloseRef)\n";
          return False;
      }
   }

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

int MESH3:: RefineTetra(TET3 *td)
{
	//   internal routine to refine a tetrahedron, loTrues for
	    
   if (td->mark) td->mark = False;         /* return if not marked */
   else return True;
   if ((td->firstSon)!=0)
     {  
       cout << "\n*** RefineTetra:??????? marked but already refined \n";
       return True;  // return if already marked 
      }

   moreRedTd++;

   if (refinementType == ShortDiag)
   {
       if ( !(RefineTdShortDiag(td)) ) { 
	   cout << "\n*** RefineTetra: tetrahedron not refined \n"; 
	   cout.flush(); abort();
       }
   }
   else if (refinementType == BeyStabilDiag) 
   {
       if ( !(RefineTdBey(td)) ) { 
	   cout << "\n*** RefineTetra: tetrahedron not refined \n"; 
	   cout.flush(); abort();
       }
   }
   else {
       cout << "\n*** RefineTetra: unknown refinement type requested\n"; 
       cout.flush(); abort();
   }
   return True;
}
//-------------------------------------------------------------------------


int MESH3:: RefineTdShortDiag(TET3 *td)
  {
	//   Refinement of a tetrahedron
	//   Method of Zang (Select the shortest diagonal in the octahedron)

    TET3   *td1, *td2, *td3, *td4, *td5, *td6, *td7, *td8;
    TR3  *tr, *tr1, *tr2, *tr3, *tr4,
         *it1, *it2, *it3, *it4, *it5, *it6, *it7, *it8,
         *t11=0,  *t12=0,  *t13=0,  *t14=0, 
         *t21=0,  *t22=0,  *t23=0,  *t24=0, 
         *t31=0,  *t32=0,  *t33=0,  *t34=0, 
         *t41=0,  *t42=0,  *t43=0,  *t44=0, 
         *t1, *t2, *t3, *t4;
    TR3  *trs[1+4];
    PT3  *p1 = td->p1,*p2 = td->p2, *p3 = td->p3, *p4 = td->p4,
         *p5, *p6, *p7, *p8, *p9, *p10; 
    EDG3 *e1 = td->e1,  *e2 = td->e2,  *e3 = td->e3, 
         *e4 = td->e4,  *e5 = td->e5,  *e6 = td->e6,
         *ie1, *ed;
    PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
    Real len59, len78, len610;
    int  k, diag59=False, diag78=False, diag610=False;
    int  depth;


     td->getTriangles(trs);
     t1=trs[1]; t2=trs[2]; t3=trs[3]; t4=trs[4];


    // Get 8 new tetrahedra, 8 new triangles, 1 new edges
    //      return False if not enough space

    td1 = getTET3();		// new TET3; 
    td2 = getTET3();
    td3 = getTET3();
    td4 = getTET3();
    td5 = getTET3();
    td6 = getTET3();
    td7 = getTET3();
    td8 = getTET3();

    noOfTetrahedra += 7;

    it1 = new TR3;		// new TR3; 
    it2 = new TR3;
    it3 = new TR3;
    it4 = new TR3;
    it5 = new TR3;
    it6 = new TR3;
    it7 = new TR3;
    it8 = new TR3;

    noOfTriangles += 8;    

    ie1 = edgAlloc.Get();		// new EDG3;

    noOfEdges += 1; 

    if ( (td8==0)||(it8==0)||(ie1==0) ) 
       {
         cout << " no memory in RefineTd \n";
         return False;
       }


   depth = (td->depth) + 1;

   ie1->depth = depth;

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

   it1->depth = depth;
   it2->depth = depth;
   it3->depth = depth;
   it4->depth = depth;
   it5->depth = depth;
   it6->depth = depth;
   it7->depth = depth;
   it8->depth = depth;

    
    //  Set the "son" triangles tik, i=1,2,3,4; k=1,2,3,4;
    //  each triangle is refined 


    // Triangle 1 
    if (!RefineTr(t1))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t1);
    tr1 = t1->firstSon;
    if (tr1==0) return False;    // No space for Refinement 

    // Triangle 2 
    if (!RefineTr(t2))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    tr2 = t2->firstSon;
    setTypeOfInnerEdges(t2);
    if (tr2==0) return False;    // No space for Refinement 


    // Triangle 3 
    if (!RefineTr(t3))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t3);
    tr3 = t3->firstSon;
    if (tr3==0) return False;    /* No space for Refinement */

    // Triangle 4 
    if (!RefineTr(t4))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t4);
    tr4 = t4->firstSon;
    if (tr4==0) return False;    /* No space for Refinement */

    p5 = e3->pm;
    p6 = e5->pm;
    p7 = e2->pm;
    p8 = e6->pm;
    p9 = e4->pm;
    p10= e1->pm;


    //  one new edge in the octahedron  

    len59  = Length(p5,p9);
    len78  = Length(p7,p8);
    len610 = Length(p6,p10);

    if ( (len59 <= len78) && (len59 <= len610) )
       {
        ie1->p1 = p5;
        ie1->p2 = p9;
        diag59 = True;
       }
       else if ( (len78 <= len59) && (len78 <= len610) )
              {
               ie1->p1 = p7;
               ie1->p2 = p8;
               diag78 = True;
              }
             else
              {
               ie1->p1 = p6;
               ie1->p2 = p10;
               diag610 = True;
              }

    ie1->pm = 0;
    //    ie1->bound = 0;   


    //  new faces         
    tr = tr1;
    for (k = 0; k<4; k++)
      {
        if (PinTr(p2,tr)) t12 = tr;
          else if (PinTr(p3,tr)) t13 = tr;
                 else if (PinTr(p4,tr)) t14 = tr;
                        else 
                         {
                          t11 = tr;
                          ed = EdinTr(tr,p6,p7);
                          it2->e1 = ed;
                          if(diag78) it5->e1 = ed;
                             else if(diag610) it7->e1 = ed;
                          ed = EdinTr(tr,p7,p9);
                          it3->e1 = ed;
                          if(diag59||diag78) it7->e1 = ed;
                          ed = EdinTr(tr,p6,p9);
                          it4->e1 = ed;
                          if(diag59) it5->e1 = ed;
                            else if(diag610) it8->e2 = ed;
                         }
        tr = tr->next;
      }

    tr = tr2;
    for (k = 0; k<4; k++)
      {
        if (PinTr(p1,tr)) t21 = tr;
          else if (PinTr(p3,tr)) t23 = tr;
                 else if (PinTr(p4,tr)) t24 = tr;
                        else 
                         {
                          t22 = tr;
                          ed = EdinTr(tr,p8,p10);
                          it1->e1 = ed;
                          if(diag78) it8->e1 = ed;
                            else if(diag610) it6->e1 = ed;
                          ed = EdinTr(tr,p8,p9);
                          it4->e2 = ed;
                          if(diag59) it6->e1 = ed;
                             else if(diag78) it7->e2 = ed;
                          ed = EdinTr(tr,p9,p10);
                          it3->e2 = ed;
                          if(diag59||diag610) it8->e1 = ed;
                         }
         tr = tr->next;
      }


    tr = tr3;
    for (k = 0; k<4; k++)
      {
        if (PinTr(p1,tr)) t31 = tr;
          else if (PinTr(p2,tr)) t32 = tr;
                 else if (PinTr(p4,tr)) t34 = tr;
                        else 
                         {
                          t33 = tr;
                          ed = EdinTr(tr,p5,p8);
                          it1->e2 = ed;
                          if(diag59||diag78) it6->e2 = ed;
                          ed = EdinTr(tr,p6,p8);
                          it4->e3 = ed;
                          if(diag78) it5->e2 = ed;
                             else if(diag610) it6->e2 = ed;
                          ed = EdinTr(tr,p6,p5);
                          it2->e2 = ed;
                          if(diag59||diag610) it5->e2 = ed;
                         }
        tr = tr->next;
      }

    
   tr = tr4;
   for (k = 0; k<4; k++)
      {
        if (PinTr(p1,tr)) t41 = tr;
          else if (PinTr(p2,tr)) t42 = tr;
                 else if (PinTr(p3,tr)) t43 = tr;
                        else 
                         {
                          t44 = tr;
                          ed = EdinTr(tr,p5,p10);
                          it1->e3 = ed;
                          if(diag59) it8->e2 = ed;
                             else if(diag610) it5->e1 = ed;
                          ed = EdinTr(tr,p7,p10);
                          it3->e3 = ed;
                          if(diag610) it7->e2 = ed;
                             else if(diag78) it8->e2 = ed;
                          ed = EdinTr(tr,p5,p7);
                          it2->e3 = ed;
                          if(diag59) it7->e2 = ed;
                            else if(diag78) it6->e1 = ed;
                        }
        tr = tr->next;
      }


    //   four new triangles on the octahedron   
    SetTP(it1,0);
    SetTP(it2,0);
    SetTP(it3,0);
    SetTP(it4,0);

    //   four new edges in the octahedron   
    it5->e3 = ie1;
    it6->e3 = ie1;
    it7->e3 = ie1;
    it8->e3 = ie1;

    SetTP(it5,0);
    SetTP(it6,0);
    SetTP(it7,0);
    SetTP(it8,0);


    //   Build the new tetrahedra.
 
    //   tetrahedra generated by chopping off the 4 corners  

    SetTdP(td1,t21,t31,t41,it1,td);
    SetTdE(td1,t21,t31,t41,it1);
    neighbourAtFace(td1,n2,t21);
    neighbourAtFace(td1,n3,t31);
    neighbourAtFace(td1,n4,t41);

    SetTdP(td2,t12,t32,t42,it2,td);
    SetTdE(td2,t12,t32,t42,it2); 
    neighbourAtFace(td2,n1,t12);
    neighbourAtFace(td2,n3,t32);
    neighbourAtFace(td2,n4,t42);

    SetTdP(td3,t13,t23,t43,it3,td);
    SetTdE(td3,t13,t23,t43,it3); 
    neighbourAtFace(td3,n1,t13);
    neighbourAtFace(td3,n2,t23);
    neighbourAtFace(td3,n4,t43);

    SetTdP(td4,t14,t24,t34,it4,td);
    SetTdE(td4,t14,t24,t34,it4); 
    neighbourAtFace(td4,n1,t14);
    neighbourAtFace(td4,n2,t24);
    neighbourAtFace(td4,n3,t34);

    //   tetrahedra generated by the inner octahedron  
    if (diag59)
      {
       SetTdP(td5,t33,it5,it4,it6,td);
       SetTdE(td5,t33,it5,it4,it6); 
       neighbourAtFace(td5,n3,t33);
       innerNeighbourAtFace(td5,it5,td6);
       innerNeighbourAtFace(td5,it4,td4);
       innerNeighbourAtFace(td5,it6,td7);
       innerNeighbourAtFace(td4,it4,td5);

       SetTdP(td6,t11,it7,it2,it5,td);
       SetTdE(td6,t11,it7,it2,it5); 
       neighbourAtFace(td6,n1,t11); 
       innerNeighbourAtFace(td6,it7,td8);
       innerNeighbourAtFace(td6,it2,td2);
       innerNeighbourAtFace(td6,it5,td5);
       innerNeighbourAtFace(td2,it2,td6);

       SetTdP(td7,t22,it8,it6,it1,td);
       SetTdE(td7,t22,it8,it6,it1); 
       neighbourAtFace(td7,n2,t22); 
       innerNeighbourAtFace(td7,it8,td8);
       innerNeighbourAtFace(td7,it6,td5);
       innerNeighbourAtFace(td7,it1,td1);
       innerNeighbourAtFace(td1,it1,td7);

       SetTdP(td8,t44,it3,it8,it7,td);
       SetTdE(td8,t44,it3,it8,it7); 
       neighbourAtFace(td8,n4,t44); 
       innerNeighbourAtFace(td8,it3,td3);
       innerNeighbourAtFace(td8,it8,td7);
       innerNeighbourAtFace(td8,it7,td6);
       innerNeighbourAtFace(td3,it3,td8);
      }
      else if (diag78)
        {
         SetTdP(td5,t33,it5,it2,it6,td);
         SetTdE(td5,t33,it5,it2,it6); 
         neighbourAtFace(td5,n3,t33); 
         innerNeighbourAtFace(td5,it5,td6);
         innerNeighbourAtFace(td5,it2,td2);
         innerNeighbourAtFace(td5,it6,td7);
         innerNeighbourAtFace(td2,it2,td5);

         SetTdP(td6,t11,it7,it4,it5,td);
         SetTdE(td6,t11,it7,it4,it5); 
         neighbourAtFace(td6,n1,t11); 
         innerNeighbourAtFace(td6,it7,td8);
         innerNeighbourAtFace(td6,it4,td4);
         innerNeighbourAtFace(td6,it5,td5);
         innerNeighbourAtFace(td4,it4,td6);

         SetTdP(td7,t44,it8,it6,it1,td);
         SetTdE(td7,t44,it8,it6,it1); 
         neighbourAtFace(td7,n4,t44); 
         innerNeighbourAtFace(td7,it8,td8);
         innerNeighbourAtFace(td7,it6,td5);
         innerNeighbourAtFace(td7,it1,td1);
         innerNeighbourAtFace(td1,it1,td7);

         SetTdP(td8,t22,it3,it8,it7,td);
         SetTdE(td8,t22,it3,it8,it7); 
         neighbourAtFace(td8,n2,t22);
         innerNeighbourAtFace(td8,it3,td3);
         innerNeighbourAtFace(td8,it8,td7);
         innerNeighbourAtFace(td8,it7,td6);
         innerNeighbourAtFace(td3,it3,td8); 
        }
        else 
          {
           SetTdP(td5,t33,it5,it1,it6,td);
           SetTdE(td5,t33,it5,it1,it6); 
           neighbourAtFace(td5,n3,t33); 
           innerNeighbourAtFace(td5,it5,td7);
           innerNeighbourAtFace(td5,it1,td1);
           innerNeighbourAtFace(td5,it6,td6);
           innerNeighbourAtFace(td1,it1,td5); 

           SetTdP(td6,t22,it6,it8,it4,td);
           SetTdE(td6,t22,it6,it8,it4); 
           neighbourAtFace(td6,n2,t22); 
           innerNeighbourAtFace(td6,it6,td5);
           innerNeighbourAtFace(td6,it8,td8);
           innerNeighbourAtFace(td6,it4,td4);
           innerNeighbourAtFace(td4,it4,td6); 

           SetTdP(td7,t44,it2,it5,it7,td);
           SetTdE(td7,t44,it2,it5,it7); 
           neighbourAtFace(td7,n4,t44);
           innerNeighbourAtFace(td7,it2,td2);
           innerNeighbourAtFace(td7,it5,td5);
           innerNeighbourAtFace(td7,it7,td8);
           innerNeighbourAtFace(td2,it2,td7); 
 
           SetTdP(td8,t11,it3,it8,it7,td);
           SetTdE(td8,t11,it3,it8,it7);
           neighbourAtFace(td8,n1,t11);
           innerNeighbourAtFace(td8,it3,td3);
           innerNeighbourAtFace(td8,it8,td6);
           innerNeighbourAtFace(td8,it7,td7);
           innerNeighbourAtFace(td3,it3,td8); 
          }


    td1->type = T_WHITE; 
    td2->type = T_WHITE;
    td3->type = T_WHITE; 
    td4->type = T_WHITE; 
    td5->type = T_WHITE; 
    td6->type = T_WHITE; 
    td7->type = T_WHITE;
    td8->type = T_WHITE; 
 
    td->firstSon = td1;
    td->type = T_RED;
 
    tetList[tetList.h]->add(td1);
    tetList[tetList.h]->add(td2);
    tetList[tetList.h]->add(td3);
    tetList[tetList.h]->add(td4);
    tetList[tetList.h]->add(td5);
    tetList[tetList.h]->add(td6);
    tetList[tetList.h]->add(td7);
    tetList[tetList.h]->add(td8);

    // destruction of the faces
    delete it1;
    delete it2;
    delete it3;
    delete it4;
    delete it5;
    delete it6;
    delete it7;
    delete it8;

    if ( !(t1->onBoundary()) ) {
       delete t11; 
       delete t12;
       delete t13;
       delete t14;
    }

    if ( !(t2->onBoundary()) ) {
       delete t21;
       delete t22;
       delete t23;
       delete t24;
    }

    if ( !(t3->onBoundary()) ) {
       delete t31;
       delete t32;
       delete t33;
       delete t34;
    }

    if ( !(t4->onBoundary()) ) {
       delete t41;
       delete t42;
       delete t43;
       delete t44;
    }

    td->returnTriangles(trs);

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

int MESH3:: faceNoInTd(TR3 *t, TET3 *tet)
{
	PT3 *p1=tet->p1, *p2=tet->p2, *p3=tet->p3, *p4=tet->p4;

	if ( PinTr(p2,t) && PinTr(p3,t) && PinTr(p4,t) ) return 1;
	if ( PinTr(p1,t) && PinTr(p3,t) && PinTr(p4,t) ) return 2;
	if ( PinTr(p1,t) && PinTr(p2,t) && PinTr(p4,t) ) return 3;
	if ( PinTr(p1,t) && PinTr(p2,t) && PinTr(p3,t) ) return 4;
	
	return 0;
}


void MESH3:: innerNeighbourAtFace(TET3* td, TR3* t, TET3* neighbourTd)
{
        int noN = faceNoInTd(t,td);

	if (neighbourTd == 0) 
	   cout << "\n\terror in innerNeighbourAtFace: neighbourTd == 0\n";

	switch (noN) {
		case 1 : 
			if (td->n1 != 0) 
			   cout << "\n\terror in innerNeighbourAtFace\n";
			td->n1 = neighbourTd; return;
		case 2 : 
			if (td->n2 != 0) 
			   cout << "\n\terror in innerNeighbourAtFace\n";
			td->n2 = neighbourTd; return;
		case 3 : 
			if (td->n3 != 0) 
			   cout << "\n\terror in innerNeighbourAtFace\n";
			td->n3 = neighbourTd; return;
		case 4 : 
			if (td->n4 != 0) 
			   cout << "\n\terror in innerNeighbourAtFace\n";
			td->n4 = neighbourTd; return;
		default: break;
   	}

	cout << "\n\n  dangerous situation in innerNeighbourAtFace\n";

	return;
}


void MESH3:: setNeighbourInTd(TET3* td, PATCH* n, int noOfNeighbour)
{
	switch (noOfNeighbour) {
		case 1 : 
			td->n1 = n; 
			return;
		case 2 : 
			td->n2 = n; 
			return;
		case 3 : 
			td->n3 = n; 
			return;
		case 4 : 
			td->n4 = n; 
			return;
		default: 
			cout << "\n\n  error in setNeighbourInTd\n"; 
			return;
		}
}


void MESH3:: neighbourAtFace(TET3* td, PATCH* N, TR3* t)
{
	TET3 *nb, *nbSon;
	int   no, no1, no2, noN;
   
        noN = faceNoInTd(t,td);

        nb = N->castToTET3();
	if ( nb == nil )  			   // N is boundary face
	{
           setNeighbourInTd(td,t,noN);	
	   switch (noN) {
		case 1 : SetTP(td,t,2,3,4); td->boundP=BOUNDARY; return;
		case 2 : SetTP(td,t,1,4,3); td->boundP=BOUNDARY; return;
		case 3 : SetTP(td,t,1,2,4); td->boundP=BOUNDARY; return;
		case 4 : SetTP(td,t,1,3,2); td->boundP=BOUNDARY; return;	   
                default: 
			cout << "\n\n  error in neighbourAtFace\n"; 
			return;
	   }
	}


	nbSon = nb->firstSon;			     // N is tetrahedron
	if (nbSon == nil)  { 
	   setNeighbourInTd(td,nb,noN);
  	   if(!nb->mark) {
	   	no = faceNoInTd(t,nb);				
	   	if (no) setNeighbourInTd(nb,td,no);		
	   }
	   return; 
	}

	int iEnd = 0;
	if      (nb->type == T_RED)      iEnd = 8;
	else if (nb->type == T_GREEN_1)  iEnd = 2;
	else if (nb->type == T_GREEN_2A) iEnd = 4;
	else if (nb->type == T_GREEN_2B) iEnd = 3;
	else if (nb->type == T_GREEN_3)  iEnd = 4;

	for (int i=1; i<=iEnd; i++)
        {
		no1 = faceNoInTd(t,nbSon);
       		if ( no1 )  
		{
		  if ( nbSon->firstSon != 0 )
		  {
			int jEnd = 0;
			if      (nbSon->type == T_RED)      jEnd = 8;
			else if (nbSon->type == T_GREEN_1)  jEnd = 2;
			else if (nbSon->type == T_GREEN_2A) jEnd = 4;
			else if (nbSon->type == T_GREEN_2B) jEnd = 3;
			else if (nbSon->type == T_GREEN_3)  jEnd = 4;

			TET3* nbSonSon = nbSon->firstSon;
			for (int j=1; j<=jEnd; j++)
			{
			   no2 = faceNoInTd(t,nbSonSon);
			   if ( no2 ) 
			   {
				setNeighbourInTd(nbSonSon,td,no2);
				setNeighbourInTd(td,nbSonSon,noN);
				return;
			   }
			   else nbSonSon = nbSonSon->next;
			}
		  }
		  setNeighbourInTd(nbSon,td,no1);
		  setNeighbourInTd(td,nbSon,noN);
	          return;
		}
            	else nbSon = nbSon->next;
	}
		   
	return;
}
//-------------------------------------------------------------------------


int MESH3:: RefineTdBey(TET3 *td)
{
    TET3  *td1, *td2, *td3, *td4, *td5, *td6, *td7, *td8;
    TR3   *tr=0, *tr1=0, *tr2=0, *tr3=0, *tr4=0,
          *it1=0, *it2=0, *it3=0, *it4=0, *it5=0, *it6=0, *it7=0, *it8,
          *t11=0, *t12=0, *t13=0, *t14=0,
          *t21=0, *t22=0, *t23=0, *t24=0,
          *t31=0, *t32=0, *t33=0, *t34=0,
          *t41=0, *t42=0, *t43=0, *t44=0,
          *t1, *t2, *t3, *t4;
    TR3  *trs[1+4];
    PT3   *p1 = td->p1,*p2 = td->p2, *p3 = td->p3, *p4 = td->p4,
          *p5, *p6, *p7, *p8, *p9, *p10; 
    EDG3  *e1 = td->e1,  *e2 = td->e2,  *e3 = td->e3, 
          *e4 = td->e4,  *e5 = td->e5,  *e6 = td->e6,
          *ie1, *ed;
    PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
    int   k;
    int   depth;


	// Regular refinement of a tetrahedron, Method of J. Bey


    td->getTriangles(trs);

    t1=trs[1]; t2=trs[2]; t3=trs[3]; t4=trs[4];

    //     Get 8 new tetrahedra, 8 new triangles, 1 new edges
    //     return False if not enough space

    td1 = getTET3();			// new TET3; 
    td2 = getTET3();	
    td3 = getTET3();	
    td4 = getTET3();	
    td5 = getTET3();	
    td6 = getTET3();	
    td7 = getTET3();	
    td8 = getTET3();	

    noOfTetrahedra += 7;

    it1 = new TR3;		// new TR3; 
    it2 = new TR3;
    it3 = new TR3;
    it4 = new TR3;
    it5 = new TR3;
    it6 = new TR3;
    it7 = new TR3;
    it8 = new TR3;

    noOfTriangles += 8;    

    ie1 = edgAlloc.Get();		// new EDG3;

    noOfEdges += 1; 

    if ((td8==0)||(it8==0)||(ie1==0)) 
       {
         cout << " no memory in RefineTd \n";
         return False;
       }


    depth = (td->depth) + 1;

    ie1->depth = depth;

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


   it1->depth = depth;
   it2->depth = depth;
   it3->depth = depth;
   it4->depth = depth;
   it5->depth = depth;
   it6->depth = depth;
   it7->depth = depth;
   it8->depth = depth;

   //	Set the "son" triangles tik, i=1,2,3,4; k=1,2,3,4;
   //	if neeccessary refine the triangles

    // Triangle 1 
    if (!RefineTr(t1))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t1);
    tr1 = t1->firstSon;
    if (tr1==0) return False;    // No space for Refinement 

    // Triangle 2 
    if (!RefineTr(t2))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    tr2 = t2->firstSon;
    setTypeOfInnerEdges(t2);
    if (tr2==0) return False;    // No space for Refinement 


    // Triangle 3 
    if (!RefineTr(t3))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t3);
    tr3 = t3->firstSon;
    if (tr3==0) return False;    /* No space for Refinement */

    // Triangle 4 
    if (!RefineTr(t4))
      {
	cout << " failed in  RefineTr called by RefineTd \n";
	return False;
      }
    setTypeOfInnerEdges(t4);
    tr4 = t4->firstSon;
    if (tr4==0) return False;    /* No space for Refinement */


    p5 = e3->pm;
    p6 = e5->pm;
    p7 = e2->pm;
    p8 = e6->pm;
    p9 = e4->pm;
    p10= e1->pm;


    //  one new edge in the octahedron       
    ie1->p1 = p6;
    ie1->p2 = p10;
    ie1->pm = 0;


    //  new faces         
    tr = tr1;
    for (k = 0; k<4; k++)
    {
	if (PinTr(p2,tr)) t12 = tr;
	else if (PinTr(p3,tr)) t13 = tr;
	     else if (PinTr(p4,tr)) t14 = tr;
		  else 
                         {
                          t11 = tr;
                          ed = EdinTr(tr,p6,p7);
                          it2->e1 = ed;
                          it7->e1 = ed;
                          ed = EdinTr(tr,p7,p9);
                          it3->e1 = ed;
                          ed = EdinTr(tr,p6,p9);
                          it4->e1 = ed;
                          it8->e2 = ed;
                         }
        tr = tr->next;
    }


    tr = tr2;
    for (k = 0; k<4; k++)
    {
	if (PinTr(p1,tr)) t21 = tr;
	else if (PinTr(p3,tr)) t23 = tr;
	     else if (PinTr(p4,tr)) t24 = tr;
		   else 
                         {
                          t22 = tr;
                          ed = EdinTr(tr,p8,p10);
                          it1->e1 = ed;
                          it6->e1 = ed;
                          ed = EdinTr(tr,p8,p9);
                          it4->e2 = ed;
                          ed = EdinTr(tr,p9,p10);
                          it3->e2 = ed;
                          it8->e1 = ed;
                         }
         tr = tr->next;
    }


    tr = tr3;
    for (k = 0; k<4; k++)
    {
	if (PinTr(p1,tr)) t31 = tr;
	else if (PinTr(p2,tr)) t32 = tr;
	     else if (PinTr(p4,tr)) t34 = tr;
		  else 
                         {
                          t33 = tr;
                          ed = EdinTr(tr,p5,p8);
                          it1->e2 = ed;
                          ed = EdinTr(tr,p6,p8);
                          it4->e3 = ed;
                          it6->e2 = ed;
                          ed = EdinTr(tr,p6,p5);
                          it2->e2 = ed;
                          it5->e2 = ed;
                         }
        tr = tr->next;
    }

    
    tr = tr4;
    for (k = 0; k<4; k++)
       {
        if (PinTr(p1,tr)) t41 = tr;
        else if (PinTr(p2,tr)) t42 = tr;
             else if (PinTr(p3,tr)) t43 = tr;
                  else 
                         {
                          t44 = tr;
                          ed = EdinTr(tr,p5,p10);
                          it1->e3 = ed;
                          it5->e1 = ed;
                          ed = EdinTr(tr,p7,p10);
                          it3->e3 = ed;
                          it7->e2 = ed;
                          ed = EdinTr(tr,p5,p7);
                          it2->e3 = ed;
                        }
        tr = tr->next;
      }



    //  four new triangles on the octahedron   
    SetTP(it1,0);
    SetTP(it2,0);
    SetTP(it3,0);
    SetTP(it4,0);

    //  four new triangles in the octahedron   
    it5->e3 = ie1;
    it6->e3 = ie1;
    it7->e3 = ie1;
    it8->e3 = ie1;

    SetTP(it5,0);
    SetTP(it6,0);
    SetTP(it7,0);
    SetTP(it8,0);


   // Build the new tetrahedra.
 
   //  tetrahedra generated by chopping off the 4 corners  

    td1->p1 = p1;  td1->p2 = p5;  td1->p3 = p10; td1->p4 = p8;
    SetTdPBey(td1,td);
    SetTdE(td1,it1,t21,t31,t41);
    neighbourAtFace(td1,n2,t21);
    neighbourAtFace(td1,n3,t31);
    neighbourAtFace(td1,n4,t41);

  
    td2->p1 = p5;  td2->p2 = p2;  td2->p3 = p7;  td2->p4 = p6;
    SetTdPBey(td2,td);
    SetTdE(td2,t12,it2,t32,t42);
    neighbourAtFace(td2,n1,t12);
    neighbourAtFace(td2,n3,t32);
    neighbourAtFace(td2,n4,t42);

    td3->p1 = p10; td3->p2 = p7;  td3->p3 = p3;  td3->p4 = p9;
    SetTdPBey(td3,td);
    SetTdE(td3,t13,t23,it3,t43);
    neighbourAtFace(td3,n1,t13);
    neighbourAtFace(td3,n2,t23);
    neighbourAtFace(td3,n4,t43);

    td4->p1 = p8;  td4->p2 = p6;  td4->p3 = p9;  td4->p4 = p4;
    SetTdPBey(td4,td);
    SetTdE(td4,t14,t24,t34,it4);
    neighbourAtFace(td4,n1,t14);
    neighbourAtFace(td4,n2,t24);
    neighbourAtFace(td4,n3,t34);

    //  tetrahedra generated by the inner octahedron  

    td5->p1 = p5;  td5->p2 = p10; td5->p3 = p8;  td5->p4 = p6;
    SetTdPBey(td5,td);
    SetTdE(td5,it6,t33,it5,it1);
    neighbourAtFace(td5,n3,t33);

    td6->p1 = p5;  td6->p2 = p10; td6->p3 = p7;  td6->p4 = p6;
    SetTdPBey(td6,td);
    SetTdE(td6,it7,it2,it5,t44);
    neighbourAtFace(td6,n4,t44);

    td7->p1 = p10; td7->p2 = p8;  td7->p3 = p6;  td7->p4 = p9;
    SetTdPBey(td7,td);
    SetTdE(td7,it4,it8,t22,it6);
    neighbourAtFace(td7,n2,t22);

    td8->p1 = p10; td8->p2 = p7;  td8->p3 = p6;  td8->p4 = p9;
    SetTdPBey(td8,td);
    SetTdE(td8,t11,it8,it3,it7);
    neighbourAtFace(td8,n1,t11);

    innerNeighbourAtFace(td1,it1,td5);
    innerNeighbourAtFace(td2,it2,td6);
    innerNeighbourAtFace(td3,it3,td8);
    innerNeighbourAtFace(td4,it4,td7);
    innerNeighbourAtFace(td5,it6,td7);
    innerNeighbourAtFace(td5,it5,td6);
    innerNeighbourAtFace(td5,it1,td1);
    innerNeighbourAtFace(td6,it7,td8);
    innerNeighbourAtFace(td6,it2,td2);
    innerNeighbourAtFace(td6,it5,td5);
    innerNeighbourAtFace(td7,it4,td4);
    innerNeighbourAtFace(td7,it8,td8);
    innerNeighbourAtFace(td7,it6,td5);
    innerNeighbourAtFace(td8,it8,td7);
    innerNeighbourAtFace(td8,it3,td3);
    innerNeighbourAtFace(td8,it7,td6);

    td1->type = T_WHITE;
    td2->type = T_WHITE;
    td3->type = T_WHITE;
    td4->type = T_WHITE;
    td5->type = T_WHITE;
    td6->type = T_WHITE;
    td7->type = T_WHITE;
    td8->type = T_WHITE;
 
    td->firstSon = td1;
    td->type = T_RED;
 
   tetList[tetList.h]->add(td1);
   tetList[tetList.h]->add(td2);
   tetList[tetList.h]->add(td3);
   tetList[tetList.h]->add(td4);
   tetList[tetList.h]->add(td5);
   tetList[tetList.h]->add(td6);
   tetList[tetList.h]->add(td7);
   tetList[tetList.h]->add(td8);
 
    // destruction of the faces
    delete it1;
    delete it2;
    delete it3;
    delete it4;
    delete it5;
    delete it6;
    delete it7;
    delete it8;

    if ( !(t1->onBoundary()) ) {
       delete t11; 
       delete t12;
       delete t13;
       delete t14;
    }

    if ( !(t2->onBoundary()) ) {
       delete t21;
       delete t22;
       delete t23;
       delete t24;
    }

    if ( !(t3->onBoundary()) ) {
       delete t31;
       delete t32;
       delete t33;
       delete t34;
    }

    if ( !(t4->onBoundary()) ) {
       delete t41;
       delete t42;
       delete t43;
       delete t44;
    }

    td->returnTriangles(trs);

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

EDG3* MESH3:: findEdgInNeighbour(PT3* p1, PT3* p2, TET3* td)
{
	TET3 *tdSon;
	EDG3 *e;
 	int iEnd;

	tdSon = td->firstSon;			
	if (tdSon == nil)  { 
            cout << "\n\n  error in findEdgInNeighbour: tdSon == nil\n" ;
	return 0; 
	}

	iEnd = 8;
	if (td->type == T_GREEN_2B) iEnd = 3;

	for (int i=1; i<=iEnd; i++)
        {
		e = EdinTet(p1,p2,tdSon);
       		if (e != 0) return e;
            	else tdSon = tdSon->next;
	}
		   
	return 0;
}
//-------------------------------------------------------------------------

void MESH3:: setTypeOfInnerEdges(TR3 *t)
{
	TR3 *tSon4;

	tSon4 = t->firstSon->next->next->next;
	tSon4->e1->type = T_RED;
	tSon4->e2->type = T_RED;
	tSon4->e3->type = T_RED;

	return;
}

int MESH3:: RefineTr(TR3 *t)
{
   TR3   *t1, *t2, *t3, *t4;
   EDG3  *e1 = t->e1,  *e2 = t->e2,  *e3 = t->e3, *ereg,
         *ie1,*ie2,*ie3, *e11, *e12, *e21, *e22, *e31, *e32;
   TET3  *nt;
   int depth = (t->depth)+1;
   Bool  nbRefined, nbFits;

   // Get 4 new triangles, 3 new edges, return False if not enough space

   t1 = new TR3;
   t2 = new TR3;
   t3 = new TR3;
   t4 = new TR3;

   // 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 = RefineEdg(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 = RefineEdg(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 = RefineEdg(e3);
	if (e31==0) return False;
    }
    else e31 = e3->firstSon;
    e32 = e31->next;
    if ((e31->p1)==(t->p2))
      { ereg = e31; e31 = e32; e32 = ereg; }



   ie1 = edgAlloc.Get();
   ie2 = edgAlloc.Get();
   ie3 = edgAlloc.Get();
   if ((t4==0)||(ie3==0)) return False;

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

   ie1->boundP = t->boundP;
   ie2->boundP = t->boundP;
   ie3->boundP = t->boundP;

   ie1->classA = t->classA;
   ie2->classA = t->classA;
   ie3->classA = t->classA;
   
   ie1->isoP   = t->isoP;
   ie2->isoP   = t->isoP;
   ie3->isoP   = t->isoP;
   


   // Build the new triangles.

   InnerEdges(t,e1,e2,e3,ie1,ie2,ie3);

   nt = t->t32;  // looking for the second neighbour

   nbFits = False; nbRefined = False;
   if ( nt != 0 )
   {
	if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
	   nbFits = True;
     	if ( nt->refined() ) nbRefined=True;
   }

   if ( nbFits && nbRefined )
   {
	EDG3 *ne1, *ne2, *ne3;
	ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
	ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
	ne3 = findEdgInNeighbour(ie3->p1,ie3->p2,nt);
	edgAlloc.Return(ie1);
	edgAlloc.Return(ie2);
	edgAlloc.Return(ie3);
        ie1=ne1; ie2=ne2; ie3=ne3;
   }
   else
   {
	noOfEdges += 3;
   	noOfTriangles += 3;

	if (depth > edgList.h)  // create new list at depth if necessary
    	{
	   DList<EDG3>* edgL = new DList<EDG3>;
	   edgList.push(edgL);
    	}
  	edgList[depth]->add(ie1); 	
	edgList[depth]->add(ie2); 	
	edgList[depth]->add(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); 
   SetTP(t2,t); 
   SetTP(t3,t); 
   SetTP(t4,t); 
   t->firstSon = t1;
   t->type = T_RED;
   t1->next = t2; t2->next = t3, t3->next = t4;
   t2->prev = t1; t3->prev = t2, t4->prev = t3;

   if (t->boundP != INTERIOR)
   {
	if (depth > trList.h)  // create new list at depth if necessary
        {
	  DList<TR3>* trL = new DList<TR3>;
	  trList.push(trL);
        }

	trList[depth]->add(t1);
	trList[depth]->add(t2);
	trList[depth]->add(t3);
	trList[depth]->add(t4);
	// noOfTriangles += 3;
   }

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


EDG3* MESH3:: RefineEdg(EDG3 *ed)
{
   EDG3 *e1, *e2;
   PT3  *pm = 0;
   int depth;

   //     Refinement of an edge.
   // --  Get two new edges; return False if not enough space

   e1 = edgAlloc.Get();
   e2 = edgAlloc.Get();
   
   if (e2==0) return 0;
   noOfEdges++;

   pm = ptAlloc.Get();
   pm->reset();   
   noOfPoints++;


   // Call the user midpoint routine and set the new fields

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

    if ( (ed->isoP == CIRCULAR) && (ed->boundP != INTERIOR) ) 
        setXYZpmCircle(pm, ed);
    else
	if ( (ed->boundP != INTERIOR) && (ed->isoP == CYLINDRICAL) )
	   setXYZpmCyl(pm, ed);
	else
	{
	    pm->x = HALF*(ed->p1->x + ed->p2->x);
	    pm->y = HALF*(ed->p1->y + ed->p2->y);
	    pm->z = HALF*(ed->p1->z + ed->p2->z);
        }

   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;
   ed->firstSon = e1; ed->pm =pm;
   e1->depth = e2->depth = ed->depth+1;

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

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

   return e1;
}
//-------------------------------------------------------------------------


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

   // set the points of a triangle, computed from the edges

   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; 
       t->isoP   = f->isoP;
       t->boundP = f->boundP;
     }
   else 
       t->boundP = INTERIOR;

   if ( (t->depth) > maxDepth )
      maxDepth = t->depth;
   return;
}
//-------------------------------------------------------------------------
void MESH3:: SetTP(TET3 *tet, TR3 *t, int noP1, int noP2, int noP3)
{
    PT3 *tetPts[1+4];
    tet->getPoints(tetPts);

    EDG3 *e1 =  t->e1, *e2 = t->e2,  *e3 = t->e3;
    PT3 *P11 = e1->p1, *P12 = e1->p2,
        *P21 = e2->p1, *P22 = e2->p2,
        *P31 = e3->p1, *P32 = e3->p2;



    t->p1 = tetPts[noP1];
    t->p2 = tetPts[noP2];
    t->p3 = tetPts[noP3];

    if(      (P11!=t->p1) && (P12!=t->p1) ) t->e1 = e1;
    else if( (P21!=t->p1) && (P22!=t->p1) ) t->e1 = e2;
    else if( (P31!=t->p1) && (P32!=t->p1) ) t->e1 = e3;

    if(      (P11!=t->p2) && (P12!=t->p2) ) t->e2 = e1;
    else if( (P21!=t->p2) && (P22!=t->p2) ) t->e2 = e2;
    else if( (P31!=t->p2) && (P32!=t->p2) ) t->e2 = e3;

    if(      (P11!=t->p3) && (P12!=t->p3) ) t->e3 = e1;
    else if( (P21!=t->p3) && (P22!=t->p3) ) t->e3 = e2;
    else if( (P31!=t->p3) && (P32!=t->p3) ) t->e3 = e3;

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

 
void TET3:: completeT(TR3 *t,TR3 *ft) const
{
/*
   PT3 *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 = ft;
   if (ft!=0)
     {
       t->depth  = (ft->depth)+1; 
       t->classA =  ft->classA; 
       t->isoP   =  ft->isoP;
       t->boundP = ft->boundP;
     }
   else 
       t->boundP = INTERIOR;

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


int MESH3:: TdVolume(TET3 *tet)
{
   Real  x1, x2, x3, y1, y2, y3, z1, z2, z3;
   Real  det;
   PT3   *P1, *P2, *P3, *P4;

   //   computes the Jacobian and the orientation of the edges 

   P1 = tet->p1;
   P2 = tet->p2;
   P3 = tet->p3;
   P4 = tet->p4;

   x1 = (P2->x) - (P1->x);
   x2 = (P2->y) - (P1->y);
   x3 = (P2->z) - (P1->z);

   y1 = (P3->x) - (P1->x);
   y2 = (P3->y) - (P1->y);
   y3 = (P3->z) - (P1->z);

   z1 = (P4->x) - (P1->x);
   z2 = (P4->y) - (P1->y);
   z3 = (P4->z) - (P1->z);

   det = x1*(y2*z3 - y3*z2) - x2*(y1*z3 -y3*z1) + x3*(y1*z2 - y2*z1);
   tet->detJ = det;	// the volume of the td: fabs(det)*SIXTH ;

   return tet->orient();
}
//-------------------------------------------------------------------------


void MESH3:: SetTdP(TET3 *td, TR3* t1, TR3* t2, TR3* t3, TR3* t4, TET3 *f)
{
    PT3  *P11 = t1->p1, *P12 = t1->p2, *P13 = t1->p3,
         *P21 = t2->p1, *P22 = t2->p2, *P23 = t2->p3,
         *P31 = t3->p1, *P32 = t3->p2, *P33 = t3->p3,
         *P41 = t4->p1, *P42 = t4->p2, *P43 = t4->p3;
    PT3  *P2, *P3;

    //   set the points of a tetrahedron td, computed from the faces;
    //   td inherits properties of its father f

    if ((P31!=P11)&&(P31!=P12)&&(P31!=P13))
              td->p1 = P31;
    else if ((P32!=P11)&&(P32!=P12)&&(P32!=P13))
              td->p1 = P32;
         else td->p1 = P33;

    if ((P41!=P21)&&(P41!=P22)&&(P41!=P23))
              td->p2 = P41;
    else if ((P42!=P21)&&(P42!=P22)&&(P42!=P23))
              td->p2 = P42;
         else td->p2 = P43;

    if ((P11!=P31)&&(P11!=P32)&&(P11!=P33))
              td->p3 = P11;
    else if ((P12!=P31)&&(P12!=P32)&&(P12!=P33))
              td->p3 = P12;
         else td->p3 = P13;

    if ((P21!=P41)&&(P21!=P42)&&(P21!=P43))
              td->p4 = P21;
    else if ((P22!=P41)&&(P22!=P42)&&(P22!=P43))
              td->p4 = P22;
         else td->p4 = P23;

    TdVolume(td);

    if (rightHand)
    {
       P2 = td->p2;
       P3 = td->p3;

       if (!(td->orient())) 
       {
            td->p2 = P3;
            td->p3 = P2;
            td->detJ = Abs(td->detJ);  // td->orient = True;
	}
   }
    
    td->n1 = nil; td->n2 = nil; td->n3 = nil; td->n4 = nil;

    td->father = f;
    if (f!=0)
    {
           td->depth =  f->depth +1; 
           td->classA = f->classA;
       }
    if (td->depth > maxDepth)  maxDepth = td->depth;

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


void MESH3:: SetTdPBey(TET3 *td,TET3 *f)
{
   //   set the points of a tetrahedron td, computed from the faces;
   //   td inherits properties of its father f

   TdVolume(td);

   td->n1 = nil; td->n2 = nil; td->n3 = nil; td->n4 = nil;

   td->father = f;
   if (f!=0)
      {
        td->depth = (f->depth)+1; 
        td->classA = f->classA;
      }
   if ( (td->depth)> maxDepth )
      maxDepth = td->depth;

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


void MESH3:: SetTPts(TET3 *td, TET3 *f)
{
   PT3  *P2, *P3;

   //   computes volume and orientation in the tetrahedron td;
   //   td inherits properties of its father f

   TdVolume(td);	// determinant

   if (rightHand)
   {
       P2 = td->p2;
       P3 = td->p3;
       
       if (!(td->orient())) 
       {
	   td->p2 = P3;
	   td->p3 = P2;
	   td->detJ = Abs(td->detJ);	// td->orient() = True;
       }
   }
   
   td->father = f;
   if (f != 0)
   { 
       td->depth = (f->depth)+1; 
       td->classA = f->classA; 
   }
   if ( (td->depth)> maxDepth )
   maxDepth = td->depth;
   return;
}
//-------------------------------------------------------------------------


EDG3* MESH3:: FindTdE(TET3 */*td*/,TR3 *t1,TR3 *t2,TR3 *t3,TR3 *t4,PT3 *p1,PT3 *p2)
{ 
   EDG3 *E11, *E12, *E13,
        *E21, *E22, *E23,
        *E31, *E32, *E33,
        *E41, *E42, *E43;

    E11 = t1->e1;
    if ( ((E11->p1==p1)&&(E11->p2==p2))||((E11->p1==p2)&&(E11->p2==p1)))
          return E11;

    E12 = t1->e2;
    if ( ((E12->p1==p1)&&(E12->p2==p2))||((E12->p1==p2)&&(E12->p2==p1)))
          return E12;

    E13 = t1->e3;
    if ( ((E13->p1==p1)&&(E13->p2==p2))||((E13->p1==p2)&&(E13->p2==p1)))
          return E13;

    E21 = t2->e1;
    if ( ((E21->p1==p1)&&(E21->p2==p2))||((E21->p1==p2)&&(E21->p2==p1)))
          return E21;

    E22 = t2->e2;
    if ( ((E22->p1==p1)&&(E22->p2==p2))||((E22->p1==p2)&&(E22->p2==p1)))
          return E22;

    E23 = t2->e3;
    if ( ((E23->p1==p1)&&(E23->p2==p2))||((E23->p1==p2)&&(E23->p2==p1)))
          return E23;

    E31 = t3->e1;
    if ( ((E31->p1==p1)&&(E31->p2==p2))||((E31->p1==p2)&&(E31->p2==p1)))
          return E31;

    E32 = t3->e2;
    if ( ((E32->p1==p1)&&(E32->p2==p2))||((E32->p1==p2)&&(E32->p2==p1)))
          return E32;

    E33 = t3->e3;
    if ( ((E33->p1==p1)&&(E33->p2==p2))||((E33->p1==p2)&&(E33->p2==p1)))
          return E33;

    E41 = t4->e1;
    if ( ((E41->p1==p1)&&(E41->p2==p2))||((E41->p1==p2)&&(E41->p2==p1)))
          return E41;

    E42 = t4->e2;
    if ( ((E42->p1==p1)&&(E42->p2==p2))||((E42->p1==p2)&&(E42->p2==p1)))
          return E42;

    E43 = t4->e3;
    if ( ((E43->p1==p1)&&(E43->p2==p2))||((E43->p1==p2)&&(E43->p2==p1)))
          return E43;

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


void MESH3:: SetTdE(TET3 *td,TR3 *t1,TR3 *t2,TR3 *t3,TR3 *t4)
{ 
   PT3  *P1 = td->p1, *P2 = td->p2, *P3 = td->p3, *P4 = td->p4;

   //   set the edges of a tetrahedron td

   td->e1 = FindTdE(td,t1,t2,t3,t4,P1,P3);
   td->e2 = FindTdE(td,t1,t2,t3,t4,P2,P3);
   td->e3 = FindTdE(td,t1,t2,t3,t4,P1,P2);
   td->e4 = FindTdE(td,t1,t2,t3,t4,P3,P4);
   td->e5 = FindTdE(td,t1,t2,t3,t4,P2,P4);
   td->e6 = FindTdE(td,t1,t2,t3,t4,P1,P4);

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


void MESH3:: InnerEdges(TR3 *t,EDG3 *e1,EDG3 *e2,EDG3 *e3,EDG3 *ie1,
			EDG3 *ie2,EDG3 *ie3)
{
   //   sets the points of the inner edges *ie1,*ie2,*ie3;
   //   Edges of the triangle t  to be refined: *e1,*e2,*e3

   ie1->p1 = e1->pm; 
   ie1->p2 = e3->pm; 
   ie1->type = T_GREEN_3; 
   ie1->classA = t->classA;	// !!!
   ie1->isoP   = t->isoP;
   ie1->boundP = t->boundP;

   ie2->p1 = e2->pm; 
   ie2->p2 = e1->pm; 
   ie2->type = T_GREEN_3; 
   ie2->classA = t->classA;
   ie2->isoP   = t->isoP;
   ie2->boundP = t->boundP;

   ie3->p1 = e3->pm; 
   ie3->p2 = e2->pm; 
   ie3->type = T_GREEN_3; 
   ie3->classA = t->classA;
   ie3->isoP   = t->isoP;
   ie3->boundP = t->boundP;

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


int MESH3:: MkGreen(TR3 *t)
{
   PT3   *p1, *p2;
   EDG3  *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew=0;
   TR3   *t1New=0, *t2New=0;
   TET3  *nt;
   int   depth;
   Bool  nbRefined, nbFits;

   // generates two green triangles as sons of t 
 
   if (t->firstSon != 0) return True;

   //   Testing if one neighbor triangle is refined and setting pi, ei
   //   corresponding to this neighbor


   if ( (e1->type == T_GREEN_1) || (e2->type == T_GREEN_1) ||
        (e3->type == T_GREEN_1) )
       cout << "Alarm in MkGreen: one edge already green \n";
   
    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 */
          }
      }

   //   Get one new edge and two new triangles

    eNew = edgAlloc.Get();

    t1New = new TR3;
    t2New = new TR3;
   

   //   Set the fields for the new objects

   eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;

   nt = t->t32;  // looking for the second neighbour

   nbFits = False; nbRefined = False;
   if ( nt != 0 )
   {
      if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
	 nbFits = True;

      if ( nt->refined() ) nbRefined=True;
   }

   if ( nbFits && nbRefined )
   {
	EDG3 *ne1;
	ne1 = findEdgInNeighbour(eNew->p1,eNew->p2,nt);
	edgAlloc.Return(eNew);
        eNew=ne1; 
   }
   else
   {
	noOfEdges++;
        noOfTriangles++;

  	eNew->type   = T_GREEN_1; 
  	eNew->classA = t->classA;
        eNew->isoP   = t->isoP;
  	eNew->boundP = t->boundP;
        depth        = t->depth+1;
  	eNew->depth  = depth;

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

   t1New->e1 = e3; t1New->e3 = eNew; t1New->type = T_GREEN_1;
   t2New->e1 = e2; t2New->e2 = eNew; t2New->type = T_GREEN_1;
   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;
      }

   SetTP(t1New,t); SetTP(t2New,t);
   t->firstSon = t1New;
   t->type = T_GREEN_1;

   // Append two new ones at their depth list
    
   depth = t1New->depth;
   
   t1New->next = t2New; t2New->prev = t1New;

   if (t->boundP != INTERIOR)
   {
	if ( depth > trList.h )  // create new list at depth if necessary
    	{
	  DList<TR3>* trL = new DList<TR3>;
	  trList.push(trL);
    	}
	trList[depth]->add(t1New);
	trList[depth]->add(t2New);
	// noOfTriangles++;
   }

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


int MESH3:: MkGreenIrr(TR3 *t)
{
   TR3   *t1=0, *t2=0, *t3=0;
   EDG3  *e1 = t->e1,  *e2 = t->e2,  *e3 = t->e3, *ereg,
         *ie1=0,*ie2=0, *e11=0, *e12=0, *e21=0, *e22=0, *e31=0, *e32=0;
   TET3  *nt;
   int   depth;
   Bool  nbFits, nbRefined, edgeChanged=False;

   //cout << "\n MkGreenIrr \n";
   //    generates three green triangles as sons of t 

   if (t->firstSon != 0) return True;

   // Get 3 new triangles, 2 new edges

    t1 = new TR3;
    t2 = new TR3;
    t3 = new TR3;

    ie1 = edgAlloc.Get();
    ie2 = edgAlloc.Get();
    

    ie1->depth = t->depth+1;
    ie2->depth = t->depth+1;
    ie1->next = ie2;
    ie2->prev = ie1;

   //  boundary value on the new triangles and edges  

   ie1->boundP   = t->boundP;
   ie1->type     = T_GREEN_2B;
   ie1->classA   = t->classA;
   ie1->isoP     = t->isoP;

   ie2->boundP   = t->boundP;
   ie2->type     = T_GREEN_2B;
   ie2->classA   = t->classA;
   ie2->isoP     = t->isoP;


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

   nt = t->t32;  // looking for the second neighbour

   nbFits = False; nbRefined = False;
   if ( nt != 0 )
   {
	if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
	   nbFits = True;
	if ( nt->refined() ) nbRefined=True;
   }

   if ((e1->firstSon)==0)         /* Edge e1 not refined */
     {
	ie1->p1 = e3->pm;
	ie1->p2 = e2->pm;
	ie2->p1 = e3->pm;
	ie2->p2 =  t->p3;

	if ( nbFits && nbRefined )
   	{
	   EDG3 *ne1, *ne2;
	   ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
	   ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
           if ( ne2 == 0 ) 
	   {
		ie2->p1 =  e2->pm;
		ie2->p2 =  t->p2;
	        ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
		edgeChanged = True;
   	   }
	   edgAlloc.Return(ie1); edgAlloc.Return(ie2);
           ie1=ne1; ie2=ne2;
   	}
  	else
   	{
	   noOfEdges += 2;
           noOfTriangles += 2;
           depth = ie1->depth;
	   if (depth > edgList.h)  // create new list at depth if necessary
    	      {
	         DList<EDG3>* edgL = new DList<EDG3>;
	         edgList.push(edgL);
    	      }
           edgList[depth]->add(ie1); 	
           edgList[depth]->add(ie2); 	
   	}

	e21 = e2->firstSon;
	e22 = e21->next;
	if ( (e21->p1) == (t->p3) )
           {ereg=e21; e21 = e22; e22 = ereg;}
	e31 = e3->firstSon;
	e32 = e31->next;
	if ( (e31->p1) == (t->p2) )
           {ereg=e31; e31 = e32; e32 = ereg;}

	if (edgeChanged) { EDG3* e = e22; e22 = e32; e32 = e;}
	t1->e1 = e31; t1->e2 = ie1; t1->e3 = e21;
	t2->e1 = e1;  t2->e2 = ie2; t2->e3 = e32;
	t3->e1 = ie1; t3->e2 = ie2; t3->e3 = e22;
     }

   else
   if ((e2->firstSon)==0)         /* Edge e2 not refined */
     {
	ie1->p1 = e1->pm;
	ie1->p2 = e3->pm;
	ie2->p1 = e1->pm;
	ie2->p2 =  t->p1;

	if ( nbFits && nbRefined )
   	{
	   EDG3 *ne1, *ne2;
	   ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
	   ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
           if ( ne2 == 0 ) 
	   {
		ie2->p1 =  e3->pm;
		ie2->p2 =  t->p3;
	        ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
		edgeChanged = True;
   	   }
	   edgAlloc.Return(ie1); edgAlloc.Return(ie2);
           ie1=ne1; ie2=ne2;
   	}
  	else
   	{
	   noOfEdges += 2;
           noOfTriangles += 2;

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

	e11 = e1->firstSon;
	e12 = e11->next;
	if ( (e11->p1) == (t->p3) )
           {ereg=e11; e11 = e12; e12 = ereg;}
	e31 = e3->firstSon;
	e32 = e31->next;
	if ( (e31->p1) == (t->p2) )
           {ereg=e31; e31 = e32; e32 = ereg;}

	if (edgeChanged) { EDG3* e = e12; e12 = e31; e31 = e;}
	t1->e1 = e32; t1->e2 = e11; t1->e3 = ie1;
	t2->e1 = e12; t2->e2 = e2;  t2->e3 = ie2;
	t3->e1 = e31; t3->e2 = ie1; t3->e3 = ie2;
     }

   else
   if ((e3->firstSon)==0)         /* Edge e3 not refined */
     {
	ie1->p1 = e2->pm;
	ie1->p2 = e1->pm;
	ie2->p1 = e2->pm;
	ie2->p2 =  t->p2;

	if ( nbFits && nbRefined )
   	{
	   EDG3 *ne1, *ne2;
	   ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
	   ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
           if ( ne2 == 0 ) 
	   {
		ie2->p1 =  e1->pm;
		ie2->p2 =  t->p1;
	        ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
		edgeChanged = True;
   	   }
	   edgAlloc.Return(ie1); edgAlloc.Return(ie2);
           ie1=ne1; ie2=ne2;
   	}
  	else
   	{
	   noOfEdges += 2;
           noOfTriangles += 2;

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

	e11 = e1->firstSon;
	e12 = e11->next;
	if ( (e11->p1) == (t->p3) )
           {ereg=e11; e11 = e12; e12 = ereg;}
	e21 = e2->firstSon;
	e22 = e21->next;
	if ( (e21->p1) == (t->p3) )
           {ereg=e21; e21 = e22; e22 = ereg;}

	if (edgeChanged) { EDG3* e = e11; e11 = e21; e21 = e;}
	t1->e1 = e12; t1->e2 = e22; t1->e3 = ie1;
	t2->e1 = ie2; t2->e2 = e21; t2->e3 = e3;
	t3->e1 = e11; t3->e2 = ie1; t3->e3 = ie2;
     }


   // Build the new triangles.
     
    t1->type = T_GREEN_2B; SetTP(t1,t); 
    t2->type = T_GREEN_2B; SetTP(t2,t); 
    t3->type = T_GREEN_2B; SetTP(t3,t); 
    t->firstSon = t1;
    t->type = T_GREEN_2B;



   // Append two new ones at their depth list
   
   depth = t1->depth;
   t1->next = t2; t2->next = t3;
   t2->prev = t1; t3->prev = t2;

   if (t->boundP != INTERIOR)
   {
	if (depth > trList.h)  		// create new list at depth if necessary
    	{
	  DList<TR3>* trL = new DList<TR3>;
	  trList.push(trL);
    	}
   	trList[depth]->add(t1);
   	trList[depth]->add(t2);
   	trList[depth]->add(t3);
	// noOfTriangles += 2;
   }

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

int MESH3:: MarkGreenTd(TET3 *td)
{
   EDG3 *eds[1+6], *ed;
   TR3  *trs[1+4];
   int i, noOfRefEdg = 0;
 
   // marks green tetrahedra to be refined

   if ( (td->type) != T_WHITE )
     cout << " Alarm in MarkGreenTd: td not white \n";  

   td->getEdges(eds);
   for (i=1; i<=6; i++)
      {
         ed = eds[i];
         if ( (ed->firstSon) != 0 ) noOfRefEdg++;
      }

   if ( noOfRefEdg == 1)
     {
       td->mark = T_GREEN_1;
       moreGreen1Td++;
     }
   else  
   if ( noOfRefEdg == 2)
     {
       td->mark = T_GREEN_2A;
       moreGreen2ATd++;

       td->getTriangles(trs);
       for ( i=1; i <= 4; i++)
          {
            if ( TestFace(trs[i]) == 2 ) 
              {
                td->mark = T_GREEN_2B;
                moreGreen2ATd--;
                moreGreen2BTd++;
              }
          }
       td->returnTriangles(trs);
     }
   else  
   if ( noOfRefEdg == 3)
     {
       td->mark = T_GREEN_3;
       moreGreen3Td++;
     }
   else td->mark = False;

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


int MESH3:: MkGreen1Td(TET3 *td)
  {
    PT3   *p1=0, *p2=0, *p3=0, *p4=0, *pm=0;
    TR3   *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
    TR3   *t1New=0, *t21=0, *t22=0, *t41=0, *t42=0;
    TR3   *trs[1+4];
    TET3  *td1,*td2;
    PATCH *n1=0, *n2=0, *n3=0, *n4=0;

    EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
    EDG3 *ed, *e1New=0, *e2New=0, *e3New=0;
    int  refEdges=0, depth = (td->depth)+1;

    // refinement of td into two green tetrahedra

    if (td->mark == T_GREEN_1) td->mark = False;
    else return True;


    td->getTriangles(trs);

      // Get  one new triangle and two tetrahedra, 
      // if not enough space: return False

   
    t1New  = new TR3;
    noOfTriangles++;

    t1New->depth = depth;
   
    td1 = getTET3();
    td2 = getTET3();
    noOfTetrahedra++;

   // Testing which edge is refined
    
    ed = e1;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
        pm = ed->pm; e3New = e5;
        t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
        n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
      }

    ed = e2;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        if (refEdges > 1 )
           {
             cout << "\n\n alarm in MkGreenI \n";
             return False;
           }
        p1=td->p3; p2=td->p1; p3=td->p2; p4=td->p4;
        pm = ed->pm; e3New = e6;
        t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
        n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
      }

    ed = e3;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        if (refEdges > 1 )
           {
             cout << "\n\n alarm in MkGreenI \n";
             return False;
           }
        p1=td->p2; p2=td->p3; p3=td->p1; p4=td->p4;
        pm = ed->pm; e3New = e4;
        t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
        n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
      }

    ed = e4;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        if (refEdges > 1 )
           {
             cout << "\n\n alarm in MkGreenI \n";
             return False;
           }
        p1=td->p4; p2=td->p1; p3=td->p3; p4=td->p2;
        pm = ed->pm; e3New = e3;
        t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
        n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
      }

    ed = e5;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        if (refEdges > 1 )
           {
             cout << "\n\n alarm in MkGreenI \n";
             return False;
           }
        p1=td->p2; p2=td->p3; p3=td->p4; p4=td->p1;
        pm = ed->pm; e3New = e1;
        t1 = trs[2]; t2 = trs[3]; t3 = trs[4]; t4 = trs[1];
        n1 = td->n2; n2 = td->n3; n3 = td->n4; n4 = td->n1;
      }

    ed = e6;
    if ( (ed->firstSon)!=0 )
      {
        refEdges +=1;
        if (refEdges > 1 )
           {
             cout << "\n\n alarm in MkGreenI \n";
             return False;
           }
        p1=td->p1; p2=td->p3; p3=td->p4; p4=td->p2;
        pm = ed->pm; e3New = e2;
        t1 = trs[1]; t2 = trs[3]; t3 = trs[4]; t4 = trs[2];
        n1 = td->n1; n2 = td->n3; n3 = td->n4; n4 = td->n2;
      }


      // Set the fields for the new objects

      MkGreen(t2);
      t = t2->firstSon;
      if ( PinTr(p1,t) ) 
         { t21 = t; t22 = t->next; }
       else
         { t22 = t; t21 = t->next; }

       e1New = EdinTr(t,p4,pm);


      MkGreen(t4);
      t = t4->firstSon;
      if ( PinTr(p1,t) ) 
         { t41 = t; t42 = t->next; }
       else
         { t42 = t; t41 = t->next; }

       e2New = EdinTr(t,p2,pm);


       t1New->e1 = e1New; t1New->e2 = e2New; t1New->e3 = e3New;
       SetTP(t1New,0); t1New->type = T_GREEN_3;


       SetTdP(td1,t21,t3,t41,t1New,td);
       SetTdE(td1,t21,t3,t41,t1New);
       neighbourAtFace(td1,n2,t21);
       neighbourAtFace(td1,n3,t3);
       neighbourAtFace(td1,n4,t41);
       td1->type = T_GREEN_1;
       
       SetTdP(td2,t22,t1,t42,t1New,td);
       SetTdE(td2,t22,t1,t42,t1New);
       neighbourAtFace(td2,n2,t22);
       neighbourAtFace(td2,n1,t1);
       neighbourAtFace(td2,n4,t42);
       td2->type = T_GREEN_1; 

       innerNeighbourAtFace(td1,t1New,td2);
       innerNeighbourAtFace(td2,t1New,td1);


       td->firstSon = td1;
       td->type = T_GREEN_1;
 

       tetList[tetList.h]->add(td1);
       tetList[tetList.h]->add(td2);


       // destruction of the faces
       delete t1New;

       if ( !(t2->onBoundary()) ) {
          delete t21; 
          delete t22;
       }

       if ( !(t4->onBoundary()) ) {
          delete t41; 
          delete t42;
       }

       td->returnTriangles(trs);


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


int MESH3:: MkGreen2ATd(TET3 *td)
  {
    PT3  *p1=0, *p2=0, *p3=0, *p4=0, *pm1=0, *pm2=0;
    TR3  *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
    TR3  *t1New=0, *t2New=0, *t3New=0, *t4New=0;
    TR3  *t11=0, *t12=0, *t21=0, *t22=0, *t31=0, *t32=0, *t41=0, *t42=0;
    TR3   *trs[1+4];
    TET3 *td1=0, *td2=0, *td3=0, *td4=0;
    PATCH *n1=0, *n2=0, *n3=0, *n4=0;
    EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
    EDG3 *ed, *edNew=0; 
    EDG3 *e1New, *e2New, *e3New, *e4New, *e5New, *e6New, *e7New, *e8New;
    int refEdges, depth = (td->depth)+1;

    // Refinement of td into four green tetrahedra

    if ( (td->mark == T_GREEN_2A) ) td->mark = False;
    else return True;

    td->getTriangles(trs);


    // Get  one new edge, four new triangles and four tetrahedra, 
    // if not enough space: return False
      

    edNew = edgAlloc.Get();
    noOfEdges += 1;

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


   t1New  = new TR3;
   t2New  = new TR3;
   t3New  = new TR3;
   t4New  = new TR3;
   noOfTriangles += 4;


   td1 = getTET3();
   td2 = getTET3();
   td3 = getTET3();
   td4 = getTET3();
   noOfTetrahedra += 3;


    //   Test which edge is refined
    
    refEdges = 0;

    ed = e1;
    if ( (ed->firstSon)!=0 )
      {
        refEdges++;
        p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
        pm1 = e1->pm; pm2 = e5->pm;
        t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
        n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
      }

    ed = e2;
    if ( (ed->firstSon)!=0 )
      {
        refEdges++;
        if ( refEdges > 1) cout << "Alarm in MkGreen2AT3 \n";
        p1=td->p3; p2=td->p1; p3=td->p2; p4=td->p4;
        pm1 = e2->pm; pm2 = e6->pm;
        t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
        n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
      }

    ed = e3;
    if ( (ed->firstSon)!=0 )
      {
        refEdges++;
        if ( refEdges > 1) cout << "Alarm in MkGreen2AT3 \n";
        p1=td->p2; p2=td->p3; p3=td->p1; p4=td->p4;
        pm1 = e3->pm; pm2 = e4->pm;
        t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
        n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
      }



    //   Set the fields for the new objects

      MkGreen(t1);
      t = t1->firstSon;
      if ( PinTr(p4,t) ) 
         { t11 = t; t12 = t->next; }
       else
         { t12 = t; t11 = t->next; }

      MkGreen(t2);
      t = t2->firstSon;
      if ( PinTr(p1,t) ) 
         { t21 = t; t22 = t->next; }
       else
         { t22 = t; t21 = t->next; }


      MkGreen(t3);
      t = t3->firstSon;
      if ( PinTr(p2,t) ) 
         { t31 = t; t32 = t->next; }
       else
         { t32 = t; t31 = t->next; }



      MkGreen(t4);
      t = t4->firstSon;
      if ( PinTr(p3,t) ) 
         { t41 = t; t42 = t->next; }
       else
         { t42 = t; t41 = t->next; }


       edNew->p1 = pm1; edNew->p2 = pm2;
       edNew->pm = 0; 
       edNew->type = T_GREEN_2A;

       e1New = EdinTr(t12,p2,pm2);
       e2New = EdinTr(t22,p4,pm1);
       e3New = EdinTr(t11,p4,pm2);
       e4New = EdinTr(t41,p2,pm1);
       e5New = EdinTr(t11,p3,pm2);
       e6New = EdinTr(t22,p3,pm1);
       e7New = EdinTr(t32,p1,pm2);
       e8New = EdinTr(t21,p1,pm1);


       t1New->e1 = edNew; t1New->e2 = e2New; t1New->e3 = e3New;
       SetTP(t1New,0); t1New->type = T_GREEN_2A;

       t2New->e1 = e1New; t2New->e2 = e4New; t2New->e3 = edNew;
       SetTP(t2New,0); t2New->type = T_GREEN_2A;

       t3New->e1 = edNew; t3New->e2 = e5New; t3New->e3 = e6New;
       SetTP(t3New,0); t3New->type = T_GREEN_2A;

       t4New->e1 = edNew; t4New->e2 = e7New; t4New->e3 = e8New;
       SetTP(t4New,0); t4New->type = T_GREEN_2A;

       SetTdP(td1,t11,t22,t1New,t3New,td);
       SetTdE(td1,t11,t22,t1New,t3New);
       neighbourAtFace(td1,n1,t11);
       neighbourAtFace(td1,n2,t22);
       td1->type = T_GREEN_2A;


       SetTdP(td2,t12,t2New,t41,t3New,td);
       SetTdE(td2,t12,t2New,t41,t3New);
       neighbourAtFace(td2,n1,t12);
       neighbourAtFace(td2,n4,t41);
       td2->type = T_GREEN_2A; 

       SetTdP(td3,t32,t21,t1New,t4New,td);
       SetTdE(td3,t32,t21,t1New,t4New);
       neighbourAtFace(td3,n3,t32);
       neighbourAtFace(td3,n2,t21);
       td3->type = T_GREEN_2A; 

       SetTdP(td4,t31,t42,t2New,t4New,td);
       SetTdE(td4,t31,t42,t2New,t4New);
       neighbourAtFace(td4,n4,t42);
       neighbourAtFace(td4,n3,t31);
       td4->type = T_GREEN_2A; 

       innerNeighbourAtFace(td1,t1New,td3);
       innerNeighbourAtFace(td1,t3New,td2);
       innerNeighbourAtFace(td2,t2New,td4);
       innerNeighbourAtFace(td2,t3New,td1);
       innerNeighbourAtFace(td3,t1New,td1);
       innerNeighbourAtFace(td3,t4New,td4);
       innerNeighbourAtFace(td4,t2New,td2);
       innerNeighbourAtFace(td4,t4New,td3);


       td->firstSon = td1;
       td->type = T_GREEN_2A;
 
       tetList[tetList.h]->add(td1);
       tetList[tetList.h]->add(td2);
       tetList[tetList.h]->add(td3);
       tetList[tetList.h]->add(td4);
  

    // destruction of the faces
    delete t1New;
    delete t2New;
    delete t3New;
    delete t4New;

    if ( !(t1->onBoundary()) ) {
       delete t11; 
       delete t12;
    }

    if ( !(t2->onBoundary()) ) {
       delete t21; 
       delete t22;
    }

    if ( !(t3->onBoundary()) ) {
       delete t31; 
       delete t32;
    }

    if ( !(t4->onBoundary()) ) {
       delete t41; 
       delete t42;
    }

    td->returnTriangles(trs);


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


int MESH3:: MkGreen2BTd(TET3 *td)
{
   PT3   *p1=0, *p2=0, *p3=0, *p4=0;
   TR3   *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
   TR3   *trs[1+4];
   TR3   *t1New=0, *t2New=0;
   TR3   *t11=0, *t12=0, *t21=0, *t22=0, *t23=0, *t31=0, *t32=0;
   TET3  *td1=0, *td2=0, *td3=0;
   PATCH *n1=0, *n2=0, *n3=0, *n4=0;

   EDG3  *e1=0, *e2=0, *e3=0, *e4=0, *e5=0, *e6=0;
   EDG3  *e1New=0, *e2New=0, *e3New=0, *e4New=0;
   int   i, depth = (td->depth)+1;

   //   Refinement of td into three green tetrahedra

   if ( (td->mark == T_GREEN_2B) ) td->mark = False;
   else return True;
  
   td->getTriangles(trs);

  
   // Get two new triangles and three tetrahedra, 
   // if not enough space: return False
     
   
   t1New  = new TR3;
   t2New  = new TR3;
   noOfTriangles += 2;

   t1New->depth = t2New->depth = depth;
   
   td1 = getTET3();
   td2 = getTET3();
   td3 = getTET3();
   noOfTetrahedra += 2;


   
   // Testing which face is green-2B refined   
   // The actual case is transformed on the td with e4, e6 refined.
         

   for (i=1; i<=4; i++)
      {
        t = trs[i];
        if ( TestFace(t) > 1 )
          { 
            MkGreenIrr(t);

            if ( i==1 )
              {
                 if ( (td->e2)->firstSon == 0)
                   {
                     t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
                     n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
                     e1 = td->e2; e2 = td->e3; e3 = td->e1;
                     e4 = td->e5; e5 = td->e6; e6 = td->e4;
                     p1 = td->p3; p2 = td->p1; p3 = td->p2; p4 = td->p4;
                   }
                 else
                 if ( (td->e4)->firstSon == 0)
                   {
                     t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
                     n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
                     e1 = td->e4; e2 = td->e1; e3 = td->e6;
                     e4 = td->e2; e5 = td->e3; e6 = td->e5;
                     p1 = td->p4; p2 = td->p1; p3 = td->p3; p4 = td->p2;
                   }
                 else
                 if ( (td->e5)->firstSon == 0)
                   {
                     t1 = trs[2]; t2 = trs[1]; t3 = trs[4]; t4 = trs[3];
                     n1 = td->n2; n2 = td->n1; n3 = td->n4; n4 = td->n3;
                     e1 = td->e5; e2 = td->e6; e3 = td->e3;
                     e4 = td->e4; e5 = td->e1; e6 = td->e2;
                     p1 = td->p2; p2 = td->p1; p3 = td->p4; p4 = td->p3;
                   }
              }

            else
            if ( i==2 )
              {
                 if ( (td->e1)->firstSon == 0)
                   {
                     t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
                     n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
                     e1 = td->e1; e2 = td->e2; e3 = td->e3;
                     e4 = td->e4; e5 = td->e5; e6 = td->e6;
                     p1 = td->p1; p2 = td->p2; p3 = td->p3; p4 = td->p4;
                   }
                 else
                 if ( (td->e4)->firstSon == 0)
                   {
                     t1 = trs[3]; t2 = trs[2]; t3 = trs[4]; t4 = trs[1];
                     n1 = td->n3; n2 = td->n2; n3 = td->n4; n4 = td->n1;
                     e1 = td->e4; e2 = td->e5; e3 = td->e2;
                     e4 = td->e6; e5 = td->e3; e6 = td->e1;
                     p1 = td->p3; p2 = td->p2; p3 = td->p4; p4 = td->p1;
                   }
                 else
                 if ( (td->e6)->firstSon == 0)
                   {
                     t1 = trs[4]; t2 = trs[2]; t3 = trs[1]; t4 = trs[3];
                     n1 = td->n4; n2 = td->n2; n3 = td->n1; n4 = td->n3;
                     e1 = td->e6; e2 = td->e3; e3 = td->e5;
                     e4 = td->e1; e5 = td->e2; e6 = td->e4;
                     p1 = td->p4; p2 = td->p2; p3 = td->p1; p4 = td->p3;
                   }
              }

            else
            if ( i==3 )
              {
                 if ( (td->e3)->firstSon == 0)
                   {
                     t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
                     n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
                     e1 = td->e3; e2 = td->e1; e3 = td->e2;
                     e4 = td->e6; e5 = td->e4; e6 = td->e5;
                     p1 = td->p2; p2 = td->p3; p3 = td->p1; p4 = td->p4;
                   }
                 else
                 if ( (td->e5)->firstSon == 0)
                   {
                     t1 = trs[4]; t2 = trs[3]; t3 = trs[2]; t4 = trs[1];
                     n1 = td->n4; n2 = td->n3; n3 = td->n2; n4 = td->n1;
                     e1 = td->e5; e2 = td->e2; e3 = td->e4;
                     e4 = td->e3; e5 = td->e1; e6 = td->e6;
                     p1 = td->p4; p2 = td->p3; p3 = td->p2; p4 = td->p1;
                   }
                 else
                 if ( (td->e6)->firstSon == 0)
                   {
                     t1 = trs[1]; t2 = trs[3]; t3 = trs[4]; t4 = trs[2];
                     n1 = td->n1; n2 = td->n3; n3 = td->n4; n4 = td->n2;
                     e1 = td->e6; e2 = td->e4; e3 = td->e1;
                     e4 = td->e5; e5 = td->e2; e6 = td->e3;
                     p1 = td->p1; p2 = td->p3; p3 = td->p4; p4 = td->p2;
                   }
              }

            else
            if ( i==4 )
              {
                 if ( (td->e1)->firstSon == 0)
                   {
                     t1 = trs[3]; t2 = trs[4]; t3 = trs[1]; t4 = trs[2];
                     n1 = td->n3; n2 = td->n4; n3 = td->n1; n4 = td->n2;
                     e1 = td->e1; e2 = td->e6; e3 = td->e4;
                     e4 = td->e3; e5 = td->e5; e6 = td->e2;
                     p1 = td->p3; p2 = td->p4; p3 = td->p1; p4 = td->p2;
                   }
                 else
                 if ( (td->e2)->firstSon == 0)
                   {
                     t1 = trs[2]; t2 = trs[4]; t3 = trs[3]; t4 = trs[1];
                     n1 = td->n2; n2 = td->n4; n3 = td->n3; n4 = td->n1;
                     e1 = td->e2; e2 = td->e4; e3 = td->e5;
                     e4 = td->e1; e5 = td->e6; e6 = td->e3;
                     p1 = td->p2; p2 = td->p4; p3 = td->p3; p4 = td->p1;
                   }
                 else
                 if ( (td->e3)->firstSon == 0)
                   {
                     t1 = trs[1]; t2 = trs[4]; t3 = trs[2]; t4 = trs[3];
                     n1 = td->n1; n2 = td->n4; n3 = td->n2; n4 = td->n3;
                     e1 = td->e3; e2 = td->e5; e3 = td->e6;
                     e4 = td->e2; e5 = td->e4; e6 = td->e1;
                     p1 = td->p1; p2 = td->p4; p3 = td->p2; p4 = td->p3;
                   }
              }

            break;
          }
      }


   // After Transform we can work on the reference td with refined e4, e6 
    
   MkGreen(t1);
   t = t1->firstSon;
   if ( PinTr(p4,t) ) 
     { t11 = t; t12 = t->next; }
   else
     { t12 = t; t11 = t->next; }

   MkGreen(t3);
   t = t3->firstSon;
   if ( PinTr(p4,t) ) 
     { t31 = t; t32 = t->next; }
   else
     { t32 = t; t31 = t->next; }

   
   // subfaces of the grenn-IIB refined face  
   t21 = t2->firstSon;
   t23 = t21->next;
   t22 = t23->next;


   //   Set the fields for the new objects

   e1New = EdinTr(t12,p2,e4->pm);
   e3New = EdinTr(t32,p2,e6->pm);
   e2New = EdinTr(t21,e4->pm,e6->pm);

   t1New->e1 = e1New; t1New->e2 = e2New; t1New->e3 = e3New;
   SetTP(t1New,0); t1New->type = T_GREEN_2B;

   SetTdP(td1,t11,t21,t1New,t31,td);
   SetTdE(td1,t11,t21,t1New,t31);
   neighbourAtFace(td1,n1,t11);
   neighbourAtFace(td1,n2,t21);
   neighbourAtFace(td1,n3,t31);
   td1->type = T_GREEN_2B;


   if ( PinTr(e6->pm,t23) )
     {
       e4New = EdinTr(t23,p3,e6->pm);
       t2New->e1 = e2; t2New->e2 = e4New; t2New->e3 = e3New;
       SetTP(t2New,0); t2New->type = T_GREEN_2B;

       SetTdP(td2,t1New,t2New,t12,t22,td);
       SetTdE(td2,t1New,t2New,t12,t22);
       neighbourAtFace(td2,n1,t12);
       neighbourAtFace(td2,n2,t22);
       td2->type = T_GREEN_2B;
 
       SetTdP(td3,t23,t32,t2New,t4,td);
       SetTdE(td3,t23,t32,t2New,t4);
       neighbourAtFace(td3,n2,t23);
       neighbourAtFace(td3,n3,t32);
       neighbourAtFace(td3,n4,t4);
       td3->type = T_GREEN_2B; 
     }
   else
     {
       e4New = EdinTr(t23,p1,e4->pm);
       t2New->e1 = e1New; t2New->e2 = e4New; t2New->e3 = e3;
       SetTP(t2New,0); t2New->type = T_GREEN_2B;

       SetTdP(td2,t1New,t2New,t32,t22,td);
       SetTdE(td2,t1New,t2New,t32,t22);
       neighbourAtFace(td2,n3,t32);
       neighbourAtFace(td2,n2,t22);
       td2->type = T_GREEN_2B;
 
       SetTdP(td3,t23,t12,t2New,t4,td);
       SetTdE(td3,t23,t12,t2New,t4);
       neighbourAtFace(td3,n2,t23);
       neighbourAtFace(td3,n1,t12);
       neighbourAtFace(td3,n4,t4);
       td3->type = T_GREEN_2B; 
     }

     innerNeighbourAtFace(td1,t1New,td2);
     innerNeighbourAtFace(td2,t1New,td1);
     innerNeighbourAtFace(td2,t2New,td3);
     innerNeighbourAtFace(td3,t2New,td2);


   td->firstSon = td1;
   td->type = T_GREEN_2B;
 
   tetList[tetList.h]->add(td1);
   tetList[tetList.h]->add(td2);
   tetList[tetList.h]->add(td3);


    // destruction of the faces
    delete t1New;
    delete t2New;

    if ( !(t1->onBoundary()) ) {
       delete t11; 
       delete t12;
    }

    if ( !(t3->onBoundary()) ) {
       delete t31; 
       delete t32;
    }

    if ( !(t2->onBoundary()) ) {
       delete t21; 
       delete t22;
       delete t23;
    }


    td->returnTriangles(trs);


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


int MESH3:: MkGreen3Td(TET3 *td)
{
   PT3  *p1, *p2, *p3, *p4;
   PT3  *p6,*p7,*p9;
   TR3  *t, *t1, *t2, *t3, *t4, *trs[1+4];
   TR3  *s=0, *s1=0, *s2=0, *s3=0, *s4=0;
   TR3  *t1New=0, *t2New=0, *t3New=0, *t4New=0, *t5New=0, *t6New=0, *t7New=0, 
        *t8New=0, *t9New=0;
   TET3 *td1,*td2,*td3,*td4;
   PATCH *n1, *n2, *n3, *n4;
   EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
   EDG3 *e1New, *e2New,*e3New;
   EDG3 *e14, *e24, *e34;
   int i, depth = (td->depth)+1;

    // Refinement of td into four green tetrahedra

   if (td->mark == T_GREEN_3) td->mark = False;
   else return True;

   td->getTriangles(trs);
   for (i=1; i<=4; i++)
      {
        t = trs[i];
        if ( TestFace(t) > 1 )
          { 
            RefineTr(t);
            break;
          }
      }
   
   
   // Get  three new edges, nine new triangles and four tetrahedra, 
   // if not enough space: return False
   
   t1New  = new TR3;
   t2New  = new TR3;
   t3New  = new TR3;
   noOfTriangles += 3;

   t1New->depth = depth; t2New->depth = depth; t3New->depth = depth;

   
   td1 = getTET3();
   td2 = getTET3();
   td3 = getTET3();
   td4 = getTET3();

   noOfTetrahedra += 3;



   // Testing if one neighbor tetrahedron is refined and setting pi, ei
   // corresponding to this neighbor
    
    t = trs[1];
    if ( (t->firstSon)!=0 && (t->type == T_RED) )
      {
        p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
        p6 = e5->pm; p7 = e2->pm; p9 = e4->pm;
        t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
        n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
      }
    else
      {
        t = trs[2];
        if ( (t->firstSon)!=0 && (t->type == T_RED) )
          {
            p1=td->p2; p2=td->p4; p3=td->p3; p4=td->p1;
            p6 = e6->pm; p7 = e4->pm; p9 = e1->pm;
	    t1 = trs[2]; t2 = trs[4]; t3 = trs[3]; t4 = trs[1];
            n1 = td->n2; n2 = td->n4; n3 = td->n3; n4 = td->n1;
          }
        else
          {
            t = trs[3];
            if ( (t->firstSon)!=0 && (t->type == T_RED) )
              {
                p1=td->p3; p2=td->p4; p3=td->p1; p4=td->p2;
                p6 = e5->pm; p7 = e6->pm; p9 = e3->pm;
	        t1 = trs[3]; t2 = trs[4]; t3 = trs[1]; t4 = trs[2];
                n1 = td->n3; n2 = td->n4; n3 = td->n1; n4 = td->n2;
              }
            else 
              {
               t = trs[4];
               if ( (t->firstSon)!=0 && (t->type == T_RED) )
                 {
                   p1=td->p4; p2=td->p1; p3=td->p3; p4=td->p2;
                   p6 = e3->pm; p7 = e1->pm; p9 = e2->pm;
	           t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
                   n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
                 }
               else return False;                 
              }
          }
      }
 

   // Set the fields for the new objects


        MkGreen(t2);
        MkGreen(t3);
        MkGreen(t4);

        s = t1->firstSon;
        for (i=0;i<4;i++)
          {
           if (PinTr(p2,s)) s2 = s; 
           else if (PinTr(p3,s)) s1 = s;
                else if (PinTr(p4,s)) s3 = s;
                     else  {s4 = s;}
           s = s->next;
          }

        s = t2->firstSon;
        e1New = EdinTr(s,p1,p9);
        if (PinTr(p4,s)) 
          { t4New = s;
            t5New = s->next; }
        else
          { t5New = s;
            t4New = s->next; }

        s = t4->firstSon;
        e2New = EdinTr(s,p1,p7);
        if (PinTr(p2,s)) 
          { t7New = s;
            t6New = s->next; }
        else
          { t6New = s;
            t7New = s->next; }

        s = t3->firstSon;
        e3New = EdinTr(s,p1,p6);
        if (PinTr(p2,s)) 
          { t8New = s;
            t9New = s->next; }
        else
          { t9New = s;
            t8New = s->next; }


        e14 = EdinTr(s1,p7,p9);
        e24 = EdinTr(s2,p6,p7);
        e34 = EdinTr(s3,p6,p9);

        t1New->e1 = e1New; t1New->e2 = e14; t1New->e3 = e2New;
        t2New->e1 = e2New; t2New->e2 = e24; t2New->e3 = e3New;
        t3New->e1 = e3New; t3New->e2 = e34; t3New->e3 = e1New;
        SetTP(t1New,0); t1New->type = T_GREEN_3;
        SetTP(t2New,0); t2New->type = T_GREEN_3;
        SetTP(t3New,0); t3New->type = T_GREEN_3;

        SetTdP(td1,s1,t1New,t5New,t6New,td);
        SetTdE(td1,s1,t1New,t5New,t6New);
        neighbourAtFace(td1,n1,s1);
        neighbourAtFace(td1,n2,t5New);
        neighbourAtFace(td1,n4,t6New);
        td1->type = T_GREEN_3;

        SetTdP(td2,s2,t2New,t7New,t8New,td);
        SetTdE(td2,s2,t2New,t7New,t8New);
        neighbourAtFace(td2,n1,s2);
        neighbourAtFace(td2,n4,t7New);
        neighbourAtFace(td2,n3,t8New);
        td2->type = T_GREEN_3; 

        SetTdP(td3,s3,t3New,t9New,t4New,td);
        SetTdE(td3,s3,t3New,t9New,t4New);
        neighbourAtFace(td3,n1,s3);
        neighbourAtFace(td3,n3,t9New);
        neighbourAtFace(td3,n2,t4New);
        td3->type = T_GREEN_3; 

        SetTdP(td4,s4,t1New,t2New,t3New,td);
        SetTdE(td4,s4,t1New,t2New,t3New);
        neighbourAtFace(td4,n1,s4);
        td4->type = T_GREEN_3; 

        innerNeighbourAtFace(td1,t1New,td4);
        innerNeighbourAtFace(td2,t2New,td4);
        innerNeighbourAtFace(td3,t3New,td4);
        innerNeighbourAtFace(td4,t1New,td1);
        innerNeighbourAtFace(td4,t2New,td2);
        innerNeighbourAtFace(td4,t3New,td3);

        td->firstSon = td1;
        td->type = T_GREEN_3;
 
   tetList[tetList.h]->add(td1);
   tetList[tetList.h]->add(td2);
   tetList[tetList.h]->add(td3);
   tetList[tetList.h]->add(td4);


    // destruction of the faces
    delete t1New;
    delete t2New;
    delete t3New;

    if ( !(t1->onBoundary()) ) {
       delete s1; 
       delete s2;
       delete s3;
       delete s4;
    }

    if ( !(t2->onBoundary()) ) {
       delete t4New; 
       delete t5New;
    }

    if ( !(t3->onBoundary()) ) {
       delete t8New; 
       delete t9New;
    }

    if ( !(t4->onBoundary()) ) {
       delete t6New; 
       delete t7New;
    }

    td->returnTriangles(trs);



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



void MESH3:: RelGreen()
{
   TR3  *t1old, *t2old, *tFather;
   EDG3 *eMid = 0, *ed, *edNext;
   int  k, depth;

   //   remove all green triangles  (on the boundary)
  
   FORALL_DOWNWARD(trList,depth)
   {
       TR3 *t = trList[depth]->first;
       while (t!=0)  
       {   
	   k = 0;
	   tFather = t->father;
	   if (tFather==0)              	{ t = t->next; continue; }
	   if ((tFather->type)!=T_GREEN_1) 	{ 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_1) { eMid = t->e1; k++;}
	   if (((t->e2)->type)==T_GREEN_1) { eMid = t->e2; k++;}
	   if (((t->e3)->type)==T_GREEN_1) { eMid = t->e3; k++;}

	   if (k!=1)
	   {
	       cout << "\n*** RemoveGreen: Triangle with impossible state k="
	       << k << "\n";
	   }
	   t = t2old->next;
	   
	   //  Remove the two green triangles; return the space       
        
	   trList[t1old->depth]->remove(t1old);
	   trList[t2old->depth]->remove(t2old);

	   delete t1old;
	   delete t2old;
           noOfTriangles--;
	   
	   tFather->firstSon = 0;
	   tFather->type = T_WHITE;
	   tFather->mark = False;

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

	   noOfEdges--;
       }
   }

   FORALL_DOWNWARD(edgList,depth)
   {
       ed = edgList[depth]->first;
       while (ed!=0)  
       {  
	   edNext = ed->next;
	   if ( (ed->type)==T_GREEN_1 || ((ed->type)==T_GREEN_3) )  
	   {  
	     edgList[ed->depth]->remove(ed);
	     edgAlloc.Return(ed);
	     noOfEdges--;
             noOfTriangles--;
	   }
	   ed = edNext;
       } 
   }

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


void MESH3:: RelGreenIrr()
{
   TR3  *t1,*t2,*t3,*tFather;
   EDG3 *e1Mid = 0, *e2Mid = 0, *ed, *edNext;
   int  k, depth;
  
   //    remove all irregular green triangles  ( on the boundary)

   FORALL_DOWNWARD(trList,depth)
   {
       TR3 *t = trList[depth]->first;
       
       // -- Loop for all triangles; ApplyT could not be used 
       
       while (t != 0)                
       {                           
	   tFather = t->father;
	   if (tFather==0)               	{ t = t->next; continue; }
	   if (tFather->type != T_GREEN_2B) 	{ t = t->next; continue; }
         
	   // --  Remove three green triangles, print internal error message
	   //     if the green edge is not found
 
	   t1 = t; t2 = t1->next; t3 = t2->next;

	   k = 0;
	   if (((t1->e1)->type)==T_GREEN_2B) { e1Mid = t1->e1; k++;}
	   if (((t1->e2)->type)==T_GREEN_2B) { e1Mid = t1->e2; k++;}
	   if (((t1->e3)->type)==T_GREEN_2B) { e1Mid = t1->e3; k++;}
	   if (k != 1)
             cout << "\n*** 1RelGreenIrr: Triangle with impossible state k="
	    	  << k << "\n";
      
	   k = 0;
	   if (((t2->e1)->type)==T_GREEN_2B) { e2Mid = t2->e1; k++;}
	   if (((t2->e2)->type)==T_GREEN_2B) { e2Mid = t2->e2; k++;}
	   if (((t2->e3)->type)==T_GREEN_2B) { e2Mid = t2->e3; k++;}
	   if (k!=1)
	       cout << "\n*** 2RelGreenIrr: Triangle with impossible state k="
	            << k << "\n";

	   t = t3->next;   			// next while-step   
 
	   // Remove the three green triangles; return the space       
	   
	   trList[t1->depth]->remove(t1);
	   trList[t2->depth]->remove(t2);
	   trList[t3->depth]->remove(t3);
	   delete t1;
	   delete t2;
	   delete t3;
	
	   noOfTriangles -= 2;
	   
	   tFather->firstSon = 0;
	   tFather->type = T_WHITE;
	   
	   edgList[e1Mid->depth]->remove(e1Mid);
	   edgList[e2Mid->depth]->remove(e2Mid);
	   edgAlloc.Return(e1Mid);
	   edgAlloc.Return(e2Mid);

	   noOfEdges -= 2;
       }
   }

   FORALL_DOWNWARD(edgList,depth)
   {
      ed = edgList[depth]->first;
       while (ed!=0)  
       {  
	   edNext = ed->next;
	   if ((ed->type)!=T_GREEN_2B)    { ed = edNext; continue; }
	   edgList[ed->depth]->remove(ed);
	   edgAlloc.Return(ed);
	   noOfEdges--;
	   noOfTriangles -= 1;

	   ed = edNext;
       } 
   }

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


void MESH3:: RelGreen1Td()
  {
    TET3  *td1old,*td2old,*tdFather;
  
    TET3 *td = tetList[refStep]->first;

    //   remove all Green-I tetrahedra

    while (td != 0)                
    {                           
        tdFather = td->father;
        if (tdFather == 0)            		{ td = td->next; continue; }
        if (tdFather->type != T_GREEN_1) 	{ td = td->next; continue; }
  
    
	//  --		two Green-I tetrahedra
	
        td1old = td; td2old = td->next;
	td = td2old->next;
	
	if 	(td1old->mark) tdFather->mark = True;     // refine father
	else if (td2old->mark) tdFather->mark = True;     // refine father
	else tdFather->mark = False; 

	//  --		correct the white neighbors
	correctWhiteNbs(td1old);
	correctWhiteNbs(td2old);

        // --  Remove the two Green-I tetrahedra from the list of tetrahedra

        tetList[refStep]->remove(td1old);
        tetList[refStep]->remove(td2old);
	returnTET3(td1old);
	returnTET3(td2old);

        tdFather->firstSon = 0;
        tdFather->type = T_WHITE;

        noOfTetrahedra--;
        noOfTriangles--;
    }

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

void MESH3:: RelGreen2ATd()
{
    TET3  *td1old,*td2old,*td3old,*td4old,*tdFather;
    EDG3 *edMid=0;
    
    TET3 *td = tetList[refStep]->first;
 
    //      remove all green-2A tetrahedra

    while (td != 0)                
    {                           
	tdFather = td->father;
	if (tdFather==0)                  { td = td->next; continue; }
	if (tdFather->type != T_GREEN_2A) { td = td->next; continue; }
	
	// Remove four green-2A tetrahedra, print internal error message
	// if the inner green-2A edge is not found
	
	td1old = td; td2old = td->next;
	td3old = td2old->next; td4old = td3old->next;
	

	td = td4old->next;
 
	if    	(td1old->mark) tdFather->mark = True;     // refine father
	else if (td2old->mark) tdFather->mark = True;     // refine father
	else if (td3old->mark) tdFather->mark = True;     // refine father
	else if (td4old->mark) tdFather->mark = True;     // refine father
	else tdFather->mark = False; 

        //   inner edge
	if ((td1old->e1)->type == T_GREEN_2A) edMid = td1old->e1;
        else
        if ((td1old->e2)->type == T_GREEN_2A) edMid = td1old->e2;
	else
	if ((td1old->e3)->type == T_GREEN_2A) edMid = td1old->e3;
	else
	if ((td1old->e4)->type == T_GREEN_2A) edMid = td1old->e4;
	else
	if ((td1old->e5)->type == T_GREEN_2A) edMid = td1old->e5;
	else
	if ((td1old->e6)->type == T_GREEN_2A) edMid = td1old->e6;

	if (edMid == 0) cout << "\n\nAlarm in RelGreen2ATd: edMid == 0\n";

	//  --		correct the white neighbours
	correctWhiteNbs(td1old);
	correctWhiteNbs(td2old);
	correctWhiteNbs(td3old);
	correctWhiteNbs(td4old);

	//  Remove the four green ones from the list of tetrahedra

	tetList[refStep]->remove(td1old);
	tetList[refStep]->remove(td2old);
	tetList[refStep]->remove(td3old);
	tetList[refStep]->remove(td4old);
	returnTET3(td1old);
	returnTET3(td2old);
	returnTET3(td3old);
	returnTET3(td4old);

	tdFather->firstSon = 0;
	tdFather->type = T_WHITE;
	
	noOfTetrahedra -= 3;
	noOfTriangles  -= 4;
	
	edgList[edMid->depth]->remove(edMid);
	edgAlloc.Return(edMid);

	noOfEdges -=1;
    }

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

void MESH3:: RelGreen2BTd()
{
    TET3  *td1old,*td2old,*td3old,*tdFather;
    
    TET3 *td = tetList[refStep]->first;

    //   remove all green-2B tetrahedra

    while (td != 0)                
    {                           
	tdFather = td->father;
	if (tdFather == 0)                { td = td->next; continue; }
        if (tdFather->type != T_GREEN_2B) { td = td->next; continue; }
  
	//  Remove three green2B tetrahedra
	
	td1old = td; td2old = td->next; td3old = td2old->next; 

	td = td3old->next;		//  next while-step 
	
	if    	(td1old->mark) tdFather->mark = True;     // refine father
	else if (td2old->mark) tdFather->mark = True;     // refine father
	else if (td3old->mark) tdFather->mark = True;     // refine father
	else tdFather->mark = False; 

	//  --		correct the white neighbours
	correctWhiteNbs(td1old);
	correctWhiteNbs(td2old);
	correctWhiteNbs(td3old);

	//  Remove the three green ones from the list of tetrahedra
	
	tetList[refStep]->remove(td1old);
	tetList[refStep]->remove(td2old);
	tetList[refStep]->remove(td3old);
	returnTET3(td1old);
	returnTET3(td2old);
	returnTET3(td3old);

	tdFather->firstSon = 0;
	tdFather->type = T_WHITE;

	noOfTetrahedra -=2;
	noOfTriangles  -= 2;
    }

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


void MESH3:: RelGreen3Td()
{
    TET3 *td1old,*td2old,*td3old,*td4old,*tdFather;
    
    TET3 *td = tetList[refStep]->first;
     
    //    remove all green-III tetrahedra

    while (td != 0)                
    {                           
	tdFather = td->father;
	if (tdFather==0)                 { td = td->next; continue; }
	if ((tdFather->type)!=T_GREEN_3) { td = td->next; continue; }
	
	//  Remove four green tetrahedra, print internal error message
	//  if the inner green faces are not found
	
	td1old = td; td2old = td->next;
	td3old = td2old->next; td4old = td3old->next;
	
	
	td = td4old->next;	// next step in while
	
	if    	(td1old->mark) tdFather->mark = True;     // refine father
	else if (td2old->mark) tdFather->mark = True;     // refine father
	else if (td3old->mark) tdFather->mark = True;     // refine father
	else if (td4old->mark) tdFather->mark = True;     // refine father
	else tdFather->mark = False; 


	//  --		correct the white neighbours
	correctWhiteNbs(td1old);
	correctWhiteNbs(td2old);
	correctWhiteNbs(td3old);
	correctWhiteNbs(td4old);

	// Substitute the old (father) tetrahedron for the four green ones
	// in the list of tetrahedra; return the space       
	
	tetList[refStep]->remove(td1old);
	tetList[refStep]->remove(td2old);
	tetList[refStep]->remove(td3old);
	tetList[refStep]->remove(td4old);
	returnTET3(td1old);
	returnTET3(td2old);
	returnTET3(td3old);
	returnTET3(td4old);

	tdFather->firstSon = 0;
	tdFather->type = T_WHITE;
	
	noOfTetrahedra -= 3;
	noOfTriangles  -= 3;
    }

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

Bool TET3:: noTdinNb(TET3* tet)
{
   TET* nbTet;

   nbTet = n1->castToTET3();
   if ( nbTet == tet)  return True; 

   nbTet = n2->castToTET3();
   if ( nbTet == tet) return True;

   nbTet = n3->castToTET3();
   if ( nbTet == tet) return True;

   nbTet = n4->castToTET3();
   if ( nbTet == tet) return True;

   return False;
}
//-------------------------------------------------------------------------

void TET3:: correctNbinTet(TET3 *tet)
{
   TET* nbTet;

   nbTet = n1->castToTET3();
   if ( nbTet == tet)  {n1 = tet->father; return; }

   nbTet = n2->castToTET3();
   if ( nbTet == tet)  {n2 = tet->father; return; }

   nbTet = n3->castToTET3();
   if ( nbTet == tet)  {n3 = tet->father; return; }

   nbTet = n4->castToTET3();
   if ( nbTet == tet)  {n4 = tet->father; return; }

// cout << "\nImpossible state in correctNbinTet\n";
}

void MESH3:: correctWhiteNbs(TET3 *tet)
{
    TET3* nbTet;
    Vector<PATCH*> nbs(4);
    int i;

    tet->getNeighbours(nbs);

    FORALL(nbs,i) 
    {
	if (nbs[i]) 
	{
	    nbTet = nbs[i]->castToTET3();
	    nbTet->correctNbinTet(tet);

           // if ( (nbTet->type != T_GREEN_1) && (nbTet->type != T_GREEN_2A) &&
		// (nbTet->type != T_GREEN_2B) && (nbTet->type != T_GREEN_3) )  
	    
	}
    }
}
//-------------------------------------------------------------------------

int MESH3:: CheckTd(TET3 *td)
{
   TET3 *nbTd;
   PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
    

   if ( n1 == 0 )
   {
   	cout << "\n\tneighbour 1==0:\n";
	td->print(cout);
	(td->p2)->print(cout); (td->p3)->print(cout); (td->p4)->print(cout);
   }

   if ( n2 == 0 )
   {
   	cout << "\n\tneighbour 2==0:\n";
	td->print(cout);
	(td->p1)->print(cout); (td->p3)->print(cout); (td->p4)->print(cout);
   }

   if ( n3 == 0 )
   {
   	cout << "\n\tneighbour 3==0:\n";
	td->print(cout);
	(td->p1)->print(cout); (td->p2)->print(cout); (td->p4)->print(cout);
   }

   if ( n4 == 0 )
   {
   	cout << "\n\tneighbour 4==0:\n";
	td->print(cout);
	(td->p1)->print(cout); (td->p2)->print(cout); (td->p3)->print(cout);
   }


   nbTd = n1->castToTET3();
   if ( nbTd ) 
	{
	   if( !nbTd->noTdinNb(td) ) {
		cout << "\n\t error in neighbor relation\n";
		td->print(cout); nbTd->print(cout);
	   }
	}

   nbTd = n2->castToTET3();
   if ( nbTd ) 
	{
	   if( !nbTd->noTdinNb(td) ) {
		cout << "\n\t error in neighbor relation\n";
		td->print(cout); nbTd->print(cout);
	   }
	}

   nbTd = n3->castToTET3();
   if ( nbTd ) 
	{
	   if( !nbTd->noTdinNb(td) ) {
		cout << "\n\t error in neighbor relation\n";
		td->print(cout); nbTd->print(cout);
	   }
	}

   nbTd = n4->castToTET3();
   if ( nbTd ) 
	{
	   if( !nbTd->noTdinNb(td) ) {
		cout << "\n\t error in neighbor relation\n";
		td->print(cout); nbTd->print(cout);
	   }
	}


   if ( n1->castToTET3() && n2->castToTET3() &&
        n3->castToTET3() && n4->castToTET3() )
   {
       if (td->boundP != INTERIOR) {
	   cout << "\n\t error in boundary definition: boundP != INTERIOR\n";
	   td->print(cout);
       }
   }
   else
   {
       if (td->boundP != BOUNDARY) {
	   cout << "\n\t error in boundary definition: boundP != BOUNDARY\n";
	   td->print(cout);
       }
   }

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


int MESH3:: consistencyCheck()
{
   PT3  *p;
   TET3 *tet;

   //int  no, noInner, noBound, noWithFather, noWithoutFather;
   //int  noWhite, noRed, noGreenI, noGreenIIA, noGreenIIB, noGreenIII;

   DListIter<TET3> iterTET(tetras());

   cout << "\n\t Begin of consistencyCheck:\n";
   
   cout << "\n\t no of points     = " << noOfPoints     << "\n";
   cout <<   "\t no of edges      = " << noOfEdges      << "\n";
   cout <<   "\t no of triangles  = " << noOfTriangles  << "\n";
   cout <<   "\t no of tetrahedra = " << noOfTetrahedra << "\n";

   int eulerC = noOfPoints - noOfEdges + noOfTriangles - noOfTetrahedra;
   cout << "\n\t Euler's characteristic = " <<  eulerC << "\n";
   if (eulerC != 1) cout << "\n\t Dangerous: Euler's characteristic not 1\n\n";

   Vector<int> PNodes(noOfPoints);

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

   while ((tet=iterTET.all()))  {
     if(!CheckTd(tet)) return False;;
   }

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


   CompAngles();

/*
   no = 0; noInner = 0; noBound = 0;
   while (p =iterPT.all())   {
	no++;
        if ( p->boundP == DIRICHLET ) noBound++;
        if ( p->boundP == INTERIOR  ) noInner++;
   }
   cout << "\n\t no of points in list  = " << no;
   cout << "\n\t no of boundary points = " << noBound;
   cout << "\n\t no of    inner points = " << noInner;


   no = 0; noInner = 0; noBound = 0; noWithFather=0; noWithoutFather=0;
   noWhite=0; noRed=0; noGreenI=0; noGreenIIA=0; noGreenIIB=0; noGreenIII=0;
   noZero=0; noGreen=0;
   while (e =iterEDG.all())  {
	no++;
        if ( e->boundP == DIRICHLET ) noBound++;
        if ( e->boundP == INTERIOR  ) noInner++;
        if ( e->type   == T_WHITE   ) noWhite++;
        if ( e->type   == T_RED     ) noRed++;
        if ( e->type   == T_GREEN_1 ) noGreenI++;
        if ( e->type   == T_GREEN_2A) noGreenIIA++;
        if ( e->type   == T_GREEN_2B) noGreenIIB++;
        if ( e->type   == T_GREEN_3 ) noGreenIII++;
        if ( e->type   == T_GREEN   ) noGreen++;
        if ( e->type   == T_ZERO    ) noZero++;

        if ( e->father != 0         ) noWithFather++;
        else if ( e->father == 0    ) noWithoutFather++;
   }

   cout << "\n\t no of edges in list        = " << no;
   cout << "\n\t no of boundary edges       = " << noBound;
   cout << "\n\t no of    inner edges       = " << noInner;
   cout << "\n\t no of edges with father    = " << noWithFather;
   cout << "\n\t no of edges without father = " << noWithoutFather;
   cout << "\n\t no of white    edges       = " << noWhite;
   cout << "\n\t no of red      edges       = " << noRed;
   cout << "\n\t no of greenI   edges       = " << noGreenI;
   cout << "\n\t no of greenIIA edges       = " << noGreenIIA;
   cout << "\n\t no of greenIIB edges       = " << noGreenIIB;
   cout << "\n\t no of greenIII edges       = " << noGreenIII;
   cout << "\n\t no of green    edges       = " << noGreen;
   cout << "\n\t no of zero     edges       = " << noZero;

   no = 0; 
   iterTET.reset();
   while (tet=iterTET.all()) no++;
   cout << "\n\t no of tets in list = " << no << "\n";

*/
   cout << "\n\t End of consistencyCheck.\n";
   cout << "\n\t ----------------------------------------------------------\n\n";

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


int MESH3:: CheckAngle(TR3 *t)
{
   PT3 *p1=t->p1, *p2=t->p2, *p3=t->p3;
   Real e1X, e1Y, e1Z;
   Real e2X, e2Y, e2Z;
   Real s12, s11, s22;
   Real alpha1, alpha2, alpha3;
   int classI;

    e1X = p2->x - p1->x;
    e1Y = p2->y - p1->y;
    e1Z = p2->z - p1->z;

    e2X = p3->x - p1->x;
    e2Y = p3->y - p1->y;
    e2Z = p3->z - p1->z;

    s12 = e1X*e2X + e1Y*e2Y + e1Z*e2Z;
    s11 = sqrt(e1X*e1X + e1Y*e1Y + e1Z*e1Z);
    s22 = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);

    alpha1 = 180.0*acos(s12/(s11*s22))/REALPI;
    if (alpha1 > alphaMax) alphaMax = alpha1;
    if (alpha1 < alphaMin) alphaMin = alpha1;

    e1X = p1->x - p2->x;
    e1Y = p1->y - p2->y;
    e1Z = p1->z - p2->z;

    e2X = p3->x - p2->x;
    e2Y = p3->y - p2->y;
    e2Z = p3->z - p2->z;

    s12 = e1X*e2X + e1Y*e2Y + e1Z*e2Z;
    s11 = sqrt(e1X*e1X + e1Y*e1Y + e1Z*e1Z);
    s22 = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);

    alpha2 = 180.0*acos(s12/(s11*s22))/REALPI;
    if (alpha2 > alphaMax) alphaMax = alpha2;
    if (alpha2 < alphaMin) alphaMin = alpha2;

    alpha3 = 180.0-alpha1-alpha2;
    if (alpha3 > alphaMax) alphaMax = alpha3;
    if (alpha3 < alphaMin) alphaMin = alpha3;
    
    classI = (int)(alpha1/5);
    alpha[classI]++;
    
    classI = (int)(alpha2/5);
    alpha[classI]++;

    classI = (int)(alpha3/5);
    alpha[classI]++;

    if ( (alpha1 < 10.) || (alpha2 < 10.) || (alpha3 < 10.) )
      {
       cout << " alphas: " << alpha1 << alpha2 << alpha3 << "\n";
       cout << " P1    : " << p1->x << p1->y << p1->z << "\n";
       cout << " P2    : " << p2->x << p2->y << p2->z << "\n";
       cout << " P3    : " << p3->x << p3->y << p3->z << "\n";
      }

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


int MESH3:: CheckCircle(TET3 *tet)
{
   PT3   *p1, *p2, *p3, *p4;
   Real  e1X, e1Y, e1Z;
   Real  e2X, e2Y, e2Z;
   Real  dyz, dxz, dxy;
   Real  len, rhoT, sigma, vT, sT, hT;
   EDG3  *eds[1+6];
   int   i;
   TR3   *t1, *t2, *t3, *t4;
   TR3   *trs[1+4];

   Real  x1, x2, x3, y1, y2, y3, z1, z2, z3;

   tet->getTriangles(trs);
   t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];

   p1 = tet->p1;
   p2 = tet->p2;
   p3 = tet->p3;
   p4 = tet->p4;

   x1 = (p2->x) - (p1->x);
   x2 = (p2->y) - (p1->y);
   x3 = (p2->z) - (p1->z);

   y1 = (p3->x) - (p1->x);
   y2 = (p3->y) - (p1->y);
   y3 = (p3->z) - (p1->z);

   z1 = (p4->x) - (p1->x);
   z2 = (p4->y) - (p1->y);
   z3 = (p4->z) - (p1->z);

   vT = SIXTH*Abs(x1*(y2*z3 - y3*z2) - x2*(y1*z3 -y3*z1) + x3*(y1*z2 - y2*z1)); 
   // = tet->volume;

   sT = ZERO; 

   p1=t1->p1; p2=t1->p2, p3=t1->p3;

   e1X = p2->x - p1->x;
   e1Y = p2->y - p1->y;
   e1Z = p2->z - p1->z;

   e2X = p3->x - p1->x;
   e2Y = p3->y - p1->y;
   e2Z = p3->z - p1->z;
   
   dyz = e1Y*e2Z - e1Z*e2Y;
   dxz = e1Z*e2X - e1X*e2Z;
   dxy = e1X*e2Y - e1Y*e2X;
   sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);

   p1=t2->p1; p2=t2->p2, p3=t2->p3;
   e1X = p2->x - p1->x;
   e1Y = p2->y - p1->y;
   e1Z = p2->z - p1->z;

   e2X = p3->x - p1->x;
   e2Y = p3->y - p1->y;
   e2Z = p3->z - p1->z;
   
   dyz = e1Y*e2Z - e1Z*e2Y;
   dxz = e1Z*e2X - e1X*e2Z;
   dxy = e1X*e2Y - e1Y*e2X;
   sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);

   p1=t3->p1; p2=t3->p2, p3=t3->p3;
   e1X = p2->x - p1->x;
   e1Y = p2->y - p1->y;
   e1Z = p2->z - p1->z;

   e2X = p3->x - p1->x;
   e2Y = p3->y - p1->y;
   e2Z = p3->z - p1->z;
   
   dyz = e1Y*e2Z - e1Z*e2Y;
   dxz = e1Z*e2X - e1X*e2Z;
   dxy = e1X*e2Y - e1Y*e2X;
   sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);

   p1=t4->p1; p2=t4->p2, p3=t4->p3;
   e1X = p2->x - p1->x;
   e1Y = p2->y - p1->y;
   e1Z = p2->z - p1->z;

   e2X = p3->x - p1->x;
   e2Y = p3->y - p1->y;
   e2Z = p3->z - p1->z;
   
   dyz = e1Y*e2Z - e1Z*e2Y;
   dxz = e1Z*e2X - e1X*e2Z;
   dxy = e1X*e2Y - e1Y*e2X;
   sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);

   sT = 0.5*sT;
   rhoT =6.0*vT/sT;

   hT =ZERO;
   tet->getEdges(eds);
   for (i=1;i<=6;i++)
     {
      e1X = (eds[i]->p1)->x;
      e1Y = (eds[i]->p1)->y;
      e1Z = (eds[i]->p1)->z;
      e2X = (eds[i]->p2)->x;
      e2Y = (eds[i]->p2)->y;
      e2Z = (eds[i]->p2)->z;

      e2X = e2X - e1X;
      e2Y = e2Y - e1Y;
      e2Z = e2Z - e1Z;

      len = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);
      if (len > hT) hT = len;
     }

   sigma = hT/rhoT;

   if (sigma > sigmaMax) sigmaMax = sigma;

   tet->returnTriangles(trs);

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


int MESH3:: CompAngles()
{
   TET3 *tet;
   TR3  *t;
   int i, noOfAngles=0;

   for (i=0;i<36;i++)
      {
       alpha[i] = 0;
      }

   alphaMax = ZERO;
   alphaMin = 90.0;
   DListIter<TR3> iterTR(triangles());
   while ((t=iterTR.all()))  CheckAngle(t);

   cout << "\n\t ----------------------------------------------------------\n";
   cout << "\t angles on boundary faces: \n";
   cout << "\n\t max. angle = " <<  alphaMax << "\n";
   cout << "\t min. angle = "   <<  alphaMin << "\n\n";

   for (i=0;i<36;i++)
      {
       cout << "\t no of angles between " << i*5 << "and " << (i+1)*5 << ":  "
            << alpha[i] << "\n";
       noOfAngles += alpha[i];
      }
   cout << "\n\t total no of angles  = " << noOfAngles << "\n";
   cout << "\n\t ----------------------------------------------------------\n";

   sigmaMax = ZERO;
   DListIter<TET3> iterTET(tetras());
   while ((tet=iterTET.all()))  CheckCircle(tet);

   cout << "\t sigmaT = " << sigmaMax << "\n"; 
   cout << "\n\t ----------------------------------------------------------\n";

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


Real MESH3:: Length(PT3 *p1, PT3 *p2)
{
   double sum = ZERO, diff;

   diff = (p1->x) - (p2->x);
   sum += diff * diff;

   diff = (p1->y) - (p2->y);
   sum += diff * diff;

   diff = (p1->z) - (p2->z);
   sum += diff * diff;
 
   return sum;
}
//-------------------------------------------------------------------------

int MESH3:: PinTr(PT3 *p, TR3 *tr)
{
   PT3 *p1 = tr->p1, *p2 = tr->p2, *p3 = tr->p3;
   if (p==p1 || p==p2 || p==p3) return True;
   else return False;
}
//-------------------------------------------------------------------------

int MESH3:: PinTet(PT3 *p, TET3 *tet)
{
   PT3 *p1 = tet->p1, *p2 = tet->p2, *p3 = tet->p3, *p4 = tet->p4;
   if (p==p1 || p==p2 || p==p3 || p==p4) return True;
   else return False;
}
//-------------------------------------------------------------------------

EDG3* MESH3:: EdinTr(TR3 *tr,PT3 *p1,PT3 *p2)
{ 
   EDG3 *e1 = tr->e1, *e2 = tr->e2, *e3 = tr->e3;
   PT3  *pt1, *pt2;
   
   pt1 = e1->p1; pt2 = e1->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e1;

   pt1 = e2->p1; pt2 = e2->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e2;

   pt1 = e3->p1; pt2 = e3->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e3;

   cout << "\nEdinTr(...): edge not found\n";
   cout.flush(); abort();
   return 0;
}
//-------------------------------------------------------------------------

EDG3* MESH3:: EdinTet(PT3 *p1,PT3 *p2,TET3 *td)
{ 
   EDG3 *e1 = td->e1, *e2 = td->e2, *e3 = td->e3,
        *e4 = td->e4, *e5 = td->e5, *e6 = td->e6;
   PT3  *pt1, *pt2;
   
   pt1 = e1->p1; pt2 = e1->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e1;

   pt1 = e2->p1; pt2 = e2->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e2;

   pt1 = e3->p1; pt2 = e3->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e3;

   pt1 = e4->p1; pt2 = e4->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e4;

   pt1 = e5->p1; pt2 = e5->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e5;

   pt1 = e6->p1; pt2 = e6->p2;
   if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
      return e6;

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

int MESH3:: TestFace(TR3 *tr)
{
   EDG3 *eds[1+3];
   int  i, noOfRefEdg = 0;
  
   tr->getEdges(eds);
   for (i=1; i<=3; i++)
      {
        if ( (eds[i]->firstSon) != 0)
          noOfRefEdg++;
      }

    return noOfRefEdg;
}
//-------------------------------------------------------------------------

int MESH3:: TestEdsInTd(TET3 *tet)
{
   EDG3 *eds[1+6];
   int  i, noOfRefEdg = 0;
  
   tet->getEdges(eds);
   for (i=1; i<=6; i++)
      {
        if ( (eds[i]->firstSon) != 0)
          noOfRefEdg++;
      }

    return noOfRefEdg;
}
//-------------------------------------------------------------------------

int MESH3:: testSons(TR3 *t, TET3 *tet)
  {
    EDG3 *e1=t->e1, *e2=t->e2, *e3=t->e3, *ed=0;
    PT3 *pm1, *pm2, *pm3;
    TET3 *tetSon;
    int noOfRefEdg = 0, i;
  
    pm1 = e1->pm; pm2 = e2->pm; pm3 = e3->pm;
    tetSon = tet->firstSon;
    for (i=1; i<=8; i++)
    {
        ed = EdinTet(pm1,pm2,tetSon);
        if (ed != 0)
           if (ed->refined())  noOfRefEdg++;
      
        ed = EdinTet(pm1,pm3,tetSon);
        if (ed != 0)
           if (ed->refined())  noOfRefEdg++;

        ed = EdinTet(pm2,pm3,tetSon);
        if (ed != 0)
           if (ed->refined())  noOfRefEdg++;

        if ( noOfRefEdg > 0 ) return noOfRefEdg;
        tetSon = tetSon->next;
    }

    return noOfRefEdg;
  }
//-------------------------------------------------------------------------


ostream& operator << (ostream& os, MESH3& t)
{
    PT3*  p;
    EDG3* e;
    TR3*  tr;
    TET3 * tet;

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

    DListIter<PT3>  iterPT (t.points());
    DListIter<TR3>  iterTR (t.triangles());
    DListIter<EDG3> iterEDG(t.edges());
    DListIter<TET3> iterTET(t.tetras());

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

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

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


void MESH3:: Info()
{
    cout << "\n Triangulation:  points   edges   tri's tetra's "
    << "  Depth=" << maxDepth << "  Steps=" << refStep
    << Form("\n%15s%8d%8d%8d%8d", " ", 
	      noOfPoints, noOfEdges, noOfTriangles, noOfTetrahedra)
    << "\n";
}
//-------------------------------------------------------------------------

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

    space += tetAlloc.MemSpace();
    space += edgAlloc.MemSpace();
    space += ptAlloc.MemSpace();
  
    space += trList[0]->first->alloc.MemSpace();

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

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

    space += EFTetAlloc.MemSpace();
    space += edgAlloc.MemSpace();
    space += ptAlloc.MemSpace();
  
    space += trList[0]->first->alloc.MemSpace();

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


void MESH3:: refinementInfo(int* noPt, int* noEdg, int* noTr, int* noTet)
{
    int blocks = varAllocator.NoOfBlocks();

    blocks += tetAlloc.NoOfBlocks();
    blocks += edgAlloc.NoOfBlocks();
    blocks += ptAlloc.NoOfBlocks();
    blocks += trList[0]->first->alloc.NoOfBlocks();

    int space = MemSpace();

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

    if (infoRefinement > 1)
    {
	cout << Form("\n%15s%8d%8d%8d%8d", "removed", 
		       noPt[1]-noPt[0], noEdg[1]-noEdg[0]
		       , noTr[1]-noTr[0], noTet[1]-noTet[0])

	<< Form("\n%15s%8d%8d%8d%8d", "generated", 
		  noPt[2]-noPt[1], noEdg[2]-noEdg[1], 
		  noTr[2]-noTr[1], noTet[2]-noTet[1])

	<< Form("\n%15s%8d%8d%8d%8d", "added",
		  noOfPoints-noPt[2], noOfEdges-noEdg[2],
		  noOfTriangles-noTr[2], noOfTetrahedra-noTet[2])
	<<  "   for green closure";
    }
    cout << Form("\n%15s%8d%8d%8d%8d   Space : %1.6f MB","total",
		 noOfPoints, noOfEdges, noOfTriangles,noOfTetrahedra, 
		   float(space)/1e6);

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

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

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

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


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

void PT3:: reset()
{
    boundP = INTERIOR;
    mark = classA = depth = NUL;
    node   = 0;
}
//-------------------------------------------------------------------------

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

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

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

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

void PT3:: readBC(BufferedParser& parser)
{
    readBCValues(&boundP, &classA, parser);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


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

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


void EDG3:: compNormalVector(Vector<Real>& /*n*/, Real* /*length*/) const
{
    notImplemented("EDG3::compNormalVector");
}
//-------------------------------------------------------------------------


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

void EDG3:: 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; 
    x[3] = L1*p1->z + L2*p2->z; 
}
//-------------------------------------------------------------------------

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

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

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


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

void TR3:: reset()
{
    p1 = p2 = p3 = 0;
    e1 = e2 = e3 = 0;
    mark = classA = depth = NUL;
    isoP   = PLANE;
    boundP = INTERIOR;
    type   = MESH::T_WHITE;
    next = prev = father = firstSon = 0;
    node = 0;
    t31 = t32 = 0;
}
//-------------------------------------------------------------------------

void TR3:: getEdges(EDG3* eds[]) const 
{ 
    eds[1] = e1;  eds[2] = e2;  eds[3] = e3; 
}
//-------------------------------------------------------------------------

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

void TR3:: 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; 
    x[3] = L1*p1->z + L2*p2->z + L3*p3->z; 
}
//-------------------------------------------------------------------------

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

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


Real TR3:: volume() const
{
    Vector<Real> u(3), v(3), w(3);

    u[1] = p2->x - p1->x;	
    u[2] = p2->y - p1->y;
    u[3] = p2->z - p1->z;

    v[1] = p3->x - p1->x;	
    v[2] = p3->y - p1->y;
    v[3] = p3->z - p1->z;

    CrossProduct(w,u,v);
    return 0.5 * Norm(w);
}
//-------------------------------------------------------------------------


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

    if (isoP > e1->isoP) e1->isoP = isoP;
    if (isoP > e2->isoP) e2->isoP = isoP;
    if (isoP > e3->isoP) e3->isoP = isoP;


    if (boundP == DIRICHLET)
    {
	e1->boundP = e2->boundP = e3->boundP = boundP;
	e1->classA = e2->classA = e3->classA = classA;

	p1->boundP = p2->boundP = p3->boundP = boundP;
	p1->classA = p2->classA = p3->classA = classA;
    }
    else if (boundP == CAUCHY || boundP == NEUMANN)
    {
	if (e1->boundP != DIRICHLET)
	{
	    e1->boundP = boundP;
	    e1->classA = classA;
	}
	if (e2->boundP != DIRICHLET) 
	{
	    e2->boundP = boundP;
	    e2->classA = classA;
	}
	if (e3->boundP != DIRICHLET) 
	{
	    e3->boundP = boundP;
	    e3->classA = classA;
	}

	if (p1->boundP != DIRICHLET)
	{
	    p1->boundP = boundP;
	    p1->classA = classA;
	}

	if (p2->boundP != DIRICHLET) 
	{
	    p2->boundP = boundP;
	    p2->classA = classA;
	}
	if (p3->boundP != DIRICHLET) 
	{
	    p3->boundP = boundP;
	    p3->classA = classA;
	}
    }
    else 
    {
	cout << "\n*** readBC: Invalid Boundary Type for Face " << boundP << "\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------


void TR3:: readHyperBC(BufferedParser& /*parser*/)
{
    if (boundP == DIRICHLET)
    {
	e1->boundP = e2->boundP = e3->boundP = boundP;
	e1->classA = e2->classA = e3->classA = classA;

	p1->boundP = p2->boundP = p3->boundP = boundP;
	p1->classA = p2->classA = p3->classA = classA;
    }
    else if (boundP == CAUCHY || boundP == NEUMANN)
    {
	if (e1->boundP != DIRICHLET)
	{
	    e1->boundP = boundP;
	    e1->classA = classA;
	}
	if (e2->boundP != DIRICHLET) 
	{
	    e2->boundP = boundP;
	    e2->classA = classA;
	}
	if (e3->boundP != DIRICHLET) 
	{
	    e3->boundP = boundP;
	    e3->classA = classA;
	}

	if (p1->boundP != DIRICHLET)
	{
	    p1->boundP = boundP;
	    p1->classA = classA;
	}
	if (p2->boundP != DIRICHLET) 
	{
	    p2->boundP = boundP;
	    p2->classA = classA;
	}
	if (p3->boundP != DIRICHLET) 
	{
	    p3->boundP = boundP;
	    p3->classA = classA;
	}
    }
    else 
    {
	cout << "\n*** readBC: Invalid Boundary Type for Face " << boundP << "\n";
	cout.flush(); abort();
    }
}

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

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

    j.J(1,2) = p2->y - p1->y;
    j.J(2,2) = p3->y - p1->y;
    j.J(3,2) = p4->y - p1->y;
    
    j.J(1,3) = p2->z - p1->z;
    j.J(2,3) = p3->z - p1->z;
    j.J(3,3) = p4->z - p1->z;
	
    j.detJ = detJ;

    //       = j.J(1,1) * (j.J(2,2)*j.J(3,3) - j.J(3,2)*j.J(2,3))
    //	     + j.J(3,1) * (j.J(1,2)*j.J(2,3) - j.J(2,2)*j.J(1,3))
    //	     + j.J(2,1) * (j.J(1,3)*j.J(3,2) - j.J(3,3)*j.J(1,2));
    // if (j.detJ <= 0.0) cout << "\n*** Jacobi-Det. <= 0 " << j.detJ
    //						<< " encountered\n";
}

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

    compJ(j);

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

void TET3:: print(ostream& os) const
{
    TET::print(os);
    os << "  type " << (int) type << "  class " << (int) classA
    << "  mark " << (int)mark << "  depth " << (int) depth << "\n";
}
//-------------------------------------------------------------------------

void TET3:: reset()
{
    mark = classA = depth = NUL;
    boundP = INTERIOR;
    detJ = -1.0;
    father = firstSon = 0;
    node = 0;
    type = MESH::T_WHITE;
    error = 0.0;

    p1 = p2 = p3 = p4 = 0;
    n1 = n2 = n3 = n4 = 0;
    e1 = e2 = e3 = e4 = e5 = e6 = 0;
}

void TET3:: setBoundary()
{
    if ( !(n1->castToTET3() && n2->castToTET3() &&
           n3->castToTET3() && n4->castToTET3()) )  boundP = BOUNDARY;
}

//-------------------------------------------------------------------------
void TET3:: getPoints(PT3* pts[]) const 
{ 
    pts[1] = p1;  pts[2] = p2;  pts[3] = p3;  pts[4] = p4;
}
//-------------------------------------------------------------------------

void TET3:: getEdges(EDG3* eds[]) const 
{ 
    eds[1] = e1;  eds[2] = e2;  eds[3] = e3; 
    eds[4] = e4;  eds[5] = e5;  eds[6] = e6; 
}
//-------------------------------------------------------------------------

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

    if (p1->node == 0 || p2->node == 0 || p3->node == 0 ||p4->node == 0)  // !!!
    {
	cout <<"\n*** TET3:: getEdgeOrientation: node == 0 !\n"; 
	cout.flush(); abort();
    }
 
    if (p1->node < p2->node) orient[1] = True;
    if (p2->node < p3->node) orient[2] = True;
    if (p3->node < p1->node) orient[3] = True;
    if (p3->node < p4->node) orient[4] = True;
    if (p2->node < p4->node) orient[5] = True;
    if (p1->node < p4->node) orient[6] = True;
}
//-------------------------------------------------------------------------


void TET3:: markTriangles(TET3* tet)
{
   TET* nbTet;

   nbTet = n1->castToTET3();
   if ( nbTet == tet) { mark++; return; }

   nbTet = n2->castToTET3();
   if ( nbTet == tet) { mark += 2; return; }

   nbTet = n3->castToTET3();
   if ( nbTet == tet) { mark += 4; return; }

   nbTet = n4->castToTET3();
   if ( nbTet == tet) { mark += 8; return; }

   cout << "\n Impossible state in markTriangles\n";
}
//-------------------------------------------------------------------------

void TET3:: getMarkedTriangles(int trsMark[])
{
   for (int i=1; i<=4; i++) trsMark[i] = 0;

   switch (mark) {
        case 0:		break;
	case 1:		trsMark[1] = 1; break;
	case 2:		trsMark[2] = 1; break;
	case 3:		trsMark[2] = 1; trsMark[1] = 1; break;
	case 4:		trsMark[3] = 1; break;
	case 5:		trsMark[3] = 1; trsMark[1] = 1; break;
	case 6:		trsMark[3] = 1; trsMark[2] = 1; break;
	case 7:		trsMark[3] = 1; trsMark[2] = 1; trsMark[1] = 1; break;
	case 8:		trsMark[4] = 1; break;
	case 9:		trsMark[4] = 1; trsMark[1] = 1; break;
	case 10:	trsMark[4] = 1; trsMark[2] = 1; break;
	case 11:	trsMark[4] = 1; trsMark[2] = 1; trsMark[1] = 1; break;
	case 12:	trsMark[4] = 1; trsMark[3] = 1; break;
	case 13:	trsMark[4] = 1; trsMark[3] = 1; trsMark[1] = 1; break;
	case 14:	trsMark[4] = 1; trsMark[3] = 1; trsMark[2] = 1; break;
	case 15:	trsMark[4] = 1; trsMark[3] = 1; trsMark[2] = 1; 
			trsMark[1] = 1; break;
	default:        cout << "\n dangerous state in getMarkedTriangles\n";
   }

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

void TET3:: getAllFaces(Vector<PATCH*>& trs) const
{ 
   TR3  *t1, *t2, *t3, *t4;
   TET3 *nt;
   void *t;

   //  t1
   if ( (nt = n1->castToTET3()) ) 
   {
       t1 = new TR3;
       t1->boundP = INTERIOR;
       t1->depth = depth;
       t1->p1 = p2;
       t1->p2 = p3;
       t1->p3 = p4;
       t1->e1 = e4;
       t1->e2 = e5;
       t1->e3 = e2;
       completeT(t1,0);
       // t1->t31 = this;    causing warnings, subst. by the next two lines 
       t       = (void*) this;
       t1->t31 = (TET3*) t;
       t1->t32 = nt;
   }
   else
   {
       t1 = n1->castToTR3(); 
       // t1->t31 = this;    causing warnings, subst. by the next two lines
       t       = (void*) this;
       t1->t31 = (TET3*) t;
       t1->t32 = 0; 
   }
   
   //  t2
   if ( (nt = n2->castToTET3()) )
   { 
       t2 = new TR3;
       t2->boundP = INTERIOR;
       t2->depth = depth;
       t2->p1 = p1;
       t2->p2 = p4;
       t2->p3 = p3;
       t2->e1 = e4;
       t2->e2 = e1;
       t2->e3 = e6;
       completeT(t2,0);     
       // t2->t31 = this;   causing warnings, subst. by the next two lines
       t       = (void*) this;
       t2->t31 = (TET3*) t;
       t2->t32 = nt;
   }
   else
   {
       t2 = n2->castToTR3(); 
       // t2->t31 = this;  causing warnings, subst. by the next two lines
       t       = (void*) this;
       t2->t31 = (TET3*) t;
       t2->t32 = 0; 
   }

   //  t3
   if ( (nt = n3->castToTET3()) ) 
   {
       t3 = new TR3;
       t3->boundP = INTERIOR;
       t3->depth = depth;
       t3->p1 = p1;
       t3->p2 = p2;
       t3->p3 = p4;
       t3->e1 = e5;
       t3->e2 = e6;
       t3->e3 = e3;
       completeT(t3,0);
       // t3->t31 = this;   causing warnings, subst. by the next two lines
       t       = (void*) this;
       t3->t31 = (TET3*) t;
       t3->t32 = nt;
   }
   else
   {
       t3 = n3->castToTR3(); 
       // t3->t31 = this;   causing warnings, subst. by the next two lines
       t       = (void*) this;
       t3->t31 = (TET3*) t;
       t3->t32 = 0; 
   }
   
   //  t4
   if ( (nt = n4->castToTET3()) ) 
   {
       t4 = new TR3;
       t4->boundP = INTERIOR;
       t4->depth = depth;
       t4->p1 = p1;
       t4->p2 = p3;
       t4->p3 = p2;
       t4->e1 = e2;
       t4->e2 = e3;
       t4->e3 = e1;
       completeT(t4,0);
       // t4->t31 = this;   causing warnings, subst. by the next two lines
       t       = (void*) this;
       t4->t31 = (TET3*) t;
       t4->t32 = nt;
   }
   else
   {
       t4 = n4->castToTR3(); 
       // t4->t31 = this;   causing warnings, subst. by the next two lines
       t       = (void*) this;
       t4->t31 = (TET3*) t;
       t4->t32 = 0; 
   }

    trs[1] = t1;  trs[2] = t2;  trs[3] = t3;  trs[4] = t4;
}
//-------------------------------------------------------------------------

void TET3:: getTriangles(TR3* trs[])
{ 
   TR3  *t1, *t2, *t3, *t4;
   TET3 *nt;

   //  t1
   if ( (nt = n1->castToTET3()) ) 
   {
       t1 = new TR3;
       t1->boundP = INTERIOR;
       t1->depth = depth;
       t1->p1 = p2;
       t1->p2 = p3;
       t1->p3 = p4;
       t1->e1 = e4;
       t1->e2 = e5;
       t1->e3 = e2;
       completeT(t1,0);
       t1->t31 = this;
       t1->t32 = nt;
   }
   else
   {
       t1 = n1->castToTR3(); 
       t1->t31 = this;
       t1->t32 = 0; 
   }
   
   //  t2
   if ( (nt = n2->castToTET3()) )
   { 
       t2 = new TR3;
       t2->boundP = INTERIOR;
       t2->depth = depth;
       t2->p1 = p1;
       t2->p2 = p4;
       t2->p3 = p3;
       t2->e1 = e4;
       t2->e2 = e1;
       t2->e3 = e6;
       completeT(t2,0);     
       t2->t31 = this;
       t2->t32 = nt;
   }
   else
   {
       t2 = n2->castToTR3(); 
       t2->t31 = this;
       t2->t32 = 0; 
   }

   //  t3
   if ( (nt = n3->castToTET3()) ) 
   {
       t3 = new TR3;
       t3->boundP = INTERIOR;
       t3->depth = depth;
       t3->p1 = p1;
       t3->p2 = p2;
       t3->p3 = p4;
       t3->e1 = e5;
       t3->e2 = e6;
       t3->e3 = e3;
       completeT(t3,0);
       t3->t31 = this;
       t3->t32 = nt;
   }
   else
   {
       t3 = n3->castToTR3(); 
       t3->t31 = this;
       t3->t32 = 0; 
   }
   
   //  t4
   if ( (nt = n4->castToTET3()) ) 
   {
       t4 = new TR3;
       t4->boundP = INTERIOR;
       t4->depth = depth;
       t4->p1 = p1;
       t4->p2 = p3;
       t4->p3 = p2;
       t4->e1 = e2;
       t4->e2 = e3;
       t4->e3 = e1;
       completeT(t4,0);
       t4->t31 = this;
       t4->t32 = nt;
   }
   else
   {
       t4 = n4->castToTR3(); 
       t4->t31 = this;
       t4->t32 = 0; 
   }

    trs[1] = t1;  trs[2] = t2;  trs[3] = t3;  trs[4] = t4;
}
//-------------------------------------------------------------------------
void TET3:: returnFaces(Vector<PATCH*>& trs) const
{
    for (int i=1; i<=4; i++) if (!(trs[i]->onBoundary())) delete trs[i];
}
//-------------------------------------------------------------------------


void TET3:: returnTriangles(TR3* trs[]) const 
{
    for (int i=1; i<=4; i++) if (!(trs[i]->onBoundary())) delete trs[i];
}
//-------------------------------------------------------------------------

void TET3:: getNeighbours(Vector<PATCH*>& neighbours) const
{
    neighbours[1] = n1->castToTET3();
    neighbours[2] = n2->castToTET3();
    neighbours[3] = n3->castToTET3();
    neighbours[4] = n4->castToTET3();
}
//-------------------------------------------------------------------------

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

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

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

// --  		 point at x -> unit coordinates

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

    compJinv(J);

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

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


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

Real TET3:: volume() const  { return oneSixth*Abs(detJ); }

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

int TET3:: NoOfSons() const
{
    if      (type == MESH::T_RED)      return 8;
    else if (type == MESH::T_GREEN_1)  return 2;
    else if (type == MESH::T_GREEN_2A) return 4;
    else if (type == MESH::T_GREEN_2B) return 3;
    else if (type == MESH::T_GREEN_3)  return 4;
    else
    {
	cout << "\n*** TET3::NoOfSons(): unknown type " << type << "\n";
	cout.flush(); abort();
	return 0;
    }
}
//-------------------------------------------------------------------------

Bool TET3:: isGreen()   const
{
    switch(type)
    {
      case MESH::T_GREEN:	
      case MESH::T_GREEN_1:	
      case MESH::T_GREEN_2A:	
      case MESH::T_GREEN_2B:	
      case MESH::T_GREEN_3:     return True;
      default: 	return False; 
    }
}

Bool TR3:: isGreen()   const
{
    switch(type)
    {
      case MESH::T_GREEN:	
      case MESH::T_GREEN_1:	
      case MESH::T_GREEN_2A:	
      case MESH::T_GREEN_2B:	
      case MESH::T_GREEN_3:     return True;
      default: 	return False; 
    }
 }

Bool EDG3:: isGreen()   const
{
    switch(type)
    {
      case MESH::T_GREEN:	
      case MESH::T_GREEN_1:	
      case MESH::T_GREEN_2A:	
      case MESH::T_GREEN_2B:	
      case MESH::T_GREEN_3:     return True;
      default: 	return False; 
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void EFTET3:: reset() 
{ 
    for (int i=0; i<4; ++i) faceNode[i] = 0; 
    TET3::reset(); 
}

void EFTET3:: getFaceNodes(Vector<int>& nodes) const
{
    int i;
    for (i=1; i<=4; ++i) nodes[i] = faceNode[i-1];
}

int& EFTET3:: getFaceNode(int face)  { return faceNode[face-1]; }

int& EFTET3:: getMyFaceNode(PATCH* patch)
{
    if (n1 == patch) return faceNode[0];
    if (n2 == patch) return faceNode[1];
    if (n3 == patch) return faceNode[2];
    if (n4 == patch) return faceNode[3];

    cout << "\n*** EFTET3:: getMyFaceNode: patch not found\n";
    cout.flush(); abort();
	
    return faceNode[0];		// dummy return
}
//-------------------------------------------------------------------------

/*-------------------------------------------------------------------------*/
void MESH3:: setXYZpmCircle(PT3* pm, const EDG3* ed)
{
    Real midPointX, midPointY, midPointZ, radius_new;
    Real hx1, hy1, hz1, hx2, hy2, hz2, radius;
    PT3 *P1 = ed->p1, *P2 = ed->p2, *p = pm;

   //  P1,P2 on surface of the ball  
     
       midPointX = 0.0;
       midPointY = 0.0;
       midPointZ = 0.0;

       hx1 = P1->x-midPointX; hy1 = P1->y-midPointY; hz1 = P1->z-midPointZ;
       radius = sqrt(hx1*hx1+hy1*hy1+hz1*hz1);

       p->x = (P1->x)+((P2->x)-(P1->x))*HALF;
       p->y = (P1->y)+((P2->y)-(P1->y))*HALF;
       p->z = (P1->z)+((P2->z)-(P1->z))*HALF;

       hx2 = p->x-midPointX; hy2 = p->y-midPointY; hz2 = p->z-midPointZ;
       radius_new = sqrt(hx2*hx2+hy2*hy2+hz2*hz2);

       p->x = midPointX + radius*hx2/radius_new;
       p->y = midPointY + radius*hy2/radius_new;
       p->z = midPointZ + radius*hz2/radius_new;

}


/*-------------------------------------------------------------------------*/
void MESH3:: setXYZpmCyl(PT3* pm, const EDG3* ed)
{
    Real radius,hx1,hy1,hx2,hy2;
    Real radius_new;
    PT3 *P1 = ed->p1, *P2 = ed->p2;

    hx1 = P1->x; hy1 = P1->y;
    radius = sqrt(hx1*hx1+hy1*hy1);

    pm->x = (P1->x)+((P2->x)-(P1->x))*HALF;
    pm->y = (P1->y)+((P2->y)-(P1->y))*HALF;
    pm->z = (P1->z)+((P2->z)-(P1->z))*HALF;

    hx2 = pm->x; hy2 = pm->y;
    radius_new = sqrt(hx2*hx2+hy2*hy2);

    pm->x = radius*hx2/radius_new;
    pm->y = radius*hy2/radius_new;

}
//-------------------------------------------------------------------------

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

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


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

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

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

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

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

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

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

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


PATCH* MESH3:: 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* MESH3:: 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* MESH3:: 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;
    }
}