//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane Del Pino

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

//  This program 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: Structured3DMesh.hpp,v 1.10 2004/05/31 18:49:32 delpinux Exp $

// This class allow 3D Structured Mesh management.

#ifndef  _STRUCTURED_3D_MESH_HPP_
#define  _STRUCTURED_3D_MESH_HPP_

#include <Mesh.hpp>

#include <string>

#include <TinyVector.hpp>

#include <Domain.hpp>
#include <CartesianHexahedron.hpp>

#include <SurfaceMeshOfQuadrangles.hpp>

#include <Index.hpp>

#include <Structured3DMeshShape.hpp>

#include <Structured3DVector.hpp>

class MeshOfHexahedra;

/*!
  \class Structured3DMesh
  This class describes Structured 3d Mesh.

  \author Stphane Del Pino
*/
class Structured3DMesh
  : public Mesh
{
public:
  typedef CartesianHexahedron ElementGeometry;

  /**
   *  If a transformation is applied the new mesh type will be
   *  MeshOfHexahedra
   */
  typedef MeshOfHexahedra Transformed;
  typedef SurfaceMeshOfQuadrangles BorderMeshType;

  enum {
    family = Mesh::volume
  };
private:
  //! The shape of the vertices' set
  const Structured3DMeshShape __verticesShape;

  //! The shape of the cells' set only the array3dshape is needed.
  const Array3DShape __cellShape;

  //! The cells set
  Vector<CartesianHexahedron> __cells;

  ReferenceCounting<BorderMeshType> __surfaceMesh; /**< The various border lists*/

  /*! Copy constructor must not be called, so it is private and generate an
    error if called. This is implemented so that the compiler does not
    generates one automagicaly.
   */
  Structured3DMesh(const Structured3DMesh& M)
  	: Mesh(M),
	  __verticesShape(M.__verticesShape),
	  __cellShape(M.__cellShape)
  {  
    ;
  }

public:
  inline bool hasSurfaceMesh() const
  {
    return (__surfaceMesh != 0);
  }

  inline BorderMeshType& surfaceMesh()
  {
    return (*__surfaceMesh);
  }

  inline const BorderMeshType& surfaceMesh() const
  {
    return (*__surfaceMesh);
  }

  //! Read-only access to the number of cells.
  inline const size_t& numberOfCells() const
  {
    return __cells.size();
  }

  typedef Mesh::T_iterator<Structured3DMesh, CartesianHexahedron> iterator;
  typedef Mesh::T_iterator<const Structured3DMesh, const CartesianHexahedron> const_iterator;

  inline Structured3DMesh::const_iterator find(const double& x,
					       const double& y,
					       const double& z) const;

  inline Structured3DMesh::const_iterator find(const TinyVector<3>& X) const
  {
    return find(X[0],X[1],X[2]);
  }

  //! Writes the mesh to the OpenDX Mesh format
  void plotOpenDX(std::ostream& fmesh,
		  const std::string& CR) const;

  //! Access to the shape of the mesh.
  const Structured3DMeshShape& shape()const
  {
    return __verticesShape;
  }

  //! Access to the Vertex (\a i,\a j,\a k) of the mesh.
  inline Vertex& vertex(const size_t i,
			const size_t j,
			const size_t k);

  //! Read-only access to the Vertex (\a i,\a j,\a k) of the mesh.
  inline const Vertex& vertex(const size_t i,
			      const size_t j,
			      const size_t k) const;

  //! Access to the Cell (\a i,\a j,\a k) of the mesh.
  inline CartesianHexahedron& cell(const size_t i,
				   const size_t j,
				   const size_t k);

  //! Read-only access to the Cell (\a i,\a j,\a k) of the mesh.
  inline const CartesianHexahedron& cell(const size_t i,
					 const size_t j,
					 const size_t k) const;

  //! Access to the Cell \a i of the mesh.
  inline CartesianHexahedron& cell(const size_t i)
  {
    return __cells[i];
  }

  //! Read-only access to the Cell \a i of the mesh.
  inline const CartesianHexahedron& cell(const size_t i) const
  {
    return __cells[i];
  }

  /*! Access to the ith Vertex.
    \warning this function should be redefined here, since it is
    overloaded.
   */
  inline Vertex& vertex(const size_t& i)
  {
    return Mesh::vertex(i);
  }

  /*! Read-only access to the ith Vertex.
    \warning this function should be redefined here, since it is
    overloaded.
   */
  inline const Vertex& vertex(const size_t& i) const
  {
    return Mesh::vertex(i);
  }

  //! Access to the Vertex indexed by \a I.
  inline Vertex& vertex(const Index& I);

  //! Read-only access to the Vertex indexed by \a I.
  inline const Vertex& vertex(const Index& I) const;

  //! Access to the Cell indexed by \a I.
  inline CartesianHexahedron& cell(const Index& I);

  //! Read-only access to the Cell indexed by \a I.
  inline const CartesianHexahedron& cell(const Index& I) const;

  /*! Returns the Index of the Cell contaning the Vertex \a V if none, returns
    error!
   */
  inline const Index cellIndex(const Vertex& V) const;

  /*! Returns the Index of the Cell contaning the point \a P if none, returns
    error!
   */
  inline const Index cellIndex(const TinyVector<3>& V) const;

  /*! Return the Index of \a V's closest Vertex in SMesh if none, return
    error!
  */
  inline const Index vertexIndex(const Vertex& V) const;

  /*! Saves the mesh into \a filename using gnuplot format.
    \warning not really used, should be removed.
   */
  void gnuplot(std::string filename);

  //! Generates references associated to a Shape
  void createReferences(const Shape& S);

  //! Returns \p true if the point \a p is inside the mesh.
  inline bool inside(const real_t& x, const real_t& y, const real_t& z) const
  {
    return ((x>=__verticesShape.a(0))&&
	    (x<=__verticesShape.b(0))&&
	    (y>=__verticesShape.a(1))&&
	    (y<=__verticesShape.b(1))&&
	    (z>=__verticesShape.a(2))&&
	    (z<=__verticesShape.b(2)));
  }

  //! Returns \p true if the point \a p is inside the mesh.
  inline bool inside(const TinyVector<3>& p) const
  {
    return this->inside(p[0], p[1], p[2]);
  }

  /*! Constructs a Structured3DMesh using \a n0 along the X-axis, \a n1 along
    the Y-axis and \a n2 along the Z-axis. The vertices \a v0 and \a v1 are
    two opposed corners of the mesh, and the edges of the mesh are parallel to
    the axis.
   */
  Structured3DMesh(const Structured3DMeshShape& s3dM);

};

