DOLFIN

Periodic boundary condition between specific vertices

Asked by Peter Kekenes-Huskey

Dear all,
I would like to apply a periodic boundary condition between specific nodes on a mesh. For instance, in the included example (based off a script in the demo directory), I would like to apply a PBC to x1=[0,0.5] and x2=[1,0.5] (along the surface). All other surface vertices are assigned a dirichlet condition. I’m able to find the vertices using the ‘inside()’ functions I’ve defined, but I don’t know how to tell the declare the mapping between x1 and x2 using ‘map()’ in the PBC function.

I’m sure there are simpler ways of doing this for a grid with aligned surface vertices, but I’d like a general solution since my meshes are very unlikely to unaligned. In practice, I might want to select a small region of vertices near x1 and x2 for this condition.

I appreciate any help on this problem!

Best,
Pete

from dolfin import *

# Periodic BCs don't work with Epetra
if parameters["linear_algebra_backend"] == "Epetra":
    print "Sorry, this demo does not work with the Epetra backend"
    import sys
    sys.exit(0)

# Create mesh and finite element
mesh = UnitSquare(2,2)
for i,c in enumerate( mesh.coordinates() ):
    print c

import numpy as np
coord1 = np.array([0.0,0.5])
coord2 = np.array([1.0,0.5])

Vv = FunctionSpace(mesh, "CG", 1)

# Source term
class Source(Expression):
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)

# Sub domain for Dirichlet boundary condition
# x free
class DirichletBoundary0(SubDomain):
    def inside(self, x, on_boundary):
      return bool( \
        np.linalg.norm( x- coord1) > DOLFIN_EPS and \
        np.linalg.norm( x- coord2) > DOLFIN_EPS and \
        on_boundary)

# Sub domain for Periodic boundary condition
class PeriodicBoundary0(SubDomain):

    def inside(self, x, on_boundary):
      result = bool( \
        (np.linalg.norm( x- coord1) < DOLFIN_EPS or \
        np.linalg.norm( x- coord2) < DOLFIN_EPS) and \
        on_boundary)
      print result
      return result

    # field component 1
    def map(self, x, y):
        #print "x y"
        print x
        y[0] = x[0] - 1.0
        y[1] = x[1]

# Create Dirichlet boundary condition
u0 = Constant(1.0)

dbc0 = DirichletBoundary0()
bc0 = DirichletBC(Vv, u0, dbc0)

pbc0 = PeriodicBoundary0()
bc1 = PeriodicBC(Vv, pbc0)

bcs = [bc0, bc1]

# Define variational problem
u = TrialFunction(Vv)
v = TestFunction(Vv)
f = Source()
a = dot(grad(u), grad(v))*dx
L = f*v*ds

# Compute solution
u = Function(Vv)
solve(a == L, u, bcs)

# Save solution to file
file = File("periodic.pvd")
file << u

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Johan Hake (johan-hake) said : 		#1

I am not expert on the inner workings of PeriodicBC. It looks to me that
the coordinate should be set correctly. Not sure that it does not work
because you only use map 1 dof?

You might want to try to follow the logic in PeriodicBC.cpp. Not sure it
gets you any wiser. I bailed out after a minute...

Johan

