// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */ 
/* License as published by the Free Software Foundation either */ 
/* version 3.0 of the License, or (at your option) any later version. */

/* This software is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/. */

#include "CollapsedMapping.hpp"

void CollapsedMapping::
mapTriNE(PointVector& tpPointsN)
{
  for(unsigned int i = 0; i<tpPointsN.rows(); i++)
  {
    tpPointsN(i,0) = 0.5 * (1.0 + tpPointsN(i,0)) * (1.0 - tpPointsN(i,1)) - 1.;
    tpPointsN(i,1) = tpPointsN(i,1);
  }
}

void CollapsedMapping::
mapTetNE(PointVector& tpPointsN)
{
  for(unsigned int i = 0; i<tpPointsN.rows(); i++)
  {
    tpPointsN(i,0) = 0.25 * (1. + tpPointsN(i,0)) * (1. - tpPointsN(i,1)) * (1. - tpPointsN(i,2)) - 1.;
    tpPointsN(i,1) = 0.5 * (1. + tpPointsN(i,1)) * (1. - tpPointsN(i,2)) - 1.;
    tpPointsN(i,2) = tpPointsN(i,2);
  }
}

void CollapsedMapping::
mapPrismNE(PointVector& tpPointsN)
{
  for(unsigned int i = 0; i<tpPointsN.rows(); i++)
  {
    tpPointsN(i,0) = 0.5 * (1. + tpPointsN(i,0)) * (1. - tpPointsN(i,1)) - 1.;
    tpPointsN(i,1) = tpPointsN(i,1);
    tpPointsN(i,2) = tpPointsN(i,2);
  }
}

void CollapsedMapping::
mapPyramidNE(PointVector& tpPointsN)
{
  for(unsigned int i = 0; i<tpPointsN.rows(); i++)
  {
    tpPointsN(i,0) = 0.5 * (1. + tpPointsN(i,0)) * (1. - tpPointsN(i,2)) - 1.;
    tpPointsN(i,1) = 0.5 * (1. + tpPointsN(i,1)) * (1. - tpPointsN(i,2)) - 1.;
    tpPointsN(i,2) = tpPointsN(i,2);
  }
}

void CollapsedMapping::
mapTriEN(PointVector& tpPointsE)
{
  for(unsigned int i = 0; i<tpPointsE.rows(); i++)
  {
    if (tpPointsE(i,1) == 1.)
      tpPointsE(i,0) = -1.;
    else
      tpPointsE(i,0) = 2. * (1.0 + tpPointsE(i,0)) / (1.0 - tpPointsE(i,1)) - 1.;
    tpPointsE(i,1) = tpPointsE(i,1);
  }
}

void CollapsedMapping::
mapPrismEN(PointVector& tpPointsE)
{
  for(unsigned int i = 0; i<tpPointsE.rows(); i++)
  {
    if (tpPointsE(i,1) == 1.)
      tpPointsE(i,0) = -1.;
    else
      tpPointsE(i,0) = 2. * (1.0 + tpPointsE(i,0)) / (1.0 - tpPointsE(i,1)) - 1.;
    tpPointsE(i,1) = tpPointsE(i,1);
    tpPointsE(i,2) = tpPointsE(i,2);
  }
}

void CollapsedMapping::
mapPyramidEN(PointVector& tpPointsE)
{
  for(unsigned int i = 0; i<tpPointsE.rows(); i++)
  {
    if (tpPointsE(i,2) == 1.)
    {
      tpPointsE(i,0) = -1.;
      tpPointsE(i,1) = -1.;
    }
    else
    {
      tpPointsE(i,0) = 2. * (1.0 + tpPointsE(i,0)) / (1.0 - tpPointsE(i,2)) - 1.;
      tpPointsE(i,1) = 2. * (1.0 + tpPointsE(i,1)) / (1.0 - tpPointsE(i,2)) - 1.;
    }
    tpPointsE(i,2) = tpPointsE(i,2);
  }
}

void CollapsedMapping::
mapTetEN(PointVector& tpPointsE)
{
  for(unsigned int i = 0; i<tpPointsE.rows(); i++)
  {
    if (tpPointsE(i,1) + tpPointsE(i,2) == 0.)
      tpPointsE(i,0) = 1.;
    else
      tpPointsE(i,0) = -2. * (1. + tpPointsE(i,0)) / (tpPointsE(i,1) + tpPointsE(i,2)) - 1.;
    if (tpPointsE(i,2) == 1.)
      tpPointsE(i,1) = -1.;
    else
      tpPointsE(i,1) = 2. * (1. + tpPointsE(i,1)) / (1. - tpPointsE(i,2)) - 1.;
    tpPointsE(i,2) = tpPointsE(i,2);
  }
}

void CollapsedMapping::
mapEN(const PointVector& points, PointVector& mappedPoints, const ElemType itemType)
{
  switch (itemType)
  {
    case TET4:
    case TET10:
    { 
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        if (points(i,1) + points(i,2) == 0.)
          mappedPoints(i,0) = 1.;
        else
          mappedPoints(i,0) = -2. * (1. + points(i,0)) / (points(i,1) + points(i,2)) - 1.;
        if (points(i,2) == 1.)
          mappedPoints(i,1) = -1.;
        else
          mappedPoints(i,1) = 2. * (1. + points(i,1)) / (1. - points(i,2)) - 1.;
        mappedPoints(i,2) = points(i,2);
      }
      return;
    } 
    case PRISM6:
    case PRISM15:
    { 
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        if (points(i,1) == 1.)
          mappedPoints(i,0) = -1.;
        else
          mappedPoints(i,0) = 2. * (1.0 + points(i,0)) / (1.0 - points(i,1)) - 1.;
        mappedPoints(i,1) = points(i,1);
        mappedPoints(i,2) = points(i,2);
      }
      return;
    } 
    case PYRAMID5:
    case PYRAMID13:
    {
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        if (points(i,2) == 1.)
        {
          mappedPoints(i,0) = -1.;
          mappedPoints(i,1) = -1.;
        }
        else
        {
          mappedPoints(i,0) = 2. * (1.0 + points(i,0)) / (1.0 - points(i,2)) - 1.;
          mappedPoints(i,1) = 2. * (1.0 + points(i,1)) / (1.0 - points(i,2)) - 1.;
        }
        mappedPoints(i,2) = points(i,2);
      }
      return;
    }
    case TRI3:
    case TRI6:
    {
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        if (points(i,1) == 1.)
          mappedPoints(i,0) = -1.;
        else
          mappedPoints(i,0) = 2. * (1.0 + points(i,0))/(1.0 - points(i,1)) - 1.;
        mappedPoints(i,1) = points(i,1);
      }
      return; 
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    case HEX8:
    case HEX20:
    case HEX27:
    case EDGE2:
    case EDGE3:
    {
      mappedPoints = points;
      return;
    }
    default:
    {
      std::cerr<< "Collapsed Mapping Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

void CollapsedMapping::
mapNE(const PointVector& points, PointVector& mappedPoints, const ElemType itemType)
{
  switch (itemType)
  {
    case TET4:
    case TET10:
    { 
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        mappedPoints(i,0) = 0.25 * (1. + points(i,0)) * (1. - points(i,1)) * (1. - points(i,2)) - 1;
        mappedPoints(i,1) = 0.5 * (1. + points(i,1)) * (1. - points(i,2)) - 1;
        mappedPoints(i,2) = points(i,2);
      }
      return;
    } 
    case PRISM6:
    case PRISM15:
    { 
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        mappedPoints(i,0) = 0.5 * (1. + points(i,0)) * (1. - points(i,1)) - 1.;
        mappedPoints(i,1) = points(i,1);
        mappedPoints(i,2) = points(i,2);
      }
      return;
    }
    case PYRAMID5:
    case PYRAMID13:
    { 
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        mappedPoints(i,0) = 0.5 * (1. + points(i,0)) * (1. - points(i,2)) - 1.;
        mappedPoints(i,1) = 0.5 * (1. + points(i,1)) * (1. - points(i,2)) - 1.;
        mappedPoints(i,2) = points(i,2);
      }
      return;
    }
    case TRI3:
    case TRI6:
    {
      for(unsigned int i = 0; i<points.rows(); i++)
      {
        mappedPoints(i,0) = 0.5 * (1.0 + points(i,0)) * (1.0 - points(i,1)) - 1;
        mappedPoints(i,1) = points(i,1);
      }
      return;
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    case HEX8:
    case HEX20:
    case HEX27:
    case EDGE2:
    case EDGE3:
    {
      mappedPoints = points;
      return;
    }
    default:
    {
      std::cerr<< "Collapsed Mapping Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}
