Mercurial > fem-fenics-eugenio
changeset 3:ce9b06cc45c7
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
| author | gedeone-octave <marco.vassallo@outlook.com> |
|---|---|
| date | Fri, 05 Jul 2013 18:07:17 +0200 |
| parents | 1974b68095fc |
| children | 00bfa5dd0dd6 |
| files | src/mesh.cc src/mesh.h |
| diffstat | 2 files changed, 352 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /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; + } +}
--- /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