# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1406912358 -7200
# Node ID ab35a8b0deefb4f83b88d0ac9f615ee4d09d84ec
# Parent  fb67b636616f88e4c48478cb8da5eb9c25967fae
Support storage of mesh markers via meshfunction

diff -r fb67b636616f -r ab35a8b0deef src/DirichletBC.cc
--- a/src/DirichletBC.cc	Wed Jul 30 21:09:52 2014 +0200
+++ b/src/DirichletBC.cc	Fri Aug 01 18:59:18 2014 +0200
@@ -18,12 +18,13 @@
 #include "boundarycondition.h"
 #include "functionspace.h"
 #include "expression.h"
+#include "meshfunction.h"
 #include "dolfin_compat.h"
 
 DEFUN_DLD (DirichletBC, args, ,
 "-*- texinfo -*-\n\
 @deftypefn {Function File} {[@var{bc}]} = \
-DirichletBC (@var{FunctionSpace}, @var{Boundary_Label},\
+DirichletBC (@var{FunctionSpace}, [@var{Boundaries},] @var{Boundary_Label},\
              @var{Function_handle}) \n\
 Specify essential boundary condition on a specific side.\n\
 The input parameters are\n\
@@ -36,6 +37,8 @@
 @var{Function handle} = ['''@'''(x, y) f1, '''@'''(x, y) f2, ...]\n\
 @item @var{Boundary_Label} is an Array which contains the label(s) of the \
 side(s) where the BC has to be applied.\n\
+@item @var{Boundaries} is an optional MeshFunction storing the markers \
+set on mesh facets. In parallel execution you must supply this argument.\n\
 @end itemize\n\
 The output @var{bc} is an object which contains the boundary conditions\n\
 @seealso{Mesh, FunctionSpace}\n\
@@ -44,7 +47,7 @@
   int nargin = args.length ();
   octave_value retval;
 
-  if (nargin < 3 || nargin > 3)
+  if (nargin < 3 || nargin > 4)
     print_usage ();
   else
     {
@@ -63,37 +66,81 @@
           mlock ();
         }
 
+      if (! meshfunction_type_loaded)
+        {
+          meshfunction::register_type ();
+          meshfunction_type_loaded = true;
+          mlock ();
+        }
+
       if (args(0).type_id () == functionspace::static_type_id ())
         {
           const functionspace & fspo = 
             static_cast<const functionspace&> (args(0).get_rep ());
           octave_fcn_handle * fh = args(1).fcn_handle_value ();
-          Array<octave_idx_type> side = args(2).array_value ();
 
-
-          if (!error_state)
+          if (nargin == 3)
             {
-              const SHARED_PTR <const dolfin::FunctionSpace>
-                & V (fspo.get_pfsp ());
+              Array<octave_idx_type> side = args(2).array_value ();
+
+              if (!error_state)
+                {
+                  const SHARED_PTR <const dolfin::FunctionSpace>
+                    & V (fspo.get_pfsp ());
 
-              octave_value_list b (3, 1);
-              octave_value_list tmp = feval (fh->function_value (), b);
-              Array<double> res = tmp(0).array_value ();
-              std::size_t l = res.length ();
+                  octave_value_list b (3, 1);
+                  octave_value_list tmp = feval (fh->function_value (), b);
+                  Array<double> res = tmp(0).array_value ();
+                  std::size_t l = res.length ();
+
+                  expression * pf;
+                  if (l > 1)
+                    pf = new expression (*fh, l);
+                  else
+                    pf = new expression (*fh);
+
+                  SHARED_PTR <const expression> f (pf);
+                  boundarycondition * pbc = new boundarycondition ();
 
-              expression * pf;
-              if (l > 1)
-                pf = new expression (*fh, l);
-              else
-                pf = new expression (*fh);
+                  for (octave_idx_type i = 0;
+                       i < side.length (); ++i)
+                      pbc->add_bc (V, f, side(i));
+                  retval = octave_value (pbc);
+                }
+            }
+          else
+            {
+              Array<octave_idx_type> side = args(3).array_value ();
+              meshfunction const & mf =
+                static_cast <meshfunction const &> (args(2).get_rep ());
+
+              if (! error_state)
+                {
+                  const SHARED_PTR <const dolfin::FunctionSpace>
+                    & V (fspo.get_pfsp ());
 
-              SHARED_PTR <const expression> f (pf);
-              boundarycondition * pbc = new boundarycondition ();
+                  octave_value_list b (3, 1);
+                  octave_value_list tmp = feval (fh->function_value (), b);
+                  Array<double> res = tmp(0).array_value ();
+                  std::size_t l = res.length ();
+
+                  expression * pf;
+                  if (l > 1)
+                    pf = new expression (*fh, l);
+                  else
+                    pf = new expression (*fh);
 
-              for (octave_idx_type i = 0;
-                   i < side.length (); ++i)
-                  pbc->add_bc (V, f, side(i));
-              retval = octave_value (pbc);
+                  SHARED_PTR <const expression> f (pf);
+                  boundarycondition * pbc = new boundarycondition ();
+
+                  SHARED_PTR <dolfin::MeshFunction <std::size_t> const> const &
+                    subdomains = mf.get_pmf ();
+
+                  for (octave_idx_type i = 0;
+                       i < side.length (); ++i)
+                      pbc->add_bc (V, f, subdomains, side(i));
+                  retval = octave_value (pbc);
+                }
             }
         }
     }
