--- a/src/Mesh.cc	Wed Jul 30 21:09:52 2014 +0200
+++ b/src/Mesh.cc	Fri Aug 01 18:59:18 2014 +0200
@@ -16,13 +16,22 @@
 */
 
 #include "mesh.h"
+#include "meshfunction.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}]} = \
+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\
@@ -31,18 +40,21 @@
 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)\
+@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}\n\
+@seealso{FunctionSpace, DirichletBC}\n\
 @end deftypefn")
 {
   int nargin = args.length ();
-  octave_value retval = 0;
+  octave_value_list retval;
 
-  if (nargin < 1 || nargin > 1)
+  if (nargin < 1 || nargin > 1 || nargout > 3)
     print_usage ();
   else
     {
@@ -60,7 +72,7 @@
             {
               std::string filename = args(0).string_value ();
               //if the filename is not valid, dolfin takes care of it
-              retval = new mesh (filename);
+              retval(0) = new mesh (filename);
             }
 
           else if (args(0).is_map () == true)
@@ -72,7 +84,26 @@
 
               if (! error_state)
                 {
-                  retval = new mesh (p, e, t);
+                  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);
+                    }
                 }
             }
 
@@ -136,14 +167,14 @@
       dolfin::MeshValueCollection<std::size_t> facet(*msh, D - 1);
       std::size_t num_side_edges = e.cols ();
 
-      unsigned const num_procs =
+      unsigned const size =
 #ifdef LATEST_DOLFIN
         dolfin::MPI::size (MPI_COMM_WORLD);
 #else
         dolfin::MPI::num_processes ();
 #endif
 
-      if (num_procs == 1)
+      if (size == 1)
         {
           if (D == 2)
             {
@@ -202,11 +233,10 @@
             }
         }
 
-
       dolfin::MeshValueCollection<std::size_t> cell (*msh, D);
       std::size_t num_cells = t.cols ();
 
-      if (num_procs == 1)
+      if (size == 1)
         {
           if (D == 2)
             {
@@ -364,3 +394,298 @@
       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, 0);
+  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(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 (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, 0);
+  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 (cell);
+  return retval;
+}