Mercurial > fem-fenics-eugenio
changeset 12:b7c74c0bdabd
New functions which get a (p,e,t) matrix from a mesh object
* Makefile.in: added the option to compile also this function
* fem_init_mesh.cc: dld function which call the mesh::get_pet method
* mesh.cc mesh.h: added the new method get_pet which return a mesh as
an octave_sclar_map. (code similar to mshm_dolfin_read)
| author | gedeone-octave <marco.vassallo@outlook.com> |
|---|---|
| date | Sat, 06 Jul 2013 00:39:34 +0200 |
| parents | 6f11638bfbac |
| children | 36882e185f02 |
| files | src/Makefile.in src/fem_init_mesh.cc src/mesh.cc src/mesh.h |
| diffstat | 4 files changed, 127 insertions(+), 48 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Makefile.in Sat Jul 06 00:31:00 2013 +0200 +++ b/src/Makefile.in Sat Jul 06 00:39:34 2013 +0200 @@ -1,6 +1,6 @@ MKOCTFILE ?= mkoctfile -OCTFILES=fem_init_env.oct fem_init_mesh.oct +OCTFILES=fem_init_env.oct fem_init_mesh.oct fem_get_mesh.oct LIBS = HAVE_DOLFIN_H = @HAVE_DOLFIN_H@ @@ -17,12 +17,18 @@ fem_init_mesh.oct: mesh.o fem_init_mesh.o $(MKOCTFILE) $(CPPFLAGS) -s mesh.o fem_init_mesh.o -o $@ $(LDFLAGS) $(LIBS) +fem_get_mesh.oct: mesh.o fem_get_mesh.o + $(MKOCTFILE) $(CPPFLAGS) -s mesh.o fem_get_mesh.o -o $@ $(LDFLAGS) $(LIBS) + mesh.o: mesh.cc mesh.h $(MKOCTFILE) $(CPPFLAGS) -c mesh.cc $(LDFLAGS) -o $@ -I. fem_init_mesh.o: fem_init_mesh.cc mesh.cc mesh.h $(MKOCTFILE) $(CPPFLAGS) -c fem_init_mesh.cc $(LDFLAGS) -o $@ -I. +fem_get_mesh.o: fem_get_mesh.cc mesh.cc mesh.h + $(MKOCTFILE) $(CPPFLAGS) -c fem_get_mesh.cc $(LDFLAGS) -o $@ -I. + fem_init_env.o: fem_init_env.cc mesh.cc mesh.h $(MKOCTFILE) $(CPPFLAGS) -c fem_init_env.cc $(LDFLAGS) -o $@ -I.
--- a/src/fem_init_mesh.cc Sat Jul 06 00:31:00 2013 +0200 +++ b/src/fem_init_mesh.cc Sat Jul 06 00:39:34 2013 +0200 @@ -17,7 +17,6 @@ #include "mesh.h" - DEFUN_DLD (fem_init_mesh, args, , "initialize a mesh from (p, e, t) or file") { int nargin = args.length (); @@ -27,34 +26,28 @@ print_usage (); else { - //octave_scalar_map a = args(1).scalar_map_value (); - //bool mesh_tl = a.contents ("mesh_tl").bool_value (); - bool mesh_tl = true; - if (!error_state) { - if (mesh_tl == false) - error("Mesh type not loaded, you should proably use fem_init()"); - else + if (args(0).is_string () == true) { - if (args(0).is_string () == true) - { - std::string filename = args(0).string_value(); - retval = new mesh (filename); - } + std::string filename = args(0).string_value (); + //if the filename is not valid, dolfin takes care of it + retval = new mesh (filename); + } + else + { + octave_scalar_map a = args(0).scalar_map_value (); + Array<double> p = a.contents ("p").matrix_value (); + Array<octave_idx_type> t = a.contents ("t").matrix_value (); + Array<octave_idx_type> e = a.contents ("e").matrix_value (); + + if (! error_state) + { + retval = new mesh (p, e, t); + } else - { - octave_scalar_map a = args(0).scalar_map_value (); - Array<double> p = a.contents ("p").matrix_value (); - Array<octave_idx_type> t = a.contents ("t").matrix_value (); - Array<octave_idx_type> e = a.contents ("e").matrix_value (); - - if (! error_state) - { - retval = new mesh (p, e, t); - } - } + error ("fem_init_mesh: the argument you provide is invalid"); } } }
--- a/src/mesh.cc Sat Jul 06 00:31:00 2013 +0200 +++ b/src/mesh.cc Sat Jul 06 00:39:34 2013 +0200 @@ -21,14 +21,100 @@ DEFINE_OCTAVE_ALLOCATOR (mesh); DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (mesh, "mesh", "mesh"); -const dolfin::Mesh & mesh::get_msh () const -{ return msh; } +const dolfin::Mesh & +mesh::get_msh () const +{ + return msh; +} + +void +mesh::set_msh (dolfin::Mesh & _msh) +{ + msh = _msh; +} + +octave_scalar_map +mesh::get_pet () const +{ + //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 ()); -void mesh::set_msh (dolfin::Mesh & _msh) -{ msh = _msh; } + 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; -mesh::mesh (Array<double>& p, Array<octave_idx_type>& e, - Array<octave_idx_type>& t) + // 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 ();
--- a/src/mesh.h Sat Jul 06 00:31:00 2013 +0200 +++ b/src/mesh.h Sat Jul 06 00:39:34 2013 +0200 @@ -27,35 +27,31 @@ public: - // Constructor - mesh () - : octave_base_value () {} + // Constructors + 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 (Array<double>& p, Array<octave_idx_type>& e, Array<octave_idx_type>& t); - - mesh (std::string _filename): octave_base_value (), msh (_filename) {} + 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() + { os << "msh : " << msh.label() << " with " << msh.num_vertices() << " vertices " << std::endl; } //destructor - ~mesh(void) { }; - + ~mesh(void) { }; bool is_defined (void) const { return true; } + //get the information from the private member + const dolfin::Mesh & get_msh (void) const; + octave_scalar_map get_pet (void) const; - //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); + // set the value of the private member + void set_msh (dolfin::Mesh & _msh); private: @@ -64,6 +60,4 @@ DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA; }; - - #endif