/*
 Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>

 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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
*/

#include "mesh.h"
#include "dolfin_compat.h"
#ifdef LATEST_DOLFIN
#include <dolfin/mesh/MeshPartitioning.h>
#endif

DEFUN_DLD (Mesh, args, ,"-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{mesh_out}]} = \
Mesh (@var{mesh_in}) \n\
Construct a mesh from file or from (p, e, t) format.\n\
The @var{mesh_in} should be either\n\
@itemize @bullet \n\
@item a string containing the name of the file where the mesh \
is stored in .xml file\n\
If the file is not a .xml file you can try to use the command\
dolfin-convert directly from the terminal. \n\
@item a a PDE-tool like structure with matrix fields (p,e,t)\
@end itemize\n\
The output @var{mesh_out} is a representation of the\
@var{mesh_in} which is compatible with fem-fenics.\n\
The easiest way for dealing with meshes is using the msh pkg. \n\
@seealso{FunctionSpace}\n\
@end deftypefn")
{
  int nargin = args.length ();
  octave_value retval = 0;

  if (nargin < 1 || nargin > 1)
    print_usage ();
  else
    {
      if (!error_state)
        {

          if (! mesh_type_loaded)
            {
              mesh::register_type ();
              mesh_type_loaded = true;
              mlock ();
            }

          if (args(0).is_string () == true)
            {
              std::string filename = args(0).string_value ();
              //if the filename is not valid, dolfin takes care of it
              retval = new mesh (filename);
            }

          else if (args(0).is_map () == true)
            {
              octave_scalar_map a = args(0).scalar_map_value ();
              Array<double> p = a.contents ("p").matrix_value ();
              Array<octave_idx_type> t = a.contents ("t").matrix_value ();
              Array<octave_idx_type> e = a.contents ("e").matrix_value ();

              if (! error_state)
                {
                  retval = new mesh (p, e, t);
                }
            }

          else
             error ("Mesh: the argument you provide is invalid");

        }
    }
  return retval;
}

