Clean boundary.py

src/bvpy/boundary_conditions/boundary.py

@@ -9,6 +9,7 @@
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
 #           Olivier Ali <olivier.ali@inria.fr>
+#           Gonzalo Revilla <gonzalo.revilla-mut@inria.fr>
 #
 #       File maintainer(s):
 #           Olivier Ali <olivier.ali@inria.fr>
@@ -19,177 +20,94 @@
 #           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
-
+import numpy as np
 import fenics as fe
-from sympy import And, Function, Or, Symbol, Not
-from sympy.parsing.sympy_parser import parse_expr
-
-from bvpy import MPI_COMM_WORLD, logger
-
-_XYZ = {'x': Symbol('x[0]'),
-        'y': Symbol('x[1]'),
-        'z': Symbol('x[2]'),
-        'all': Symbol('on_boundary'),
-        'near': Function('near')}
 
+from bvpy import MPI_COMM_WORLD
 
-class Boundary(object):
-    """Defines the boundary of a domain.
-
-    Parameters
-    ----------
-    expr : str
-        An expression used to compute the boundary.
 
+class Boundary:
+    """Allows to identify subdomains of the boundary given a parametric expression.
 
     Attributes
     ----------
-    _string : str
-        The sting passed to the class
-    _expr : str
-        The parsed string for fenics.
-    _domain : SubDomain
-        The subdomain fenics defining the boundary.
-    _kwargs : dict
-        Optional argument.
-
-    Examples
-    --------
-    >>> Boundary("near(x+a, 2)", a=0) & Boundary("all")
+    string : str
+        Parametric expression that defines the boundary subdomain.
+    kwargs : dict
+        Parameters (if any) of the parametric expression.
+    subdomain : fe.CompiledSubDomain
+        Fenics object that represents the current subdomain.
+
+    Methods
+    -------
+    mark_mesh_function()
+        Given a fenics mesh function, mark the elements of the current boundary subdomain.
 
     """
     def __init__(self, expr, **kwargs):
-        """Constructor of the Boundary class.
-
+        """
         Parameters
         ----------
         expr : str
-            An expression used to compute the boundary.
+            An expression used to compute the subdomain.
+        **kwargs
+            Expression parameters that will be sent to the fenics object. (See next example)
 
-        Yields
-        -------
-        Boundary
-            A boundary instance.
+        Examples
+        --------
+        >>> b1 = Boundary("near(x, 0)") & Boundary("near(x+a, 2)", a=0)
+        >>> b2 = Boundary("all")
+        >>> b3 = Boundary("near(x, 0)") | Boundary("near(x, 1)")
+        >>> b4 = ~Boundary("near(x, 1)")
 
         """
-        self._string = expr
-        self._expr = parse_expr(expr, local_dict=_XYZ, evaluate=False)
-        self._domain = fe.CompiledSubDomain(str(self._expr), **kwargs,
-                                            mpi_comm=MPI_COMM_WORLD)
-        self._kwargs = kwargs
+        self.__string = expr
+        replacements = {'x': 'x[0]', 'y': 'x[1]', 'z': 'x[2]', 'all': 'on_boundary', '~': '!'}
+        for old, new in replacements.items():
+            expr = expr.replace(old, new)
+        self.__kwargs = kwargs
 