On 06/16/2012 03:31 AM, Peter Kekenes-Huskey wrote:
> New question #200584 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/200584
>
> Dear all,
> I would like to apply a periodic boundary condition between specific nodes on a mesh. For instance, in the included example (based off a script in the demo directory), I would like to apply a PBC to x1=[0,0.5] and x2=[1,0.5] (along the surface). All other surface vertices are assigned a dirichlet condition. I’m able to find the vertices using the ‘inside()’ functions I’ve defined, but I don’t know how to tell the declare the mapping between x1 and x2 using ‘map()’ in the PBC function.
>
> I’m sure there are simpler ways of doing this for a grid with aligned surface vertices, but I’d like a general solution since my meshes are very unlikely to unaligned. In practice, I might want to select a small region of vertices near x1 and x2 for this condition.
>
> I appreciate any help on this problem!
>
> Best,
> Pete
>
>
> from dolfin import *
>
> # Periodic BCs don't work with Epetra
> if parameters["linear_algebra_backend"] == "Epetra":
> print "Sorry, this demo does not work with the Epetra backend"
> import sys
> sys.exit(0)
>
> # Create mesh and finite element
> mesh = UnitSquare(2,2)
> for i,c in enumerate( mesh.coordinates() ):
> print c
>
> import numpy as np
> coord1 = np.array([0.0,0.5])
> coord2 = np.array([1.0,0.5])
>
>
> Vv = FunctionSpace(mesh, "CG", 1)
>
> # Source term
> class Source(Expression):
> def eval(self, values, x):
> dx = x[0] - 0.5
> dy = x[1] - 0.5
> values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)
>
> # Sub domain for Dirichlet boundary condition
> # x free
> class DirichletBoundary0(SubDomain):
> def inside(self, x, on_boundary):
> return bool( \
> np.linalg.norm( x- coord1)> DOLFIN_EPS and \
> np.linalg.norm( x- coord2)> DOLFIN_EPS and \
> on_boundary)
>
>
> # Sub domain for Periodic boundary condition
> class PeriodicBoundary0(SubDomain):
>
> def inside(self, x, on_boundary):
> result = bool( \
> (np.linalg.norm( x- coord1)< DOLFIN_EPS or \
> np.linalg.norm( x- coord2)< DOLFIN_EPS) and \
> on_boundary)
> print result
> return result
>
> # field component 1
> def map(self, x, y):
> #print "x y"
> print x
> y[0] = x[0] - 1.0
> y[1] = x[1]
>
> # Create Dirichlet boundary condition
> u0 = Constant(1.0)
>
> dbc0 = DirichletBoundary0()
> bc0 = DirichletBC(Vv, u0, dbc0)
>
> pbc0 = PeriodicBoundary0()
> bc1 = PeriodicBC(Vv, pbc0)
>
> bcs = [bc0, bc1]
>
> # Define variational problem
> u = TrialFunction(Vv)
> v = TestFunction(Vv)
> f = Source()
> a = dot(grad(u), grad(v))*dx
> L = f*v*ds
>
> # Compute solution
> u = Function(Vv)
> solve(a == L, u, bcs)
>
> # Save solution to file
> file = File("periodic.pvd")
> file<< u
>
>

Revision history for this message
Peter Kekenes-Huskey (huskeypm) said : 		#2

Thanks Johan. I think it might take me a little while to understand exactly how the cpp file is doing things, though extract_dof_pairs() seems to be the important function to consider. It too calls a 'sub_domain->map' function, which is presumably the one I would need to set correctly. Alternatively, it is easiest just to enforce the vertices to align on opposing faces at the mesh generation stage?

cheers
pete

On Jun 20, 2012, at 2:55 PM, Johan Hake wrote:

