# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1407432096 -7200
# Node ID a61fc34334cac05740dd69568610061ee17f0563
# Parent  68cae2998775700ace6da041406df28c2707c860
Use meshfunction to mark subdomains

diff -r 68cae2998775 -r a61fc34334ca inst/BilinearForm.m
--- a/inst/BilinearForm.m	Thu Aug 07 14:56:35 2014 +0200
+++ b/inst/BilinearForm.m	Thu Aug 07 19:21:36 2014 +0200
@@ -29,9 +29,10 @@
 ## The optional arguments are the @var{coefficient_1}, @var{coefficient_2} 
 ## which specify the parameters for the BilinearForm as stated in the ufl file.
 ## They can be either a Constant, a Function or an Expression.
+## You can also optionally provide a MeshFunction marking subdomains.
 ##
 ## @seealso{import_ufl_BilinearForm, import_ufl_Problem, FunctionSpace, 
-## LinearForm, ResidualForm}
+## LinearForm, ResidualForm, MeshFunction}
 ## @end deftypefn
 function a = BilinearForm (name, U, V, varargin)
 
diff -r 68cae2998775 -r a61fc34334ca inst/LinearForm.m
--- a/inst/LinearForm.m	Thu Aug 07 14:56:35 2014 +0200
+++ b/inst/LinearForm.m	Thu Aug 07 19:21:36 2014 +0200
@@ -30,9 +30,10 @@
 ## which specify the parameters for the LinearForm with the same name which
 ## was used in the ufl file.
 ## They can be either a Constant, a Function or an Expression.
+## You can also optionally provide a MeshFunction marking subdomains.
 ##
 ## @seealso{import_ufl_LinearForm, import_ufl_Problem, BilinearForm, 
-## ResidualForm, BilinearForm}
+## ResidualForm, MeshFunction}
 ## @end deftypefn
 
 function a = LinearForm (name, V, varargin)
diff -r 68cae2998775 -r a61fc34334ca inst/private/generate_lhs.m
--- a/inst/private/generate_lhs.m	Thu Aug 07 14:56:35 2014 +0200
+++ b/inst/private/generate_lhs.m	Thu Aug 07 19:21:36 2014 +0200
@@ -25,6 +25,7 @@
 #include <fem-fenics/coefficient.h>\n\
 #include <fem-fenics/function.h>\n\
 #include <fem-fenics/functionspace.h>\n\
+#include <fem-fenics/meshfunction.h>\n\
 #include <fem-fenics/dolfin_compat.h>\n\
 \n\
 DEFUN_DLD (@@UFL_NAME@@_LinearForm, args, , "" b = fem_lhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF)"")\n\
@@ -75,6 +76,13 @@
                   mlock ();\n\
                 }\n\
 \n\
+              if (! meshfunction_type_loaded)\n\
+                {\n\
+                  meshfunction::register_type ();\n\
+                  meshfunction_type_loaded = true;\n\
+                  mlock ();\n\
+                }\n\
+\n\
               for (std::size_t i = 1; i < nargin; ++i)\n\
                 {\n\
                   if (args(i).type_id () == coefficient::static_type_id ())\n\
@@ -98,6 +106,22 @@
                       L.set_coefficient (n, pfun);\n\
                       ++nc;\n\
                     }\n\
+\n\
+                  if (args(i).type_id () == meshfunction::static_type_id ())\n\
+                    {\n\
+                      meshfunction const & mf_arg\n\
+                        = static_cast <meshfunction const &> (args(i).get_rep ());\n\
+\n\
+                      std::string type = mf_arg.get_str ();\n\
+                      dolfin::MeshFunction <std::size_t> const &\n\
+                        mf = mf_arg.get_mf ();\n\
+                      if (type == ""dx"")\n\
+                        L.dx = mf;\n\
+                      else if (type == ""ds"")\n\
+                        L.ds = mf;\n\
+                      else if (type == ""dS"")\n\
+                        L.dS = mf;\n\
+                    }\n\
                  }\n\
 \n\
               if (nc != ncoef)\n\
diff -r 68cae2998775 -r a61fc34334ca inst/private/generate_rhs.m
--- a/inst/private/generate_rhs.m	Thu Aug 07 14:56:35 2014 +0200
+++ b/inst/private/generate_rhs.m	Thu Aug 07 19:21:36 2014 +0200
@@ -25,6 +25,7 @@
 #include <fem-fenics/coefficient.h>\n\
 #include <fem-fenics/function.h>\n\
 #include <fem-fenics/functionspace.h>\n\
+#include <fem-fenics/meshfunction.h>\n\
 #include <fem-fenics/dolfin_compat.h>\n\
 \n\
 DEFUN_DLD (@@UFL_NAME@@_BilinearForm, args, , ""A = fem_rhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF)"")\n\
@@ -79,6 +80,13 @@
                   mlock ();\n\
                 }\n\
 \n\
+              if (! meshfunction_type_loaded)\n\
+                {\n\
+                  meshfunction::register_type ();\n\
+                  meshfunction_type_loaded = true;\n\
+                  mlock ();\n\
+                }\n\
+\n\
               for (std::size_t i = 2; i < nargin; ++i)\n\
                 {\n\
                   if (args(i).type_id () == coefficient::static_type_id ())\n\
@@ -102,6 +110,22 @@
                       a.set_coefficient (n, pfun);\n\
                       ++nc;\n\
                     }\n\
