Surface subdomain fails to mark CellFunction on a 3D mesh

Question

Hi everyone!

I'm currently trying to impose periodic BC on a 3D mesh. Basically I want the solutions of my PDE to be equal on opposite normal-to-midplan faces of a 3D plate.

The mesh is defined with GMSH, with extra attention to make the corresponding boundary meshes to be matching (enforced in GMSH through Periodic Surface command).

Basically I import the mesh (shown in figure) and define my periodic BC so as to identify the bottom and left sides in the inside procedure:

mesh = Mesh('Omega3D_GMSH.xml');
class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
    return bool((near(x[0],-0.0075) or near(x[1],-0.0075))\
                 and not ((near(x[0],0.0075) and near(x[1],-0.0075))\
                 or (near(x[0],-0.0075) and near(x[1],0.0075))) and on_boundary)

def map(self, x, y):
    if near(x[0],00075) and near(x[1],0.0075):
        y[0] = x[0] - 0.015
        y[1] = x[1] - 0.015
        y[2] = x[2]
    elif near(x[0],00075): 
        y[0] = x[0] - 0.015
        y[1] = x[1]
        y[2] = x[2]
    elif near(x[1],0.0075):
        y[0] = x[0]
        y[1] = x[1] - 0.015
        y[2] = x[2]
    else:
        y[0]=1.0
        y[1]=1.0
        y[2]=1.0

Then I create a CellFunction on a BoundaryMesh in order to mark the interior hole:

bmesh = BoundaryMesh(mesh, "exterior");
mapping = bmesh.entity_map(2);
part_of_boundary = CellFunction("size_t", bmesh, 0);
for cell in cells(bmesh):
    if ((abs(Facet(mesh, mapping[cell.index()]).normal().x()) > 0) or \
       (abs(Facet(mesh, mapping[cell.index()]).normal().y()) > 0)) and \
       (abs(Facet(mesh, mapping[cell.index()]).distance(Point(0.0,0.0,0.0))) < 0.0025):
           part_of_boundary[cell]=1

PeriodicBoundary().mark(part_of_boundary,2);

ds = ds[part_of_boundary];

plot(part_of_boundary);

And to be sure that my Periodic BC is set as I want I marked the CellFunction on the corresponding part and plotted the resulting Subdomain map, but it doesn't work at all, and the solution is indeed non periodic.

As you can see the bottom and left sides are not marked through the command and I can't figure out why. Thank you for your answers,

Regards,

MFV.

Answer 1

The condition in inside is wrong. If you change it with return True you'll see that the marking works correctly, so it looks like the conditional never evaluates to True. You can verify this by doing, in inside,

b = bool(...)
if b:
    print "working!"
return b

If I'm right, you should never see the line printed.

Comments

Thank you for your answer, indeed the condition is wrong and after further testing it looks like it is the on_boundary statement that makes it incorrect. In fact it seems that the on_boundary method simply doesn't work on the mesh.

Deleting the on_boundary condition in the inside function of the PeriodicBoundary() SubDomain class makes it work correctly.

Thanks again for your help,

Regards,

MFV