mesh::mesh (Array<double>& p,
            Array<octave_idx_type>& e,
            Array<octave_idx_type>& t)
{
  std::size_t D = p.rows ();
  if (D < 2 || D > 3)
    error ("Mesh constructor: only 2D or 3D meshes are supported");
  else
    {
      dolfin::MeshEditor editor;
      SHARED_PTR <dolfin::Mesh> msh (new dolfin::Mesh ());
      editor.open (*msh, D, D);
      editor.init_vertices (p.cols ());
      editor.init_cells (t.cols ());

      if (D == 2)
        {
          for (uint i = 0; i < p.cols (); ++i)
            editor.add_vertex (i,
                               p.xelem (0, i),
                               p.xelem (1, i));

          for (uint i = 0; i < t.cols (); ++i)
            editor.add_cell (i,
                             t.xelem (0, i) - 1,
                             t.xelem (1, i) - 1,
                             t.xelem (2, i) - 1);
        }

      if (D == 3)
        {
          for (uint i = 0; i < p.cols (); ++i)
            editor.add_vertex (i,
                               p.xelem (0, i),
                               p.xelem (1, i),
                               p.xelem (2, i));

          for (uint i = 0; i < t.cols (); ++i)
            editor.add_cell (i,
                             t.xelem (0, i) - 1,
                             t.xelem (1, i) - 1,
                             t.xelem (2, i) - 1,
                             t.xelem (3, i) - 1);
        }

      editor.close ();

      // store information associated with e
      msh->init (D - 1);
      dolfin::MeshValueCollection<std::size_t> facet(*msh, D - 1);
      std::size_t num_side_edges = e.cols ();

      unsigned const num_procs =
#ifdef LATEST_DOLFIN
        dolfin::MPI::size (MPI_COMM_WORLD);
#else
        dolfin::MPI::num_processes ();
#endif

      if (num_procs == 1)
        {
          if (D == 2)
            {
              for (uint i = 0; i < num_side_edges; ++i)
                {
                  dolfin::Vertex v (*msh, e.xelem (0, i) - 1);
                  for (dolfin::FacetIterator f (v); ! f.end (); ++f)
                    {
                      if ((*f).entities(0)[0] == e.xelem (0, i) - 1
                          && (*f).entities(0)[1] == e.xelem (1, i) - 1
                          || (*f).entities(0)[0] == e.xelem (1, i) - 1
                          && (*f).entities(0)[1] == e.xelem (0, i) - 1)
                        {
                          std::pair <std::size_t, std::size_t>
                            idxvl ((*f).index (), e.xelem (4, i));
                          msh->domains ().set_marker (idxvl, D - 1);
                          break;
                        }
                    }
                }
            }

          if (D == 3)
            {
              for (uint i = 0; i < num_side_edges; ++i)
                {
                  dolfin::Vertex v (*msh, e.xelem (0, i) - 1);
                  for (dolfin::FacetIterator f (v); ! f.end (); ++f)
                    {
                      if ((*f).entities(0)[0] == e(0, i) - 1
                          && (*f).entities(0)[1] == e.xelem (1, i) - 1
                          && (*f).entities(0)[2] == e.xelem (2, i) - 1
                          || (*f).entities(0)[0] == e.xelem (0, i) - 1
                          && (*f).entities(0)[1] == e.xelem (2, i) - 1
                          && (*f).entities(0)[2] == e.xelem (1, i) - 1
                          || (*f).entities(0)[0] == e.xelem (1, i) - 1
                          && (*f).entities(0)[1] == e.xelem (0, i) - 1
                          && (*f).entities(0)[2] == e.xelem (2, i) - 1
                          || (*f).entities(0)[0] == e.xelem (1, i) - 1
                          && (*f).entities(0)[1] == e.xelem (2, i) - 1
                          && (*f).entities(0)[2] == e.xelem (0, i) - 1
                          || (*f).entities(0)[0] == e.xelem (2, i) - 1
                          && (*f).entities(0)[1] == e.xelem (0, i) - 1
                          && (*f).entities(0)[2] == e.xelem (1, i) - 1
                          || (*f).entities(0)[0] == e.xelem (2, i) - 1
                          && (*f).entities(0)[1] == e.xelem (1, i) - 1
                          && (*f).entities(0)[2] == e.xelem (0, i) - 1)
                        {
                          std::pair <std::size_t, std::size_t> 
                            idxvl ((*f).index (), e.xelem (9, i));
                          msh->domains ().set_marker (idxvl, D - 1);
                          break;
                        }
                    }
                }
            }
        }


      dolfin::MeshValueCollection<std::size_t> cell (*msh, D);
      std::size_t num_cells = t.cols ();

      if (num_procs == 1)
        {
          if (D == 2)
            {
              for (uint i = 0; i < num_cells; ++i)
                {
                  dolfin::Vertex v (*msh, t.xelem (0, i) - 1);
                  for (dolfin::CellIterator f (v); ! f.end (); ++f)
                    {
                      if ((*f).entities(0)[0] == t.xelem (0, i) - 1
                          && (*f).entities(0)[1] == t.xelem (1, i) - 1
                          && (*f).entities(0)[2] == t.xelem (2, i) - 1
                          || (*f).entities(0)[0] == t.xelem (0, i) - 1
                          && (*f).entities(0)[1] == t.xelem (2, i) - 1
                          && (*f).entities(0)[2] == t.xelem (1, i) - 1
                          || (*f).entities(0)[0] == t.xelem (1, i) - 1
                          && (*f).entities(0)[1] == t.xelem (0, i) - 1
                          && (*f).entities(0)[2] == t.xelem (2, i) - 1
                          || (*f).entities(0)[0] == t.xelem (1, i) - 1
                          && (*f).entities(0)[1] == t.xelem (2, i) - 1
                          && (*f).entities(0)[2] == t.xelem (0, i) - 1
                          || (*f).entities(0)[0] == t.xelem (2, i) - 1
                          && (*f).entities(0)[1] == t.xelem (0, i) - 1
                          && (*f).entities(0)[2] == t.xelem (1, i) - 1
                          || (*f).entities(0)[0] == t.xelem (2, i) - 1
                          && (*f).entities(0)[1] == t.xelem (1, i) - 1
                          && (*f).entities(0)[2] == t.xelem (0, i) - 1)
                        {
                          std::pair <std::size_t, std::size_t>
                            idxvl ((*f).index (), t.xelem (3, i));
                          msh->domains ().set_marker (idxvl, D);
                          break;
                        }
                    }
                }
            }

          if (D == 3)
            {
              for (uint i = 0; i < num_cells; ++i)
                {
                  dolfin::Vertex v (*msh, t.xelem (0, i) - 1);
                  for (dolfin::CellIterator f (v); ! f.end (); ++f)
                    {
                      if ((*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1
                           || (*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (0, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1

                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1
                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (1, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1

                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1
                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (3, i) - 1
                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (3, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1
                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1
                           || (*f).entities(0)[0] == t.xelem (2, i) - 1
                           && (*f).entities(0)[1] == t.xelem (3, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1

                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (0, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1
                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (2, i) - 1
                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (1, i) - 1
                           && (*f).entities(0)[2] == t.xelem (2, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1
                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (0, i) - 1
                           && (*f).entities(0)[3] == t.xelem (1, i) - 1
                           || (*f).entities(0)[0] == t.xelem (3, i) - 1
                           && (*f).entities(0)[1] == t.xelem (2, i) - 1
                           && (*f).entities(0)[2] == t.xelem (1, i) - 1
                           && (*f).entities(0)[3] == t.xelem (0, i) - 1)
                        {
                          std::pair <std::size_t, std::size_t>
                            idxvl ((*f).index (), t.xelem (4, i));
                          msh->domains ().set_marker (idxvl, D);
                          break;
                        }
                    }
                }
            }
        }

      dolfin::MeshPartitioning::build_distributed_mesh (*msh);
      pmsh = msh;
    }
}