esys-escript

Periodical Boundary Condition

Asked by Bastien Grandet

Hi,

I'm trying to use periodical BCs.

In the simple script below, are those PBCs defined the way they should?

This script is a slightly modified version of linearElastic.py I found with the source code. This script reproduces a simple shear problem.

Cheers,
Bastien Grandet

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys import finley
from esys.weipa import saveVTK

press0=1.
lamb=1.
nu=0.3

# this sets the hook tensor:
def setHookTensor(w,l,n):
   C=Tensor4(0.,w)
   for i in range(w.getDim()):
     for j in range(w.getDim()):
          C[i,i,j,j]+=l
          C[j,i,j,i]+=n
          C[j,i,i,j]+=n

   return C

# generate mesh: here 10x20 mesh of order 2
domain=finley.Rectangle(10,10,1,l0=1.0,l1=1.0,periodic0=True)

# get handel to nodes and elements:
e=Function(domain)
fe=FunctionOnBoundary(domain)
n=ContinuousFunction(domain)
v=Vector(0. , ContinuousFunction(domain))
#
# set a mask msk of type vector which is one for nodes and components set be a constraint:
#
msk=whereZero(n.getX()[1])*[1.,1.]+whereZero(1-n.getX()[1])*[1.,1.] #whereZero(n.getX()[0])*[0.,1.]+
#
#Values of the displacement BC
#
v=v+whereZero(1-n.getX()[1])*[0.01,0.]

# assemble the linear system:
mypde=LinearPDE(domain)
mypde.setValue(A=setHookTensor(e,lamb,nu),q=msk,r=v)
mypde.setSymmetryOn()
mypde.getSolverOptions().setVerbosityOn()
mypde.getSolverOptions().setPreconditioner(mypde.getSolverOptions().AMG)

# solve for the displacements:
u_d=mypde.getSolution()

# get the gradient and calculate the stress:
g=grad(u_d)
stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g))

#Update geometry
domain.setX(domain.getX()+u_d)

# write the hydrostatic pressure:
saveVTK("result.vtu",displacement=u_d,pressure=trace(stress)/domain.getDim(),boundary=msk)

Question information

Language:
English Edit question
Status:
Answered
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Lutz Gross (l-gross) said : 		#1

For periodic boundary conditions interpolation of Data with the ContinuousFunction() attribute to Data with the Solution() attribute produces undefined results for data points where values are shared due to periodic boundary conditions (in your example these are the points x with x[0]=0 and x[0]=1). If Data with the ContinuousFunction() are not periodic there are two different values available for a pair of points with identical values due to periodic boundary conditions. In this case the interpolation process needs to enforce periodicity and makes an (unpredictable) choice between the two available values.

In your program the statement "domain.getX()+u_d" to set new coordinates for the domain causes the problem.
The automatic interpolation within binary operations performs an interpolation of the node coordinates (with ContinuousFunction() attribute) to periodic Data as outlined above.

To avoid the problem manual interpolation should be inserted:

    domain.setX( domain.getX()+interpolate(u_d, ContinuousFunction(domain)) )

(or change the order of arguments in the addition). We will consider a revise of the interpolation policy for this case to avoid the problem.

Revision history for this message
Lutz Gross (l-gross) said : 		#2

see https://mantis.esscc.uq.edu.au/view.php?id=644

Can you help with this problem?

Provide an answer of your own, or ask Bastien Grandet for more information if necessary.

To post a message you must log in.