# HG changeset patch
# User gedeone-octave <marcovass89@hotmail.it>
# Date 1375533153 -7200
# Node ID 3e49ef16d74a9617f222af94c3f307f166bff486
# Parent  670a5d91c397204e6cc3263e6ecf2e93ba3a24f1
The methods of the mesh class are implemented in the file where they are needed.

* fem_get_mesh.cc: implement the method for extracting the matrix.
* fem_init_mesh.cc: implement the constructor
* This is done to avoid conflict of multiple definition functions.

diff -r 670a5d91c397 -r 3e49ef16d74a src/fem_get_mesh.cc
--- a/src/fem_get_mesh.cc	Sat Aug 03 14:29:04 2013 +0200
+++ b/src/fem_get_mesh.cc	Sat Aug 03 14:32:33 2013 +0200
@@ -53,3 +53,83 @@
     }
   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;
+
+}
diff -r 670a5d91c397 -r 3e49ef16d74a src/fem_init_mesh.cc
--- a/src/fem_init_mesh.cc	Sat Aug 03 14:29:04 2013 +0200
+++ b/src/fem_init_mesh.cc	Sat Aug 03 14:32:33 2013 +0200
@@ -80,3 +80,262 @@
     }
   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;
+      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;
+    }
+}
diff -r 670a5d91c397 -r 3e49ef16d74a src/mesh.cc
--- a/src/mesh.cc	Sat Aug 03 14:29:04 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,359 +0,0 @@
-/*
- 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;
-    }
-}