/*
 $Id: triang.cc,v 1.4 1996/10/11 15:15:34 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "triang.h"
#include "triangA.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;

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

Jacobian:: Jacobian() : J(3,3), Jinv(3,3) { }

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


void PATCH:: notImplemented(const char* s) const
{
    cout << "\n*** PATCH:: " << s << " not implemented\n";  cout.flush(); abort();
}
//-------------------------------------------------------------------------

int PATCH:: MLNode(int /*depth*/) const { notImplemented("MLNode"); return 0; }

PATCH* PATCH:: getMidPoint() const { notImplemented("getMidPoint"); return 0; }	

static const Vector<PT*>    pDummyForCompiler(1);
static const Vector<EDG*>   eDummyForCompiler(1);
const Vector<PT*>&  PATCH:: getPoints() const { notImplemented("getPoints()");
						return pDummyForCompiler; }
const Vector<EDG*>& PATCH:: getEdges()  const { notImplemented("getEdges()");
						return eDummyForCompiler; }

static const Vector<PATCH*> fDummyForCompiler(1);
const Vector<PATCH*>& PATCH:: getFaces() const { notImplemented("getFaces()");
						return fDummyForCompiler; }

void PATCH:: getAllFaces (Vector<PATCH*>& /*f*/) const  { notImplemented("getAllFaces"); }
void PATCH:: returnFaces (Vector<PATCH*>& /*f*/) const  { return; }

void PATCH:: getNeighbours(Vector<PATCH*>& /*neighbours*/)  const
					 { notImplemented("getNeighbours");}
void PATCH:: getEdgeOrientation(Vector<Bool>& /*orient*/)  const
					 { notImplemented("getEdgeOrientation");}



int PATCH:: noOfPoints() const { notImplemented("noOfPoints"); return 0; } 
int PATCH:: noOfEdges()  const { notImplemented("noOfEdges");  return 0; }
int PATCH:: noOfFaces()  const { notImplemented("noOfFaces");  return 0; }

int PATCH:: noOfNeighbours()  const 
			{ notImplemented("noOfNeighboursFaces"); return 0;}


float PATCH:: getError() const  { notImplemented("getError"); return 0; }
void  PATCH:: setError(float /*x*/) { notImplemented("setError"); }

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

void PATCH:: getFaceNodes(Vector<int>& /*nodes*/) const 
{
    notImplemented("getFaceNodes");
}
//-------------------------------------------------------------------------

int& PATCH:: getFaceNode(int /*face*/)			// for EFTET3
{
    static int dummy;
    notImplemented("getFaceNode");
    return dummy;
}

int& PATCH:: getMyFaceNode(PATCH* /*patch*/)		// for EFTET3
{
    static int dummy;
    notImplemented("getMyFaceNode");
    return dummy;
}
//-------------------------------------------------------------------------

void PATCH:: realCoordinates(const Vector<Real>& /*u*/, Vector<Real>& /*x*/) 	const
{
    notImplemented("realCoordinates(p1)");
}

void PATCH:: unitCoordinates(const Vector<Real>& /*u*/, Vector<Real>& /*x*/) 	const
{
    notImplemented("unitCoordinates(p1)");
}
//-------------------------------------------------------------------------

void PATCH:: centerOfGravity(Vector<Real>& /*xG*/) const
{
    notImplemented("centerOfGravity");
}
//-------------------------------------------------------------------------

void PATCH:: getNodeCoord(Vector<Real>& /*x*/) const
{
    notImplemented("getCoordinates(p1)");
}
//-------------------------------------------------------------------------

