/* $Id: triang3.h,v 1.3 1996/10/11 15:15:49 bzferdma Exp $ */

#ifndef TRIANG3_H
#define TRIANG3_H

#include "triangA.h"
#include "alloc.h"


class PT3;
class EDG3;
class TR3; 
class TET3;

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


class PT3 : public PT
{
  public:

    char  boundP, mark, depth, classA;
    PT3   *next, *prev;
    Real  x,y,z;

    PT3() { reset(); }
    virtual void reset();

    virtual PATCH* Next() { return next; }
    virtual void setMark() { mark = 1; }
    virtual void unMark()  { mark = 0; }
    virtual Bool marked() const { return mark; }

    virtual int spaceDim() const { return 3; }
    virtual int Class() const { return classA; }
    virtual int Depth()    const { return depth; }

    virtual Bool isDirichlet() const { return boundP==DIRICHLET; }

    virtual void getNodeCoord(Vector<Real>& v) const;
    virtual void getCoordinates(Vector<Real>& v) const;
    virtual void print(ostream& os) const;

    void readBC(BufferedParser& parser);
};
//-------------------------------------------------------------------------


class EDG3 : public EDG
{
  private:
    static Vector<PT*>  p;

  public:

    char  boundP, type, depth, mark, classA, isoP;
    PT3  *p1, *p2, *pm;
    EDG3 *next, *prev, *father, *firstSon;

    EDG3() { reset(); }
    virtual void reset();

    virtual PATCH* Next() { return next; }
    virtual void setMark() { mark = 1; }
    virtual void unMark()  { mark = 0; }
    virtual Bool marked() const { return mark; }

    virtual int spaceDim() const { return 3; }

    virtual int Class() const { return classA; }
    virtual int Depth() const { return depth; }

    virtual PATCH* getMidPoint() const { return pm; }

    virtual PATCH* Father() const { return father; }

    virtual Bool isDirichlet() const { return boundP==DIRICHLET; }
    virtual Bool refined() const { return Bool(firstSon); }
    virtual Bool isGreen() const;

    virtual void compNormalVector(Vector<Real>& normal, Real* length) const;

    virtual void realCoordinates(const Vector<Real>& u, Vector<Real>& x) const;

    virtual void getNodeCoord(Vector<Real>& v) const;
    virtual void getCoordinates(Vector<Real>& p1, Vector<Real>& p2) const;

    virtual const Vector<PT*>& getPoints() const { p[1]=p1; p[2]=p2; return p; }

    virtual Real hMax() const;

    virtual void print(ostream& os) const;
};
//-------------------------------------------------------------------------


class TR3 : public TR
{
  private:
    static Vector<PT*>    p;
    static Vector<EDG*>   e;
    static Vector<PATCH*> f;

  public:

    PT3		*p1,*p2,*p3;
    EDG3	*e1,*e2,*e3;
    TR3		*next, *prev, *father, *firstSon;
    TET3  	*t31, *t32;
    char	mark, type, depth, boundP, classA, isoP;

    static StaticAllocator<TR3> alloc;

    void* operator new(size_t /*size*/) { return alloc.Get(); }
    void  operator delete(void* tr) { alloc.Return((TR3*) tr); }

    TR3() { reset(); }
    virtual ~TR3() {  }

    virtual void reset();

    virtual PATCH* Next() { return next; }
    virtual int spaceDim() const { return 3; }

    virtual void setMark() { mark = 1; }
    virtual void unMark()  { mark = 0; }
    virtual Bool marked() const { return mark; }

    virtual int Class() const { return classA; }
    virtual Bool refined() const { return Bool(firstSon); }
    virtual Bool isGreen()   const;
    virtual int Depth()    const { return depth; }

    virtual PATCH* Father() const { return father; }

    virtual Bool onBoundary() const  { return (t31==0 || t32==0); }
    virtual Bool isDirichlet() const { return boundP==DIRICHLET; }
    virtual Bool isNeumann() const   { return boundP==NEUMANN; }
    virtual Bool isCauchy() const    { return boundP==CAUCHY; }

    virtual Real volume() const;

    virtual void realCoordinates(const Vector<Real>& u, Vector<Real>& x) const;
    virtual void getCoordinates(Vector<Real>& p1, Vector<Real>& p2,
				Vector<Real>& p3) const;