//! Access to the Vertex (\a i,\a j,\a k) of the mesh.
inline Vertex& Structured3DMesh::vertex(const size_t i,
					const size_t j,
					const size_t k)
{
  return (*__verticesSet)[__verticesShape(i, j, k)];
}

//! Read-only access to the Vertex (\a i,\a j,\a k) of the mesh.
inline const Vertex& Structured3DMesh::vertex(const size_t i,
					      const size_t j,
					      const size_t k) const
{
  return (*__verticesSet)[__verticesShape(i, j, k)];
}

//! Access to the Cell (\a i, \a j, \a k).
inline CartesianHexahedron& Structured3DMesh::cell(const size_t i,
						   const size_t j,
						   const size_t k)
{
  return (__cells[__cellShape(i, j, k)]);
}

//! Read-only access to the Cell (\a i, \a j, \a k).
inline const CartesianHexahedron& Structured3DMesh::cell(const size_t i,
							 const size_t j,
							 const size_t k) const
{
  return (__cells[__cellShape(i, j, k)]);
}
//! Access to the Vertex indexed by \a I.
inline Vertex& Structured3DMesh::vertex(const Index& I)
{
  return (*__verticesSet)[__verticesShape(I[0], I[1], I[2])];
}

//! Read-only access to the Vertex indexed by \a I.
inline const Vertex& Structured3DMesh::vertex(const Index& I) const
{
  return (*__verticesSet)[__verticesShape(I[0], I[1], I[2])];
}

