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

#include "problem.h"
#include "physics.h"

#include "linsystem.h"
#include "precondsg.h"
#include "precondmg.h"

#include "triang.h"
#include "elements.h"
#include "int.h"

#include "dirichlet.h"
#include "materials.h"
#include "adapt.h"
#include "nodeco.h"

#include "cmdpars.h"
extern CmdPars  Cmd;

extern ostream* infoFile;

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


Problem:: Problem (char* fileName0, symmetryType symmetry0, int spaceDim0,
		   int nComp0) : 

	symmetry(symmetry0), spaceDim(spaceDim0), nComp(nComp0),
	Fct(0,0), Energy(0,0), Error(0,0)
{ 
    fileName	  = 0; 

    material	  = 0;
    dirichletBCs  =  dirichletBCsDLY = 0;
    precond	  = 0;
    Ab 		  = 0;
    element	  = 0;
    interface 	  = 0;
    errorEstimator= 0;
    statistic     = new Statistic();

    fileName = new char[strlen(fileName0)+10];
    strcpy(fileName, fileName0);


    	   // inquire the problem-specific parameters:

    biSection = Cmd.isTrue("biSection");

    infoSolution = True;  	Cmd.get("infoSolution", &infoSolution);

    plotSolution   = False;	Cmd.get("plotSolution", &plotSolution);
    postScript     = False;	Cmd.get("postScript", &postScript);
    writeSolution  = False;	Cmd.get("writeSolution",&writeSolution);
    compareSolution  = False;	Cmd.get("compareSolution",&compareSolution);

    timeAssembler    = 0;  Cmd.get("timeAssembler", &timeAssembler);
    accTimeAssembler = 0;  Cmd.get("accTimeAssembler", &accTimeAssembler);

    maxRefSteps =  50;	Cmd.get("maxRefSteps", &maxRefSteps);
    minRefSteps =   0;	Cmd.get("minRefSteps", &minRefSteps);

    maxRefDepth = 50;   Cmd.get("maxRefDepth", &maxRefDepth);

    globConvTest = estimRelError;
}
//-------------------------------------------------------------------------


Problem:: ~Problem() 
{ 
    delete fileName; 

    delete precond;
    delete material;
    delete dirichletBCs; delete dirichletBCsDLY;
    delete Ab; 

    delete element;
    delete interface;

    delete errorEstimator;
}
//-------------------------------------------------------------------------

// NonLinearity* Problem:: getNonLinearity() const { return 0; }

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

void Problem:: newLinSystem()
{
    if (precond      == 0) missingObject("newLinSystem","precond");
    if (dirichletBCs == 0) missingObject("newLinSystem","dirichletBCs");

    delete Ab;
    Ab = new LinSystem(spaceDim, symmetry, precond, dirichletBCs);
}
//-------------------------------------------------------------------------

void Problem:: newPreconditioner()
{
    delete precond;

    if      (Cmd.isSet("precond","jacobi")) precond = new Jacobi(); 
    else if (Cmd.isSet("precond","sgs"))    precond = new SGS();  
    else if (Cmd.isSet("precond","ssor"))   precond = new SSOR();
    else if (Cmd.isSet("precond","TRssor")) precond = new TrSSOR();
    else if (Cmd.isSet("precond","ILU"))    precond = new ILU();

    else if (Cmd.isSet("precond","MGsgs"))  precond = new MGSGS();
    else if (Cmd.isSet("precond","MGssor")) precond = new MGSSOR();
    else if (Cmd.isSet("precond","MGjac"))  precond = new MGJacobi();
    else if (Cmd.isSet("precond","MGCG"))   precond = new MGCG();

    else if (Cmd.isSet("precond","HB"))     precond = new HB();
    else if (Cmd.isSet("precond","AMGjac")) precond = new AMGJacobi();
    else if (Cmd.isSet("precond","AMGsgs")) precond = new AMGSGS();

    else if (Cmd.isSet("precond","MLJacobi")) precond = new MLJacobi();
    else if (Cmd.isSet("precond","MLsgs"))    precond = new MLSGS();
    else if (Cmd.isSet("precond","MLssor"))   precond = new MLSSOR();

    else if (Cmd.isSet("precond","AMLJacobi"))	precond = new AMLJacobi();
    else if (Cmd.isSet("precond","BPX"))	precond = new AMLJacobi();
    else if (Cmd.isSet("precond","AMLsgs"))    	precond = new AMLSGS();

    else if (Cmd.isSet("precond","none"))   precond = new DummyPreconditioner();
    else MissingParameter("Preconditioner");
}
//-------------------------------------------------------------------------

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

