/*
 $Id: triang1.cc,v 1.5 1996/11/06 12:26:19 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "triang1.h"
#include "utils.h"
#include "numerics.h"
#include "mzibutil.h"
#if NOKASKPLOT == 0
#include "feplot1.h"
#endif

#include "cmdpars.h"
extern CmdPars Cmd;

Vector<PT*>      EDG1::p(2);
Vector<EDG*>     EDG1::e(1);
Vector<PATCH*>   EDG1::f(2);


MESH1:: MESH1(const char* inFileName, Bool readMesh)
		 
      : MESH(inFileName), ptList(0,25), edgList(0,25), 
	ptAlloc(ElementsInBlock), edgAlloc(ElementsInBlock)
{
    noOfPoints = 0;
    noOfEdges  = 0;

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

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

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


void MESH1:: readTriangulation(const char* fileName0)
{
    int  i, maxIndex;

    DList<PT1>*   ptL = new DList<PT1>;    ptList.push(ptL);
    DList<EDG1>* edgL = new DList<EDG1>;  edgList.push(edgL);

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


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

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

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

    readPoints(parser, ptAdr);


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


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


    EDG1* e;
    DListIter<EDG1> iter(edges());
    while ((e=iter.all())) e->setBoundary();


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

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

    while (parser.getValue(&ptNo) == parser.ok)  
    {
	if (ptNo > ptAdr.high()) pointIndexError(ptNo);
	
	p = ptAlloc.Get();			//new PT1;

	ptAdr[ptNo] = p;
	parser.getValue(&p->x);
	p->node = ++count;			// used in Find...

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

void MESH1:: readElements(BufferedParser& parser, Vector<PT1*>& ptAdr)
{
    int  pNo1, pNo2, classA;
    EDG1 *e;

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

	e = getEDG1();			// new EDG1;
	e->p1 = ptAdr[pNo1];
	e->p2 = ptAdr[pNo2];
	setRightHand(e);

	parser.getValue(&classA);
	e->classA = classA;
	SetEP(e);

	edgList[0]->add(e);
	++noOfEdges;
    }
}
//-------------------------------------------------------------------------

void MESH1:: writeTriangulation(const char* outFileText, const Vector<Real>& u)
									const 
{
   int  i;
   int  node;
   EDG1 *ed;
   PT1  *p;
   FILE *outFile;

   Vector<int> PNodes(noOfPoints);

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


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

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


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

   iterPts.reset();
   while ((p=iterPts.all()))  {
      if (!p->onBoundary()) continue;
      if (p->isDirichlet())
                  fprintf(outFile, "%d d %d\n",
                  p->node, p->classA );
      else if (p->isNeumann())
                  fprintf(outFile, "%d n %d\n",
                  p->node, p->classA );
      else if (p->isCauchy())
                  fprintf(outFile, "%d c %d\n",
                  p->node, p->classA );
   }

   fprintf(outFile, "END\n");


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


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

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

void MESH1:: setRightHand(EDG1* e)
{
    if (e->p2->x < e->p1->x)
    {
	PT1* p = e->p2;
	e->p2 = e->p1;
	e->p1 = p;
    }
}
//-------------------------------------------------------------------------

void MESH1:: readBoundary(BufferedParser& parser, Vector<PT1*>& ptAdr) 
{
    int pNo;
    PT1 *p;

    while (parser.getValue(&pNo) == parser.ok)	
    {
	p = ptAdr[pNo];
	if (!p->onBoundary()) boundaryWarning(pNo);
	p->readBC(parser);
    }
}
//-------------------------------------------------------------------------


void MESH1:: Refine()
{
    EDG1* e;
    int noPt[2], noEdg[2];
    Timer timer, accTimer;

    noPt[0]  = noOfPoints;
    noEdg[0] = noOfEdges;
    
    
    DList<EDG1>* edgL = new DList<EDG1>;  
    edgList.push(edgL);

    DListIter<EDG1> iter(edges());
    while ((e=iter.all()))  { if (e->mark) RefineEdge(e); }

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

    ++refStep;

    if (infoRefinement) refinementInfo (noPt, noEdg);
    timeInfo(timer, accTimer);

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


void MESH1:: Flatten()
{
    int noReleasedEdges=0;
    int noPoints=0, noEdges=0;
    int k;

    EDG1* ed;
    PT1*  p;
    DList<EDG1>* edgL = new DList<EDG1>;
    DList<PT1>*  ptL  = new DList<PT1>;

    Stack<EDG1*> edgToRetain(25), edgToRelease(25);
    Stack<PT1*> ptToRetain(25);

    DListIter<EDG1> iterEdg(edges());
    DListIter<PT1>  iterPt(points());

    while ((ed=iterEdg.allOfHistory()))
    {
	if (!ed->refined())
	{
	    edgToRetain.push(ed);
	    noEdges++;
	}
	else
	{
	    edgToRelease.push(ed);
	    noReleasedEdges++;
	}
    }

    while ((p=iterPt.allOfHistory()))
    {
	ptToRetain.push(p);
	noPoints++;
    }

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

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

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

    for (k=edgToRetain.l; k<=edgToRetain.h; k++)
    {
	ed = edgToRetain[k];
	ed->father = nil;
	ed->depth = 0;
	edgL->add(ed);
    }

    while (!edgToRelease.empty())
    {
	ed = edgToRelease.pop();
	edgAlloc.Return(ed);
    }

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

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


Bool MESH1:: RefineEdge(EDG1 *ed)
{
    EDG1 *e1, *e2;
    PT1  *pm;

    e1 = getEDG1();
    e2 = getEDG1();	
    ++noOfEdges;

    pm = ptAlloc.Get();			// new PT1;  assert(pm != 0);	
    ++noOfPoints;

    ed->pm = pm;
    pm->depth  = ed->depth + 1;
    pm->x = 0.5*(ed->p1->x + ed->p2->x);

    e1->p1 = ed->p1; 
    e1->p2 = pm;

    e2->p1 = pm; 
    e2->p2 = ed->p2;

    e1->father = e2->father = ed; 
    e1->type   = e2->type   = ed->type;
    e1->classA = e2->classA = ed->classA;
    e1->depth  = e2->depth  = ed->depth + 1;

    ed->firstSon = e1; 
    DelEP(ed);
    SetEP(e1);
    SetEP(e2);

    edgList[edgList.h]->add(e1); 	
    edgList[edgList.h]->add(e2); 	

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

    if (e1->depth > maxDepth) maxDepth = e1->depth;

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

// --		 SetEP(e):  set edge e on its points


void MESH1:: SetEP(EDG1 *e) 
{
    PT1 *p1, *p2;;
    int k = 0;

    p1 = e->p1;
    if      (p1->e1 == 0) { p1->e1 = e; ++k; }
    else if (p1->e2 == 0) { p1->e2 = e; ++k; }

    p2 = e->p2;
    if      (p2->e1 == 0) { p2->e1 = e; ++k; }
    else if (p2->e2 == 0) { p2->e2 = e; ++k; }

    if (k != 2)
    {
	cout << "\n*** SetEP: Edge with impossible state, k = " << k << "\n";
        printf("\n e->p1: %p",e->p1);
        printf("\n e->p2: %p",e->p2);
	e->print(cout);
	cout.flush(); abort();
    }

    if (p1->onBoundary() || p2->onBoundary()) e->boundP = BOUNDARY;
}
//-------------------------------------------------------------------------

// --		 DelEP(e):  delete edge e on its points


void MESH1:: DelEP(EDG1 *e) 
{
    PT1* p;
    int k = 0;

    p = e->p1;
    if      (p->e1 == e) { p->e1 = 0; ++k; }
    else if (p->e2 == e) { p->e2 = 0; ++k; }

    p = e->p2;
    if      (p->e1 == e) { p->e1 = 0; ++k; }
    else if (p->e2 == e) { p->e2 = 0; ++k; }

    if (k != 2)
    {
	cout << "\n*** DelEP: impossible state, k = " << k << ")\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

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

int MESH1:: MemSpace()
{
    int space = varAllocator.MemSpace();
    space += edgAlloc.MemSpace();
    space += ptAlloc.MemSpace();
    return space;
}
//-------------------------------------------------------------------------

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

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

    int space = MemSpace();

    /*
      FORALL (edgList,depth)
      {
      space += edgList[depth]->noOfElements * sizeof(EDG1);
      if (depth <= ptList.h) 
      space += ptList[depth]->noOfElements * sizeof(PT1);
      }
      */

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

    if (infoRefinement > 1)
    {
	cout << Form("\n%15s%10d%10d", "generated", 
		       noPt[1]-noPt[0], noEdg[1]-noEdg[0]);
    }

    cout << Form("\n%15s%10d%10d   Space : %1.6f MB","total",
		   noOfPoints, noOfEdges, float(space)/1e6);

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

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

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

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