> Your question #200584 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200584
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> I am not expert on the inner workings of PeriodicBC. It looks to me that
> the coordinate should be set correctly. Not sure that it does not work
> because you only use map 1 dof?
>
> You might want to try to follow the logic in PeriodicBC.cpp. Not sure it
> gets you any wiser. I bailed out after a minute...
>
> Johan
>
> On 06/16/2012 03:31 AM, Peter Kekenes-Huskey wrote:
>> New question #200584 on DOLFIN:
>> https://answers.launchpad.net/dolfin/+question/200584
>>
>> Dear all,
>> I would like to apply a periodic boundary condition between specific nodes on a mesh. For instance, in the included example (based off a script in the demo directory), I would like to apply a PBC to x1=[0,0.5] and x2=[1,0.5] (along the surface). All other surface vertices are assigned a dirichlet condition. I’m able to find the vertices using the ‘inside()’ functions I’ve defined, but I don’t know how to tell the declare the mapping between x1 and x2 using ‘map()’ in the PBC function.
>>
>> I’m sure there are simpler ways of doing this for a grid with aligned surface vertices, but I’d like a general solution since my meshes are very unlikely to unaligned. In practice, I might want to select a small region of vertices near x1 and x2 for this condition.
>>
>> I appreciate any help on this problem!
>>
>> Best,
>> Pete
>>
>>
>> from dolfin import *
>>
>> # Periodic BCs don't work with Epetra
>> if parameters["linear_algebra_backend"] == "Epetra":
>> print "Sorry, this demo does not work with the Epetra backend"
>> import sys
>> sys.exit(0)
>>
>> # Create mesh and finite element
>> mesh = UnitSquare(2,2)
>> for i,c in enumerate( mesh.coordinates() ):
>> print c
>>
>> import numpy as np
>> coord1 = np.array([0.0,0.5])
>> coord2 = np.array([1.0,0.5])
>>
>>
>> Vv = FunctionSpace(mesh, "CG", 1)
>>
>> # Source term
>> class Source(Expression):
>> def eval(self, values, x):
>> dx = x[0] - 0.5
>> dy = x[1] - 0.5
>> values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)
>>
>> # Sub domain for Dirichlet boundary condition
>> # x free
>> class DirichletBoundary0(SubDomain):
>> def inside(self, x, on_boundary):
>> return bool( \
>> np.linalg.norm( x- coord1)> DOLFIN_EPS and \
>> np.linalg.norm( x- coord2)> DOLFIN_EPS and \
>> on_boundary)
>>
>>
>> # Sub domain for Periodic boundary condition
>> class PeriodicBoundary0(SubDomain):
>>
>> def inside(self, x, on_boundary):
>> result = bool( \
>> (np.linalg.norm( x- coord1)< DOLFIN_EPS or \
>> np.linalg.norm( x- coord2)< DOLFIN_EPS) and \
>> on_boundary)
>> print result
>> return result
>>
>> # field component 1
>> def map(self, x, y):
>> #print "x y"
>> print x
>> y[0] = x[0] - 1.0
>> y[1] = x[1]
>>
>> # Create Dirichlet boundary condition
>> u0 = Constant(1.0)
>>
>> dbc0 = DirichletBoundary0()
>> bc0 = DirichletBC(Vv, u0, dbc0)
>>
>> pbc0 = PeriodicBoundary0()
>> bc1 = PeriodicBC(Vv, pbc0)
>>
>> bcs = [bc0, bc1]
>>
>> # Define variational problem
>> u = TrialFunction(Vv)
>> v = TestFunction(Vv)
>> f = Source()
>> a = dot(grad(u), grad(v))*dx
>> L = f*v*ds
>>
>> # Compute solution
>> u = Function(Vv)
>> solve(a == L, u, bcs)
>>
>> # Save solution to file
>> file = File("periodic.pvd")
>> file<< u
>>
>>
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/dolfin/+question/200584/+confirm?answer_id=0
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/dolfin/+question/200584
>
> You received this question notification because you asked the question.

Revision history for this message
Johan Hake (johan-hake) said : 		#3

If you remove the check for coord2 in the inside method the mapping works.

   (np.linalg.norm(x-coord1) < DOLFIN_EPS) and on_boundary

This method should give you a set of coordinates with associated dofs.
Then map is applied on all coordinates to check if the mapped coordinate
is inside the domain. In a way you are trying to map back the
coordinates. If the mapped coordinate is inside the domain you have a
match.

So to let inside only return true for the first domain is a good start!

Johan