void PATCH:: getCoordinates(Vector<Real>& /*p1*/) 	const
{
    notImplemented("getCoordinates(p1)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/) const
{
    notImplemented("getCoordinates(p1,p2)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/,
			   Vector<Real>& /*p3*/) 	const
{
    notImplemented("getCoordinates(p1,p2,p3)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/, Vector<Real>& /*p3*/,
			   Vector<Real>& /*p4*/) const
{
    notImplemented("getCoordinates(p1,p2,p3,p4)");
}
//-------------------------------------------------------------------------

void PATCH:: compJ(Jacobian& /*j*/, Bool /*checkDetJ*/)    const { notImplemented("compJ"); }
void PATCH:: compJinv(Jacobian& /*j*/) const { notImplemented("compJinv"); }

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

Real PATCH:: volume() const    { notImplemented("volume"); return 0; }

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

Bool PATCH:: onBoundary() const { notImplemented("onBoundary");  return 0; }
Bool PATCH:: isDirichlet()const { notImplemented("isDirichlet"); return 0; }
Bool PATCH:: isNeumann()  const { notImplemented("isNeumann");   return 0; }
Bool PATCH:: isCauchy()   const { notImplemented("isCauchy");    return 0; }
Bool PATCH:: isGreen()    const { notImplemented("isGreen");     return 0; }
Bool PATCH:: isVolatile() const { notImplemented("isVolatile");  return 0; }

PATCH* PATCH:: Father()   const { notImplemented("Father");   return 0; }
PATCH* PATCH:: FirstSon() const { notImplemented("FirstSon"); return 0; }
PATCH* PATCH:: Partner() const { notImplemented("Partner"); return 0; }
void   PATCH:: setPartner(PATCH* /*p*/) { notImplemented("setPartner"); }

int PATCH:: NoOfSons() const { notImplemented("NoOfSons"); return 0; }

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


Real PATCH:: hMax() const
{
    int i;
    Real h, maxH = 0.0;

    const Vector<EDG*>& edges = getEdges();
    
    FORALL(edges,i)
    {
	h = edges[i]->hMax();
	maxH = Max(h,maxH);
    }
    return maxH;
}
//-------------------------------------------------------------------------

void PATCH:: setMark() 	    { notImplemented("setMark"); }
void PATCH:: unMark()	    { notImplemented("unMark"); }
Bool PATCH:: marked() const  { notImplemented("marked"); return 0; }

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

const PATCH* PATCH:: findPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/) const
{
    notImplemented("findPatch"); return 0;
}
Bool PATCH:: inPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/) const
{
    notImplemented("inPatch");  return 0;
}
 
PATCH* PATCH:: findPartner(const Vector<Real>& /*xG*/, int /*maxDepth*/)
{
    notImplemented("findPartner"); return 0;
}
//-------------------------------------------------------------------------


void PATCH:: readBCValues(char* boundP, char* classA, BufferedParser& parser)
{
    char BCType;

    parser.getValue(&BCType);
    BCType = toupper(BCType);

    //if (BCType == 'I') { *boundP = INTERIOR; 	*classA = 0; }
    // else
    if 	    (BCType == 'D') *boundP = DIRICHLET;
    else if (BCType == 'N') *boundP = NEUMANN;
    else if (BCType == 'C') *boundP = CAUCHY;
    else {
	    cout << "\n*** readBC: Unknown Boundary Type " << BCType << "\n";
	    cout.flush(); abort();
	}

    int cl;
    parser.getValue(&cl);	
    if (cl <= 0)
	{
	    cout << "\n*** Boundary condition class must be > 0 (read: " 
	         << cl << ")\n";
	    cout.flush(); abort();
	}
    *classA = cl;			// int -> char
}


void PATCH:: readBCValues(char* boundP, char* classA, char* isoP, BufferedParser& parser)
{
    char BCType;

    parser.getValue(&BCType);
    BCType = toupper(BCType);
    //if (BCType == 'I') { *boundP = INTERIOR; 	*classA = 0; }
    // else
    if 	    (BCType == 'D') *boundP = DIRICHLET;
    else if (BCType == 'N') *boundP = NEUMANN;
    else if (BCType == 'C') *boundP = CAUCHY;
    else {
	    cout << "\n*** readBC: Unknown Boundary Type " << BCType << "\n";
	    cout.flush(); abort();
	}

    int cl;
    parser.getValue(&cl);	
    if (cl <= 0)
	{
	    cout << "\n*** Boundary condition class must be > 0 (read: " 
	         << cl << ")\n";
	    cout.flush(); abort();
	}
    *classA = cl;			// int -> char


    char is;
    if(parser.getSpecValue(&is) == parser.ok)
    {
	is = toupper(is);

	if 	(is == 'S') *isoP   = CIRCULAR;
	else if (is == 'C') *isoP   = CYLINDRICAL;
	else                *isoP   = PLANE;
    }	
    else *isoP   = PLANE;
}

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


//--------------------------  base class MESH  ----------------------------

 
MESH:: MESH(const char* inFileName)  :  varAllocator(ElementsInBlock*100),
					elemIter(0), ptIter(0), edgIter(0)
{
    fileName = new char[strlen(inFileName)+1];
    strcpy(fileName, inFileName);

    maxDepth = 0;
    refStep  = 0;

    innerBoundary   = False;	if (Cmd.isTrue("innerBoundary")) innerBoundary=True;

    infoRefinement = 0;     Cmd.get("infoRefinement", &infoRefinement); 
    timeRefinement = 0;     Cmd.get("timeRefinement", &timeRefinement);
    accTimeRefinement = 0;  Cmd.get("accTimeRefinement", &accTimeRefinement);
    accTime = 0.0;
}
//-------------------------------------------------------------------------

MESH:: ~MESH() { delete fileName; }

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

void MESH:: notImplemented(const char* s) const
{
    cout << "\n*** MESH::" << s << " not implemented\n";  cout.flush(); abort();
}
//-------------------------------------------------------------------------

void MESH:: Flatten() { notImplemented("Flatten"); }

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

void MESH:: writeTriangulation(const char* /*outFileText*/,const Vector<Real>& /*u*/) 
									   const
{
    notImplemented("writeTriangulation");
}
//-------------------------------------------------------------------------


void MESH:: missingItem(const char* item, const char* fileName0)
{
    cout << "\n*** File " << fileName0 << " without " << item << "\n";
    cout.flush(); abort(); 
}
//-------------------------------------------------------------------------

void MESH:: pointIndexError(int no)
{
    cout << "\n*** readTriangulation: point index too high: " << no << "\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void MESH:: checkMaxIndex(int maxIndex)
{
    if (maxIndex < 0) {
	cout << "\n*** readTriangulation: maxIndex must be >= 0 (read: "
	     << maxIndex << ")\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

void MESH:: boundaryWarning(int pNo1, int pNo2, int pNo3)
{
    cout << "\n*** Warning --- readTriangulation: boundary element specified"
         << "\n\tthat is not on exterior boundary; points: " << pNo1;

    if (pNo2) cout << pNo2;
    if (pNo3) cout << pNo3;
    cout << "\n";
}
//-------------------------------------------------------------------------

void MESH:: timeInfo(Timer& timer, Timer& accTimer)
{
    if (timeRefinement)  { cout << "\n\tGrid Refinement: "; timer.cpu(); }
    if (accTimeRefinement) 
    {
	accTime += accTimer.cpu(False);
	cout << Form("\tAccumulated time: %1.2f sec.\n", accTime); 
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MESH:: resetEdgIter()
{
    if (edgIter != 0)
    {
	cout << "\n*** resetEdgIter: edgIter != 0"
	     << " \t (previous operation not completed?)\n";
	cout.flush(); abort();
    }
    edgIterLevel = 0;
}
//-------------------------------------------------------------------------


void MESH:: resetPtIter()
{
    if (ptIter != 0)
    {
	cout << "\n*** resetPtIter: ptIter != 0"
	     << " \t (previous operation not completed?)\n";
	cout.flush(); abort();
    }
    ptIterLevel = 0;
}
//-------------------------------------------------------------------------


const PATCH* MESH:: findPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/, 
			      const PATCH* /*newPatch*/) const
{
    notImplemented("findPatch");  return 0;
}
//-------------------------------------------------------------------------

void MESH:: patchNotFound() const
{
    cout << "\n*** MESH: patch not found\n";     cout.flush(); abort();
}
//-------------------------------------------------------------------------


// --			set the partners of depth 0

void MESH:: setInitialPartners(MESH* prevMesh)
{
    PATCH* patch, *partner;

    if (prevMesh==0) return;

    prevMesh->resetElemIter();
    resetElemIter();
    
    while ((patch=elemIterAll()))
    {
	partner = prevMesh->elemIterAllOfHistory();
	patch->setPartner(partner);
    }

    // --	  dummy operation on rest of elements:

    while ((partner=prevMesh->elemIterAllOfHistory())) { } 
}
//-------------------------------------------------------------------------


void MESH:: setPartners()
{
    Timer  timer;
    int    k, noOfSons;
    PATCH* fPartner, *p, *son, *father, *fPartnerSon;

    const int lastStep = True;
    resetElemIter(lastStep);


    while ((p = elemIterAllOfHistory()))
    {
	if (p->Partner() == 0)
	{
	    father  = p->Father();
	    fPartner = father->Partner();
	    son = p;
	    noOfSons = father->NoOfSons();

	    if (father->isGreen() || fPartner->isGreen() || 
		!fPartner->refined() || fPartner->Depth() < father->Depth())
	    {
		for (k=1; k<=noOfSons; ++k)
		{
		    son->setPartner(fPartner);
		    son = son->Next();
		}
	    }	
	    else
	    {
		fPartnerSon = fPartner->FirstSon();

		for (k=1; k<=noOfSons; ++k)
		{
		    son->setPartner(fPartnerSon);
		    son = son->Next();
		    fPartnerSon = fPartnerSon->Next();
		}
	    }
	}
    }

    if (Cmd.isTrue("TimeTransport")) 
      { cout << "\n\tFind Partners: "; timer.cpu(); }
}
//-------------------------------------------------------------------------


void MESH:: newNodeStack(PT* p, int depth, int targetDepth)
{
    static const int nodeStackSize = sizeof(NodeStack);

    if (depth >= targetDepth) 
    {
	cout << "\n***MESH:: newNodeStack: depth >= targetDepth (" 
	     << depth << ">=" << targetDepth << ")\n";
	cout.flush(); abort();
    }
    p->nodeStack = (NodeStack*) varAllocator.Get(nodeStackSize);
    p->nodeStack->init(&varAllocator, depth, targetDepth);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void NodeStack:: init(VarSizeAllocator* varAlloc0, int depth, int targetDepth)
{
    static const int intSize = sizeof(int);

    varAlloc = varAlloc0; 

    l   = depth;
    h   = depth-1;
    top = targetDepth;

    v  = (int*) varAlloc->Get(intSize*(top-l+1));
    v -= l;
}
//-------------------------------------------------------------------------

void NodeStack:: push(int node)
{
    if (h==top) extend();
    v[++h] = node;
}
//-------------------------------------------------------------------------

void NodeStack:: extend()
{
    static const int intSize = sizeof(int);

    top += top;
    int* vnew = (int*) varAlloc->Get(intSize*(top-l+1));

    vnew -= l;
    for (int i=l; i<=h; ++i) vnew[i] = v[i];

    // delete v;   	was allocated from varAlloc
    v = vnew;
}
//-------------------------------------------------------------------------

void NodeStack:: checkBounds(int i) const
{
    if (i<l || i>h) 
    {
	cout << "\n*** Vector index " << i << " out of range\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

void NodeStack:: print() const
{
    int i;
    for (i=l; i<=h; ++i)  cout << v[i] << " ";
    cout << "\n";
}