/*
 $Id: adaptnl.cc,v 1.2 1996/10/04 15:06:22 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "adaptnl.h"

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

#include "problemnl.h"
#include "nonlin.h"
#include "dirichletA.h" // !!!

#include "elements.h"
#include "triang.h"
#include "int.h"
#include "linsystem.h"
#include "precondnl.h"
#include "sysmatml.h"
#include "family.h"

#include "cmdpars.h"
extern CmdPars Cmd;

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

RK1:: RK1(Element* elementDLY0, Element* elementLag0, Interface* interface0)  

    : DLY(elementDLY0, interface0), lagrangeElementDLY(elementLag0), 
      doRefine(1)
{ }
//-------------------------------------------------------------------------

RK1:: ~RK1() { delete lagrangeElementDLY; }

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


Bool RK1:: estimateError(Problem& problem)
{
    if (infoErrorEstimator) cout << "\n\tError Est.: RK1";

    distributeError(problem);
    enforcedPMediaRefinement(problem);

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


void RK1:: distributeError(Problem& problem)
{
    int i;
    Vector<Num> error(1), uQ(1), ADiag(1);

    solveQuadDefectProblem(problem, error, uQ, ADiag);
    
    interfaceDLY->solToHB(error);

    EnergyError = 0.0;
    for (i=problem.u.high()+1; i<=error.high(); ++i) 
    {
	error[i] = error[i] * error[i] * ADiag[i];	// error on edge
	EnergyError += Abs(error[i]);                   // sum up global Error
    }
    EnergyError *= 0.5;

    distributeEdgeErrors(problem, error);
}
//-------------------------------------------------------------------------


void RK1:: solveQuadDefectProblem(Problem& problem, Vector<Num>& error, 
				  Vector<Num>& uQ, Vector<Num>& ADiag)
{
    int i, j;
    Num sum;
    PATCH* patch;

    MESH* mesh = problem.Mesh();
    Vector<Num>&  u = problem.u;

    const int edNodeM1 = u.high();
    int   edNode   = edNodeM1;

    interfaceDLY->setHighOrderNodes(&edNode);
    
    error.resize(edNode);
    ADiag.resize(edNode);
    uQ.resize(edNode);

    FORALL(error,i) error[i] = ADiag[i] = 0.0;


    Vector<Num> uPrevOnNewMesh(1);
    TransientProblem* trP = 0;

    if ((trP = problem.castToTransientProblem()))
    {
	uPrevOnNewMesh.resize(edNodeM1+1, edNode);
	trP->solutionToNewMesh(uPrevOnNewMesh, interfaceDLY);
    }

    const int low  = problem.element->NoOfNodes()+1;
    const int high = elementDLY->NoOfNodes();

    Matrix<Num> AElem(high, high);
    Vector<Num> bElem(high);

    Vector<int> globalNodes(high);
    Jacobian Jac;

    Vector<Real> LNorm(high);
    Vector<Real> lagrangeAreaWeights(error.high());
    FORALL(lagrangeAreaWeights,i) lagrangeAreaWeights[i] = 0.0;

    Matrix<int> FullDiag(high,high), LowDiag(high,high);

    FORALL_ROWS(FullDiag,i)
    FORALL_COLUMNS(FullDiag,j) FullDiag(i,j) = LowDiag(i,j) = False;

    for (i=1; i<=high; ++i) FullDiag(i,i) = True;
    for (i=1; i< low; ++i)  LowDiag(i,i) = True;


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

	problem.assemble(*elementDLY, *patch, Jac, AElem, bElem, 
			 elementDLY->DLYPattern, True);

	interfaceDLY->getGlobalNodes(patch, globalNodes);  

	for (i=low; i<=high; ++i)     			// sum up diagonal AQQ
	{				     		//   and right hand side
	    if (globalNodes[i] == 0) continue;
	    ADiag[globalNodes[i]] += AElem(i,i);
	    error[globalNodes[i]] += bElem[i];

	    sum = 0.0;
	    for (j=1; j<low; ++j)  	// subtract 'coupling terms' from rhs
	    {
		if (globalNodes[j])  sum += AElem(i,j) * u[globalNodes[j]];
	    }
	    error[globalNodes[i]] -= sum;
	}


	problem.assemble(*lagrangeElementDLY, *patch, Jac, AElem, bElem, 
			 &LowDiag, True);


	for (i=1; i<low; ++i)  // sum up order 1 diagonal terms of Lagrange-Elem.
	{ 
	    if (globalNodes[i])  ADiag[globalNodes[i]] += AElem(i,i);
	}

	FORALL(LNorm,i) LNorm[i] = 0.0;

	lagrangeElementDLY->assembleLNorm(*patch, LNorm, Jac, &FullDiag);

	for (i=1; i<=high; ++i) 
	  if (globalNodes[i]) lagrangeAreaWeights[globalNodes[i]] += LNorm[i];
    }
    

    if (trP)	     // transient problem: contrib. of previous step to rhs
    {
	const NonLinearity* nonLin = problem.getNonLinearity();

	FORALL(uPrevOnNewMesh,i) 
	 error[i] += lagrangeAreaWeights[i] * 
	 				nonLin->HValue(real(uPrevOnNewMesh[i]));
    }

    problem.Ab->Residual(error,u); 
    interfaceDLY->rhsToNB(error);     // residual from hier. to nodal Basis

    FORALL(uQ,i) uQ[i] = 0.0;
    FORALL(u,i)  uQ[i] = u[i];

    interfaceDLY->solToNB(uQ);     // interpolate solution from lin. to
                               	   //  quad nodal Basis

    nonLinCorrection(error, uQ, ADiag, lagrangeAreaWeights, problem);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void RK1:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ, 
			    const Vector<Num>& ADiag, 
			    const Vector<Real>& lagrangeAreaWeight,
			    Problem& problem)
{
    int  i;
    Num  res, a;
    Real uppDefO, lowDefO;		// dummies
    Bool critical;
    const int dim = error.high();
    
    Vector<Real> extUppO(dim), extLowO(dim);

    const NonLinearity* nonLin = problem.getNonLinearity();
    nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);


    // --	    one non-linear Jacobi-step:

    for (i=1; i<=dim; ++i)
    {
	res = error[i];
	nonLin->defectCorrection(error[i], res, uQ[i], ADiag[i], a, uppDefO, 
				 lowDefO, critical, extUppO[i], extLowO[i], 
				 lagrangeAreaWeight[i]);
    }
}
//-------------------------------------------------------------------------


void RK1:: enforcedPMediaRefinement(Problem& problem)
{
    int i;
    PATCH* patch;

    MESH* mesh = problem.Mesh();

    doRefine.resize(mesh->noOfElements());
    FORALL(doRefine,i)  doRefine[i] = False;


    NonLinearity* NonLin = problem.getNonLinearity();

    if (NonLin->castToPorousMedia())
    {
	// -- force elements at the 'frontal boundary' to be refined:
	
	Vector<Num>& u = problem.u;

	Vector<int> nodes(problem.element->NoOfNodes());
	
	Num uMin = u[1];
	Num uMax = u[1];
	FORALL(u,i) 
	{
	    if (u[i] < uMin) uMin = u[i];
	    if (u[i] > uMax) uMax = u[i];
	}
	
	Num upperBound = 1e-4;  // upperBound=2.525e-5
	upperBound *= (uMax - uMin);
	
	int no = 0;
	mesh->resetElemIter();
	while ((patch = mesh->elemIterAll()))
	{
	    ++no;
	    problem.interface->getGlobalNodes(patch, nodes);
	    
	    uMin = 2.0*upperBound;
	    Bool allNodesCritical = True;
	    
	    FORALL(nodes,i)
	    if (nodes[i]) 
	    { 
		if (!NonLin->isCritical(nodes[i])) allNodesCritical = False;
		if (u[nodes[i]] < uMin) uMin = u[nodes[i]]; 
	    }
	    
	    if (uMin < upperBound) doRefine[no] = True;
	    
	    if (allNodesCritical)  doRefine[no] = False;  
	}
	
	
	// -- 	extend forced refinement pattern:
	
	int noOfNeighbPatches=0;
	
	no = 0;
	mesh->resetElemIter();
	while ((patch = mesh->elemIterAll()))
	{
	    ++no;
	    if (doRefine[no]) patch->setMark();
	    else	      patch->unMark();
	    
	    noOfNeighbPatches = patch->noOfNeighbours();
	}
	
	Vector<PATCH*> neighb(noOfNeighbPatches);
	
	no = 0;
	mesh->resetElemIter();
	while ((patch = mesh->elemIterAll()))
	{
	    ++no;
	    if (!doRefine[no])
	    {
		patch->getNeighbours(neighb);
		FORALL(neighb,i) 
		{
		    if (neighb[i])
		      if (neighb[i]->marked()) 
		      { 
			  doRefine[no] = True; 
			  break;
		      }
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


RKA:: RKA(Element* elementDLY0, Element* elementLag0, Interface* interface0)  

	: RK1(elementDLY0, elementLag0, interface0)  
{ }
//-------------------------------------------------------------------------


Bool RKA:: estimateError(Problem& problem)
{
    if (infoErrorEstimator) cout << "\n\tError Est.: RKA";

    distributeError(problem);
    enforcedPMediaRefinement(problem);

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


void RKA:: distributeError(Problem& problem)
{
    int i;
    Vector<Num> error(1), uQ(1);

    MLSparseMatrix AQ(problem.Ab->symmetry, problem.SpaceDim());

    solveQuadProblem(problem, error, uQ, AQ);

    AQ.Mult(uQ,error); 
    EnergyError = 0.5 * Abs(cdot(uQ,error));

    interfaceDLY->solToHB(error);

    for (i=problem.u.high()+1; i<=error.high(); ++i) 
    {
	error[i] = error[i] * error[i] * AQ.Diag(i);	// error on edge
    }

    distributeEdgeErrors(problem, error);

    //cout << "\n-- EnergyErr " << EnergyError << "  EdgeErr " << EdgeError 
    //     << "  d=" << 100*Abs(EdgeError-EnergyError)/EnergyError << " %\n";
}
//-------------------------------------------------------------------------


void RKA:: solveQuadProblem(Problem& problem, Vector<Num>& error, 
			    Vector<Num>& uQ,  MLSparseMatrix& AQ)
{
    int i;
    PATCH* patch;

    MESH* mesh = problem.Mesh();
    const Vector<Num>&  u = problem.u;

    const int edNodeM1 = u.high();
    int   edNode   = edNodeM1;

    interfaceDLY->setHighOrderNodes(&edNode);

    uQ.resize(edNode);
    error.resize(edNode);
    FORALL(error,i) error[i] = 0.0;

    Vector<Num> uPrevOnNewMesh(1);
    TransientProblem* trP = 0;
    if ((trP = problem.castToTransientProblem()))
    {
	uPrevOnNewMesh.resize(1,edNode);
	trP->solutionToNewMesh(uPrevOnNewMesh, interfaceDLY);
    }


    // -- 	full assembling procedure for new lin. system:


    AQ.extend(edNode);

    LinSystem AbQ(problem.spaceDim, AQ.symmetry, problem.precond, 
		  problem.dirichletBCsDLY);

    AbQ.update(&AQ);
    AbQ.reset();


    const int high = lagrangeElementDLY->NoOfNodes();

    Matrix<Num> AElem(high,high);
    Vector<Num> bElem(high);

    Vector<int> globalNodes(high);
    Jacobian Jac;

    Vector<Real> LNorm(high);
    Vector<Real> lagrangeAreaWeights(error.high());
    FORALL(lagrangeAreaWeights,i) lagrangeAreaWeights[i] = 0.0;


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

	problem.assemble(*lagrangeElementDLY,*patch,Jac,AElem,bElem,0,True);
						// 0: use default symm. pattern

	interfaceDLY->getGlobalNodes(patch, globalNodes);  
	AbQ.store(AElem, bElem, globalNodes);


	FORALL(LNorm,i) LNorm[i] = 0.0;
	lagrangeElementDLY->assembleLNorm(*patch,LNorm,Jac,0);

	for (i=1; i<=high; ++i) 
	  if (globalNodes[i]) lagrangeAreaWeights[globalNodes[i]] += LNorm[i];
    }
    
    Real time = 0.0;

    if (trP)	     // transient problem: contrib. of previous step to rhs
    {
	const NonLinearity* nonLin = problem.getNonLinearity();

	FORALL(uPrevOnNewMesh,i) 
	  AbQ.b[i] += lagrangeAreaWeights[i] * nonLin->HValue(uPrevOnNewMesh[i]);

	time = trP->time;
    }


    interfaceDLY->updateDirichletBCs(time);
    AbQ.setDirichletBCs(*problem.dirichletBCsDLY);

    // -- end of assembling proc.


    FORALL(uQ,i) uQ[i] = 0.0;
    FORALL(u,i)  uQ[i] = u[i];

    interfaceDLY->solToNB(uQ);     // interpolate solution from lin. to
                                   //  quad nodal Basis
    AbQ.Residual(error,uQ);

    nonLinCorrection(error, uQ, AQ, lagrangeAreaWeights, problem);
}
//-------------------------------------------------------------------------


void RKA:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ, 
			    MLSparseMatrix& AQ,
			    const Vector<Real>& lagrangeAreaWeight,
			    Problem& problem)
{
    int  i;
    Num  eNew, res, a;
    Real uppDefO, lowDefO;	
    Bool critical;
    const int quadDim = error.high();
    
    Vector<Real> extLowO(quadDim), extUppO(quadDim);

    const NonLinearity* nonLin = problem.getNonLinearity();
    nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);


    if (Cmd.isSet("RKSmooth","GS"))	//  non-linear Gauss-Seidel steps:
    {
	AQ.symmetricExpansion();
	
	Vector<Num> r(quadDim);
	FORALL(r,i) { r[i] = error[i];  error[i] = 0.0; }
	
	for (i=1; i<=quadDim; ++i)
	{
	    if (problem.dirichletBCsDLY->isSet(i))  continue;
	    
	    res = AQ.MultRow(i,error);
	    res = r[i] - res;
	    
	    nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	}
	
	for (i=quadDim; i>=1; --i)
	{
	    if (problem.dirichletBCsDLY->isSet(i)) continue;
	    
	    res = AQ.MultRow(i,error);
	    res = r[i] - res;
	    
	    nonLin->defectCorrection(eNew, res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	    error[i] += eNew;
	}
    }
    else				     // non-linear Jacobi-step:
    {
	for (i=1; i<=quadDim; ++i)
	{
	    if (problem.dirichletBCsDLY->isSet(i))  continue;
	    
	    res = error[i];
	    nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	} 
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


RK2:: RK2(Element* elementDLY0, Element* elementLag0, Interface* interface0)  

	: RKA(elementDLY0, elementLag0, interface0)  
{ }
//-------------------------------------------------------------------------


Bool RK2:: estimateError(Problem& problem)
{
    if (infoErrorEstimator) cout << "\n\tError Est.: RK2";

    if (problem.precond->mode() != multiLevel)
    {
	cout << "\n*** RKC:: Preconditioner must be of multi-level-type:\n"
	     << "NonLinearMLSG or TrcNonLinearMLSG\n";
	cout.flush(); abort();
    }

    distributeError(problem);
    enforcedPMediaRefinement(problem);

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


void RK2:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ, 
			    MLSparseMatrix& AQ,
			    const Vector<Real>& lagrangeAreaWeight,
			    Problem& problem)
{
    int  i;
    Num  eNew, res, a;
    Real uppDefO, lowDefO;	
    Bool critical;
    const int quadDim = error.high();
    
    Vector<Real> extLowO(quadDim), extUppO(quadDim);

    const NonLinearity* nonLin = problem.getNonLinearity();
    nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);


    const int linDim = nonLin->lowObstacle.high();

    Vector<Num> r(quadDim), Diag0(quadDim), aux(quadDim);

    FORALL(r,i) 
    { 
	r[i] = error[i];  
	error[i] = 0.0; 
	Diag0[i] = AQ.Diag(i);
    }

    NonLinearMLGS* precond = (NonLinearMLGS*) problem.precond;

    Vector<Real>& uppDefObstacle = *precond->uppDefObstacle.Top();
    Vector<Real>& lowDefObstacle = *precond->lowDefObstacle.Top();

    FORALL(uppDefObstacle,i)  
    {
	uppDefObstacle[i] = lowDefObstacle[i] = 0.0;
	precond->setUnCritical(i);
    }


    if (Cmd.isSet("RKSmooth","GS"))	//  non-linear Gauss-Seidel steps:
    {
	AQ.symmetricExpansion();
	
	for (i=1; i<=quadDim; ++i)
	{
	    if (problem.dirichletBCsDLY->isSet(i)) continue;
	    
	    res = AQ.MultRow(i,error);
	    res = r[i] - res;
	    nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	}
	
	for (i=quadDim; i>=1; --i)
	{
	    if (problem.dirichletBCsDLY->isSet(i)) continue;
	    
	    res = AQ.MultRow(i,error);
	    res = r[i] - res;
	    nonLin->defectCorrection(eNew, res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	    r[i] -= res;
	    error[i]   += eNew;
	    AQ.Diag(i) += a;
	    
	    if (i <= linDim) 
	    {
		uppDefObstacle[i] = uppDefO - real(error[i]);
		lowDefObstacle[i] = lowDefO - real(error[i]);
	    
		if (critical) precond->setCritical(i);
		else	  precond->setUnCritical(i);
	    } 
	}
    }
    else				     // non-linear Jacobi-step:
    {
	for (i=1; i<=quadDim; ++i)
	{
	    if (problem.dirichletBCsDLY->isSet(i)) continue;
	    
	    res = r[i];
	    nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a, 
				     uppDefO, lowDefO, critical, extUppO[i], 
				     extLowO[i], lagrangeAreaWeight[i]);
	    AQ.Diag(i) += a;
	    r[i] -= res;
	    
	    if (i <= linDim) 
	    {
		uppDefObstacle[i] = uppDefO - real(error[i]);
		lowDefObstacle[i] = lowDefO - real(error[i]);
	    
		if (critical) precond->setCritical(i);
		else	      precond->setUnCritical(i);
	    } 
	} 
    }


    // --  restrict Matrix and residual to hierarchical linear basis:


    Generation gen(quadDim);
    interfaceDLY->setHBGeneration(gen);		  // for Galerkin-Restriction

    MLMatrix* ALin = precond->Al.Top()->castToMLMatrix();
    ALin->reset();
    ALin->markNodes();
    ALin->GalerkinRestriction(AQ, gen);

    AQ.Mult(aux,error);
    FORALL(r,i) aux[i] = r[i] - aux[i];		// the new residual
    interfaceDLY->rhsToHB(aux);

    const int maxLevel = precond->maxLevel; 
    Vector<Num> rLin(linDim);
    FORALL(rLin,i) 
    {
	if (precond->isCritical(i,maxLevel))     rLin[i] = 0.0;
	else if (problem.dirichletBCs->isSet(i)) rLin[i] = 0.0;
	else  rLin[i] = aux[i];
    }


    // initialize the multi-level preconditioner and call one V-Cycle:


    const int nPreSmooth0  = precond->nPreSmooth;
    const int nPostSmooth0 = precond->nPostSmooth;
    precond->nPreSmooth  = 1;
    precond->nPostSmooth = 1;

    Vector<Num> b(1); 					// a dummy
    precond->initialize(ALin, problem.u, b);

    Vector<Num>& eMG = precond->eMG;

    precond->el.push(&eMG);
    precond->MGCycle(eMG, *ALin, rLin, precond->maxLevel);
    precond->el.pop();

    FORALL(aux,i) aux[i] = 0.0;
    FORALL(eMG,i) 
    {
	if (!precond->isCritical(i, maxLevel) &&
	    !problem.dirichletBCs->isSet(i))  aux[i] = eMG[i];
    }
    interfaceDLY->solToNB(aux);

    FORALL(error,i) error[i] += aux[i];


    // -- 	    restore preconditioner settings:

    precond->nPreSmooth  = nPreSmooth0;
    precond->nPostSmooth = nPostSmooth0;

    FORALL(precond->Al,i) precond->Al[i]->removeSymmetricExpansion();

    FORALL(Diag0,i)  AQ.Diag(i) = Diag0[i];		// restore diagonal
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------