# Doubly-periodic BC syntax

Asked by Andrew McRae on 2012-11-15

I am trying to solve a problem on a doubly-periodic domain. However, the answer is clearly wrong, as it has differing values at "equivalent" nodes at the corners. This strongly suggests that the 4 corners are not being contracted to a single degree of freedom. I have no idea whether my program is wrong, or whether this is a bug, so I'm asking whether my code below is correct.

from dolfin import *
set_log_level(ERROR)
mesh = UnitSquare(1, 1)
E = FunctionSpace(mesh, 'CG', 2)
S = FunctionSpace(mesh, 'BDM', 1)

class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ( (x[0]<DOLFIN_EPS and x[1]<(1-DOLFIN_EPS)) or (x[1] < DOLFIN_EPS and x[0]<(1-DOLFIN_EPS)))
def map(self, x, y):
if (x[0] > (1 - DOLFIN_EPS)) and (x[1] > (1 - DOLFIN_EPS)):
y[0] = x[0] - 1
y[1] = x[1] - 1
elif x[0] > (1 - DOLFIN_EPS):
y[0] = x[0] - 1
y[1] = x[1]
else:
y[0] = x[0]
y[1] = x[1] - 1
bcEp = PeriodicBC(E, PeriodicBoundary())

u2_expr = Expression(("1 + x[0]*(1-x[0])","x[1]*(1-x[1])")) # curl-free
u2 = Function(interpolate(u2_expr,S))

v = TrialFunction(E)
w = TestFunction(E)
om2 = Function(E)
solve(inner(v,w)*dx==-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx,om2,bcEp)
print om2.vector().array()

As of the current dev version, the node numbering for P2 elements in UnitSquare(1,1) appears to be:
8 7 4
6 3 1
5 0 2

The difference between the values at 1 and 6 is 0, but the difference between the 0-7 pair is O(10^-16). None of the corner pairs match. The 5-8 pair has a difference of O(10^-16), while all the other pairs have O(1) differences.

In the UnitSquare(2,2) case, the numbering seems to be
24 23 19 12 11
22 21 18 10 08
20 17 16 06 03
15 13 07 05 01
14 09 04 00 02

The 0-12 and 4-19 pairs weren't numerically identical, with errors at the 10^-16 / 10^-17 scale. Again, none of the corner pairs match, while the horizontal pairs 8-22, 3-20, 1-15 have differences that are numerically 0.

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:

## This question was reopened

 Mikael Mortensen (mikael-mortensen) said on 2012-11-15: #1

Hi,

I may be wrong, but I don't think it will work out of the box with two periodic directions in one PeriodicBC object? I know there is nothing in there right now that handles those corners. I used to work with two periodic directions using two PeriodicBC objects, but you cannot do that anymore either.

The handling of periodic boundaries are about to change. If you are in a hurry you can try this branch here: lp:/~mikael-mortensen/dolfin/periodic_to_dofmap, where the periodic boundaries are incorporated directly into the dofmaps. It works with two and three periodic directions as well, but the code is under review and not mature nor supported yet.

Best regards

Mikael

Den Nov 15, 2012 kl. 3:11 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
> Description changed to:
> I am trying to solve a problem on a doubly-periodic domain. However,
> the answer is clearly wrong, as it has differing values at "equivalent"
> nodes at the corners. This strongly suggests that the 4 corners are not
> being contracted to a single degree of freedom. I have no idea whether
> my program is wrong, or whether this is a bug, so I'm asking whether my
> code below is correct.
>
> from dolfin import *
> set_log_level(ERROR)
> mesh = UnitSquare(1, 1)
> E = FunctionSpace(mesh, 'CG', 2)
> S = FunctionSpace(mesh, 'BDM', 1)
>
> class PeriodicBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and ( (x[0]<DOLFIN_EPS and x[1]<(1-DOLFIN_EPS)) or (x[1] < DOLFIN_EPS and x[0]<(1-DOLFIN_EPS)))
> def map(self, x, y):
> if (x[0] > (1 - DOLFIN_EPS)) and (x[1] > (1 - DOLFIN_EPS)):
> y[0] = x[0] - 1
> y[1] = x[1] - 1
> elif x[0] > (1 - DOLFIN_EPS):
> y[0] = x[0] - 1
> y[1] = x[1]
> else:
> y[0] = x[0]
> y[1] = x[1] - 1
> bcEp = PeriodicBC(E, PeriodicBoundary())
>
> u2_expr = Expression(("1 + x[0]*(1-x[0])","x[1]*(1-x[1])")) # curl-free
> u2 = Function(interpolate(u2_expr,S))
>
> v = TrialFunction(E)
> w = TestFunction(E)
> om2 = Function(E)
> solve(inner(v,w)*dx==-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx,om2,bcEp)
> print om2.vector().array()
>
> As of the current dev version, the node numbering for P2 elements in UnitSquare(1,1) appears to be:
> 8 7 4
> 6 3 1
> 5 0 2
>
> The difference between the values at 1 and 6 is 0, but the difference
> between the 0-7 pair is O(10^-16). None of the corner pairs match. The
> 5-8 pair has a difference of O(10^-16), while all the other pairs have
> O(1) differences.
>
> In the UnitSquare(2,2) case, the numbering seems to be
> 24 23 19 12 11
> 22 21 18 10 08
> 20 17 16 06 03
> 15 13 07 05 01
> 14 09 04 00 02
>
> The 0-12 and 4-19 pairs weren't numerically identical, with errors at
> the 10^-16 / 10^-17 scale. Again, none of the corner pairs match, while
> the horizontal pairs 8-22, 3-20, 1-15 have differences that are
> numerically 0.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2012-11-15: #2