    virtual const Vector<PT*>& getPoints() const
				{ p[1]=p1; p[2]=p2; p[3]=p3; return p; }
    virtual const Vector<EDG*>& getEdges() const
				{ e[1]=e1; e[2]=e2; e[3]=e3; return e; }
    virtual const Vector<PATCH*>& getFaces() const
				{ f[1]=e1; f[2]=e2; f[3]=e3; return f; }

    virtual void getEdges(EDG3* eds[]) const;

    virtual void print(ostream& os) const;

    void readBC(BufferedParser& parser);
    void readHyperBC(BufferedParser& parser);

  protected:

    virtual TR3*  castToTR3() { return this; }
};
//-------------------------------------------------------------------------


class TET3 : public TET
{
  private:
    static Vector<PT*>    p;
    static Vector<EDG*>   e;
    static Vector<PATCH*> f;

  public:

    Real	detJ;
    PT3  	*p1, *p2, *p3, *p4;
    EDG3  	*e1, *e2, *e3, *e4, *e5, *e6;
    PATCH   	*n1, *n2, *n3, *n4;
    TET3   	*next, *prev, *father, *firstSon;
    float 	error;
    char        type, depth, mark, classA, boundP;

    TET3() { reset(); }
    virtual void reset();

    virtual PATCH* Next() { return next; }
    virtual int  spaceDim() const { return 3; }
    virtual Real volume() const;
    virtual Bool orient() const { return detJ >= 0.0; }

    virtual void setMark() { mark = 1; }
    virtual void unMark()  { mark = 0; }
    virtual Bool marked() const { return mark; }

    virtual int  Class()   const { return classA; }
    virtual int  Depth()   const { return depth; }
    virtual Bool refined() const { return Bool(firstSon); }
    virtual Bool isGreen()   const;

    virtual void setBoundary();
    virtual Bool  onBoundary() const { return boundP==BOUNDARY; }

    virtual PATCH* Father()   const { return father; }
    virtual PATCH* FirstSon() const { return firstSon; }
    virtual int NoOfSons() const;

    virtual void getPoints(PT3* pts[]) const; 

    virtual const Vector<PT*>& getPoints() const
				{ p[1]=p1; p[2]=p2; p[3]=p3; p[4]=p4; return p; }
    virtual const Vector<EDG*>& getEdges() const
				{ e[1]=e1; e[2]=e2; e[3]=e3; e[4]=e4; e[5]=e5; e[6]=e6; return e; }

    virtual const Vector<PATCH*>& getFaces() const
		 { f[1]=n1->castToTR3(); f[2]=n2->castToTR3(); 
		   f[3]=n3->castToTR3(); f[4]=n4->castToTR3(); return f; }

    virtual void getAllFaces(Vector<PATCH*>& trs)  const;  
    virtual void returnFaces(Vector<PATCH*>& trs)  const;  

    virtual void getEdgeOrientation(Vector<Bool>& orient)  const;

    virtual void getNeighbours(Vector<PATCH*>& neighbours) const;

    virtual void realCoordinates(const Vector<Real>& u, Vector<Real>& x) const;
    virtual void unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const;

    virtual void getCoordinates(Vector<Real>& p1, Vector<Real>& p2,
				Vector<Real>& p3, Vector<Real>& p4) const;
    virtual void compJ(Jacobian& J, Bool checkDetJ=True)  const;
    virtual void compJinv(Jacobian& J) const;

    virtual void getEdges(EDG3* eds[]) const;
    virtual void getTriangles(TR3* trs[]);  
    virtual void returnTriangles(TR3* trs[]) const;  

    void getMarkedTriangles(int trsMark[]);  
    void markTriangles(TET3* tet);  
    Bool noTdinNb(TET3* tet);  
    void correctNbinTet(TET3* tet);  

    virtual const PATCH* findPatch(const Vector<Real>& x, Vector<Real>& xUnit) 
    									const;
    virtual Bool inPatch(const Vector<Real>& x, Vector<Real>& xUnit) const;

    virtual void print(ostream& os) const;

    virtual float getError() const { return error; }
    virtual void  setError(float x) { error = x; }

  protected:

    virtual TET3* castToTET3() { return this; }
    void completeT(TR3* t, TR3* f) const;
};

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


class EFTET3 : public TET3
{
  public:

    int  faceNode[4];

    EFTET3() { reset(); }
    virtual void reset();

