diff --git a/src/bvpy/boundary_conditions/boundary.py b/src/bvpy/boundary_conditions/boundary.py
index 395cf123c1ebd89d6f281987ef1ebdd02891579f..8a8979a35c2e0b4f4a93cfb664cecc6f356cb7f6 100644
--- a/src/bvpy/boundary_conditions/boundary.py
+++ b/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")
diff --git a/src/bvpy/boundary_conditions/dirichlet.py b/src/bvpy/boundary_conditions/dirichlet.py
index b1d018d0e78e1831a70f32c5e422f7b85451d45a..1220fa39b7c6eb22dad36b3580cd969f3cc813e6 100644
--- a/src/bvpy/boundary_conditions/dirichlet.py
+++ b/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)
diff --git a/src/bvpy/boundary_conditions/neumann.py b/src/bvpy/boundary_conditions/neumann.py
index a6b4c3c7e9e887e0f05161ee007baeb3addd45a5..db07aaaa52e89c68c0303132ada0c207d6381ed2 100644
--- a/src/bvpy/boundary_conditions/neumann.py
+++ b/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:
diff --git a/test/test_boundary.py b/test/test_boundary.py
index 5725177e9fdda963ce025404c923124c7415e1a1..361517202f5773558c2b16b74389f73f98cc7784 100644
--- a/test/test_boundary.py
+++ b/test/test_boundary.py
@@ -1,4 +1,5 @@
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]))