//=========================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//=========================================================================

#ifndef __smoab_model_FaceSets_h
#define __smoab_model_FaceSets_h

#include "CellSets.h"
#include <map>
#include <set>

namespace smoab {
namespace model {

const smoab::EntityHandle VOID_REGION=0;

//----------------------------------------------------------------------------
class FaceCellSet : public CellSet
{
public:
  FaceCellSet(int id, smoab::EntityHandle p,const smoab::Range& c):
    CellSet(p,c),
    ID(id)
    {}

  int faceId() const { return ID; }
  void overrideFaceId(int i) { ID = i; } //USE AT YOUR OWN RISK

private:
  int ID;
};


//----------------------------------------------------------------------------
typedef std::vector<FaceCellSet> FaceCellSets;


//----------------------------------------------------------------------------
//class that store the regions that a faces is adjacent too
struct FacesAdjRegions
{
  FacesAdjRegions(int f, smoab::EntityHandle r0, smoab::EntityHandle r1):
    FaceId(f),
    Region0(r0),
    Region1(r1)
    {
    if (r0 > r1)
      {
      std::swap(this->Region0,this->Region1);
      }
    }

  bool operator<(const FacesAdjRegions& other) const
    {
    return (this->FaceId < other.FaceId);
    }

  smoab::EntityHandle otherId(smoab::EntityHandle other) const
    {
    if(other == Region0)
      {
      return Region1;
      }
    return Region0;
    }