    virtual void getFaceNodes(Vector<int>& nodes) const;
    virtual int& getFaceNode(int face);
    virtual int& getMyFaceNode(PATCH* patch);
};
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


class MESH3 : public MESH
{

  public:

    MESH3(const char* inFileName, Bool readMesh=True);
    ~MESH3();

    virtual MESH3* castToMESH3() { return this; }

    void writeTriangulation(const char* outFileText,const Vector<Real>& u) const;
    void writeZIBFormat   (const char* outFileText,const Vector<Real>& u) const;
    void writeGrapeFormat (const char* outFileText,const Vector<Real>& u) const;
    void writeAVSFormat   (const char* outFileText,const Vector<Real>& u) const;
    void writeHyperFormat (const char* outFileText,const Vector<Real>& u) const;

    virtual void Refine();

    virtual int noOfElements() { return noOfTetrahedra; }
    virtual int spaceDim() const { return 3; }

    const Stack<DList<PT3>*>&  points()    const { return ptList; }
    const Stack<DList<EDG3>*>& edges ()    const { return edgList; }
    const Stack<DList<TR3>*>&  triangles() const { return trList; }
    const Stack<DList<TET3>*>& tetras()    const { return tetList; }

    virtual void Info();
    virtual int MemSpace();


  protected:

    int	    noOfPoints, noOfEdges, noOfTriangles, noOfTetrahedra;

    enum   { ElementsInBlock = 250 };

    enum RefinementType { ShortDiag, BeyStabilDiag };
    RefinementType	refinementType;

    Bool    rightHand;

    int     moreRedTd, moreGreen1Td, moreGreen2ATd, moreGreen2BTd, moreGreen3Td;
    int     alpha[36];
    Real    alphaMax, alphaMin, sigmaMax;

    Stack<DList<PT3>*>  ptList;   
    Stack<DList<TR3>*>  trList;
    Stack<DList<EDG3>*> edgList;
    Stack<DList<TET3>*> tetList;

    Allocator<PT3>   ptAlloc;
    Allocator<EDG3>  edgAlloc;
    Allocator<TET3>  tetAlloc;

    virtual TET3* getTET3() 	        { return tetAlloc.Get(); }
    virtual void  returnTET3(TET3* tet) { tetAlloc.Return(tet); }

    virtual const PATCH* findPatch(const Vector<Real>& x, Vector<Real>& xUnit,
				   const PATCH* newPatch)  const;
    

    // --		element iterator

    virtual void  resetElemIter(Bool lastStep=False);
    virtual PATCH* elemIterAll();
    virtual PATCH* elemIterAllOfHistory();

    virtual PATCH* ptIterAll();

    virtual PATCH* edgIterAll();
    virtual PATCH* edgIterAllOfHistory();

    
    void readTriangulation(const char* fileName);
    void readZIBFormat    (const char* fileName);
    void readOldZIBFormat (const char* fileName);
    void readHyperFormat  (const char* fileName);
    void readGrapeFormat  (const char* fileName);

    void readPoints  (BufferedParser& parser, Vector<PT3*>& ptAdr);
    void readElements(BufferedParser& parser, Vector<PT3*>& ptAdr);
    void readBoundary(BufferedParser& parser, Vector<PT3*>& ptAdr, 
		      Vector<Stack<TR3*>*>& trsOfPts);

    void readHyperPoints  (BufferedParser& parser, int pointsSection, 
			   Vector<PT3*>& ptAdr);
    void readHyperElements(BufferedParser& parser, int tetrahedronsSection, 
			   Vector<PT3*>& ptAdr);
    void readHyperBoundary(BufferedParser& parser, 
			   int trianglesSection, int triangleDataSection,
			   Vector<PT3*>& ptAdr, Vector<Stack<TR3*>*>& trsOfPts);
    void readHyperClass(BufferedParser& parser, int tetrahedronDataSection);
    
    int   RefineTetra(TET3* t);
    int   MakeGreen();
    void  RemoveGreen();
    void  refinementInfo(int* noPt,int* noEdg,int* noTr,int* noTet);

    EDG3* FindEdgByPts(PT3*p1,PT3 *p2, Vector<Stack<EDG3*>*>& edsOfPts);
    TR3*  FindTrByEds(PT3 *p1,PT3 *p2,PT3 *p3,EDG3 *e1,EDG3 *e2,EDG3 *e3,
		      Vector<Stack<TR3*>*>& trsOfPts);
    int   completeTetrahedron(TET3* td, Vector<Stack<EDG3*>*>& edsOfPts,
			      Vector<Stack<TR3*>*>& trsOfPts);
    Bool  pointOnEdge(PT3* p, EDG3* ed);
    Bool  edgeOnTri(EDG3* ed, TR3* t);

