Is it possible to specify more than one periodic boundary condition?

Question

Before Mikael's replacement of PeriodicBC, it was obvious how to apply more than one periodic boundary condition. I'm not so sure it's obvious now. I tried the natural approach:

diff --git a/demo/undocumented/periodic/python/demo_periodic.py b/demo/undocumented/periodic/python/demo_periodic.py
index bab6b55..7d1f2f3 100644
--- a/demo/undocumented/periodic/python/demo_periodic.py
+++ b/demo/undocumented/periodic/python/demo_periodic.py
@@ -44,7 +44,7 @@ class DirichletBoundary(SubDomain):
         return bool((x[1] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS)) and on_boundary)

 # Sub domain for Periodic boundary condition
-class PeriodicBoundary(SubDomain):
+class PeriodicBoundaryX(SubDomain):

     # Left boundary is "target domain" G
     def inside(self, x, on_boundary):
@@ -55,12 +55,24 @@ class PeriodicBoundary(SubDomain):
         y[0] = x[0] - 1.0
         y[1] = x[1]

+class PeriodicBoundaryY(SubDomain):
+
+    # Bottom boundary is "target domain" G
+    def inside(self, x, on_boundary):
+        return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary)
+
+    # Map right boundary (H) to left boundary (G)
+    def map(self, x, y):
+        y[0] = x[0]
+        y[1] = x[1] - 1.0
+
 # Create periodic boundary condition
-pbc = PeriodicBoundary()
+pbcx = PeriodicBoundaryX()
+pbcy = PeriodicBoundaryY()

 # Create mesh and finite element
 mesh = UnitSquareMesh(32, 32)
-V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)
+V = FunctionSpace(mesh, "CG", 1, constrained_domain=[pbcx, pbcy])


 # Create Dirichlet boundary condition

but it errors on the creation of the FunctionSpace:

[pef@caoimhe:/data/pfarrell/dolfin/dolfin/demo/undocumented/periodic/python] (master) $ python demo_periodic.py 
Traceback (most recent call last):
  File "demo_periodic.py", line 75, in <module>
    V = FunctionSpace(mesh, "CG", 1, constrained_domain=[pbcx, pbcy])
  File "/usr/lib/python2.7/dist-packages/dolfin/functions/functionspace.py", line 390, in __init__
    FunctionSpaceBase.__init__(self, mesh, element, constrained_domain)
  File "/usr/lib/python2.7/dist-packages/dolfin/functions/functionspace.py", line 76, in __init__
    "Illegal argument, not a subdomain: " + str(constrained_domain))
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/common.py", line 2182, in dolfin_error
    return _common.dolfin_error(*args)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create function space.
*** Reason:  Illegal argument, not a subdomain: [<__main__.PeriodicBoundaryX; proxy of <Swig Object of type 'boost::shared_ptr< dolfin::SubDomain > *' at 0x1dba060> >, <__main__.PeriodicBoundaryY; proxy of <Swig Object of type 'boost::shared_ptr< dolfin::SubDomain > *' at 0x1dba0f0> >].
*** Where:   This error was encountered inside functionspace.py.
*** Process: 0
*** -------------------------------------------------------------------------

What's the DOLFINic way to express this?

Answer 1

It's still possible with two - or even three periodic directions. You just have to place it all in the same map. I use this for two periodic directions (x and z) in a rectangular channel flow problem:

Lx = 2.*pi
Ly = 2.
Lz = pi
Nx = 257
Ny = 193
Nz = 193
mesh = BoxMesh(0., -Ly/2., -Lz/2., Lx, Ly/2., Lz/2., Nx, Ny, Nz)

class PeriodicDomain(SubDomain):

    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool((near(x[0], 0) or near(x[2], -Lz/2.)) and 
            (not ((near(x[0], Lx) and near(x[2], -Lz/2.)) or 
                  (near(x[0], 0) and near(x[2], Lz/2.)))) and on_boundary)

    def map(self, x, y):
        if near(x[0], Lx) and near(x[2], Lz/2.):
            y[0] = x[0] - Lx
            y[1] = x[1] 
            y[2] = x[2] - Lz
        elif near(x[0], Lx):
            y[0] = x[0] - Lx
            y[1] = x[1]
            y[2] = x[2]
        elif near(x[2], Lz/2.):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - Lz
        else:
            y[0] = -1000
            y[1] = -1000
            y[2] = -1000

See if you can figure out the logic:-) The final else just needs to map somewhere off the boundary, doesn't need to be 1000. I agree two subdomains is more natural.

For the 2D demo this should work:

# Sub domain for Periodic boundary condition
class PeriodicBoundary(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 bool((near(x[0], 0) or near(x[1], 0)) and 
                (not ((near(x[0], 0) and near(x[1], 1)) or 
                        (near(x[0], 1) and near(x[1], 0)))) and on_boundary)

    def map(self, x, y):
        if near(x[0], 1) and near(x[1], 1):
            y[0] = x[0] - 1.
            y[1] = x[1] - 1.
        elif 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.

Comments

Hi Mikael,

Thanks for clearing that up, I understand it.

I think the ability to separate them cleanly was useful -- I think it's unfortunate/nonobvious that you have to pack all of the separate periodic boundaries into the one object.

Could FunctionSpace be made to take in a list of constrained_domains?

Or maybe we could compose the SubDomains somehow before passing it to the FunctionSpace constructor?

---

I think the ability to separate them cleanly was useful -- I think it's unfortunate/nonobvious that you have to pack all of the separate periodic boundaries into the one object.

I agree and have argued in this favor. I think the single map is too complex, because it needs to handle all the corners. With two or three SubDomains containing one single periodic direction each the handling of corners could be taken care of under the hood.

Then again, if you have these recepies for maps available it's not really a big problem. A bigger problem as I see it with the current design is the fact that periodic maps of mesh entities now are recomputed for every new FunctionSpace.

Mikael

---

Hi Mikael, may I ask why it is necessary to map some points to -1000 in the first example of your answer but not in the second? Generally, for which mesh nodes will the map() function be called, and where should each vertex be mapped to (depending on its position on the boundary or inside the mesh)? My understanding from your examples is that only the vertices on the boundaries that are to be identified need to be mapped, but it would be good to have that confirmed. Many thanks!