Forwarding a question that was erroneously bumped from the list.
On Jan 12, 2011, at 4:34 PM, Changsoo Kim wrote:
> Is there an option to assign a periodic boundary condition on a mesh created
> from Gmsh? I want to impose a periodic BC along x-direction (only), I know we
> can use "PeriodicGrid2DLeftRight", but not on the mesh from Gmsh. Any
> thoughts?
In principle, you can call mesh._connectFaces(faces0, faces1) on any mesh.
faces0 and faces1 are lists of face IDs. The problem is getting those lists of
IDs and making sure that they're in the right order.
To get the list of IDs, you can do things like
X, Y = mesh.getFaceCenters()
facesLeft = mesh.getExteriorFaces() & numerix.isclose(X, Xleft, atol=1e-10)
facesRight = mesh.getExteriorFaces() & numerix.isclose(X, Xright, atol=1e-10)
faces0 = numerix.nonzero(facesLeft)[0]
faces1 = numerix.nonzero(facesRight)[0]
There is no guarantee, though, that faces0 are lined up with faces1, or even
that the faces are the same size (they must match up exactly).
You might be able to do something with sorting the IDs on the corresponding Y
values, obtained with Y[faces0]
but you still need to be sure that gmsh is generating faces that can match.
Here's a script that illustrates it. This requires the tip of trunk (r4101)
because I had to relax a check on matching orientations.
The banding that appears is a viewer artifact.
from fipy import *
# make triangulated mesh in a square domain
mesh = Gmsh2D('''
cellSize = 0.2;
Point(1) = {0, 0, 0, cellSize};
Point(2) = {10., 0, 0, cellSize};
Point(3) = {10., 10., 0, cellSize};
Point(4) = {0., 10., 0, cellSize};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
''' % locals())
x, y = mesh.getCellCenters()
X, Y = mesh.getFaceCenters()
# get the left and right faces
facesLeft = mesh.getExteriorFaces() & numerix.isclose(X, 0., atol=1e-10)
facesRight = mesh.getExteriorFaces() & numerix.isclose(X, 10., atol=1e-10)
# get the IDs of the faces
faces0 = numerix.nonzero(facesLeft)[0]
faces1 = numerix.nonzero(facesRight)[0]
# sort the IDs of the faces in ascending Y
faces0 = faces0[argsort(Y[faces0])]
faces1 = faces1[argsort(Y[faces1])]
# connect the faces
mesh._connectFaces(faces0=faces0, faces1=faces1)
var = CellVariable(mesh=mesh)
# make a circular seed
var.setValue(((x - 5.)**2 + (y - 5.)**2 <= 2**2) * 1.)
viewer = Viewer(vars=var, datamin=0., datamax=1.)
viewer.plot()
# calculate the maximum timestep
dx = numerix.sqrt(mesh.getCellVolumes().min())
speed = 1.0
cfl = 0.5
dt = cfl * dx / speed
# advect the circle
eq = TransientTerm() + VanLeerConvectionTerm(coeff = [[speed], [0.]])
for step in range(1000.):
eq.solve(var, dt=dt)
viewer.plot()
