# HG changeset patch
# User gedeone-octave <marco.vassallo@outlook.com>
# Date 1373040437 -7200
# Node ID ce9b06cc45c77e175d94abd5eb96a665ba52fc4c
# Parent  1974b68095fc1cd1d2d5dbe81f5ea267693e6033
New class for dealing with mesh

* mesh.h: header file: our class derive from dolfin::mesh and octave_base_value
* mesh.cc: we only need to implement the constructor

diff -r 1974b68095fc -r ce9b06cc45c7 src/mesh.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mesh.cc	Fri Jul 05 18:07:17 2013 +0200
@@ -0,0 +1,283 @@
+/*
+ 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>
+
+DEFINE_OCTAVE_ALLOCATOR (mesh);
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (mesh, "mesh", "mesh");
+
+const dolfin::Mesh & mesh::get_msh () const
+{ return msh; } 
+
+void mesh::set_msh (dolfin::Mesh & _msh)
+{ msh = _msh; }
+
+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;
+      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;
+    }
+}
diff -r 1974b68095fc -r ce9b06cc45c7 src/mesh.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mesh.h	Fri Jul 05 18:07:17 2013 +0200
@@ -0,0 +1,69 @@
+/*
+ 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/>.
+*/
+
+#ifndef _MESH_OCTAVE_
+#define _MESH_OCTAVE_
+
+#include <dolfin.h>
+#include <octave/oct.h>
+#include <octave/oct-map.h>
+
+class mesh : public octave_base_value, public dolfin::Mesh
+{
+
+ public:
+
+  // Constructor
+  mesh ()
+    : octave_base_value () {}
+
+  mesh (const dolfin::Mesh& _msh) : octave_base_value (), msh (_msh) {}
+
+  mesh (Array<double>& p, Array<octave_idx_type>& e, 
+                    Array<octave_idx_type>& t);
+
+
+  mesh (std::string _filename): octave_base_value (), msh (_filename) {} 
+
+  //a method for printing
+  void print (std::ostream& os, bool pr_as_read_syntax = false) const
+  { os << "msh : " << msh.label() << " with " << msh.num_vertices() 
+       << " vertices " << std::endl; }
+
+  //destructor
+  ~mesh(void) {  };
+
+
+  bool is_defined (void) const { return true; }
+
+
+  //get the information from the private member
+  const dolfin::Mesh  & get_msh (void) const;
+
+  // set the value of private members
+  void set_msh   (dolfin::Mesh & _msh);
+ 
+ private:
+
+  dolfin::Mesh msh;
+  DECLARE_OCTAVE_ALLOCATOR;
+  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA;
+
+};
+
+
+#endif