void Problem:: missingObject(const char* function, const char* object) const
{
    cout << "\n*** " << function << ": " << object << " missing\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void Problem:: missingMaterialTerm(const char* name, int classVal) const
{
    cout << "\n*** assemble: " << name << " value for material class " 
         << classVal << " missing\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

Bool Problem:: globalConvergenceTest(Real globalPrecision, int step)
{
    switch (globConvTest)
    {
      case estimRelError: return estimRelErrorConvTest(globalPrecision, step);
      case estimAbsError: return estimAbsErrorConvTest(globalPrecision, step);

      //case extrapolError: return extrapolatedErrorConvTest(globalPrecision, 
      //						       step);
      default: cout << "\n*** Problem: no valid convergence Test specified\n";
			  cout.flush(); abort(); return 0;
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Problem:: contractVector(Vector<Num>& newVec, Vector<int>& oldPos, 
			      Vector<Num>& oldVec)
{
    int i;
    if (newVec.high() < oldVec.high())
      { cout << "\nProblem.contract: inconsistent vector sizes \n"; cout.flush(); abort(); }

    FORALL(newVec,i)  { if (oldPos[i]) oldVec[oldPos[i]] = newVec[i]; }
}

void Problem:: expandVector(Vector<Num>& oldVec, Vector<int>& oldPos, 
			    Vector<Num>& newVec, Num fill)
{
    int i;
    if (newVec.high() < oldVec.high())
      { cout << "\nProblem.expand: inconsistent vector sizes \n"; cout.flush(); abort(); }

    FORALL(newVec,i)  
    { 
	if (oldPos[i]) newVec[i] = oldVec[oldPos[i]];
	else 	       newVec[i] = fill; 
    }
}
//-------------------------------------------------------------------------


Bool Problem:: estimRelErrorConvTest(Real globalPrecision, int step)
{
/* 
    if (equal(Energy[step],0.0)) return False;

    Real relError = Error[step]/Energy[step];

    if (Cmd.isTrue("writeCpuTime")) *infoFile << "\t " << relError << "\n";

    if (relError < globalPrecision) return True;
    else return False;
*/

    Bool ConvOk =   (Error[step] < globalPrecision * Energy[step])||
                    (Error[step]==0); 
   
    if (Cmd.isTrue("writeCpuTime")) 
      { if (equal(Energy[step],0.0)) 
	  {
	    if (ConvOk) *infoFile << "\t < " << globalPrecision << "\n";    
	    else  *infoFile << "\t  " << machMax(Real(0.0)) << "\n"; 
	  }
	else
	  {
	    Real relError = Error[step]/Energy[step]; 
	    *infoFile << "\t " << relError << "\n";
	  }
      }
    return ConvOk;
}
//-------------------------------------------------------------------------


Bool Problem:: estimAbsErrorConvTest(Real globalPrecision, int step)
{
    if ( (Error[step] < globalPrecision) || (Error[step]==0) ) 
         return True;
    else return False;
}
//-------------------------------------------------------------------------


Bool Problem:: extrapolatedErrorConvTest(Real globalPrecision, int step)
{
    Real ratio, expo, extrapolatedError = machMax(Real(0));
    static int dimM1;
    Bool status = False;

    if (step > 0)
    {
	extrapolatedError = Abs(Energy[step]-Energy[step-1]);

	expo = 2.0/Real(spaceDim);
	ratio = Real(dimM1)/Real(u.size());  
	ratio = pow(ratio,expo);

	extrapolatedError *= ratio/(1.0-ratio); 

	if ( (extrapolatedError < globalPrecision*Energy[step]) ||
	     (extrapolatedError ==0) )  status = True;
    }

    if (Cmd.isTrue("writeCpuTime")) 
    	*infoFile << "\t " << quot(extrapolatedError,Energy[step]) << "\n";

    dimM1 = u.size();
    return status;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Problem:: SolutionInfo(int step)
{
    static int i, dimM1;
    static Timer accTimer;
    static Real accTime = 0.0;
    Real error, extrapolatedError=0.0, ratio=0.0, expo;

    if (step > 0)
    {
	extrapolatedError = 
	  Abs(quot(Energy[step]-Energy[step-1],Energy[step]));

	expo = 2.0/Real(spaceDim);	
	ratio = Real(dimM1)/Real(u.size());  
	ratio = pow(ratio,expo);

	extrapolatedError *= 100.0 * ratio/(1.0-ratio);	

	ratio = quot(Error[step],Error[step-1]);
    }	

    if (!equal(Energy[step],0.0)) error = 100.0*Error.Top()/Energy.Top();
    else 			  error = Error.Top();
    
    cout << Form("\n %7s %5s %6s %11s","RefStep", "level", "nodes", "Energy");

    if (!equal(Energy[step],0.0))  cout << Form(" %8s","Error(%)");
    else        		  cout << Form(" %8s","Error");

    if (step > 0) cout << Form(" %7s %8s", "ratio", "dE(%)");

    cout << Form(" %9s", "cpu(sec)");
    // cout << Form(" %9s", "mem(MB)");


    cout << Form("\n %7d %5d %6d %11.6g", step, interface->MaxDepth(), 
		   interface->Dim(), Energy.Top());
    cout << Form(" %8.3g", error);

    if (step > 0)
    {
	cout << Form(" %7.2g %8.3g", ratio,
		       quot(100.0*(Energy[step]-Energy[step-1]),Energy[step]));
    }


    accTime = accTimer.cpu(False);
    cout << Form(" %9.2f", accTime);


    //int space = mesh->MemSpace();  //  + Ab->MemSpace() + precond->MemSpace();
    //cout << Form(" %9.6f", float(space)/1e6);

    /*
       cout << "\n Step " << step << ": E=" << Form("%1.5g", Energy[step]);
       if (!equal(Energy[step],0.0))
       cout << " Err=" << Form("%1.3g %%", 100.*Error[step]/Energy[step]);
       cout << Form(" (ext=%1.3g)", 100*extrapolatedError);
       cout << Form(" ratio(Err)=%1.2f",(Error[step] / Error[step-1]));
       cout << Form(" dE(%1d)=%1.3g", step-1, 
       100*(Energy[step]-Energy[step-1])/Energy[step]) << " %";
       */

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

    cout << "\n ";  for (i=1; i<=76; ++i) cout << '-';  cout << "\n";

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

    if (Cmd.isTrue("EnergyNorm")) cout << "\n E-Norm:  " << EnergyNorm(u) <<"\n";
    if (Cmd.isTrue("L2Norm"))     cout << "\n L2-Norm: " << L2Norm(u) << "\n";

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


    dimM1 = u.size();

    while ( Pause() == picture )
	interface->post (u, step, fileName);    
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Problem:: compEnergy(Vector<Num>& sol, Real* energy, Num* fct)
{
    Bool print=True;
    Ab->compEnergy(sol, energy, fct, print);
}
//-------------------------------------------------------------------------


Real Problem:: L2Norm(const Vector<Num>& x) const 
{
    int  dim;
    Real L2 = 0.0;
    PATCH* patch;

    dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<int> globalNodes(dim); 
    Jacobian Jac;
    
    Mesh()->resetElemIter();

    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem);

	element->assembleL2Norm(*patch, AElem, Jac, 0);
	interface->getGlobalNodes(patch, globalNodes);  

	L2 += energyProd(AElem, x, globalNodes);
    }

    // L2 = sqrt(2.0*L2);		// annihilate factor 1/2 in energy
    return L2;
}
//-------------------------------------------------------------------------


Real Problem:: L2MassNorm(const Vector<Num>& x) const 
{
    int  dim;
    Real L2 = 0.0;
    PATCH* patch;

    dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<int> globalNodes(dim); 
    Jacobian Jac;
    
    Mesh()->resetElemIter();

    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem);

	element->assembleMass(*patch, AElem, Jac, 0);
	interface->getGlobalNodes(patch, globalNodes);  

	L2 += energyProd(AElem, x, globalNodes);
    }

    // L2 = sqrt(2.0*L2);		// annihilate factor 1/2 in energy
    return L2;
}
//-------------------------------------------------------------------------


Real Problem:: EllipticNorm(const Vector<Num>& x) const  
{
    int  dim;
    Real E = 0.0;
    PATCH* patch;

    dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<Num> bElem(dim);
    Vector<int> globalNodes(dim); 
    Jacobian Jac;
    
    Mesh()->resetElemIter();

    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem,bElem);

	element->assembleEllip(*patch, AElem, Jac, 0);
	interface->getGlobalNodes(patch, globalNodes);  
	
	E += energyProd(AElem, x, globalNodes);
    }
    return E;
}
//-------------------------------------------------------------------------


