/*
 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 "meshfunction.h"
#include "dolfin_compat.h"
#ifdef LATEST_DOLFIN
#include <dolfin/mesh/MeshPartitioning.h>
#endif

octave_value
compute_facet_markers (dolfin::Mesh const &, Array <octave_idx_type> const &,
                       std::size_t const);
octave_value
compute_cell_markers (dolfin::Mesh const &, Array <octave_idx_type> const &,
                      std::size_t const);

DEFUN_DLD (Mesh, args, nargout, "-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{mesh_out}, \
@var{facet_markers}, @var{cell_markers}]} = \
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 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\
Optional outputs @var{facet_markers} and @var{cell_markers} \
are MeshFunctions holding the markers extracted from the PDE-tool \
representation. \n\
The easiest way for dealing with meshes is using the msh pkg. \n\
@seealso{FunctionSpace, DirichletBC}\n\
@end deftypefn")
{
  int nargin = args.length ();
  octave_value_list retval;

  if (nargin < 1 || nargin > 1 || nargout > 3)
    { 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(0) = 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(0) = new mesh (p, e, t);

                  if (nargout >= 2)
                    {
                      if (! meshfunction_type_loaded)
                        {
                          meshfunction::register_type ();
                          meshfunction_type_loaded = true;
                          mlock ();
                        }

                      mesh const & msh_ov =
                        static_cast <mesh const &> (retval(0).get_rep ());
                      dolfin::Mesh const & msh = msh_ov.get_msh ();
                      std::size_t const D = p.rows ();

                      retval(1) = compute_facet_markers (msh, e, D);
                      if (nargout == 3)
                        { retval(2) = compute_cell_markers (msh, t, D); }
                    }
                }
            }

          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 size =
#ifdef LATEST_DOLFIN
        dolfin::MPI::size (MPI_COMM_WORLD);
#else
        dolfin::MPI::num_processes ();
#endif

      if (size == 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 (size == 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;
    }
}

octave_value
compute_facet_markers (dolfin::Mesh const & _msh,
                       Array <octave_idx_type> const & e,
                       std::size_t const D)
{
  dolfin::MeshFunction <std::size_t> facet (_msh, D - 1,
                                            std::numeric_limits <std::size_t>
                                            ::max ());
  std::size_t const num_side_edges = e.cols ();
  std::vector <std::size_t> const & global_vertices =
    _msh.topology ().global_indices (0);
  std::vector <std::size_t> const & global_facets =
    _msh.topology ().global_indices (D - 1);

  if (D == 2)
    {
      for (std::size_t i = 0; i < num_side_edges; ++i)
        {
          std::size_t local_vertex = 0;
          std::size_t const e_index = e.xelem (0, i) - 1;
          while (local_vertex < global_vertices.size () &&
                 e_index != global_vertices[local_vertex])
            { ++local_vertex; }

          if (local_vertex < global_vertices.size ())
            {
              dolfin::Vertex v (_msh, local_vertex);
              for (dolfin::FacetIterator f (v); ! f.end (); ++f)
                {
                  std::size_t const & vertex0 =
                    global_vertices[f->entities(0)[0]];
                  std::size_t const & vertex1 =
                    global_vertices[f->entities(0)[1]];

                  if (vertex0 == e.xelem (0, i) - 1
                      && vertex1 == e.xelem (1, i) - 1
                      || vertex0 == e.xelem (1, i) - 1
                      && vertex1 == e.xelem (0, i) - 1)
                    {
                      facet[*f] = e.xelem (4, i);
                      break;
                    }
                }
            }
        }
    }

  if (D == 3)
    {
      for (std::size_t i = 0; i < num_side_edges; ++i)
        {
          std::size_t local_vertex = 0;
          std::size_t const e_index = e.xelem (0, i) - 1;
          while (local_vertex < global_vertices.size () &&
                 e_index != global_vertices[local_vertex])
            { ++local_vertex; }

          if (local_vertex < global_vertices.size ())
            {
              dolfin::Vertex v (_msh, local_vertex);
              for (dolfin::FacetIterator f (v); ! f.end (); ++f)
                {
                  std::size_t const & vertex0 =
                    global_vertices[f->entities(0)[0]];
                  std::size_t const & vertex1 =
                    global_vertices[f->entities(0)[1]];
                  std::size_t const & vertex2 =
                    global_vertices[f->entities(0)[2]];

                  if (vertex0 == e.xelem (0, i) - 1
                      && vertex1 == e.xelem (1, i) - 1
                      && vertex2 == e.xelem (2, i) - 1
                      || vertex0 == e.xelem (0, i) - 1
                      && vertex1 == e.xelem (2, i) - 1
                      && vertex2 == e.xelem (1, i) - 1
                      || vertex0 == e.xelem (1, i) - 1
                      && vertex1 == e.xelem (0, i) - 1
                      && vertex2 == e.xelem (2, i) - 1
                      || vertex0 == e.xelem (1, i) - 1
                      && vertex1 == e.xelem (2, i) - 1
                      && vertex2 == e.xelem (0, i) - 1
                      || vertex0 == e.xelem (2, i) - 1
                      && vertex1 == e.xelem (0, i) - 1
                      && vertex2 == e.xelem (1, i) - 1
                      || vertex0 == e.xelem (2, i) - 1
                      && vertex1 == e.xelem (1, i) - 1
                      && vertex2 == e.xelem (0, i) - 1)
                    {
                      facet[*f] = e.xelem (9, i);
                      break;
                    }
                }
            }
        }
    }

  octave_value retval = new meshfunction ("ds", facet);
  return retval;
}

octave_value
compute_cell_markers (dolfin::Mesh const & _msh,
                      Array <octave_idx_type> const & t,
                      std::size_t const D)
{
  dolfin::MeshFunction <std::size_t> cell (_msh, D,
                                           std::numeric_limits <std::size_t>
                                           ::max ());
  std::size_t const num_cells = t.cols ();
  std::vector <std::size_t> const & global_vertices =
    _msh.topology ().global_indices (0);
  std::vector <std::size_t> const & global_cells =
    _msh.topology ().global_indices (D);

  if (D == 2)
    {
      for (std::size_t i = 0; i < num_cells; ++i)
        {
          std::size_t local_vertex = 0;
          std::size_t const t_index = t.xelem (0, i) - 1;
          while (local_vertex < global_vertices.size () &&
                 t_index != global_vertices[local_vertex])
            { ++local_vertex; }

          if (local_vertex < global_vertices.size ())
            {
              dolfin::Vertex v (_msh, local_vertex);
              for (dolfin::CellIterator f (v); ! f.end (); ++f)
                {
                  std::size_t const & vertex0 =
                    global_vertices[f->entities(0)[0]];
                  std::size_t const & vertex1 =
                    global_vertices[f->entities(0)[1]];
                  std::size_t const & vertex2 =
                    global_vertices[f->entities(0)[2]];

                  if (vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (0, i) - 1)
                    {
                      cell[*f] = t.xelem (3, i);
                      break;
                    }
                }
            }
        }
    }

  if (D == 3)
    {
      for (std::size_t i = 0; i < num_cells; ++i)
        {
          std::size_t local_vertex = 0;
          std::size_t const t_index = t.xelem (0, i) - 1;
          while (local_vertex < global_vertices.size () &&
                 t_index != global_vertices[local_vertex])
            { ++local_vertex; }

          if (local_vertex < global_vertices.size ())
            {
              dolfin::Vertex v (_msh, local_vertex);
              for (dolfin::CellIterator f (v); ! f.end (); ++f)
                {
                  std::size_t const & vertex0 =
                    global_vertices[f->entities(0)[0]];
                  std::size_t const & vertex1 =
                    global_vertices[f->entities(0)[1]];
                  std::size_t const & vertex2 =
                    global_vertices[f->entities(0)[2]];
                  std::size_t const & vertex3 =
                    global_vertices[f->entities(0)[3]];

                  if (vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (0, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (1, i) - 1

                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (0, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (1, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (0, i) - 1

                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (3, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (3, i) - 1
                      && vertex3 == t.xelem (0, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (2, i) - 1
                      && vertex1 == t.xelem (3, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (0, i) - 1

                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (0, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (2, i) - 1
                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (1, i) - 1
                      && vertex2 == t.xelem (2, i) - 1
                      && vertex3 == t.xelem (0, i) - 1
                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (0, i) - 1
                      && vertex3 == t.xelem (1, i) - 1
                      || vertex0 == t.xelem (3, i) - 1
                      && vertex1 == t.xelem (2, i) - 1
                      && vertex2 == t.xelem (1, i) - 1
                      && vertex3 == t.xelem (0, i) - 1)
                    {
                      cell[*f] = t.xelem (4, i);
                      break;
                    }
                }
            }
        }
    }

  octave_value retval = new meshfunction ("dx", cell);
  return retval;
}