-    def __call__(self):
-        """Fetches the fenics object defining the boundary.
+        self.__subdomain = fe.CompiledSubDomain(expr, **kwargs, mpi_comm=MPI_COMM_WORLD)
 
-        Returns
-        -------
-        Domain
-            Description of returned object.
+    @property
+    def string(self):
+        return self.__string
 
-        """
-        return self._domain
-
-    def __or__(self, other):
-        """Intersects two boundary domains together.
-
-        Parameters
-        ----------
-        other : type
-            Description of parameter `other`.
+    @property
+    def kwargs(self):
+        return self.__kwargs
 
-        Returns
-        -------
-        type
-            Description of returned object.
+    @property
+    def subdomain(self):
+        return self.__subdomain
 
+    def __or__(self, other):
+        """Creates a boundary instance: The addition of the 2 given boundaries objects.
         """
-        new = Boundary('near(x, 0)')
-        s1 = Symbol(str(self._expr))
-        s2 = Symbol(str(other._expr))
-        new._expr = Or(s1, s2)
-        new._string = self._string+' | '+other._string
-        new._kwargs = {**self._kwargs, **other._kwargs}
-        str_expr = str(new._expr).replace('&', '&&')\
-                                 .replace('|', '||')\
-                                 .replace('~', '!')
-        new._domain = fe.CompiledSubDomain(str_expr, **new._kwargs,
-                                           mpi_comm=MPI_COMM_WORLD)
-        return new
+        or_string = '(' + self.string+' || '+other.string+')'
+        return Boundary(or_string, **{**self.kwargs, **other.kwargs})
 
     def __and__(self, other):
-        """Adds two boundary domains together.
-
-        Parameters
-        ----------
-        other : type
-            Description of parameter `other`.
-
-        Returns
-        -------
-        type
-            Description of returned object.
-
+        """Creates a boundary instance: The intersection of the 2 given boundaries objects.
         """
-        new = Boundary('near(x, 0)')
-        s1 = Symbol(str(self._expr))
-        s2 = Symbol(str(other._expr))
-        new._expr = And(s1, s2)
-        new._string = self._string+' & '+other._string
-        new._kwargs = {**self._kwargs, **other._kwargs}
-        str_expr = str(new._expr).replace('&', '&&')\
-                                 .replace('|', '||')\
-                                 .replace('~', '!')
-        logger.debug('expression in and: %s ', str_expr)
-        new._domain = fe.CompiledSubDomain(str_expr, **new._kwargs,
-                                           mpi_comm=MPI_COMM_WORLD)
-
-        return new
+        and_string = '(' + self.string+' && '+other.string+')'
+        return Boundary(and_string, **{**self.kwargs, **other.kwargs})
 
     def __invert__(self):
-        """Inverts a boundary domain.
-
-        Parameters
-        ----------
-        other : type
-            Description of parameter `other`.
-        Returns
-        -------
-        type
-            Description of returned object.
+        """Creates a boundary instance: the inverse of the current one.
         """
-        new = Boundary('near(x, 0)')
-        s1 = Symbol(str(self._expr))
-        new._expr = Not(s1)
-        new._string = '~ ('+self._string+')'
-        new._kwargs = self._kwargs
-        str_expr = str(new._expr).replace('&', '&&')\
-                                 .replace('|', '||')\
-                                 .replace('~', '!')
-        logger.debug('expression in and: %s ', str_expr)
-        new._domain = fe.CompiledSubDomain(str_expr, **new._kwargs,
-                                           mpi_comm=MPI_COMM_WORLD)
-
-        return new
-
-    def markMeshFunction(self, valueTrue,
-                         mesh=None, valueFalse=0, meshfunc=None,
-                         dim_topology=0, type='double'):
-        """Short summary.
+        invert_string = '~('+self.string+')'
+        return Boundary(invert_string, **self.kwargs)
+
+    def mark_mesh_function(self, mf :fe.MeshFunction, value_true: int=1):
+        """Given a fenics mesh function, mark the elements of the current boundary subdomain.
 
         Parameters
         ----------
-        valueTrue : type
-            Description of parameter `valueTrue`.
-        mesh : type
-            Description of parameter `mesh` (the default is None).
-        valueFalse : type
-            Description of parameter `valueFalse` (the default is 0).
-        meshfunc : type
-            Description of parameter `meshfunc` (the default is None).
-        dim_topology : type
-            Description of parameter `dim_topology` (the default is 0).
-        type : type
-            Description of parameter `type` (the default is 'double').
+        mf : fe.Function
+            The fenics mesh function that will be marked.
+        value_true : int
+            Label that will have the elements of the current boundary subdomain.
 
         Returns
         -------
@@ -197,11 +115,12 @@ class Boundary(object):
             Description of returned object.
 
         """
