/*
 $Id: int.cc,v 1.6 1997/07/11 09:46:11 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "int.h"

#include "utils.h"
#include "family.h"
#include "connect.h"

#include "triang.h"
#if NOKASKPLOT == 0
#include "feplot.h"
#endif
#include "elements.h"

#include "linsystem.h"
#include "precond.h"

#include "cmdpars.h"
extern CmdPars Cmd;


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

#if NOKASKPLOT == 0
Stack<FEPlot*> Interface:: FEPlotStack(25);
FEPlot*        Interface:: fePlot = 0;
#endif
//-------------------------------------------------------------------------


Interface:: Interface(MESH* mesh0, const Element* element0, 
		      DirichletBCs* dirichletBCs0, LinSystem* Ab0,
		      Preconditioner* precond0, int spaceDim0, int nComp0)

	 : mesh(mesh0), element(element0), dirichletBCs(dirichletBCs0),
  	   Ab(Ab0), precond(precond0), 
	   familyTree(0),
	   noOfNodes(-1,DefaultStackSize), pointNodes(-1,DefaultStackSize),
	   spaceDim(spaceDim0), nComp(nComp0)
{

    noOfNodes.push(0);			// no node on level '-1'
    pointNodes.push(0);			// no node on level '-1'

    timeUpdate = 0;	 	Cmd.get("timeUpdate",   &timeUpdate);
    accTime    = 0;	 	Cmd.get("accTimeUpdate",&accTime);
    accTimeUpdate = 0.0;

    plotSize = 0.5;  		Cmd.get("plotSize", &plotSize);
    plotKeep = False;  		Cmd.get("plotKeep", &plotKeep);

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

Interface:: ~Interface() 
{ 
    delete familyTree;
#if NOKASKPLOT == 0
    if (plotKeep) FEPlotStack.push(fePlot);
    // else  delete fePlot;    
    // fePlot = 0;
#endif
}

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

void Interface:: notImplemented(const char* name) const
{
    cout << "\n*** class Interface: function " << name << " not implemented\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void Interface:: interpolateSolution(Vector<Num>& /*u*/) const
		       		      { notImplemented("interpolateSolution"); }
void Interface:: updateMGFamilyTree() { notImplemented("updateMGFamilyTree"); }
void Interface:: getGlobalMLNodes(const PATCH* /*t*/, Vector<int>& /*globalNodes*/,
				  int /*depth*/) const
				      { notImplemented("getGlobalMLNodes"); }
void Interface:: setMLNodeNumbers(int /*maxDepth*/, int /*maxDepthM1*/, 
				  int /*targetDepth*/)
		       		      { notImplemented("setMLNodeNumbers"); }
void Interface:: updateMLFamilyTree(int /*maxDepth*/, int /*maxDepthM1*/)
		       		      { notImplemented("updateMLFamilyTree"); }
void Interface:: interpolateMLSolution(Vector<Num>& /*u*/,
				       int /*maxDepth*/, int /*maxDepthM1*/) const
		       		      { notImplemented("interpolateMLSolution");}

void Interface:: setHighOrderNodes(int* /*edNode*/)    
					{ notImplemented("setHighOrderNodes"); }

void Interface:: solToNB(Vector<Num>& /*u*/) { notImplemented("solToNB"); }
void Interface:: rhsToNB(Vector<Num>& /*b*/) { notImplemented("rhsToNB"); }
void Interface:: solToHB(Vector<Num>& /*u*/) { notImplemented("solToHB"); }
void Interface:: rhsToHB(Vector<Num>& /*b*/) { notImplemented("rhsToHB"); }

void Interface:: setHBGeneration(Generation& /*gen*/) 
					{ notImplemented("setHBGeneration"); }
 
void Interface:: getNodeCoordinates(NodeCoordinates& /*nc*/)  const
					{ notImplemented("getNodeCoordinates"); }
//-------------------------------------------------------------------------

void Interface:: print()
{
    cout << "\nInterface:\n";
    if (familyTree) familyTree->print();
}
//-------------------------------------------------------------------------

int Interface:: MaxDepth() const { return mesh->MaxDepth(); }

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

