I am curious to see what you find does the trick.

As for my coastline problem, the good news is that I can now import the geometry.  The bad news is  that when I solve the weak for of the PDE it does not impose the Dirichlet BCs.  

At the moment I am trying

bc = DirichletBC(z.sub(0), 0.0, "on_boundary")

I remember having a similar problem a while back and I needed to specify the name of the boundary and use that to impose the condition.

What is the current practice for meshes that we import?

Cheers, Francis

Yes, that did the trick.  I used your idea in many different tests that I have done and basically added the second of the following two lines and that seemed to do the trick.

Thank you!

...

Plane Surface(1) = {1};

Physical Surface(2) = {1};