#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: meshPeriodic.c,v 1.00 2001/09/25 06:48:57 knepley Exp $";
#endif

/*
     Defines the interface to functions on a periodic mesh
*/

#include "src/mesh/meshimpl.h"    /*I "mesh.h" I*/

#undef  __FUNCT__
#define __FUNCT__ "MeshIsPeriodic"
/*@
  MeshIsPeriodic - This function returns the periodicity of the mesh.

  Not collective

  Input Parameter:
. mesh - The mesh

  Output Parameter:
. per  - The flag for periodicity

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshIsPeriodicDimension(), MeshSetPeriodicDimension(), MeshPeriodicX(), MeshPeriodicRelativeX(), MeshPeriodicDiffX()
@*/
int MeshIsPeriodic(Mesh mesh, PetscTruth *per)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(per);
  *per = mesh->isPeriodic;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshIsPeriodicDimension"
/*@
  MeshIsPeriodicDimension - This function returns the periodicity of a given dimension the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- dir  - The coordinate direction

  Output Parameter:
. per  - The flag for periodicity

  Level: intermediate

.keywords: mesh, coordinates, periodic. dimension
.seealso: MeshIsPeriodic(), MeshSetPeriodicDimension(), MeshPeriodicX(), MeshPeriodicRelativeX(), MeshPeriodicDiffX()
@*/
int MeshIsPeriodicDimension(Mesh mesh, int dir, PetscTruth *per)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(per);
  if ((dir < 0) || (dir >= mesh->dim)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Direction %d does not appear in this mesh", dir)
  *per = mesh->isPeriodicDim[dir];
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetPeriodicDimension"
/*@
  MeshSetPeriodicDimension - This function sets the periodicity of a given dimension the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
. dir  - The coordinate direction
- per  - The flag for periodicity

  Level: intermediate

.keywords: mesh, coordinates, periodic. dimension
.seealso: MeshIsPeriodic(), MeshIsPeriodicDimension(), MeshPeriodicX(), MeshPeriodicRelativeX(), MeshPeriodicDiffX()
@*/
int MeshSetPeriodicDimension(Mesh mesh, int dir, PetscTruth per)
{
  int d;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if ((dir < 0) || (dir >= mesh->dim)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Direction %d does not appear in this mesh", dir)
  mesh->isPeriodicDim[dir] = per;
  for(mesh->isPeriodic = PETSC_FALSE, d = 0; d < mesh->dim; d++) {
    if (mesh->isPeriodicDim[d] == PETSC_TRUE) mesh->isPeriodic = PETSC_TRUE;
  }
  PetscFunctionReturn(0);
}

/* If optimized, these functions are replaced by macros in $PETSC_DIR/include/mesh.h */
#ifndef MESH_PERIODIC_OPTIMIZED

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicX"
/*@C
  MeshPeriodicX - Returns the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- x    - The original coordinate

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicY(), MeshPeriodicZ(), MeshPeriodicRelativeX(), MeshPeriodicDiffX()
@*/
double MeshPeriodicX(Mesh mesh, double x)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[0] == PETSC_FALSE)
    PetscFunctionReturn(x);
  if (x - PETSC_MACHINE_EPSILON > mesh->endX) {
    if (x - mesh->sizeX > mesh->endX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", x);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(x - mesh->sizeX);
  } else if (x + PETSC_MACHINE_EPSILON < mesh->startX) {
    if (x + mesh->sizeX < mesh->startX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", x);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(x + mesh->sizeX);
  } else {
    PetscFunctionReturn(x);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicY"
/*@C
  MeshPeriodicY - Returns the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- y    - The original coordinate

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicX(), MeshPeriodicZ(), MeshPeriodicRelativeY(), MeshPeriodicDiffY()
@*/
double MeshPeriodicY(Mesh mesh, double y)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[1] == PETSC_FALSE)
    PetscFunctionReturn(y);
  if (y - PETSC_MACHINE_EPSILON > mesh->endY) {
    if (y - mesh->sizeY > mesh->endY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", y);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(y - mesh->sizeY);
  } else if (y + PETSC_MACHINE_EPSILON < mesh->startY) {
    if (y + mesh->sizeY < mesh->startY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", y);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(y + mesh->sizeY);
  } else {
    PetscFunctionReturn(y);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicZ"
/*@C
  MeshPeriodicZ - Returns the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- z    - The original coordinate

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicX(), MeshPeriodicY(), MeshPeriodicRelativeZ(), MeshPeriodicDiffZ()
@*/
double MeshPeriodicZ(Mesh mesh, double z)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[2] == PETSC_FALSE)
    PetscFunctionReturn(z);
  if (z - PETSC_MACHINE_EPSILON > mesh->endZ) {
    if (z - mesh->sizeZ > mesh->endZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", z);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(z - mesh->sizeZ);
  } else if (z + PETSC_MACHINE_EPSILON < mesh->startZ) {
    if (z + mesh->sizeZ < mesh->startZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one\n", z);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(z + mesh->sizeZ);
  } else {
    PetscFunctionReturn(z);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicRelativeX"
/*@C
  MeshPeriodicRelativeX - Returns the value modulo the period of the mesh, relative
  to a given coordinate.

  Not collective

  Input Parameters:
+ mesh - The mesh
. x    - The original coordinate
- xO   - The origin defining the period

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

  Note:
  This function is normally used to have a geometric object reside in a single
  period of the mesh. For instance, a triangle that needs to be drawn.

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicRelativeY(), MeshPeriodicRelativeZ(), MeshPeriodicX(), MeshPeriodicDiffX()
@*/
double MeshPeriodicRelativeX(Mesh mesh, double x, double xO)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[0] == PETSC_FALSE)
    PetscFunctionReturn(x);
  if (x - PETSC_MACHINE_EPSILON > xO + 0.5*mesh->sizeX) {
    if (x - mesh->sizeX > xO + mesh->sizeX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", x, xO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(x - mesh->sizeX);
  } else if (x + PETSC_MACHINE_EPSILON < xO - 0.5*mesh->sizeX) {
    if (x + mesh->sizeX < xO - mesh->sizeX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", x, xO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(x + mesh->sizeX);
  } else {
    PetscFunctionReturn(x);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicRelativeY"
/*@C
  MeshPeriodicRelativeY - Returns the value modulo the period of the mesh, relative
  to a given coordinate.

  Not collective

  Input Parameters:
+ mesh - The mesh
. y    - The original coordinate
- yO   - The origin defining the period

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

  Note:
  This function is normally used to have a geometric object reside in a single
  period of the mesh. For instance, a triangle that needs to be drawn.

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicRelativeX(), MeshPeriodicRelativeZ(), MeshPeriodicY(), MeshPeriodicDiffY()
@*/
double MeshPeriodicRelativeY(Mesh mesh, double y, double yO)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[1] == PETSC_FALSE)
    PetscFunctionReturn(y);
  if (y - PETSC_MACHINE_EPSILON > yO + 0.5*mesh->sizeY) {
    if (y - mesh->sizeY > yO + mesh->sizeY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", y, yO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(y - mesh->sizeY);
  } else if (y + PETSC_MACHINE_EPSILON < yO - 0.5*mesh->sizeY) {
    if (y + mesh->sizeY < yO - mesh->sizeY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", y, yO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(y + mesh->sizeY);
  } else {
    PetscFunctionReturn(y);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicRelativeZ"
/*@C
  MeshPeriodicRelativeZ - Returns the value modulo the period of the mesh, relative
  to a given coordinate.

  Not collective

  Input Parameters:
+ mesh - The mesh
. z    - The original coordinate
- zO   - The origin defining the period

  Output Parameter:
. ret  - The normalized coordinate

  Level: intermediate

  Note:
  This function is normally used to have a geometric object reside in a single
  period of the mesh. For instance, a triangle that needs to be drawn.

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicRelativeX(), MeshPeriodicRelativeY(), MeshPeriodicZ(), MeshPeriodicDiffZ()
@*/
double MeshPeriodicRelativeZ(Mesh mesh, double z, double zO)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[2] == PETSC_FALSE)
    PetscFunctionReturn(z);
  if (z - PETSC_MACHINE_EPSILON > zO + 0.5*mesh->sizeZ) {
    if (z - mesh->sizeZ > zO + mesh->sizeZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", z, zO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(z - mesh->sizeZ);
  } else if (z + PETSC_MACHINE_EPSILON < zO - 0.5*mesh->sizeZ) {
    if (z + mesh->sizeZ < zO - mesh->sizeZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Coordinate %g has winding number greater than one about %g\n", z, zO);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(z + mesh->sizeZ);
  } else {
    PetscFunctionReturn(z);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicDiffX"
/*@C
  MeshPeriodicDiffX - Takes the difference of two x-coordinates in the domain, and returns
  the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- diff - The original difference

  Output Parameter:
. ret  - The normalized difference

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicDiffY(), MeshPeriodicDiffZ(), MeshPeriodicX(), MeshPeriodicRelativeX()
@*/
double MeshPeriodicDiffX(Mesh mesh, double diff)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[0] == PETSC_FALSE)
    PetscFunctionReturn(diff);
  if (diff - PETSC_MACHINE_EPSILON > 0.5*mesh->sizeX) {
    if (diff - mesh->sizeX > 0.5*mesh->sizeX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff - mesh->sizeX);
  } else if (diff + PETSC_MACHINE_EPSILON < -0.5*mesh->sizeX) {
    if (diff + mesh->sizeX < -0.5*mesh->sizeX) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff + mesh->sizeX);
  } else {
    PetscFunctionReturn(diff);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicDiffY"
/*@C
  MeshPeriodicDiffY - Takes the difference of two y-coordinates in the domain, and returns
  the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- diff - The original difference

  Output Parameter:
. ret  - The normalized difference

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicDiffX(), MeshPeriodicDiffZ(), MeshPeriodicY(), MeshPeriodicRelativeY()
@*/
double MeshPeriodicDiffY(Mesh mesh, double diff)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[1] == PETSC_FALSE)
    PetscFunctionReturn(diff);
  if (diff - PETSC_MACHINE_EPSILON > 0.5*mesh->sizeY) {
    if (diff - mesh->sizeY > 0.5*mesh->sizeY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff - mesh->sizeY);
  } else if (diff + PETSC_MACHINE_EPSILON < -0.5*mesh->sizeY) {
    if (diff + mesh->sizeY < -0.5*mesh->sizeY) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff + mesh->sizeY);
  } else {
    PetscFunctionReturn(diff);
  }
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPeriodicDiffZ"
/*@C
  MeshPeriodicDiffZ - Takes the difference of two z-coordinates in the domain, and returns
  the value modulo the period of the mesh.

  Not collective

  Input Parameters:
+ mesh - The mesh
- diff - The original difference

  Output Parameter:
. ret  - The normalized difference

  Level: intermediate

.keywords: mesh, coordinates, periodic
.seealso: MeshPeriodicDiffX(), MeshPeriodicDiffY(), MeshPeriodicZ(), MeshPeriodicRelativeZ()
@*/
double MeshPeriodicDiffZ(Mesh mesh, double diff)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->isPeriodicDim[2] == PETSC_FALSE)
    PetscFunctionReturn(diff);
  if (diff - PETSC_MACHINE_EPSILON > 0.5*mesh->sizeZ) {
    if (diff - mesh->sizeZ > 0.5*mesh->sizeZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff - mesh->sizeZ);
  } else if (diff + PETSC_MACHINE_EPSILON < -0.5*mesh->sizeZ) {
    if (diff + mesh->sizeZ < -0.5*mesh->sizeZ) {
      PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference %g has winding number greater than one\n", diff);
      PetscFunctionReturn(0.0);
    }
    PetscFunctionReturn(diff + mesh->sizeZ);
  } else {
    PetscFunctionReturn(diff);
  }
}

#endif /* MESH_PERIODIC_OPTIMIZED */