diff -r fb67b636616f -r ab35a8b0deef src/Makefile.in
--- a/src/Makefile.in	Wed Jul 30 21:09:52 2014 +0200
+++ b/src/Makefile.in	Fri Aug 01 18:59:18 2014 +0200
@@ -56,7 +56,7 @@
 Mesh.oct: Mesh.o
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) Mesh.o -o $@ $(LIBS)
 
-Mesh.o: Mesh.cc mesh.h dolfin_compat.h
+Mesh.o: Mesh.cc mesh.h dolfin_compat.h meshfunction.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c Mesh.cc -o $@
 
 fem_get_mesh.oct: fem_get_mesh.o
@@ -68,7 +68,7 @@
 DirichletBC.oct: DirichletBC.o
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) DirichletBC.o -o $@ $(LIBS)
 
-DirichletBC.o: DirichletBC.cc functionspace.h boundarycondition.h expression.h dolfin_compat.h
+DirichletBC.o: DirichletBC.cc functionspace.h boundarycondition.h expression.h dolfin_compat.h meshfunction.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c DirichletBC.cc -o $@
 
 Expression.oct: Expression.o
diff -r fb67b636616f -r ab35a8b0deef src/Mesh.cc
--- 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;
+}
diff -r fb67b636616f -r ab35a8b0deef src/boundarycondition.h
--- a/src/boundarycondition.h	Wed Jul 30 21:09:52 2014 +0200
+++ b/src/boundarycondition.h	Fri Aug 01 18:59:18 2014 +0200
@@ -60,6 +60,17 @@
       bcu.push_back(p);
     }
 
+  void
+  add_bc (SHARED_PTR <dolfin::FunctionSpace const> const & V,
+          SHARED_PTR <dolfin::GenericFunction const> f,
+          SHARED_PTR <dolfin::MeshFunction <std::size_t> const> const & sd,
+          std::size_t n)
+    {
+      SHARED_PTR <dolfin::DirichletBC const> 
+        p (new dolfin::DirichletBC (V, f, sd, n));
+      bcu.push_back(p);
+    }
+
  private:
 
   std::vector<SHARED_PTR 
diff -r fb67b636616f -r ab35a8b0deef src/meshfunction.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/meshfunction.h	Fri Aug 01 18:59:18 2014 +0200
@@ -0,0 +1,60 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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/>.
+*/
+
+#ifndef _MESHFUNCTION_OCTAVE_
+#define _MESHFUNCTION_OCTAVE_
+
+#include <dolfin.h>
+#include <octave/oct.h>
+#include "dolfin_compat.h"
+
+class meshfunction : public octave_base_value
+{
+  public:
+
+  meshfunction (void)
+    : octave_base_value () {}
+
+  meshfunction (dolfin::MeshFunction <std::size_t> const & _mf)
+    : octave_base_value (), pmf (new dolfin::MeshFunction <std::size_t> (_mf)) {}
+
+  bool
+  is_defined (void) const
+    { return true; }
+
+  dolfin::MeshFunction <std::size_t> const &
+  get_mf (void) const
+    { return *pmf; }
+
+  SHARED_PTR <dolfin::MeshFunction <std::size_t> const> const &
+  get_pmf (void) const
+    { return pmf; }
+
+  private:
+
+  SHARED_PTR <dolfin::MeshFunction <std::size_t> const> pmf;
+
+  DECLARE_OCTAVE_ALLOCATOR;
+  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA;
+};
+
+static bool meshfunction_type_loaded = false;
+
+DEFINE_OCTAVE_ALLOCATOR (meshfunction);
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (meshfunction, "meshfunction", "meshfunction");
+
+#endif