/*
 Copyright (C) 2013 Marco Vassallo

 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 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>

octave_scalar_map
mesh::get_pet () const
{
  const dolfin::Mesh & msh (*pmsh);
  //p matrix
  uint D = msh.topology ().dim ();
  std::size_t num_v = msh.num_vertices ();
  Matrix p (D, num_v);
  std::copy (msh.coordinates ().begin (),
             msh.coordinates ().end (),
             p.fortran_vec ());

  // e has 7 rows in 2d, 10 rows in 3d
  msh.init (D - 1, D);
  std::size_t num_f = msh.num_facets ();
  dim_vector dims;
  dims.resize (2);
  dims(0) = D == 2 ? 7 : 10;
  dims(1) = num_f;
  Array<octave_idx_type> e (dims, 0);
  octave_idx_type *evec = e.fortran_vec ();
  uint D2 = D * D;
  octave_idx_type l = 0, m = 0;

  dolfin::MeshFunction <std::size_t> facet_domains;
  if (! msh.domains ().is_empty ())
      if (msh.domains ().num_marked (D-1) != 0)
        facet_domains = * (msh.domains ().facet_domains ());

  for (dolfin::FacetIterator f (msh); ! f.end (); ++f)
    {
      if ((*f).exterior () == true)
        {
          l = 0;
          for (dolfin::VertexIterator v (*f); ! v.end (); ++v, ++l)
            e.xelem (l, m) = (*v).index () + 1;

          if (! facet_domains.empty ())
            e.xelem (D2, m) = facet_domains[*f];

          ++m;
        }
    }

  dims(1) = m;
  e.resize (dims);

  for (octave_idx_type j = e.rows () - 2;
       j < e.numel () - 2; j += e.rows ())
    evec[j] = 1;

  // t matrix
  dims(0) = D + 2;
  dims(1) = msh.num_cells ();
  Array<octave_idx_type> t (dims, 1);
  std::vector<unsigned int> my_cells = msh.cells ();
  std::size_t n = 0;

  dolfin::MeshFunction<std::size_t> cell_domains;
  if (! msh.domains ().is_empty ())
      if (msh.domains ().num_marked (D) != 0)
        cell_domains = * (msh.domains ().cell_domains ());

  for (octave_idx_type j = 0; j < t.cols (); ++j)
    {
      for (octave_idx_type i = 0; i < D + 1; ++i, ++n)
        t.xelem (i, j) += my_cells[n];

       if (! cell_domains.empty ())
         t.xelem (D + 1, j) = cell_domains[j];
    }

  octave_scalar_map a;
  a.setfield ("p", p);
  a.setfield ("e", e);
  a.setfield ("t", t);
  return a;

}


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;
      dolfin::Mesh msh;
      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 (D - 1);
      std::size_t num_side_edges = e.cols ();

      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)
                    {
                      facet.set_value ((*f).index (), e.xelem (4, i), msh);
                      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)
                    {
                      facet.set_value ((*f).index (), e.xelem (9, i), msh);
                      break;
                    }
                }
            }
        }

      *(msh.domains ().markers (D - 1)) = facet;

      // store information associated with t
      dolfin::MeshValueCollection<std::size_t> cell (D);
      std::size_t num_cells = t.cols ();

      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)
                    {
                      cell.set_value ((*f).index (), t.xelem (3, i), msh);
                      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)
                    {
                      cell.set_value ((*f).index (), t.xelem (4, i), msh);
                      break;
                    }
                }
            }
        }

      *(msh.domains ().markers (D)) = cell;

     boost::shared_ptr<const dolfin::Mesh> my_m (new dolfin::Mesh (msh)); 
     pmsh = my_m;
    }
}