void Interface:: timeInfo(Timer& timer, Timer& accTimer)
{
    if (timeUpdate) { cout << "\n\tUpdate: "; timer.cpu();}
    if (accTime) 
    {
	accTimeUpdate += accTimer.cpu(False);
	cout << Form("\tAccumulated time: %1.2f seconds\n", accTimeUpdate); 
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Interface:: updatePrecond()
{
    precond->update(A, familyTree, dirichletBCs);

    // --		      infos:

    if (Cmd.isTrue("infoAb"))  Ab->memSpaceInfo();
    if (Cmd.isTrue("printAb")) 
    	{ cout << "Ab assembled with DBc:\n"; Ab->print(); }
    // if (Cmd.isTrue("printMatLabFormat"))  Ab->A->printMatLabFormat();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Interface:: plot(Vector<Num>& u, int step, const char* fileName, 
		      int timeStep)
{
    if (Cmd.isTrue("graphics"))
    { 
        if (spaceDim == 3)
    	{
	   cout << "\n*** No graphic for spaceDim 3 available\n";
	   return;
    	}

	if (Cmd.isTrue("plotPointNodes")) 
	{
	    keepOrClearPlot(fileName, step, timeStep);
#if NOKASKPLOT == 0
	    fePlot->plotPointNodes(); 
#endif
	    Pause(); 
	}

	if (Cmd.isTrue("plotEdgeNodes")) 
	{
	    keepOrClearPlot(fileName, step, timeStep);
#if NOKASKPLOT == 0
	    fePlot->plotEdgeNodes(); 
#endif
	    Pause(); 
	}

	if (Cmd.isTrue("plotTriangleNodes") && spaceDim==2) 
	{
	    keepOrClearPlot(fileName, step, timeStep);
#if NOKASKPLOT == 0
	    fePlot->plotTriangleNodes();
#endif
	    Pause(); 
	}	    

	keepOrClearPlot(fileName, step, timeStep);
#if NOKASKPLOT == 0
	plotSolution(u, fePlot);
#endif
    }
}
//-------------------------------------------------------------------------


void Interface:: post(Vector<Num>& u, int /*step*/, const char* fileName, 
		      int timeStep)
{
 #if NOKASKPLOT == 0
   static int noPostScript = 0;
    FEPlot* fePost=0;

    char s[10], psFile[256];
    strcpy (psFile, fileName);
    
    if (strchr(psFile,'.'))  strcpy(strchr(psFile,'.'), "--");
    if (timeStep > 0) sprintf(strstr(psFile,"--"), "-%1d--", timeStep);
    
    
    // --		      interactive

    // cout  << "Enter 'y' for mesh plot: "; cout.flush(); fgets(s,10,stdin);
    cout  << "Postscript picture of mesh plot: ";

    s[0] = 'y';
    if (strchr(s,'y') || strchr(s,'Y')) 
    {
	strcpy(strstr(psFile,"--") + 2, "mesh.ps");
	completePostFileName(psFile, &noPostScript);
	
	fePost = newFEPlot(PS, psFile);
	fePost->plotElements(); 

	if (Cmd.isTrue("plotBoundary"))    fePost->plotBoundary();
	delete fePost;
    }
    
    // cout  << "Enter 'y' for solution plot: "; cout.flush(); fgets(s,10,stdin);
    cout  << "Postscript picture of solution plot: ";
    s[0] = 'y';
    if (strchr(s,'y') || strchr(s,'Y')) 
    {
	strcpy(strstr(psFile,"--") + 2, "sol.ps");
	completePostFileName(psFile, &noPostScript);
	
	fePost = newFEPlot(PS, psFile);
	plotSolution(u, fePost); 
	delete fePost;
    }
#endif
}
//-------------------------------------------------------------------------

   
// 			non-interactive:


void Interface:: autoPost(Vector<Num>& u, int /*step*/, const char* fileName, 
			  int timeStep)
{
#if NOKASKPLOT == 0
    static int noPostScript = 0;
    FEPlot* fePost=0;
    
    char psFile[256];
    strcpy (psFile, fileName);
    
    if (strchr(psFile,'.'))  strcpy(strchr(psFile,'.'), "--");
    if (timeStep > 0) sprintf(strstr(psFile,"--"), "-%1d--", timeStep);
    
    
    strcpy(strstr(psFile,"--") + 2, "mesh.ps");
    completePostFileName(psFile, &noPostScript);
    
    fePost = newFEPlot(PS, psFile);
    fePost->plotElements(); 
    delete fePost; 
    
    strcpy(strstr(psFile,"--") + 2, "sol.ps");
    completePostFileName(psFile, &noPostScript);
    
    fePost = newFEPlot(PS, psFile);
    plotSolution(u, fePost); 
    delete fePost;
#endif
}
//-------------------------------------------------------------------------


void Interface:: keepOrClearPlot(const char* fileName, int step, int timeStep)
{
#if NOKASKPLOT == 0
    char caption[MaxLineLength];
    
    if (plotKeep) sprintf(caption, "%s:  Ref.Step %1d", fileName, step);
    else 	  sprintf(caption, "%s:  Ref.Step varying", fileName);

    if (timeStep >= 0) 
    {
	sprintf(caption+strlen(caption), "  --  Time Step");
	if (plotKeep) sprintf(caption+strlen(caption), " %1d", timeStep);
    }

    if  (fePlot == 0) fePlot = newFEPlot(SCREEN, caption, plotSize);
    else if (plotKeep) 
    {
	FEPlotStack.push(fePlot);
	fePlot = newFEPlot(SCREEN, caption, plotSize);
    }
    else 
    {
	fePlot->updateMesh(mesh);
	fePlot->clear();
    }
#endif
}
//-------------------------------------------------------------------------


void Interface:: setSolutionVector(Vector<Real>& u1, const Vector<Num>& u)
{
    int i, node;

    if      (Cmd.isTrue("plotImag"))  FORALL(u1,i) u1[i] = ::imag(u[i]);
    else if (Cmd.isTrue("plotPhase")) FORALL(u1,i) u1[i] = ::arg (u[i]);
    else if (Cmd.isTrue("plotabs"))   FORALL(u1,i) u1[i] = ::Abs (u[i]);
    else		              FORALL(u1,i) u1[i] = ::real(u[i]);
 
    if (nComp > 1)	  	// manipulate u1 to trick plot-routines
    {
	int n = 1;	  Cmd.get("PlotComponent", &n);

	if (n<1 || n>nComp)
	  { cout << "\n*** Invalid PlotComponent: " << n << "\n"; cout.flush(); abort(); }

	for (node=1; n<=u1.high(); n+=nComp)
	{
	    for (i=1; i<=nComp; ++i) u1[node++] = u1[n];
	}
    }
}
//-------------------------------------------------------------------------


void Interface:: plotSolution(Vector<Num>& u, FEPlot* localFEPlot)
{
#if NOKASKPLOT == 0
    Vector<Real> u1(u.high());

    setSolutionVector(u1, u);

    if (spaceDim==1) 
    {
	if (Cmd.isTrue("plotElements"))    localFEPlot->plotElements();
	if (Cmd.isTrue("plotIsoLines"))    localFEPlot->plotSolution(u1);
	return;
    }

    if (spaceDim==2)  
    {
	if (Cmd.isTrue("plot3D"))  localFEPlot->plot3D(u1);
        else
	{
	 if (Cmd.isTrue("plotElements"))    localFEPlot->plotElements();
	 if (Cmd.isTrue("plotBoundary"))    localFEPlot->plotBoundary();
	 if (Cmd.isTrue("plotIsoLines"))    localFEPlot->plotSolution(u1);
	}
        return;
    }
#endif
}
//-------------------------------------------------------------------------

void Interface:: write(Vector<Num>& u, int /*step*/, const char* fileName) 
{ 
    Vector<Real> u1(u.high());

    setSolutionVector(u1, u);
    mesh->writeTriangulation(fileName, u1);
}

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


void Interface:: completePostFileName(char* psFile, int* noPostScript)
{
    FILE* fp;

    while (1)
    {
	++(*noPostScript);
	sprintf(strstr(psFile,".ps")+3, ".%1d", *noPostScript);

	fp = fopen(psFile,"r");
	if (fp) fclose(fp);
	else break;
    }
}
//-------------------------------------------------------------------------

/* void Interface:: compressPointNodes()
   {
   int    node;
   PATCH* p;
   
   mesh->resetPtIter();
   while (p=mesh->ptIterAll()) 
   {
   node = p->getNode();
   MCNode.baseNode(node, nComp);
   p->setNode(node);
   }
   }
   void Interface:: expandPointNodes()
   {
   int    node;
   PATCH* p;
   
   mesh->resetPtIter();
   while (p=mesh->ptIterAll()) 
   {
   node = p->getNode();
   node = MCNode.node(node, 1, nComp);
   p->setNode(node);
   }
}
*/
//-------------------------------------------------------------------------


FEPlot* Interface:: newFEPlot(int plotType, char* caption, float size) 
{ 
#ifndef NOASKPLOT
    if (spaceDim == 3)
    {
	cout << "\n*** FEPlot for spaceDim 3 not available\n";
	return 0;
    }

    return mesh->newFEPlot(plotType, caption, size);
#else
    return 0;
#endif
}