-        if meshfunc is not None:
-            self._domain.mark(meshfunc, valueTrue) # this mark the subdomain meshfunc using the condition
-        else:
-            meshfunc = fe.MeshFunction(type, mesh, dim_topology)
-            meshfunc.set_all(valueFalse)
-            self._domain.mark(meshfunc, valueTrue)
-
-        return meshfunc
+        if np.any(mf.array() == value_true):
+            raise ValueError(f'The mesh function already contains some elements with the label {value_true}')
+
+        self.subdomain.mark(mf, value_true)
+
+        if np.all(mf.array() != value_true):
+            raise ValueError(f'The mesh function has not marked any vertex')
+
+        # TODO : assert that all the marked vertex are on the boundary. Maybe use fe.BoundaryMesh(mesh, "exterior")

src/bvpy/boundary_conditions/dirichlet.py

@@ -86,16 +86,16 @@ class ConstantDirichlet(object):
         """
         if self._subspace is None:
             dir = fe.DirichletBC(functionSpace,
-                                 fe.Constant(self._val), self._boundary(),
+                                 fe.Constant(self._val), self._boundary.subdomain,
                                  method=self._method)
         else:
             dir = fe.DirichletBC(functionSpace.sub(self._subspace),
                                  fe.Constant(self._val),
-                                 self._boundary(), method=self._method)
+                                 self._boundary.subdomain, method=self._method)
             
         if not dir.get_boundary_values():
             logger.warning('No nodes are marked for this domain: '
-                           + str(self._boundary._string))
+                           + str(self._boundary.string))
 
         return dir
 
@@ -109,7 +109,7 @@ class ConstantDirichlet(object):
 
         """
         return "<ConstantDirichlet object, location: "\
-               + str(self._boundary._string) + ", value: " + str(self._val)+">"
+               + str(self._boundary.string) + ", value: " + str(self._val)+">"
 
 
 class VariableDirichlet(object):
@@ -193,17 +193,17 @@ class VariableDirichlet(object):
 
         if self._subspace is None:
             dir = fe.DirichletBC(functionSpace, expr,
-                                 self._boundary(),
+                                 self._boundary.subdomain,
                                  method=self._method)
         else:
             dir = fe.DirichletBC(functionSpace.sub(self._subspace),
                                  expr,
-                                 self._boundary(),
+                                 self._boundary.subdomain,
                                  method=self._method)
 
         if not dir.get_boundary_values():
             logger.warning('No nodes are marked for this domain: '
-                           + str(self._boundary._string))
+                           + str(self._boundary.string))
 
         return dir
 
@@ -217,7 +217,7 @@ class VariableDirichlet(object):
 
         """
         return "<VariableDirichlet object, location: "\
-               + str(self._boundary._string) + ", value: "\
+               + str(self._boundary.string) + ", value: "\
                + str(self._val) + ">"
 
 
@@ -345,11 +345,11 @@ class NormalDirichlet(object):
 
         """
         return "<NormalDirichlet object, location: "\
-               + str(self._boundary._string) + ", value: " + str(self._val)+">"
+               + str(self._boundary.string) + ", value: " + str(self._val)+">"
 
     def apply(self, functionSpace):
         return fe.DirichletBC(functionSpace,
                               boundary_normal(functionSpace.mesh(),
                                               scale=self._val),
-                              self._boundary(),
+                              self._boundary.subdomain,
                               method=self._method)

src/bvpy/boundary_conditions/neumann.py

@@ -69,11 +69,11 @@ class NeumannInterface(object):
         raise NotImplementedError
 
     def markMeshFunction(self, mf):
-        self._boundary.markMeshFunction(self._id, meshfunc=mf)
+        self._boundary.mark_mesh_function(mf, self._id)
 
     def __repr__(self):
         return "<" + self.__class__.__name__ + " object, location: " +\
