Add option for flat Domain + constrain functionSpace with periodicity

example/poisson.py

-from bvpy.templates.vforms import PoissonForm
-from bvpy.templates.domains import HemiSphere
-from bvpy.boundary_conditions import dirichlet
-from bvpy import BVP
-
-# bvp = BVP()
-
-hsphere = HemiSphere(radius=1, resolution=0.1)
-# bvp.set_domain(domain)
-
-poisson = PoissonForm(source=10)
-poisson.add_neumann_term(100)
-# bvp.set_vform(vform)
-
-# boundary = Boundary('all')
-diri = dirichlet(0, boundary='near(x, 0)')
-# bvp.add_boundary_condition(diri)
-
-bvp = BVP(domain=hsphere, vform=poisson, bc=diri)
-bvp.solve()

example/poisson_rectangle.py

-from bvpy.templates.vforms import PoissonForm
-from bvpy.templates.domains import Rectangle
-from bvpy.boundary_conditions import dirichlet
-from bvpy import BVP
-
-# bvp = BVP()
-
-rectangle = Rectangle(resolution=0.05)
-# bvp.set_domain(domain)
-
-poisson = PoissonForm(source=0)
-poisson.add_neumann_term(100, boundary='near(y,1)')
-# bvp.set_vform(vform)
-
-# boundary = Boundary('all')
-diri = dirichlet(0, boundary='near(x,0)')
-# bvp.add_boundary_condition(diri)
-
-bvp = BVP(domain=rectangle, vform=poisson, bc=diri)
-bvp.solve()
-
-from bvpy.utils.visu import plot
-plot(bvp.solution)

src/bvpy/boundary_conditions/periodic.py

+import fenics as fe
+
+
+class SquarePeriodicBoundaryLR(fe.SubDomain):
+
+    # Left boundary is "target domain" G
+    def inside(self, x, on_boundary):
+        return on_boundary and fe.near(x[0], 0)
+
+    # Map right boundary (H) to left boundary (G)
+    def map(self, x, y, z):
+        y[0] = x[0] - 1.0
+        y[1] = x[1]
+
+
+# Sub domain for Periodic boundary condition
+class SquarePeriodicBoundaryBT(fe.SubDomain):
+
+    # Bottom boundary is "target domain" G
+    def inside(self, x, on_boundary):
+        return on_boundary and fe.near(x[1], 0)
+
+    # Map bottom boundary (H) to top boundary (G)
+    def map(self, x, y):
+        y[0] = x[0]
+        y[1] = x[1] - 1.0
+
+
+# Sub domain for Periodic boundary condition
+class SquarePeriodicBoundary(fe.SubDomain):
+
+    # Left boundary is "target domain" G
+    def inside(self, x, on_boundary):
+        # return True if on left or bottom boundary
+        # AND NOT on one of the two corners (0, 1) and (1, 0)
+        return (fe.near(x[0], 0)
+                or fe.near(x[1], 0)) and (not ((fe.near(x[0], 0)
+                                          and fe.near(x[1], 1))
+                                          or (fe.near(x[0], 1)
+                                              and fe.near(x[1],
+                                                          0)))) and on_boundary
+
+    def map(self, x, y):
+        if fe.near(x[0], 1) and fe.near(x[1], 1):
+            y[0] = x[0] - 1.
+            y[1] = x[1] - 1.
+        elif fe.near(x[0], 1):
+            y[0] = x[0] - 1.
+            y[1] = x[1]
+        else:   # near(x[1], 1)
+            y[0] = x[0]
+            y[1] = x[1] - 1.

src/bvpy/bvp.py

@@ -76,7 +76,13 @@ class BVP(object):
         else:
             logger.error('no elements found in vform:{form}'.format(form=vform))
 
-        self.functionSpace = fe.FunctionSpace(self.domain.mesh, reduce(operator.mul, elems))
+        if self.vform._constrain is not None:
+            self.functionSpace = fe.FunctionSpace(self.domain.mesh,
+                                                  reduce(operator.mul, elems),
+                                                  constrained_domain=self.vform._constrain)
+        else:
+            self.functionSpace = fe.FunctionSpace(self.domain.mesh,
+                                                  reduce(operator.mul, elems))
 
         self.trialFunction = fe.TrialFunction(self.functionSpace)
         self.testFunction = fe.TestFunction(self.functionSpace)

src/bvpy/domains/abstract.py

@@ -139,7 +139,8 @@ class AbstractDomain(ABC):
                  cell_type='triangle',
                  resolution=1,
                  verbosity=False,
-                 algorithm='Delaunay'):
+                 algorithm='Delaunay',
+                 dim=3):
         """Generates an instance of the AbstractDomain class.
 
         Parameters
@@ -182,6 +183,7 @@ class AbstractDomain(ABC):
         self.set_resolution(resolution)
         self.set_meshing_algorithm(algorithm)
 
+        self.dim = dim
         self.mesh = None
         self.cdata = None
         self.bdata = None
@@ -408,8 +410,8 @@ class AbstractDomain(ABC):
         """
         path = self._dir_name + '/mesh.msh'
         mesh_io = meshio.read(path, file_format='gmsh')
-        # if ignore_z:
-        #     mesh.points = mesh.points[:, :2]
+
+        mesh_io.points = mesh_io.points[:, :self.dim]
         # mesh.prune()
 
         mesh = fe.Mesh()

src/bvpy/vforms/abstract.py

@@ -49,6 +49,7 @@ class AbstractVform(ABC):
         self.dx = fe.dx
 
         self._bterm = []
+        self._constrain = None
 
     @property
     def elements(self):
@@ -82,6 +83,8 @@ class AbstractVform(ABC):
         elif isinstance(val, (str, FunctionType)):
             self._bterm.append((create_expression(val, degree=degree), boundary))
 
+    def constrain(self, constrain):
+        self._constrain = constrain()
 
 
     @abstractmethod