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

#include "adapt.h"

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

#include "problem.h"
#include "problemtr.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;

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


ErrorEstimator:: ErrorEstimator() 
{

    elementError = 0;

    minRefinementRatio = 0.05;  Cmd.get("minRefRatio", &minRefinementRatio);

    infoErrorEstimator = 0;  Cmd.get("infoErrorEstimator", &infoErrorEstimator);
    timeErrorEstimator = 0;  Cmd.get("timeErrorEstimator", &timeErrorEstimator);
    accTime = 0;  	     Cmd.get("accTimeErrorEstimator", &accTime);


    refStrategy = uniform;
    if      (Cmd.isSet("refStrategy","maxValue")) refStrategy=maxValue;
    else if (Cmd.isSet("refStrategy","extraPol")) refStrategy=extrapolation;
    else if (Cmd.isSet("refStrategy","uniform"))  refStrategy=uniform;
    else if (Cmd.isSet("refStrategy","random"))   refStrategy=random;
    else  MissingParameter("refStrategy");
}
//-------------------------------------------------------------------------

ErrorEstimator:: ~ErrorEstimator() { }

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


Bool ErrorEstimator:: estimateError(Problem& /*problem*/) 
{
    if (infoErrorEstimator) cout << "\n    No Error Estimator";

    EnergyError = machMax(Real(0.0));
    return False;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


DLY:: DLY(Element* elementDLY0, Interface* interfaceDLY0)  

    : ErrorEstimator(), elementDLY(elementDLY0), interfaceDLY(interfaceDLY0)
{ }
//-------------------------------------------------------------------------

DLY:: ~DLY()
{
    delete elementDLY;
    delete interfaceDLY;
}
//-------------------------------------------------------------------------


EFDLY:: EFDLY(Element* elementDLY0, Interface* interface0)  

    : DLY(elementDLY0, interface0)  
{ }
//-------------------------------------------------------------------------

ResidualError:: ResidualError() : ErrorEstimator()  { }

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

QuadStdTriangleError:: QuadStdTriangleError() : ErrorEstimator() { }

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


void ErrorEstimator:: adapt(Problem& problem, Real* energyError, 
			    Real reqAbsGlobalPrec, Real* linENorm)
{
    int i, markedElements = 0;
    static Real accTimeError = 0.0;
    Timer timer, accTimer;

    elementError = new Vector<Real>(problem.Mesh()->noOfElements());
    FORALL(*elementError,i) (*elementError)[i] = 0.0;


    if (estimateError(problem) == True)
      setRefinementFlags(*problem.Mesh(), reqAbsGlobalPrec, &markedElements);
	
    if (markedElements == 0) 
      UniformRefinement(*problem.Mesh(), reqAbsGlobalPrec, &markedElements);


    delete elementError;
    elementError = 0;

    *energyError = EnergyError;


    if (linENorm != 0) *linENorm = LinENorm;

    //if (Cmd.isTrue("L2Norm")) cout << "\tL2-Error: " << L2Error << "\n";

    if (timeErrorEstimator)  
	{ cout << "\n\tError Estimation: "; timer.cpu(); cout.flush(); }
    if (accTime)
    {
	accTimeError += accTimer.cpu(False);
	cout << Form("\tAccumulated time:  %1.2f sec.\n", accTimeError); 
    }
}
//-------------------------------------------------------------------------


void ErrorEstimator:: setRefinementFlags(MESH& mesh, Real reqGlobalPrec, 
					 int* markedElems)
{
    switch (refStrategy)
    {
      case maxValue:      MaxValue(mesh, reqGlobalPrec, markedElems); 
			  break;
      case random: 	  RandomRefinement(mesh, reqGlobalPrec, markedElems); 
			  break;
      case extrapolation: ExtrapolRefinement(mesh, reqGlobalPrec, markedElems);
			  break;
      case uniform:  	  UniformRefinement(mesh, reqGlobalPrec, markedElems);
			  break;
    }
}
//-------------------------------------------------------------------------


void ErrorEstimator:: MaxValue(MESH& mesh, Real reqAbsGlobalPrec, int* marked)
{
    int i, loops=0;
    int markedElements;
    PATCH* patch;

    Real maxError = 0.0, limit;
    
    FORALL (*elementError, i)  
        if ((*elementError)[i] > maxError)  maxError = (*elementError)[i];


    Real limitFactor = 0.25;  	Cmd.get("limitFactor", &limitFactor);
    limit = limitFactor*maxError;
    if (!(EnergyError == 0))
    limit = Max(limit, reqAbsGlobalPrec/(EnergyError)*sqrt(limit*maxError));

    do
    {
	markedElements = 0;
	mesh.resetElemIter();
	i = 0;

	while ((patch = mesh.elemIterAll()))
	{
	    ++i;
	    if ( (*elementError)[i] > limit || DoRefine(i) )  
	    {
                patch->setMark();
		++markedElements;
	    }
	    else patch->unMark();
	}    
	limit *= 0.9;
	if (loops < 100) ++loops;
	else 	         break;
    }
    while (markedElements < minRefinementRatio * mesh.noOfElements());

    *marked = markedElements;

    if (infoErrorEstimator) 
    {
	cout << "\t -- Refinement: MaxValue\n";
	cout << "\tMarked Elements: " << markedElements 
	<< " (" << int(100*markedElements/mesh.noOfElements()) <<  " %), "
	<< " Limit reduction loops: " << loops <<"\n";
    }
}
//-------------------------------------------------------------------------


void ErrorEstimator:: ExtrapolRefinement(MESH& mesh, Real reqAbsGlobalPrec,
					 int* marked)
{
    int i, loops=0;
    int markedElements;
    PATCH* elem, *elemFather;

    Real err, extrapolErr, fatherErr, limit = 0.0, maxError = 0.0;

    mesh.resetElemIter();
    i = 0;
    while ((elem = mesh.elemIterAll()))
    {
	err = (*elementError)[++i];
	elem->setError(err);
	elemFather = elem->Father();

	if (elemFather != 0)
	{
            fatherErr = elemFather->getError();
            if (equal(fatherErr,0.0)) extrapolErr = 0.0;
            else 		      extrapolErr = err*err/fatherErr;
	}
	else extrapolErr = 0.0;   // err
	
	if (extrapolErr > err)   extrapolErr = err;   // rather unlikely
	if (extrapolErr > limit) limit = extrapolErr; 

        if (err > maxError)  maxError = err;
    }

    if (mesh.maxDepth==0) limit=0.5*maxError; 

    Real limitFactor = 0.5;  	Cmd.get("limitFactor", &limitFactor);
    limit *= limitFactor;		// 0.5 is empirical

    if (!(EnergyError==0.0))
	limit = Max(limit, reqAbsGlobalPrec/(EnergyError)*sqrt(limit*maxError));

    do
    {
	markedElements = 0;
	mesh.resetElemIter();
	i = 0;
	while ((elem = mesh.elemIterAll()))
	{
	    ++i;
	    if ((*elementError)[i] > limit  || DoRefine(i))  
	    {
		elem->setMark();
		++markedElements;
	    }
	    else elem->unMark();
	}    
	limit *= 0.9;
	if (loops < 100) ++loops;
	else 	         break;
    }
    while (markedElements < minRefinementRatio * mesh.noOfElements());

    *marked = markedElements;

    if (infoErrorEstimator) 
    {
	cout << "\t -- Refinement: Extrapolation\n";
	cout << "\tMarked Elements: " << markedElements 
	<< " (" << int(100*markedElements/mesh.noOfElements()) <<  " %), "
	<< " Limit reduction loops: " << loops <<"\n";
    }
}
//-------------------------------------------------------------------------


void ErrorEstimator:: UniformRefinement(MESH& mesh, Real /*reqAbsGlobalPrec*/,
					int* marked)
{
    PATCH* patch;

    mesh.resetElemIter();
    while ((patch = mesh.elemIterAll()))  patch->setMark(); 

    *marked = mesh.noOfElements();

    if (infoErrorEstimator) 
    {
	cout << "\t -- Refinement: uniform\n";
	cout << "\tMarked Elements: " << mesh.noOfElements() << "\n";
    }
}
//-------------------------------------------------------------------------


void ErrorEstimator:: RandomRefinement(MESH& mesh, Real /*reqAbsGlobalPrec*/, 
				       int* marked)
{
    int i;
    PATCH* patch;

    mesh.resetElemIter();
    i = 0;

    while ((patch = mesh.elemIterAll()))
    {
	++i;
	if (i%10 == 1)  patch->setMark();
	else 		patch->unMark();
    }

    *marked = i;

    if (infoErrorEstimator) 
    {
	cout << "\t -- Refinement: random\n";
	cout << "\tMarked Elements: "<< float(i)/10<<"\n";
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool ResidualError:: estimateError(Problem& /*problem*/)
{
    
    cout << "\n\tError Est.: ResidualError not implemented\n" ; 
    cout.flush(); abort();
    /*
      int i, depth, noTri, noTri1, noTri2, ip;
      TR*  tr;
      
      if (infoErrorEstimator) cout << "\n\tError Est.: Residual Error";
      
      MESH& mesh    = *problem.Mesh();
      Element& 	   element = *problem.element;
      Interface& interface = *problem.interface;
      Vector<Num>&   u       =  problem.u;
      
      
      noTri = 0;
      mesh.setAllInnerNodes(&noTri);  		// Number triangles:
      
      int spaceDim = element.spaceDim;
      
      EDG* ed;
      Vector<Num> gradElem(spaceDim);			// local gradient
      Vector<Num> uElem(element.noOfNodes);		// local solution
      Vector(int) node (element.noOfNodes);
      
      Matrix<Num>   grad  (spaceDim, mesh.noOfTriangles);
      Vector<Real> mat(mesh.noOfTriangles);
      
      
      //  calculate gradients:
      
      Jacobian Jac;
      ip = 1;
      FORALL (mesh.trList, depth)
      for (TrListProc tp1(mesh.trList[depth],all,&tr); tr; tr=tp1.All())
      {
      interface.getGlobalNodes(*tr, node);
      FORALL(uElem,i) uElem[i] = u[node[i]];
      
      Jac.compJinv(*tr);
      element.gradient(*tr, Jac, uElem, gradElem, ip);
      
      mat[tr->node] = problem.material->E(tr->classA);
      for (i=1; i<=spaceDim; ++i)  
      grad(i,tr->node) = mat[tr->node] * gradElem[i];
      // cout << "\n--- gradient: " << gradElem[1] << "  " << gradElem[2];
      }
      
      
      Num   flux1, flux2, js;	
      Real normal[4], length;
      
      FORALL (mesh.edgList, depth)
      for (EdgListProc ep(mesh.edgList[depth],all,&ed); ed; ed=ep.All())
      {
      if (ed->boundP == DIRICHLET) continue;
      
      ed->compNormalVector(normal, &length);
      
      noTri2 = 0;
      if (ed->t1 == 0)  noTri1 = ed->t2->node;
      else  
      {
      noTri1 = ed->t1->node;
      if (ed->t2) noTri2 = ed->t2->node;
      }
      
      flux1 = flux2 = 0.0;
      
      for (i=1; i<=spaceDim; ++i) flux1 += normal[i] * grad(i,noTri1);  
      
      if (noTri2) 
      for (i=1; i<=spaceDim; ++i)  flux2 += normal[i] * grad(i,noTri2);
      
      js = flux1-flux2;		
      js = js*js * length*length / 24.;	
      
      if (noTri1 && noTri2) js *= 0.5;
      
      (*elementError)[noTri1] += abs(js / mat[noTri1]); 
      if (noTri2) 
      (*elementError)[noTri2] += abs(js / mat[noTri2]);
      }
      
      EnergyError = 0.0;
      FORALL (*elementError, i) EnergyError += (*elementError)[i];
      */
    return False;
}
//-------------------------------------------------------------------------


		// returns RELATIVE error in per cent !


Bool QuadStdTriangleError:: estimateError(Problem& /*problem*/)
{

    cout << "\n\tError Est.: QuadStdTriangleError not implemented\n" ; 
    cout.flush(); abort();

    /*
      int i, j, k, depth;
      TR* t;
      Num TotalLinEnergy=0.0, LinEnergy, TotalQuadEnergy=0.0, QuadEnergy;
      
      
      MESH& mesh    = *problem.Mesh();
      Element& 	   element = *problem.element;
      Interface& interface = *problem.interface;
      Vector<Num>&   u 	   =  problem.u;
      
      if (infoErrorEstimator) cout << "\n\tError Est.: QuadStdTriangleError";
      
      
      int noTri = 0;
      mesh.setAllInnerNodes(&noTri);	  // Number triangles
      
      
      // assemble elements and calculate energies:
      
      
      int dimElem = element.noOfNodes;
      Matrix<Num> AElem(dimElem, dimElem);
      Vector<Num> Au(dimElem), uElem(dimElem);
      Vector(int) globalNodes(element.noOfNodes);
      Jacobian Jac;
      
      
      FORALL (mesh.trList, depth)
      for (TrListProc trp2(mesh.trList[depth],all,&t); t; t=trp2.All())
      {
      
      Jac.compJinv(*t);
      
      element.initAb(AElem);
      interface.getGlobalNodes(*t, globalNodes);  
      
      problem.assembleErrorEstimator(*t, Jac, AElem);
      
      if (problem.symmetry == sym)			// set upper triangle
      {
      FORALL_ROWS(AElem,i)
      for (j=1; j<i; ++j)  AElem(j,i) = AElem(i,j);
      }
      
      FORALL(Au,i) Au[i] = 0.0;
      FORALL(uElem,i) {
      if (globalNodes[i]) uElem[i] = u[globalNodes[i]];
      else		uElem[i] = 0.0;
      }
      
      for (i=1; i<=3; ++i)		    	// compute linear Au = A*u
      {
      if (globalNodes[i]==0) continue;
      for (j=1; j<=3; ++j)  Au[i] += AElem(i,j) * uElem[j]; 
      }
      
      LinEnergy = 0.0;
      for (i=1; i<=3; ++i) LinEnergy += conj(uElem[i]) * Au[i];
      LinEnergy *= 0.5;
      
      
      for (i=1; i<=3; ++i)	     		// compute quadratic Au = A*u
      {
      if (globalNodes[i]==0) continue;
      for (j=4; j<=6; ++j)  Au[i] += AElem(i,j) * uElem[j]; 
      }
      for (i=4; i<=6; ++i)
      {
      if (globalNodes[i]==0) continue;
      for (j=1; j<=6; ++j)  Au[i] += AElem(i,j) * uElem[j]; 
      }	
      
      QuadEnergy = 0.0;
      for (i=1; i<=6; ++i) QuadEnergy += Au[i] * uElem[i];
      QuadEnergy *= 0.5;
      
      TotalLinEnergy  += LinEnergy;
      TotalQuadEnergy += QuadEnergy;
      
      (*elementError)[t->node] = 0.25 * abs(real(LinEnergy-QuadEnergy));
      }
      
      EnergyError = 0.25 * (TotalLinEnergy-TotalQuadEnergy);
      
      EnergyError = abs(100 * EnergyError / TotalQuadEnergy);
      LinENorm = abs(TotalLinEnergy);
      
      
      if (infoErrorEstimator)
      {
      cout << "\n\nTotalQuadEnergy = " << TotalQuadEnergy 
      <<   "  TotalLinEnergy = " << TotalLinEnergy
      << "  Error = " << EnergyError << " %\n";
      }
      
      // cout << "\n elementError: " << *elementError <<"\n";
  */
    return False;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool DLY:: estimateError(Problem& problem)
{
    if (infoErrorEstimator)
    {
	if (Cmd.isSet("ErrorEstim","ModDLY")) cout << "\n\tError Est.: ModDLY";
	else 				    cout << "\n\tError Est.: DLY";
    }
    Vector<Num> error(1);


    Bool status = solveQuadDefectProblem(problem, error);

    if (status == True)	distributeEdgeErrors(problem, error);	

    return status;
}   
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool DLY:: solveQuadDefectProblem(Problem& problem, Vector<Num>& b)
{
    int i, j;
    Num sum;
    PATCH* patch;

    const Vector<Num>&  u = problem.u;

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

    interfaceDLY->setHighOrderNodes(&edNode);

    if (edNode == edNodeM1) 			// only dirichlet-edges !
    {
	EnergyError = 0.0;
	return False;
    }

    b.resize(edNodeM1+1,edNode);
    Vector<Num> ADiag(edNodeM1+1,edNode);

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


    // -- 	    local assembling for DLY:

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

    if ((trP = problem.castToTransientProblem()))
    {
	if (trP->transientMode)
	{
	    uPrevOnNewMesh.resize(b.low(),b.high());
	    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;


    problem.Mesh()->resetElemIter();

    while ((patch = problem.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 and right hand side
	{
	    if (globalNodes[i] == 0) continue;

	    ADiag[globalNodes[i]] += AElem(i,i);
	    b    [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]];
	    }
	    b[globalNodes[i]] -= sum;
	}


	if (trP)	     // transient problem (implicit Euler assumed)
	{
	  if (trP->transientMode)
	  {
	    for (i=low; i<=high; ++i)  AElem(i,i) = 0.0;

	    elementDLY->assembleLumpedP(*patch, AElem, Jac,  
					   elementDLY->DLYPattern);
	    for (i=low; i<=high; ++i)  
	      if (globalNodes[i])
	        b[globalNodes[i]] += AElem(i,i) * uPrevOnNewMesh[globalNodes[i]];
	  }
        }
    }

    FORALL(b,i) b[i] = conj(b[i])* b[i] / ADiag[i];	// error on edge

    EnergyError = 0.0;
    FORALL(b,i) EnergyError += Abs(b[i]);		// sum up global Error
    EnergyError *= 0.5;

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


void DLY:: distributeEdgeErrors(Problem& problem, Vector<Num>& error)
{
    int i, n, node, no=0;
    PATCH* patch;

    MESH* mesh = problem.Mesh();
    mesh->resetElemIter();

    const int nComp = problem.NComp();


    if (problem.SpaceDim() == 1)
    {
	while ((patch = mesh->elemIterAll()))
	{
            ++no;
	    if ((node=patch->getNode()))
	      for (n=1; n<=nComp; ++n) (*elementError)[no] += Abs(error[node++]);
	}
    }
    else
    {
	while ((patch = mesh->elemIterAll()))
	{
	    ++no;
	    const Vector<EDG*>& edges = patch->getEdges();

	    FORALL(edges,i)
	    {
		if (edges[i])
		{
		    if ((node = edges[i]->getNode()))
		    {
			for (n=1; n<=nComp; ++n) 
			  (*elementError)[no] += Abs(error[node++]);
		    } 
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool EFDLY:: estimateError(Problem& problem)
{
    int i, node;
    PATCH* patch;

    if (infoErrorEstimator) cout << "\n\tError Est.: EF-DLY";

    if (problem.spaceDim != 3) 
    {
	cout << "\n*** Estimator EFDLY: dimension not 3\n\n";  cout.flush(); abort();
    }

    Vector<Num> error(1);

    Bool status = solveQuadDefectProblem(problem, error);

    if (status == True)
    {
	// --		distribute errors to elements
	int noOfFaces = problem.element->NoOfFaces();
	Vector<int> faceNodes(noOfFaces);
	
	
	problem.Mesh()->resetElemIter();
	int no = 0;
	
	while ((patch = problem.Mesh()->elemIterAll()))
	{
	    ++no;
	    const Vector<EDG*>& edges = patch->getEdges();

	    FORALL(edges,i)
	    {
		if (edges[i])
		{
		    node = edges[i]->getNode();
		    if (node) (*elementError)[no] += Abs(error[node]);
						//  support: ~ 6 tetrahedra)
		} 
	    }
	    
	    patch->getFaceNodes(faceNodes);
	    FORALL(faceNodes,i)
	    {
		node = faceNodes[i];
		if (node) (*elementError)[no] += 3.0 * Abs(error[node]); 
					// different support: 2 tetrahedra)
	    } 
	}
    }

    return status;
}
//-------------------------------------------------------------------------