ostream& operator << (ostream& os, const MESH1& t)
{
    PT1*  p;
    EDG1* e;

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

    DListIter<PT1>  iterPT (t.points());
    DListIter<EDG1> iterEDG(t.edges());

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

    os << "\nRefinement step " << t.refStep;
    os << "  --   points: " << t.noOfPoints;
    os << "  edges: " << t.noOfEdges;
    os << "\n";
    return os;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

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

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

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

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

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

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


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

void EDG1:: compJ(Jacobian& j, Bool /*checkDetJ*/)  const
{
    j.J(1,1) = p2->x - p1->x;
    j.detJ = j.J(1,1);
}

void EDG1:: compJinv(Jacobian& j) const
{
    compJ(j);
    j.Jinv(1,1) = 1./j.J(1,1);
}
//-------------------------------------------------------------------------

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

void EDG1:: reset()  
{ 
    mark = classA = type = depth = NUL;
    boundP = INTERIOR;
    p1 = p2 = pm = 0;
    node = 0; 
    father = firstSon = 0;
    error = 0.0;
}

void EDG1:: setBoundary()
{
    if (p1->onBoundary() || p2->onBoundary()) boundP = BOUNDARY;
}
//-------------------------------------------------------------------------



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

    if (p2->e1 != this) neighbours[2] = p2->e1;
    else	        neighbours[2] = p2->e2;
}
//-------------------------------------------------------------------------

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

void EDG1:: 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; 
}
//-------------------------------------------------------------------------

// --  		 point at x -> unit coordinates

void EDG1:: unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const
{
    u[1] = (x[1] - p1->x) / (p2->x - p1->x);
}
//-------------------------------------------------------------------------

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

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

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

Real EDG1:: hMax() const { return Abs(p2->x - p1->x); }

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


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

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


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

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

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

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

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

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

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

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


PATCH* MESH1:: 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* MESH1:: 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* MESH1:: 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;
    }
}