  int FaceId;
  smoab::EntityHandle Region0;
  smoab::EntityHandle Region1;
};
//----------------------------------------------------------------------------
smoab::model::FaceCellSets findFaceSets(smoab::CellSets shells,
                                 smoab::CellSets boundaries,
                                 std::set<smoab::model::FacesAdjRegions>& faceMaps)
{
  typedef smoab::CellSets::iterator iterator;
  typedef smoab::model::FaceCellSets::iterator faceIterator;
  typedef std::set<smoab::model::FacesAdjRegions>::const_iterator FaceAdjIterator;

  //we need to properly label each unique face in shells
  //we do this by intersecting each shell with each other shell
  //to find shell on shell contact, and than we intersect each
  //resulting shell with the boundary conditions
  //the end result of these intersections will be the new modelfaces
  int faceId = 1;
  smoab::model::FaceCellSets shellFaces;

  //first intersect each shell with each other shell
  std::set<smoab::model::FacesAdjRegions> shellFaceContacts;
  for(iterator i=shells.begin();i!= shells.end();++i)
    {
    //copy the cells so we can add a face that represents
    //all the cells of the region that aren't shared with another region
    int numCells = i->cells().size(); //size() on range is slow, so cache it
    for(iterator j = i+1;
        j != shells.end() && numCells > 0;
        ++j)
      {
      //intersect i and j to make a new face
      smoab::Range intersection = smoab::intersect(i->cells(),j->cells());
      if(!intersection.empty())
        {
        //don't want to increment faceId when the intersection is empty
        smoab::model::FaceCellSet face(faceId,i->entity(),intersection);
        shellFaces.push_back(face);
        i->erase(intersection);
        j->erase(intersection);
        numCells -= intersection.size();

        //add this to the face map
        smoab::model::FacesAdjRegions faceInfo(faceId,i->entity(),j->entity());
        shellFaceContacts.insert(faceInfo);

        ++faceId;
        }
      }
    //if all the cells for shell i are used, don't add a new
    //empty face
    if(numCells > 0)
      {
      smoab::model::FaceCellSet face(faceId,i->entity(),i->cells());
      shellFaces.push_back(face);

      //add this to the face map
      smoab::model::FacesAdjRegions faceInfo(faceId,VOID_REGION,i->entity());
      shellFaceContacts.insert(faceInfo);

      ++faceId;
      }
    }

  //now we have all the faces that match shell on shell contact
  //we know process all the new faces to see if they intersect
  //with any boundary sets. A boundary set can span multiple
  //shells so we want to process it as a second loop
  if(boundaries.size()==0)
    {
    //no boundaries we are done.
    faceMaps = shellFaceContacts;
    return shellFaces;
    }

  //store the end before we start adding boundary faces, which
  //we don't need to check agianst other boundaries
  faceId = 0; //reset the faced id

  //store in a new face set, expanding the current one causes incorrect results
  smoab::model::FaceCellSets faces;
  for(faceIterator i=shellFaces.begin();i != shellFaces.end();++i)
    {
    //determine from the shell faces if the new face we are creating
    //is bounded by two regions or just one
    smoab::model::FacesAdjRegions idToSearchFor(i->faceId(),-3,-2);
    FaceAdjIterator adjRegions = shellFaceContacts.find(idToSearchFor);
    smoab::EntityHandle otherRegionId = adjRegions->otherId(i->entity());

    int numCells = i->cells().size(); //size() on range is slow, so cache it
    for(iterator j=boundaries.begin();j != boundaries.end(); ++j)
      {
      smoab::Range intersect = smoab::intersect(i->cells(),j->cells());
      if(!intersect.empty())
          {
          //don't want to increment faceId when the intersection is empty
          smoab::model::FaceCellSet face(faceId,j->entity(),intersect);
          faces.push_back(face);
          i->erase(intersect);
          numCells -= intersect.size();
          smoab::model::FacesAdjRegions faceInfo(faceId,i->entity(),otherRegionId);
          faceMaps.insert(faceInfo);

          ++faceId;
          }
      }
    if(numCells > 0)
      {
      smoab::model::FaceCellSet face(faceId,i->entity(),i->cells());
      faces.push_back(face);
      smoab::model::FacesAdjRegions faceInfo(faceId,i->entity(),otherRegionId);
      faceMaps.insert(faceInfo);

      ++faceId;
      }
    }
  return faces;
}

//----------------------------------------------------------------------------
template<typename T>
std::vector<T> faceIdsPerCell(const smoab::model::FaceCellSets& faces)
{
  typedef smoab::model::FaceCellSets::const_iterator iterator;
  typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
  typedef smoab::Range::const_iterator RangeIterator;

  //find all the cells that are in the faceCellSet, and than map
  //the proper face id to that relative position, here comes lower_bounds!
  std::vector<smoab::EntityHandle> searchableCells;
  smoab::Range faceRange = smoab::getAllCells(faces);
  smoab::RangeToVector(faceRange,searchableCells);

  //faceIds will be the resulting array
  std::vector<T> faceIds(searchableCells.size());

  //construct the start and end iterators for the lower bounds call

  EntityHandleIterator s_begin = searchableCells.begin();
  EntityHandleIterator s_end = searchableCells.end();

  //search the face cell sets
  for(iterator i=faces.begin(); i!=faces.end(); ++i)
    {
    T value = static_cast<T>(i->faceId());
    const smoab::Range& entitiesCells = i->cells();
    for(RangeIterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
      {
      EntityHandleIterator result = std::lower_bound(s_begin,
                                                     s_end,
                                                     *j);
      std::size_t newId = std::distance(s_begin,result);
      faceIds[newId] = value;
      }
    }
  return faceIds;
}

//----------------------------------------------------------------------------
template<typename T>
std::map<smoab::EntityHandle,T> findRegions(
    const std::set<smoab::model::FacesAdjRegions>& regionFaceAdj,
    const smoab::Interface& interface)
  {
  //lets talk about what this is doing.
  //When we are looking to determine regions we need to walk the geometric
  //parent child relationship. We are basically trying to find
  //for each region the 3d geometric parent. that will become the region id
  //for that cell

  smoab::GeomTag threeDTag(3);
  smoab::EntityHandle root = interface.getRoot();
  smoab::Range realParents = interface.findEntitiesWithTag(threeDTag,root);

  std::set<smoab::EntityHandle> childIds;
  typedef typename std::set<smoab::model::FacesAdjRegions>::const_iterator iterator;
  for(iterator i=regionFaceAdj.begin(); i!=regionFaceAdj.end();++i)
    {
    childIds.insert(i->Region0);
    childIds.insert(i->Region1);
    }

  T newParentId = T();

  typedef typename std::map<smoab::EntityHandle,T>::iterator parentIterator;
  std::map<smoab::EntityHandle,T> foundParents;
  std::map<smoab::EntityHandle,T> returnValues;
  typedef typename std::set<smoab::EntityHandle>::const_iterator set_iterator;
  for(set_iterator i=childIds.begin(); i != childIds.end(); ++i)
    {
    typedef typename smoab::Range::const_iterator range_iterator;

    /*
     *  IF A REGION IS SET TO VOID_REGION WE NEED TO PUSH THAT VALUE DOWN
     *  AS THE MATERIAL, SINCE THE MOAB DEFAULT TAG VALUE WILL
     *  BE CONSIDIERED A REGION, AND WE WANT TO SAY IT BOUNDS THE
     *  VOID REGION
     *
     */
    if(*i == VOID_REGION)
      {
      //we represents void with -1, while moab represents it with zero
      returnValues[VOID_REGION]= -1;
      continue;
      }

    smoab::Range parents = interface.findParentEntities(*i);
    T currentParentId = newParentId;
    if(!parents.empty())
      {
      //we presume that the intersection will only return a single value
      smoab::Range s = smoab::intersect(parents,realParents);
      const smoab::EntityHandle& p = s.front();

      //check if we have this parent in the map already
      parentIterator fp = foundParents.find(p);
      if(fp == foundParents.end())
        {
        //this parent id has never been seen before. Add it to the found parent
        //map incase another child is owned by it.
        std::pair<smoab::EntityHandle,T> parentMatch(p,currentParentId);
        foundParents.insert(parentMatch);
        ++newParentId;
        }
      else
        {
        currentParentId = fp->second;
        }
      }
    else
      {
      //no match so we increment the parent id
      ++newParentId;
      }
    std::pair<smoab::EntityHandle,T> match(*i,currentParentId);
    returnValues.insert(match);
    }

  return returnValues;
  }

//----------------------------------------------------------------------------
template<typename T>
std::map<smoab::EntityHandle,T> findDomains(const smoab::CellSets& parentSets,
    const smoab::Tag *domainTag,
    const smoab::Interface& interface)
  {

  //the basic model of attack is as follow.
  //for each entity in the parent sets, find the tag value and convert it to
  //type T

  std::map<smoab::EntityHandle,T> returnValues;

  typedef smoab::CellSets::const_iterator iterator;
  for(iterator i=parentSets.begin(); i != parentSets.end(); ++i)
    {
    std::pair<smoab::EntityHandle,T> key_value;
    key_value.first = i->entity();
    key_value.second = interface.getTagData<T>(*domainTag,i->entity());
    returnValues.insert(key_value);
    }
  return returnValues;
  }

 } } //smoab::model

#endif
