#include "PollutionError.h"
#include "mpi.h"

namespace AMDiS
{

  PollutionErrorKonvex::BoundaryTraverseStack::BoundaryTraverseStack(Mesh *mesh, 
								     PartitionStatus ps_)
  {
    TraverseStack myStack;
    ElInfo *swap = 
      myStack.traverseFirst(mesh,-1, 
			    Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_NEIGH);
    ps = ps_;

    WorldVector<double>* centerWorld;
    DimVec<double> center = DimVec<double>(2, DEFAULT_VALUE, 1.0 / 3.0);
    while (swap != NULL) {
      if (proofConditions(swap)) {
	// if conditions are true, push the midpoint to the end of the list. a better 
	// algorithm should use the boundary edge
	centerWorld = new WorldVector<double>();
	swap->coordToWorld(center, *centerWorld);
	push_back(centerWorld);
      }
      swap = myStack.traverseNext(swap);
    }
  }
	
  bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo) 
  {
    if (elInfo == NULL) 
      return false;

    Element *curEl = elInfo->getElement();
    PartitionElementData *ped = 
      static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));

    TEST_EXIT_DBG(ped != NULL)("No PartitionElementData found");    

    if (ped != NULL && ped->getPartitionStatus() == OVERLAP) {
      //check neighbour
      int nrNeigh = curEl->getGeo(NEIGH);
      for (int i = 0; i < nrNeigh; i++) {
	Element* curNeigh = elInfo->getNeighbour(i);
	if (curNeigh != NULL) { 		    
	  ped = 
	    static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
	  if (ped != NULL && ped->getPartitionStatus() == ps)
	    return true;
	}
      }
    }
    return false;
  }

  PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r)
    : ResidualEstimator(name,r),
      center(2,DEFAULT_VALUE, 1.0 / 3.0)
  {
    FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
    GET_PARAMETER(0, name + "->C1p", "%f", &C1p);
    GET_PARAMETER(0, name + "->C2p", "%f", &C2p);   
  }

  void PollutionErrorKonvex::init(double timestep)  
  {
    FUNCNAME("PollutionErrorKonvex::init");
    ResidualEstimator::init(timestep);

    // set the distance (ie. width of the overlapped area)
    // approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
    // retrieve all elements at the boundary overlapped -> out
    
    bts = new BoundaryTraverseStack(mesh, IN);	
  }

  void PollutionErrorKonvex::estimateElement(ElInfo *elInfo)  
  {
    FUNCNAME("PollutionErrorKonvex::estimateElement()");

    // switch between inner and outer
    // inner part --> use H_1Norm
    // outer part --> use L_2Norm

    PartitionElementData* ped = 
      static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
    TEST_EXIT_DBG(ped != NULL)("No partition Elementdata found");

    PartitionStatus ps = ped->getPartitionStatus();
    TEST_EXIT_DBG(ps != UNDEFINED)("Elementstatus undefined");

    // save old values
    double C0 = ResidualEstimator::C0;
    double C1 = ResidualEstimator::C1;

    if (ps == OUT) {
      //outer part
      ResidualEstimator::norm = L2_NORM;		
      ResidualEstimator::C0 *= C2p;

      //compute distance
      double d = getDistance(elInfo);
      ResidualEstimator::C1 *= C2p / d;
      
    } else {
      //inner part
      ResidualEstimator::norm = H1_NORM;
      ResidualEstimator::C0 *= C1p;
      ResidualEstimator::C1 *= C1p;
    }

    ResidualEstimator::estimateElement(elInfo);
    //set old values
    ResidualEstimator::C0 = C0;
    ResidualEstimator::C1 = C1;
  }

  double PollutionErrorKonvex::getDistance(ElInfo* element)  
  {
    if (bts==NULL)
      MSG("bts is null");
    TEST_EXIT_DBG(bts!=NULL)("no bts");
    double distance = 0.0;
    double curDistance = 1.0;
    //get center in world coord
    WorldVector<double> centerWorld;
    element->coordToWorld(center,centerWorld);

    //start traversal over all Midpoints at the boundary overlap <-> in      
    for (BoundaryTraverseStack::iterator myIt = bts->begin(); 
	 myIt != bts->end(); myIt++) {
      curDistance = 0;	
      for (int i = 0;i < dim;i++)
	curDistance += ((*(*myIt))[i]-centerWorld[i]) * ((*(*myIt))[i] - centerWorld[i]);

      distance = max(curDistance, distance);
    }
    return distance;
  }
  
  void PollutionErrorKonvex::exit(bool output = true) 
  {
    ResidualEstimator::exit(output);
    delete bts;
  }
}