/*
 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"

DEFUN_DLD (fem_get_mesh, args, ,"-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{mesh}]} = \
fem_get_mesh (@var{fem_mesh}) \n\
Return a (p, e, t) representation of @var{fem_mesh}\n\
The @var{mesh_to_read} should be an object created with \
fem_init_mesh().\n\
The output @var{mesh} is a PDE-tool like structure\n\
with matrix fields (p,e,t).\n\
@seealso{fem_init_mesh}\n\
@end deftypefn")
{
  int nargin = args.length ();
  octave_value retval;

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

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

      if (args(0).type_id () == mesh::static_type_id ())
        {
          const mesh & msh = static_cast<const mesh&> (args(0).get_rep ());
          retval = octave_value (msh.get_pet ());
        }
      else
        error ("fem_get_mesh: the argument is not a mesh type");
    }
  return retval;
}

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;

}