//! Access to the Cell indexed by \a I.
inline CartesianHexahedron& Structured3DMesh::cell(const Index& I)
{
  return (__cells[__cellShape(I[0], I[1], I[2])]);
}

//! Read-only access to the Cell indexed by \a I.
inline const CartesianHexahedron& Structured3DMesh::cell(const Index& I) const
{
  return (__cells[__cellShape(I[0], I[1], I[2])]);
}

/*! Returns the Index of the Cell contaning the point \a P if none, returns
  error!
*/
inline const Index Structured3DMesh::cellIndex(const TinyVector<3>& V) const
{
  const real_t& x = V[0];
  const real_t& y = V[1];
  const real_t& z = V[2];

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int(std::floor((x - x0)/__verticesShape.hx()));
  int j = int(std::floor((y - y0)/__verticesShape.hy()));
  int k = int(std::floor((z - z0)/__verticesShape.hz()));

  bool foundCell = inside(x,y,z);

  if (!foundCell) {
    i = (i<0) ? 0 : i;
    j = (j<0) ? 0 : j;
    k = (k<0) ? 0 : k;
  }

  int nx_1 = __cellShape.nx() - 1;
  int ny_1 = __cellShape.ny() - 1;
  int nz_1 = __cellShape.nz() - 1;

  i = (i>nx_1) ? nx_1 : i;
  j = (j>ny_1) ? ny_1 : j;
  k = (k>nz_1) ? nz_1 : k;

  Index I(i,j,k);
  I.isGood() = foundCell;

  return I;
}

/*! Return the Index of \a V's closest Vertex in SMesh if none, return
  error!
*/
inline const Index Structured3DMesh::vertexIndex(const Vertex& V) const
{
  const real_t& x = V.x();
  const real_t& y = V.y();
  const real_t& z = V.z();

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int((x - x0)/__verticesShape.hx()+0.5);
  int j = int((y - y0)/__verticesShape.hy()+0.5);
  int k = int((z - z0)/__verticesShape.hz()+0.5);

  bool foundCell = true;
  if (((i<0)||(i>static_cast<int>(__verticesShape.nx()-1)))
      &&((j<0)||(j>static_cast<int>(__verticesShape.ny()-1)))
      &&((k<0)||(k>static_cast<int>(__verticesShape.nz()-1)))) {
	foundCell = false;
  }

  Index I(i,j,k);
  I.isGood() = foundCell;

  return I;
}

inline Structured3DMesh::const_iterator 
Structured3DMesh::find(const double& x,
		       const double& y,
		       const double& z) const
{
  bool foundCell = inside(x,y,z);

  if (!foundCell)
    return Structured3DMesh::const_iterator(*this,
					    Structured3DMesh::const_iterator::End);

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int(std::floor((x - x0)/__verticesShape.hx()));
  int j = int(std::floor((y - y0)/__verticesShape.hy()));
  int k = int(std::floor((z - z0)/__verticesShape.hz()));

  {
    int nx_1 = __cellShape.nx() - 1;
    int ny_1 = __cellShape.ny() - 1;
    int nz_1 = __cellShape.nz() - 1;

    // this transformation is to operate if the (x,y,z) is a point on
    // one of the faces of the box.

    i = (i>nx_1) ? nx_1 : i;
    j = (j>ny_1) ? ny_1 : j;
    k = (k>nz_1) ? nz_1 : k;
  }

  return Structured3DMesh::const_iterator(*this,
					  __cellShape(i, j, k));
}

#endif // _STRUCTURED_3D_MESH_HPP_