Mercurial > fem-fenics-eugenio
changeset 79:3e49ef16d74a
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.
| author | gedeone-octave <marcovass89@hotmail.it> |
|---|---|
| date | Sat, 03 Aug 2013 14:32:33 +0200 |
| parents | 670a5d91c397 |
| children | 16ccfaf1476a |
| files | src/fem_get_mesh.cc src/fem_init_mesh.cc src/mesh.cc |
| diffstat | 3 files changed, 339 insertions(+), 359 deletions(-) [+] |
line wrap: on
line diff
--- 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; + +}
--- 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; + } +}
--- 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; - } -}