Real Problem:: EnergyNorm(const Vector<Num>& x) 
{
    int  dim;
    Real E = 0.0;
    PATCH* patch;

    dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<Num> bElem(dim);
    Vector<int> globalNodes(dim); 
    Jacobian Jac;
    
    Mesh()->resetElemIter();

    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem,bElem);

	assemble(*element, *patch, Jac, AElem, bElem, 0, True);
	interface->getGlobalNodes(patch, globalNodes);  
	
	E += energyProd(AElem, x, globalNodes);
    }
    return E;
}
//-------------------------------------------------------------------------

Real Problem:: energyProd(Matrix<Num>& AElem, const Vector<Num>& sol, 
			  const Vector<int>& globalNodes) const 
{
    int i, j, row, col;
    Num sum;
    Num E = 0.0;
    
    const int dim = globalNodes.high();

    if (symmetry == sym)				// upper triangle
    {
	for (i=1; i<=dim; ++i)
	for (j=1; j<i;    ++j) AElem(j,i) = AElem(i,j);
    }
    
    for (i=1; i<=dim; ++i)
    {
	if ((row=globalNodes[i]))
	{
	    sum = 0.0;
	    for (j=1; j<=dim; ++j)
	      if ((col=globalNodes[j])) sum += AElem(i,j) * sol[col];

	    E += conj(sol[row]) * sum; 
	}
    }
    return 0.5*Abs(E);

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


Bool Problem:: solve(Real relGlobalPrecision)
{
    Bool absPrecision=False;
    if (Cmd.isTrue("absPrecision")) absPrecision=True;

    return staticSolution(relGlobalPrecision, absPrecision);
}
//-------------------------------------------------------------------------
void Problem:: compare(Vector<Num>& v,  Real time)
{
    if (material->trueSolKnown())
      {
	Real maxError=0.0, error=0.0;
	Num trueSol=0.0;
	int i;
	const int dim      = interface->Dim();
//	const spaceDim = interface->SpaceDim();

	Vector<Real> x(spaceDim);
	NodeCoordinates nc(spaceDim, 1, dim);

	interface->getNodeCoordinates(nc);

	FORALL(v,i)
	  {
	    nc.getCoordinate(i,x);
	    trueSol = material->trueSolInPoint(x,time);
	    error = Abs(v[i] - trueSol);
	    if (error > maxError) maxError = error;
	  }

	cout << "\n\tCompare with true solution: Max. error at mesh points=" 
	  << maxError << "\n";
      }
    else cout << "\n Can not compare, since true solution is not known\n";
}
//-------------------------------------------------------------------------


Bool Problem:: staticSolution(Real globalPrecision, Bool absPrecision, Real time)
{
    int i, step;
    Num fct;
    Real error=0.0, energy=1.0, relPrecision;

    if (absPrecision) globConvTest = estimAbsError;    
    else	      globConvTest = estimRelError;   

    Energy.resize(0,DefaultStackSize);
    Error.resize (0,DefaultStackSize);
    Fct.resize   (0,DefaultStackSize);

    newPreconditioner();
    newLinSystem();
    newMesh();
    newInterface();
    newErrorEstimator();

    interface->updateLevel0(u);			

    for (step=0; True; ++step)
    {
	if (Cmd.isTrue("uReset"))  FORALL(u,i) u[i] = 0.0;

	assembleGlobal(time);

       // cout << "\n\tNo of Dirichlet points = " 
       //      << dirichletBCs->noOfConstraints() << "\n";

	interface->updatePrecond(); 

	if (absPrecision && !(energy==0)) relPrecision = globalPrecision/energy;
	else		  relPrecision = globalPrecision;

	if (statistic->active)
	{
	    statistic->ZD_IntWrite(statistic->idRefStep,step);
	    statistic->ZD_IntWrite(statistic->idDepth,interface->MaxDepth());
	    statistic->ZD_IntWrite(statistic->idNodes,interface->Dim());
	}


        Ab->statistic = statistic;

	Ab->solve(u, energy, error, relPrecision, step);

	compEnergy(u, &energy, &fct); 
 
	Energy.push(energy);
	Fct.push(fct);  

	if (plotSolution)       interface->plot (u, step, fileName);
	if (postScript)         interface->post (u, step, fileName);
	if (writeSolution)      interface->write(u, step, fileName);
	if (compareSolution)    compare(u,time);

	if (absPrecision)
	  errorEstimator->adapt(*this, &error, globalPrecision);
	else 
	  errorEstimator->adapt(*this, &error, relPrecision*energy);
	Error.push(error);
	
	if (statistic->active)
	{
	    statistic->ZD_RealWrite(statistic->idDiscErr,error);
	}

	if (infoSolution) SolutionInfo(step);

	if (step >= minRefSteps)
	{
	    if (globalConvergenceTest(globalPrecision, step)) return True;

	    if (step == maxRefSteps) 
	    { 
		cout << "\n* maxRefStep " << maxRefSteps << " reached\n";
		return False; 
	    }
	    if (interface->MaxDepth() > maxRefDepth)
	    {
		cout << "\n* maxRefDepth " << maxRefDepth << " exceeded\n";
		return False; 
	    }
	}

	interface->refine(u);
    }
}
//-------------------------------------------------------------------------


void Problem:: assembleGlobal(Real time)
{
    static Real accTime = 0.0;
    Timer timer, accTimer;

    PATCH* patch;
    const int dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<Num> bElem(dim);
    Vector<int> nodes(dim);
    Jacobian Jac;

    Ab->reset();

    Mesh()->resetElemIter();
    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem,bElem);
	  

	assemble(*element, *patch, Jac, AElem, bElem, 0);
  

	interface->getGlobalNodes(patch, nodes);
	Ab->store(AElem, bElem, nodes);
    }

    interface->updateDirichletBCs(time);
    Ab->setDirichletBCs(*dirichletBCs);


    // -- 	 		infos:

    if (timeAssembler)  { cout << "\n\tGlobal Assembling: "; timer.cpu(); }
    if (accTimeAssembler)
    {
	accTime += accTimer.cpu(False);
	cout << Form("\tAccumulated time:  %1.2f sec.\n", accTime); 
    }
}
//-------------------------------------------------------------------------


void Problem:: resolveMesh(Real lambda0)
{
    Timer timer;
    int    nElem;
    PATCH* patch;
    Real   lambda, h, eps, mu;

    lambda0 *= 0.999;

    while(1)
    {
	nElem = 0;
	Mesh()->resetElemIter();

	while ((patch = Mesh()->elemIterAll()))
	{
	    mu  = 1/material->E(patch->Class());
	    eps =   material->M(patch->Class());
	    
	    h = patch->hMax();
	    lambda = lambda0 / sqrt(eps*mu);
	    
	    if (h > lambda) { patch->setMark(); ++nElem; }
	    else  patch->unMark();
	}

	if (nElem == 0) break;

	Mesh()->Refine();
    } 

    if (spaceDim != 3) Mesh()->Flatten();   // !!!

    if (timeAssembler)  { cout << "\n\tResolve Mesh: "; timer.cpu(); }
}