+\n\
+                  if (args(i).type_id () == meshfunction::static_type_id ())\n\
+                    {\n\
+                      meshfunction const & mf_arg\n\
+                        = static_cast <meshfunction const &> (args(i).get_rep ());\n\
+\n\
+                      std::string type = mf_arg.get_str ();\n\
+                      dolfin::MeshFunction <std::size_t> const &\n\
+                        mf = mf_arg.get_mf ();\n\
+                      if (type == ""dx"")\n\
+                        a.dx = mf;\n\
+                      else if (type == ""ds"")\n\
+                        a.ds = mf;\n\
+                      else if (type == ""dS"")\n\
+                        a.dS = mf;\n\
+                    }\n\
                  }\n\
 \n\
               if (nc != ncoef)\n\
diff -r 68cae2998775 -r a61fc34334ca src/Mesh.cc
--- a/src/Mesh.cc	Thu Aug 07 14:56:35 2014 +0200
+++ b/src/Mesh.cc	Thu Aug 07 19:21:36 2014 +0200
@@ -491,7 +491,7 @@
         }
     }
 
-  octave_value retval = new meshfunction (facet);
+  octave_value retval = new meshfunction ("ds", facet);
   return retval;
 }
 
@@ -690,6 +690,6 @@
         }
     }
 
-  octave_value retval = new meshfunction (cell);
+  octave_value retval = new meshfunction ("dx", cell);
   return retval;
 }
diff -r 68cae2998775 -r a61fc34334ca src/MeshFunction.cc
--- a/src/MeshFunction.cc	Thu Aug 07 14:56:35 2014 +0200
+++ b/src/MeshFunction.cc	Thu Aug 07 19:21:36 2014 +0200
@@ -23,18 +23,24 @@
 DEFUN_DLD (MeshFunction, args, nargout, "-*- texinfo -*-\n\
 @deftypefn {Function File} {[@var{meshfunc}]} = \
 MeshFunction (@var{Mesh}, @var{filename})\n\
+@deftypefnx {Function File} {[@var{meshfunc}]} = \
+MeshFunction (@var{type}, @var{Mesh}, @var{filename})\n\
 Initialize a meshfunction with the values contained in a XDMF 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 \
 facets where Dirichlet boundary conditions are to be applied.\n\
-@seealso{DirichletBC, save}\n\
+Optionally, you can set a label identifying the @var{type} of markers, \
+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\
 @end deftypefn")
 {
 
   int nargin = args.length ();
   octave_value retval;
   
-  if (nargin < 2 || nargin > 2 || nargout > 1)
+  if (nargin < 2 || nargin > 3 || nargout > 1)
     print_usage ();
   else
     {
@@ -52,12 +58,21 @@
           mlock ();
         }
 
-      if (args(0).type_id () == mesh::static_type_id ()
-          && args(1).is_string ())
+      octave_idx_type offset = 0;
+      std::string label;
+
+      if (args(0).is_string ())
         {
-          std::string filename = args(1).string_value ();
+          label = args(0).string_value ();
+          offset = 1;
+        }
+
+      if (args(0 + offset).type_id () == mesh::static_type_id ()
+          && args(1 + offset).is_string ())
+        {
+          std::string filename = args(1 + offset).string_value ();
           mesh const & msh_arg =
-            static_cast<mesh const &> (args(0).get_rep ());
+            static_cast<mesh const &> (args(0 + offset).get_rep ());
 
           if (!error_state)
             {
@@ -66,7 +81,7 @@
               filename += ".xdmf";
 
               try
-                { retval = new meshfunction (pmsh, filename); }
+                { retval = new meshfunction (label, pmsh, filename); }
               catch (std::runtime_error &)
                 { error ("error reading file"); }
             }
diff -r 68cae2998775 -r a61fc34334ca src/meshfunction.h
--- a/src/meshfunction.h	Thu Aug 07 14:56:35 2014 +0200
+++ b/src/meshfunction.h	Thu Aug 07 19:21:36 2014 +0200
@@ -29,13 +29,17 @@
   meshfunction (void)
     : octave_base_value () {}
 
-  meshfunction (dolfin::MeshFunction <std::size_t> const & _mf)
-    : octave_base_value (), pmf (new dolfin::MeshFunction <std::size_t> (_mf)) {}
+  meshfunction (std::string const & label,
+                dolfin::MeshFunction <std::size_t> const & _mf)
+    : octave_base_value (), pmf (new dolfin::MeshFunction <std::size_t> (_mf)),
+      str (label) {}
 
-  meshfunction (SHARED_PTR <dolfin::Mesh const> mesh,
+  meshfunction (std::string const & label,
+                SHARED_PTR <dolfin::Mesh const> mesh,
                 std::string const & filename)
     : octave_base_value (),
-      pmf (new dolfin::MeshFunction <std::size_t> (mesh, filename)) {}
+      pmf (new dolfin::MeshFunction <std::size_t> (mesh, filename)),
+      str (label) {}
 
   bool
   is_defined (void) const
@@ -53,9 +57,14 @@
   get_pmf (void) const
     { return pmf; }
 
+  std::string const &
+  get_str (void) const
+    { return str; }
+
   private:
 
   SHARED_PTR <dolfin::MeshFunction <std::size_t> const> pmf;
+  std::string str;
 
   DECLARE_OCTAVE_ALLOCATOR;
   DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA;