//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane 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.cpp,v 1.7 2004/05/31 18:49:32 delpinux Exp $


// This class allow to Manipulate 3D Scalar Variable.

#include <fstream>
#include <cmath>

#include <Hexahedron.hpp>
#include <Structured3DMesh.hpp>
#include <Scene.hpp>

//! Writes the mesh to the OpenDX Mesh format
void Structured3DMesh::plotOpenDX(std::ostream& fmsh, const std::string& CR) const
{
  fferr(2) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
  std::exit(1);
}

/*!
  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::Structured3DMesh(const Structured3DMeshShape& SMShape)
  : Mesh(Mesh::cartesianHexahedraMesh,
	 Mesh::volume,
	 SMShape.numberOfVertices()),
    __verticesShape(SMShape),
    __cellShape(__verticesShape.shape()-1),
    __cells(SMShape.numberOfCells())
{
  // w[ijk][01] variable respresent the coefficient needed to 
  // compute vertices position between a and b.
  // for example wj1 is the weight computed for the jth point in the
  // 'y' direction for v1

  const size_t nx = __verticesShape.nx();
  const size_t ny = __verticesShape.ny();
  const size_t nz = __verticesShape.nz();

  const size_t nx_1 = __cellShape.nx();
  const size_t ny_1 = __cellShape.ny();
  const size_t nz_1 = __cellShape.nz();

  for (size_t i=0; i<nx; i++) {
    real_t wi0 = static_cast<real_t>(i) / static_cast<real_t>(nx_1);
    real_t wi1 = static_cast<real_t>(nx_1-i) / static_cast<real_t>(nx_1);;
    for (size_t j=0; j<ny; j++) {
      real_t wj0 = static_cast<real_t>(j) / static_cast<real_t>(ny_1);
      real_t wj1 = static_cast<real_t>(ny_1-j) / static_cast<real_t>(ny_1);
      for (size_t k=0; k<nz; k++) {
	real_t wk0 = static_cast<real_t>(k) / static_cast<real_t>(nz_1);
	real_t wk1 = static_cast<real_t>(nz_1-k)/static_cast<real_t>(nz_1);
	Vertex& v = vertex(i,j,k);
	v.x() = wi0 * __verticesShape.b(0) + wi1 * __verticesShape.a(0);
	v.y() = wj0 * __verticesShape.b(1) + wj1 * __verticesShape.a(1);
	v.z() = wk0 * __verticesShape.b(2) + wk1 * __verticesShape.a(2);
      }
    }
  }

  __edges.reserve(nx_1*ny*nz + nx*ny_1*nz + nx*ny*nz_1);
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny; j++) {
      for (size_t k=0; k<nz; k++) {
	__edges.push_back(Edge(vertex(i,j,k), vertex(i+1,j,k)));
      }
    }
  }

  for (size_t i=0; i<nx; i++) {
    for (size_t j=0; j<ny_1; j++) {
      for (size_t k=0; k<nz; k++) {
	__edges.push_back(Edge(vertex(i,j,k), vertex(i,j+1,k)));
      }
    }
  }

  for (size_t i=0; i<nx; i++) {
    for (size_t j=0; j<ny; j++) {
      for (size_t k=0; k<nz_1; k++) {
	__edges.push_back(Edge(vertex(i,j,k), vertex(i,j,k+1)));
      }
    }
  }

  // Generating Mesh's cells.
  // As Cells generation needs Vertices's references
  // This must be done after.

  size_t n=0;
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      for (size_t k=0; k<nz_1; k++) {
	__cells[n] = CartesianHexahedron(vertex(  i,  j,  k),
					 vertex(i+1,  j,  k),
					 vertex(i+1,j+1,  k),
					 vertex(  i,j+1,  k),
					 vertex(  i,  j,k+1),
					 vertex(i+1,  j,k+1),
					 vertex(i+1,j+1,k+1),
					 vertex(  i,j+1,k+1));
	n++;
      }
    }
  }

  /**
   * Now generates Faces meshes
   * 
   */
  Vector<Quadrangle>* pQuadrangles
    = new Vector<Quadrangle>(2*ny_1*nz_1+2*nx_1*nz_1+2*nx_1*ny_1);
  Vector<Quadrangle>& quadrangles = *pQuadrangles;

  size_t currentCell = 0;
  /// x=xmin
  for (size_t j=0; j<ny_1; j++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(0,  j,  k);
      const Vertex& V1 = vertex(0,  j,k+1);
      const Vertex& V2 = vertex(0,j+1,k+1);
      const Vertex& V3 = vertex(0,j+1,  k);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 0);
      Q.setMother(&(cell(0,j,k)));

      currentCell++;
    }
  }

  /// x=xmax
  for (size_t j=0; j<ny_1; j++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(nx_1,  j,  k);
      const Vertex& V1 = vertex(nx_1,j+1,  k);
      const Vertex& V2 = vertex(nx_1,j+1,k+1);
      const Vertex& V3 = vertex(nx_1,  j,k+1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 1);
      Q.setMother(&(cell(nx_1-1,j,k)));

      currentCell++;
    }
  }

  /// y=ymin
  for (size_t i=0; i<nx_1; i++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(  i,0,  k);
      const Vertex& V1 = vertex(i+1,0,  k);
      const Vertex& V2 = vertex(i+1,0,k+1);
      const Vertex& V3 = vertex(  i,0,k+1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 2);
      Q.setMother(&(cell(i,0,k)));

      currentCell++;
    }
  }

  /// y=ymax
  for (size_t i=0; i<nx_1; i++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(  i,ny_1,  k);
      const Vertex& V1 = vertex(  i,ny_1,k+1);
      const Vertex& V2 = vertex(i+1,ny_1,k+1);
      const Vertex& V3 = vertex(i+1,ny_1,  k);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 3);
      Q.setMother(&(cell(i,ny_1-1,k)));

      currentCell++;
    }
  }

  /// z=zmin
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      const Vertex& V0 = vertex(  i,  j,0);
      const Vertex& V1 = vertex(  i,j+1,0);
      const Vertex& V2 = vertex(i+1,j+1,0);
      const Vertex& V3 = vertex(i+1,  j,0);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 4);
      Q.setMother(&(cell(i,j,0)));

      currentCell++;
    }
  }


  /// z=zmax
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      const Vertex& V0 = vertex(  i,  j,nz_1);
      const Vertex& V1 = vertex(i+1,  j,nz_1);
      const Vertex& V2 = vertex(i+1,j+1,nz_1);
      const Vertex& V3 = vertex(  i,j+1,nz_1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 5);
      Q.setMother(&(cell(i,j,nz_1-1)));
      currentCell++;
    }
  }
  __surfaceMesh = new SurfaceMeshOfQuadrangles(__verticesSet, pQuadrangles);
  (*__surfaceMesh).setBackgroundMesh(this);
}

//! Generates references associated to a POVRef ...
void Structured3DMesh::createReferences(const Shape& S)
{
  for (size_t i=0; i<this->numberOfVertices(); ++i) {
    Vertex& V = vertex(i);

    V.reference() = S.inside(V);
  }
}

//! Saves the mesh into \a filename using gnuplot format.
void Structured3DMesh::gnuplot(std::string filename)
{
  std::ofstream gout(filename.c_str());

  const size_t nx = __verticesShape.nx();
  const size_t ny = __verticesShape.ny();
  const size_t nz = __verticesShape.nz();

  const size_t nx_1 = nx - 1;
  const size_t ny_1 = ny - 1;
  const size_t nz_1 = nz - 1;
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      for (size_t k=0; k<nz_1; k++) {

	const Vertex& v0 = vertex(i,j,k);
	const Vertex& vi = vertex(i+1,j,k);
	const Vertex& vj = vertex(i,j+1,k);
	const Vertex& vk = vertex(i,j,k+1);

	gout << v0 << "\n";
 	gout << vi << "\n\n\n";
 	gout << v0 << "\n";
 	gout << vj << "\n\n\n";
 	gout << v0 << "\n";
 	gout << vk << "\n\n\n";
      }
    }
  }
  

  for (size_t j=0; j<ny_1; j++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& v0 = vertex(nx_1,j,k);
      const Vertex& vj = vertex(nx_1,j+1,k);
      const Vertex& vk = vertex(nx_1,j,k+1);
      gout << v0 << "\n";
      gout << vj << "\n\n\n";
      gout << v0 << "\n";
      gout << vk << "\n\n\n";
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& v0 = vertex(i,  ny_1, k);
      const Vertex& vi = vertex(i+1,ny_1,k);
      const Vertex& vk = vertex(i,  ny_1,k+1);
      gout << v0 << "\n";
      gout << vi << "\n\n\n";
      gout << v0 << "\n";
      gout << vk << "\n\n\n";
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      const Vertex& v0 = vertex(i,  j,  nz_1);
      const Vertex& vi = vertex(i+1,j,  nz_1);
      const Vertex& vj = vertex(i,  j+1,nz_1);
      gout << v0 << "\n";
      gout << vi << "\n\n\n";
      gout << v0 << "\n";
      gout << vj << "\n\n\n";
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    const Vertex& v0 = vertex(i,  ny_1,  nz_1);
    const Vertex& vi = vertex(i+1,ny_1,  nz_1);
    gout << v0 << "\n";
    gout << vi << "\n\n\n";
  }

  for (size_t j=0; j<ny_1; j++) {
    const Vertex& v0 = vertex(nx_1,  j,  nz_1);
    const Vertex& vj = vertex(nx_1,j+1,  nz_1);
    gout << v0 << "\n";
    gout << vj << "\n\n\n";
  }

  for (size_t k=0; k<nz_1; k++) {
    const Vertex& v0 = vertex(nx_1, ny_1,k);
    const Vertex& vk = vertex(nx_1, ny_1,k+1);
    gout << v0 << "\n";
    gout << vk << "\n\n\n";
  }
}