-               str(self._boundary._string) + ", value: " + str(self._val) + ">"
+               str(self._boundary.string) + ", value: " + str(self._val) + ">"
 
 
 class ConstantNeumann(NeumannInterface):
@@ -186,7 +186,7 @@ class NormalNeumann(NeumannInterface):
 
         """
         return "<NormalNeumann object, location: "\
-               + str(self._boundary._string) + ", value: " + str(self._val)+">"
+               + str(self._boundary.string) + ", value: " + str(self._val)+">"
 
     def apply(self, functionSpace, mf=None):
         if mf is not None:

test/test_boundary.py

 import unittest
+import numpy as np
 from bvpy.boundary_conditions.boundary import Boundary
 import fenics as fe
 
@@ -7,37 +8,67 @@ class TestBoundary(unittest.TestCase):
     def setUp(self):
         """Initiates the test.
         """
+        self.mesh = fe.UnitSquareMesh(2, 2)
+        self.mf = fe.MeshFunction("size_t", self.mesh, 0)  # 0 means defined on the mesh vertices
 
     def tearDown(self):
         """Concludes and closes the test.
         """
 
+    def assert_mf(self, computed, expected):
+        self.assertTrue(np.array_equal(computed, expected))
+
+    def test_near(self):
+        b1 = Boundary("near(x, 0)")
+        b1.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([1, 0, 0, 1, 0, 0, 1, 0, 0]))
+
+    def test_all(self):
+        b1 = Boundary("all")
+        b1.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([1, 1, 1, 1, 0, 1, 1, 1, 1]))
+
+    def test_kwargs(self):
+        b1 = Boundary("near(x, a)", a=1)
+        b1.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([0, 0, 1, 0, 0, 1, 0, 0, 1]))
+
     def test_and(self):
         b1 = Boundary("near(x, 0)")
         b2 = Boundary("near(y, 0)")
+
         b3 = b1 & b2
-        self.assertEqual(b3._string, "near(x, 0) & near(y, 0)")
+        b3.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([1, 0, 0, 0, 0, 0, 0, 0, 0]))
 
     def test_or(self):
         b1 = Boundary("near(x, 0)")
         b2 = Boundary("near(y, 0)")
+
         b3 = b1 | b2
-        self.assertEqual(b3._string, "near(x, 0) | near(y, 0)")
+        b3.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([1, 1, 1, 1, 0, 0, 1, 0, 0]))
 
     def test_invert(self):
         b1 = Boundary("near(x, 0)")
         b2 = ~ b1
-        self.assertEqual(b2._string, "~ (near(x, 0))")
+        b2.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([0, 1, 1, 0, 1, 1, 0, 1, 1]))
 
-    def test_mark(self):
+    def test_and_or(self):
         b1 = Boundary("near(x, 0)")
-        mesh = fe.UnitSquareMesh(2, 2)
-        mf = fe.MeshFunction("size_t", mesh, 0, 0)
-        b1.markMeshFunction(1, meshfunc=mf)
-        self.assertEqual(mf[0], 1)
-        self.assertEqual(mf[3], 1)
-        self.assertEqual(mf[6], 1)
-        mf2 = b1.markMeshFunction(1, type="size_t", mesh=mesh)
-        self.assertEqual(mf2[0], 1)
-        self.assertEqual(mf2[3], 1)
-        self.assertEqual(mf2[6], 1)
+        b2 = Boundary("near(y, 0)")
+        b3 = Boundary("near(x, 1)")
+
+        b4 = (b1 & ~b2) | b3
+        b4.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([0, 0, 1, 1, 0, 1, 1, 0, 1]))
+
+    def test_or_and(self):
+        b1 = Boundary("near(x, 0)")
+        b2 = Boundary("near(y, 0)")
+        b3 = Boundary("near(x, 1)")
+
+        b4 = b1 & (b2 | b3)
+        b4.mark_mesh_function(self.mf, 1)
+        self.assert_mf(self.mf.array(), np.array([1, 0, 0, 0, 0, 0, 0, 0, 0]))