# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1407850970 -7200
# Node ID 53039ac903687f1904c170ba3593ffc98a6d08f0
# Parent  6b37560b7cbbfb1e0713093a182f110033745ecb
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

diff -r 6b37560b7cbb -r 53039ac90368 INDEX
--- 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
diff -r 6b37560b7cbb -r 53039ac90368 src/Makefile.in
--- 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
 
diff -r 6b37560b7cbb -r 53039ac90368 src/MeshFunction.cc
--- 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");
diff -r 6b37560b7cbb -r 53039ac90368 src/SubDomain.cc
--- /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
diff -r 6b37560b7cbb -r 53039ac90368 src/mark.cc
--- /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;
+}
diff -r 6b37560b7cbb -r 53039ac90368 src/meshfunction.h
--- 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; }
diff -r 6b37560b7cbb -r 53039ac90368 src/subdomain.h
--- /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
diff -r 6b37560b7cbb -r 53039ac90368 src/subdomain_rep.h
--- /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