Mercurial > fem-fenics-eugenio
changeset 267:53039ac90368
Mark meshfunction using subdomain
* src/SubDomain.cc: return subdomain to use for marking
* src/MeshFunction.cc: return a user defined meshfunction
* src/mark.cc: with the information from subdomain, mark a meshfunction
author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
---|---|
date | Tue, 12 Aug 2014 15:42:50 +0200 |
parents | 6b37560b7cbb |
children | 61830a4f9ab9 |
files | INDEX src/Makefile.in src/MeshFunction.cc src/SubDomain.cc src/mark.cc src/meshfunction.h src/subdomain.h src/subdomain_rep.h |
diffstat | 8 files changed, 318 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/INDEX Fri Aug 08 09:52:59 2014 +0200 +++ b/INDEX Tue Aug 12 15:42:50 2014 +0200 @@ -10,6 +10,8 @@ Mesh FunctionSpace SubSpace + SubDomain + mark Problem variables Constant Expression
--- a/src/Makefile.in Fri Aug 08 09:52:59 2014 +0200 +++ b/src/Makefile.in Tue Aug 12 15:42:50 2014 +0200 @@ -45,7 +45,9 @@ interpolate.oct \ is_master_node.oct \ barrier.oct \ - MeshFunction.oct + MeshFunction.oct \ + mark.oct \ + SubDomain.oct all: $(OCTFILES) @@ -191,6 +193,18 @@ barrier.o: barrier.cc dolfin_compat.h CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -o $@ -c $< +SubDomain.oct: SubDomain.o + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) $< -o $@ $(LIBS) + +SubDomain.o: SubDomain.cc subdomain.h + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@ + +mark.oct: mark.o + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) $< -o $@ $(LIBS) + +mark.o: mark.cc subdomain.h meshfunction.h dolfin_compat.h + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@ + clean: $(RM) *.o core octave-core *.oct *~ *.xml *.a
--- a/src/MeshFunction.cc Fri Aug 08 09:52:59 2014 +0200 +++ b/src/MeshFunction.cc Tue Aug 12 15:42:50 2014 +0200 @@ -25,6 +25,8 @@ MeshFunction (@var{Mesh}, @var{filename})\n\ @deftypefnx {Function File} {[@var{meshfunc}]} = \ MeshFunction (@var{type}, @var{Mesh}, @var{filename})\n\ +@deftypefnx {Function File} {[@var{meshfunc}]} = \ +MeshFunction (@var{type}, @var{Mesh}, @var{dim}, @var{value}) \n\ Initialize a meshfunction with the values contained in a file.\n\ The output @var{meshfunc} is an object which contains a representation of the \ meshfunction in @var{filename} which can be used to mark subdomains or \ @@ -33,14 +35,17 @@ in order to use @var{meshfunc} to mark subdomains in BilinearForm or \ LinearForm. Three labels are recognized: ""dx"" for cell markers, \ ""ds"" for exterior facets and ""dS"" for interior facets.\n\ -@seealso{DirichletBC, save, BilinearForm, LinearForm}\n\ +With the third invocation listed above it is possible to initialize \ +a meshfunction of topological dimension @var{dim} on @var{Mesh}, setting all \ +of its values to @var{value}, which is an integer.\n\ +@seealso{DirichletBC, save, BilinearForm, LinearForm, mark}\n\ @end deftypefn") { int nargin = args.length (); octave_value retval; - if (nargin < 2 || nargin > 3 || nargout > 1) + if (nargin < 2 || nargin > 4 || nargout > 1) print_usage (); else { @@ -67,23 +72,41 @@ offset = 1; } - if (args(0 + offset).type_id () == mesh::static_type_id () - && args(1 + offset).is_string ()) + if (args(0 + offset).type_id () == mesh::static_type_id ()) { - std::string filename = args(1 + offset).string_value (); mesh const & msh_arg = static_cast<mesh const &> (args(0 + offset).get_rep ()); - if (!error_state) + if (args(1 + offset).is_string ()) { - SHARED_PTR <dolfin::Mesh const> const & - pmsh = msh_arg.get_pmsh (); + std::string filename = args(1 + offset).string_value (); + if (!error_state) + { + SHARED_PTR <dolfin::Mesh const> const & + pmsh = msh_arg.get_pmsh (); - try - { retval = new meshfunction (label, pmsh, filename); } - catch (std::runtime_error &) - { error ("error reading file"); } + try + { retval = new meshfunction (label, pmsh, filename); } + catch (std::runtime_error &) + { error ("error reading file"); } + } } + else if (args(1 + offset).is_real_scalar () && + args(2 + offset).is_real_scalar ()) + { + std::size_t dim = args(1 + offset).ulong_value (); + std::size_t value = args(2 + offset).ulong_value (); + + if (!error_state) + { + SHARED_PTR <dolfin::Mesh const> const & + pmsh = msh_arg.get_pmsh (); + + retval = new meshfunction (label, pmsh, dim, value); + } + } + else + error ("invalid input arguments"); } else error ("invalid input arguments");
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SubDomain.cc Tue Aug 12 15:42:50 2014 +0200 @@ -0,0 +1,56 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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 3 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 <dolfin.h> +#include "subdomain.h" + +DEFUN_DLD (SubDomain, args, nargout, +"-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{subdomain}]} = \ +SubDomain (@var{inside}, @var{boundary})\n\ +Initialize a subdomain with a function handle, @var{inside}, and a binary \ +flag. When the latter is true, @var{subdomain} will only contain boundary \ +entities. @var{subdomain} can be used to mark a meshfunction. \n\ +@seealso{mark, MeshFunction}\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + if (nargin < 2 || nargin > 2 || nargout > 1) + print_usage (); + else + { + if (args(0).is_function_handle () && args(1).is_bool_scalar ()) + { + if (! subdomain_type_loaded) + { + subdomain::register_type (); + subdomain_type_loaded = true; + mlock (); + } + + octave_fcn_handle * pfh = args(0).fcn_handle_value (); + bool on_boundary = args(1).bool_value (); + retval = new subdomain (*pfh, on_boundary); + } + else + error ("invalid input arguments"); + } + + return retval; +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mark.cc Tue Aug 12 15:42:50 2014 +0200 @@ -0,0 +1,77 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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 3 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 <dolfin.h> +#include "meshfunction.h" +#include "subdomain.h" +#include "dolfin_compat.h" + +DEFUN_DLD (mark, args, nargout, +"-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{markers}]} = \ +mark (@var{subdomain}, @var{meshfunction}, @var{number})\n\ +Mark @var{meshfunction} with @var{number} on the entities contained in \ +@var{subdomain}. The output, @var{markers}, is a copy of @var{meshfunction} \ +where entities inside @var{subdomain} are marked.\n\ +@seealso{MeshFunction, SubDomain}\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + if (nargin < 3 || nargin > 3 ||nargout > 1) + print_usage (); + else + { + if (! meshfunction_type_loaded) + { + meshfunction::register_type (); + meshfunction_type_loaded = true; + mlock (); + } + if (! subdomain_type_loaded) + { + subdomain::register_type (); + subdomain_type_loaded = true; + mlock (); + } + + if (args(0).type_id () == subdomain::static_type_id () && + args(1).type_id () == meshfunction::static_type_id () && + args(2).is_real_scalar ()) + { + std::size_t number = args(2).ulong_value (); + subdomain const & sd = static_cast <subdomain const &> + (args(0).get_rep ()); + meshfunction const & mf = static_cast <meshfunction const &> + (args(1).get_rep ()); + + if (! error_state) + { + dolfin::MeshFunction <std::size_t> out = mf.get_mf (); + SHARED_PTR <subdomain_rep const> const & sdrep = sd.get_psd (); + + sdrep->mark (out, number); + retval = new meshfunction (mf.get_str (), out); + } + } + else + error ("invalid input arguments"); + } + + return retval; +}
--- a/src/meshfunction.h Fri Aug 08 09:52:59 2014 +0200 +++ b/src/meshfunction.h Tue Aug 12 15:42:50 2014 +0200 @@ -41,6 +41,13 @@ pmf (new dolfin::MeshFunction <std::size_t> (mesh, filename)), str (label) {} + meshfunction (std::string const & label, + SHARED_PTR <dolfin::Mesh const> mesh, + std::size_t dim, std::size_t value) + : octave_base_value (), + pmf (new dolfin::MeshFunction <std::size_t> (mesh, dim, value)), + str (label) {} + bool is_defined (void) const { return true; }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/subdomain.h Tue Aug 12 15:42:50 2014 +0200 @@ -0,0 +1,65 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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 3 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 _SUBDOMAIN_OCTAVE_ +#define _SUBDOMAIN_OCTAVE_ + +#include "subdomain_rep.h" +#include "dolfin_compat.h" + +class subdomain : public octave_base_value +{ + public: + + subdomain () + : octave_base_value () {} + + subdomain (octave_fcn_handle const & fh, bool on_boundary) + : octave_base_value (), rep (new subdomain_rep (fh, on_boundary)) {} + + ~subdomain (void) {} + + void + print (std::ostream& os, bool pr_as_read_syntax = false) const + { + os << "SubDomain" << std::endl; + } + + bool + is_defined (void) const + { return true; } + + SHARED_PTR <subdomain_rep const> const & + get_psd (void) const + { return rep; } + + private: + + SHARED_PTR <subdomain_rep const> rep; + + DECLARE_OCTAVE_ALLOCATOR; + DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA; +}; + +static bool subdomain_type_loaded = false; + +DEFINE_OCTAVE_ALLOCATOR (subdomain); +DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (subdomain, + "subdomain", + "subdomain"); + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/subdomain_rep.h Tue Aug 12 15:42:50 2014 +0200 @@ -0,0 +1,61 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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 3 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 _SUBDOMAIN_REP_OCTAVE_ +#define _SUBDOMAIN_REP_OCTAVE_ + +#include <dolfin.h> +#include <octave/oct.h> +#include <octave/oct-map.h> +#include <octave/ov-fcn-handle.h> +#include <octave/ov-fcn.h> +#include <octave/parse.h> +#include <octave/octave.h> + +class subdomain_rep : public dolfin::SubDomain +{ + public: + + subdomain_rep (void) + : dolfin::SubDomain (), pfh (NULL) {} + + subdomain_rep (octave_fcn_handle const & fh, bool on_boundary) + : dolfin::SubDomain (), pfh (new octave_fcn_handle (fh)), + check_boundary (on_boundary) {} + + ~subdomain_rep (void) { delete pfh; } + + bool + inside (dolfin::Array<double> const & x, bool on_boundary) const + { + octave_value_list b; + b.resize (x.size ()); + for (std::size_t i = 0; i < x.size (); ++i) + b(i) = x[i]; + octave_value_list tmp = feval (pfh->function_value (), b); + bool retval = tmp(0).bool_value (); + if (check_boundary) + retval = retval && on_boundary; + return retval; + } + + private: + + octave_fcn_handle * pfh; + bool check_boundary; +}; +#endif