    EDG3* FindEdgByPts(PT3*p1,PT3 *p2);
    TR3*  FindTrByEds(PT3 *p1,PT3 *p2,PT3 *p3,EDG3 *e1,EDG3 *e2,EDG3 *e3);
    int   CompleteTd(TET3* td);

    int   CloseRedTet(TET3* tet);
    int   RefineTdShortDiag(TET3* td);
    int   RefineTdBey(TET3* td);
    int   RefineTr(TR3* t);
    EDG3* RefineEdg(EDG3* ed);
    int   TdVolume(TET3* tet);
    void  SetTP(TR3* t,TR3* f);
    void  SetTP(TET3 *tet,TR3* t, int noP1, int noP2, int noP3); 
    void  SetTdP(TET3* td,TR3* t1,TR3* t2,TR3* t3,TR3* t4,TET3* f);
    void  SetTdPBey(TET3* td,TET3* f);
    void  SetTPts(TET3* td,TET3* f);
    EDG3* FindTdE(TET3* td,TR3* t1,TR3* t2,TR3* t3,TR3* t4,PT3* p1,PT3* p2);
    void  SetTdE(TET3 *td,TR3* t1,TR3* t2,TR3* t3,TR3* t4);
    void  InnerEdges(TR3* t,EDG3* e1,EDG3* e2,EDG3* e3,EDG3* ie1,EDG3* ie2,
		     EDG3* ie3);
    int   MkGreen(TR3* t);
    int   MkGreenIrr(TR3* t);
    int   MarkGreenTd(TET3* td);
    int   MkGreen1Td(TET3* td);
    int   MkGreen2ATd(TET3* td);
    int   MkGreen2BTd(TET3* td);
    int   MkGreen3Td(TET3* td);
    void  RelGreen();
    void  RelGreenIrr();
    void  RelGreen1Td();
    void  RelGreen2ATd();
    void  RelGreen2BTd();
    void  RelGreen3Td();

    virtual int consistencyCheck();
    int   CheckTd(TET3 *td);
    int   CheckAngle(TR3* t);
    int   CheckCircle(TET3* t);
    int   CompAngles();

    int   testSons(TR3 *tr, TET3 *tet);
    Real  Length(PT3 *p1,PT3* p2);
    int   PinTr(PT3* p,TR3* tr);
    int   PinTet(PT3* p,TET3* tet);
    EDG3* EdinTr(TR3* tr,PT3* p1,PT3* p2);
    EDG3* EdinTet(PT3* p1,PT3* p2,TET3* td);
    int   TestFace(TR3* tr);
    int   TestEdsInTd(TET3* tet);

    int    faceNoInTd(TR3* t, TET3* tet);
    void   neighbourAtFace(TET3* td, PATCH* N, TR3* t);
    void   innerNeighbourAtFace(TET3* td, TR3* t, TET3* neighbourTd);
    void   setNeighbourInTd(TET3* td, PATCH* n, int noOfNeighbour);
    EDG3*  findEdgInNeighbour(PT3* p1, PT3* p2, TET3* n);

    void   setNeighbsOfTetra(TR3 *t);
    void   makeFacesOfTet(TET3 *td, TR3** t1, TR3** t2, TR3** t3, TR3** t4);
    void   correctWhiteNbs(TET3 *tet);

    void   setTypeOfInnerEdges(TR3 *t);

    void setXYZpmCircle(PT3* pm, const EDG3* e);
    void setXYZpmCyl(PT3* pm, const EDG3* e);

    friend ostream& operator << (ostream& os, MESH3& t);
};
//-------------------------------------------------------------------------


class EFMESH3 : public MESH3
{
  public:

    EFMESH3(const char* inFileName, Bool readMesh=True);
    virtual ~EFMESH3() { }

  protected:

    Allocator<EFTET3>  EFTetAlloc;

    virtual TET3* getTET3() 	        { return (TET3*) EFTetAlloc.Get(); }
    virtual void  returnTET3(TET3* tet) { EFTetAlloc.Return((EFTET3*) tet); }

    int MemSpace();
};


#endif