On 06/21/2012 12:11 AM, Peter Kekenes-Huskey wrote:
> Question #200584 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200584
>
> Status: Answered => Open
>
> Peter Kekenes-Huskey is still having a problem:
> Thanks Johan. I think it might take me a little while to understand
> exactly how the cpp file is doing things, though extract_dof_pairs()
> seems to be the important function to consider. It too calls a
> 'sub_domain->map' function, which is presumably the one I would need to
> set correctly. Alternatively, it is easiest just to enforce the vertices
> to align on opposing faces at the mesh generation stage?
>
>
> cheers
> pete
>
> On Jun 20, 2012, at 2:55 PM, Johan Hake wrote:
>
>> Your question #200584 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/200584
>>
>> Status: Open => Answered
>>
>> Johan Hake proposed the following answer:
>> I am not expert on the inner workings of PeriodicBC. It looks to me that
>> the coordinate should be set correctly. Not sure that it does not work
>> because you only use map 1 dof?
>>
>> You might want to try to follow the logic in PeriodicBC.cpp. Not sure it
>> gets you any wiser. I bailed out after a minute...
>>
>> Johan
>>
>> On 06/16/2012 03:31 AM, Peter Kekenes-Huskey wrote:
>>> New question #200584 on DOLFIN:
>>> https://answers.launchpad.net/dolfin/+question/200584
>>>
>>> Dear all,
>>> I would like to apply a periodic boundary condition between specific nodes on a mesh. For instance, in the included example (based off a script in the demo directory), I would like to apply a PBC to x1=[0,0.5] and x2=[1,0.5] (along the surface). All other surface vertices are assigned a dirichlet condition. I’m able to find the vertices using the ‘inside()’ functions I’ve defined, but I don’t know how to tell the declare the mapping between x1 and x2 using ‘map()’ in the PBC function.
>>>
>>> I’m sure there are simpler ways of doing this for a grid with aligned surface vertices, but I’d like a general solution since my meshes are very unlikely to unaligned. In practice, I might want to select a small region of vertices near x1 and x2 for this condition.
>>>
>>> I appreciate any help on this problem!
>>>
>>> Best,
>>> Pete
>>>
>>>
>>> from dolfin import *
>>>
>>> # Periodic BCs don't work with Epetra
>>> if parameters["linear_algebra_backend"] == "Epetra":
>>> print "Sorry, this demo does not work with the Epetra backend"
>>> import sys
>>> sys.exit(0)
>>>
>>> # Create mesh and finite element
>>> mesh = UnitSquare(2,2)
>>> for i,c in enumerate( mesh.coordinates() ):
>>> print c
>>>
>>> import numpy as np
>>> coord1 = np.array([0.0,0.5])
>>> coord2 = np.array([1.0,0.5])
>>>
>>>
>>> Vv = FunctionSpace(mesh, "CG", 1)
>>>
>>> # Source term
>>> class Source(Expression):
>>> def eval(self, values, x):
>>> dx = x[0] - 0.5
>>> dy = x[1] - 0.5
>>> values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)
>>>
>>> # Sub domain for Dirichlet boundary condition
>>> # x free
>>> class DirichletBoundary0(SubDomain):
>>> def inside(self, x, on_boundary):
>>> return bool( \
>>> np.linalg.norm( x- coord1)> DOLFIN_EPS and \
>>> np.linalg.norm( x- coord2)> DOLFIN_EPS and \
>>> on_boundary)
>>>
>>>
>>> # Sub domain for Periodic boundary condition
>>> class PeriodicBoundary0(SubDomain):
>>>
>>> def inside(self, x, on_boundary):
>>> result = bool( \
>>> (np.linalg.norm( x- coord1)< DOLFIN_EPS or \
>>> np.linalg.norm( x- coord2)< DOLFIN_EPS) and \
>>> on_boundary)
>>> print result
>>> return result
>>>
>>> # field component 1
>>> def map(self, x, y):
>>> #print "x y"
>>> print x
>>> y[0] = x[0] - 1.0
>>> y[1] = x[1]
>>>
>>> # Create Dirichlet boundary condition
>>> u0 = Constant(1.0)
>>>
>>> dbc0 = DirichletBoundary0()
>>> bc0 = DirichletBC(Vv, u0, dbc0)
>>>
>>> pbc0 = PeriodicBoundary0()
>>> bc1 = PeriodicBC(Vv, pbc0)
>>>
>>> bcs = [bc0, bc1]
>>>
>>> # Define variational problem
>>> u = TrialFunction(Vv)
>>> v = TestFunction(Vv)
>>> f = Source()
>>> a = dot(grad(u), grad(v))*dx
>>> L = f*v*ds
>>>
>>> # Compute solution
>>> u = Function(Vv)
>>> solve(a == L, u, bcs)
>>>
>>> # Save solution to file
>>> file = File("periodic.pvd")
>>> file<< u
>>>
>>>
>>
>> --
>> If this answers your question, please go to the following page to let us
>> know that it is solved:
>> https://answers.launchpad.net/dolfin/+question/200584/+confirm?answer_id=0
>>
>> If you still need help, you can reply to this email or go to the
>> following page to enter your feedback:
>> https://answers.launchpad.net/dolfin/+question/200584
>>
>> You received this question notification because you asked the question.
>

Can you help with this problem?

Provide an answer of your own, or ask Peter Kekenes-Huskey for more information if necessary.

To post a message you must log in.