Ah, thank you. Do you have an idea when these changes will be merged into the trunk?

I previously used multiple PeriodicBC objects in the stable version, and I noticed it complained when I switched to the dev version.

 Andrew McRae (andymc) said on 2012-11-15: #3

Thanks Mikael Mortensen, that solved my question.

 Mikael Mortensen (mikael-mortensen) said on 2012-11-15: #4

Den Nov 15, 2012 kl. 3:50 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
> Andrew McRae posted a new comment:
> Ah, thank you. Do you have an idea when these changes will be merged
> into the trunk?

Hopefully sooner than later. Meanwhile you can do like this:

from dolfin import *
set_log_level(ERROR)
mesh = UnitSquare(1, 1)

class PeriodicBoundaryX(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.)

def map(self, x, y):
y[0] = x[0] - 1
y[1] = x[1]

class PeriodicBoundaryY(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.)

def map(self, x, y):
y[0] = x[0]
y[1] = x[1] - 1.

left = PeriodicBoundaryX()
bottom = PeriodicBoundaryY()

E = FunctionSpace(mesh, 'CG', 2)
S = FunctionSpace(mesh, 'BDM', 1)

bcEp = [PeriodicBC(E, left), PeriodicBC(E, bottom)]

u2_expr = Expression(("1 + x[0]*(1-x[0])","x[1]*(1-x[1])")) # curl-free
u2 = Function(interpolate(u2_expr,S))

v = TrialFunction(E)
w = TestFunction(E)
om2 = Function(E)
#solve(inner(v,w)*dx==-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx,om2)
A = assemble(inner(v,w)*dx)
b = assemble(-inner(as_vector([w.dx(1),-w.dx(0)]),u2)*dx)
[bc.apply(A, b) for bc in bcEp]
solve(A, om2.vector(), b)

print om2.vector().array()

Mikael

>
> I previously used multiple PeriodicBC objects in the stable version, and
> I noticed it complained when I switched to the dev version.

>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2012-11-15: #5

This is in the stable, or in your branch? (Or both?)

 Mikael Mortensen (mikael-mortensen) said on 2012-11-15: #6

This should work with the latest in trunk, not the stable one.

Mikael

Den Nov 15, 2012 kl. 5:01 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
> Andrew McRae posted a new comment:
> This is in the stable, or in your branch? (Or both?)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2012-11-15: #7

Amazing, thank you!

(I didn't notice that you've split up the solve, so at first I got the same 'cannot have multiple periodic BC' error)

 Andrew McRae (andymc) said on 2013-03-07: #8

Hi Mikael,

Now that your work is in trunk, what is the correct way to do doubly-periodic?

 Mikael Mortensen (mikael-mortensen) said on 2013-03-07: #9

Hi,

Try this for a doubly periodic domain for a rectangle mesh [-pi, pi]x[-pi, pi]

class PeriodicDomain(SubDomain):
"""Doubly periodic for rectangle [-pi, pi] x [-pi, pi]"""
def inside(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], -pi) or near(x[1], -pi)) and
(not ((near(x[0], -pi) and near(x[1], pi)) or
(near(x[0], pi) and near(x[1], -pi)))) and on_boundary)

def periodic_map(x, y):
if near(x[0], pi) and near(x[1], pi):
y[0] = x[0] - 2.0*pi
y[1] = x[1] - 2.0*pi
elif near(x[0], pi):
y[0] = x[0] - 2.0*pi
y[1] = x[1]
elif near(x[1], pi):
y[0] = x[0]
y[1] = x[1] - 2.0*pi
else:
y[0] = 1e8
y[1] = 1e8

Use it like in the periodic demo.

Best regards

Mikael

Den Mar 7, 2013 kl. 1:21 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
> Status: Solved => Open
>
> Andrew McRae is still having a problem:
> Hi Mikael,
>
> Now that your work is in trunk, what is the correct way to do doubly-
> periodic?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2013-03-07: #10

Thanks. I have it working for the CG space, but if I try to use the BDM space, I get

Traceback (most recent call last):
File "solntest.py", line 9, in <module>
S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 390, in __init__
FunctionSpaceBase.__init__(self, mesh, element, constrained_domain)
File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 89, in __init__
dolfin_dofmap = cpp.DofMap(ufc_dofmap, mesh, constrained_domain)
File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 668, in __init__
_fem.DofMap_swiginit(self,_fem.new_DofMap(*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
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to complete call to function build_ufc().
*** Reason: Assertion slave_master_entities->find(0) != slave_master_entities->end() failed.
*** Where: This error was encountered inside /home/atm112/fenics/src/dolfin/current/dolfin/fem/DofMapBuilder.cpp (line 426).
*** Process: 0
*** -------------------------------------------------------------------------

Is the new PeriodicBC not yet implemented for some elements?

 Mikael Mortensen (mikael-mortensen) said on 2013-03-07: #11

Hi,

Could you give me a copy of the script you're running? I have no experience with BDM, but I know DG elements do not work with periodic yet. Did BDM use to work with the old PeriodicBC?

Mikael

Den Mar 7, 2013 kl. 2:20 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
> Andrew McRae posted a new comment:
> Thanks. I have it working for the CG space, but if I try to use the BDM
> space, I get
>
> Traceback (most recent call last):
> File "solntest.py", line 9, in <module>
> S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 390, in __init__
> FunctionSpaceBase.__init__(self, mesh, element, constrained_domain)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 89, in __init__
> dolfin_dofmap = cpp.DofMap(ufc_dofmap, mesh, constrained_domain)
> File "/home/atm112/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 668, in __init__
> _fem.DofMap_swiginit(self,_fem.new_DofMap(*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
> ***
> ***
> *** Remember to include the error message listed below and, if possible,
> *** include a *minimal* running example to reproduce the error.
> ***
> *** -------------------------------------------------------------------------
> *** Error: Unable to complete call to function build_ufc().
> *** Reason: Assertion slave_master_entities->find(0) != slave_master_entities->end() failed.
> *** Where: This error was encountered inside /home/atm112/fenics/src/dolfin/current/dolfin/fem/DofMapBuilder.cpp (line 426).
> *** Process: 0
> *** -------------------------------------------------------------------------
>
> Is the new PeriodicBC not yet implemented for some elements?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2013-03-07: #12

Yes, BDM worked with the old PeriodicBC.

from dolfin import *
mesh = UnitSquareMesh(5, 5)
class PeriodicBoundaryX(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
def map(self, x, y):
y[0] = x[0] - 1
y[1] = x[1]
pbc = PeriodicBoundaryX()
E = FunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
#S = VectorFunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)

 Mikael Mortensen (mikael-mortensen) said on 2013-03-07: #13

Ok, thanks. Looks like there's an assertion in DofMapBuilder that shouldn't be there. At least when I remove it it seems to work.

I'll make a patch as soon as I know for sure.

Mikael

Den Mar 7, 2013 kl. 3:21 PM skrev Andrew McRae:

> Question #214267 on DOLFIN changed:
>
>
> Andrew McRae is still having a problem:
> Yes, BDM worked with the old PeriodicBC.
>
> from dolfin import *
> mesh = UnitSquareMesh(5, 5)
> class PeriodicBoundaryX(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and near(x[0], 0)
> def map(self, x, y):
> y[0] = x[0] - 1
> y[1] = x[1]
> pbc = PeriodicBoundaryX()
> E = FunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
> S = FunctionSpace(mesh, 'BDM', 1, constrained_domain=pbc)
> #S = VectorFunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Andrew McRae (andymc) said on 2013-03-07: #14

Thanks.

 Andrew McRae (andymc) said on 2013-03-10: #15

By the way, is there any reason you included an 'else' clause in the code you suggested?

def periodic_map(x, y):
if near(x[0], pi) and near(x[1], pi):
y[0] = x[0] - 2.0*pi
y[1] = x[1] - 2.0*pi
elif near(x[0], pi):
y[0] = x[0] - 2.0*pi
y[1] = x[1]
elif near(x[1], pi):
y[0] = x[0]
y[1] = x[1] - 2.0*pi
else:
y[0] = 1e8
y[1] = 1e8

 Mikael Mortensen (mikael-mortensen) said on 2013-03-10: #16

On Sunday, 10 March 2013, Andrew McRae wrote:

> Question #214267 on DOLFIN changed:
>
> Andrew McRae posted a new comment:
> By the way, is there any reason you included an 'else' clause in the
> code you suggested?
>
>
Yes, to make sure all points off boundaries are mapped off the slave
boundary. Perhaps we should try to get rid of it though? With the algorithm
used to find matching entities the map is called for entities on and off
the boundaries. So when none of the first three if statements are true y
will simply not change its value. If it happens that y was on the
slave boundary in the previous check then inside(y) will still be true. A
remedy I see now could be to set default y value in C++ before the map is
called. Then the last else should not be needed. Will look into it.

Mikael

> def periodic_map(x, y):
> if near(x[0], pi) and near(x[1], pi):
> y[0] = x[0] - 2.0*pi
> y[1] = x[1] - 2.0*pi
> elif near(x[0], pi):
> y[0] = x[0] - 2.0*pi
> y[1] = x[1]
> elif near(x[1], pi):
> y[0] = x[0]
> y[1] = x[1] - 2.0*pi
> else:
> y[0] = 1e8
> y[1] = 1e8
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

 Andrew McRae (andymc) said on 2013-03-10: #17

It seems to work no matter what value it's set to (including 0, 0.5, 1... which I'd otherwise expect to cause issues!). But if it's removed altogether then things break, as I just found out...