Thank you. This is the .geo file that I am using is copied below. First I tried splines and that seemed to have a line that goes outside of the domain whereas I have now switched to lines. It's a very, very simple set of coastlines but if this works I can do more complicated coastlines. Unfortunately, when I try to compile it I do get a warning about orienting normals. The .msh file looks okay but clearly there is something not quite right. But it still does not important in firedrake though. fpoulin@domlt32:~/Research/2016-Christine/OceanBasins$ gmsh -2 test3.geo Info : Running 'gmsh -2 test3.geo' [Gmsh 2.10.1, 1 node, max. 1 thread] Info : Started on Thu Nov 2 15:49:09 2017 Info : Reading 'test3.geo'... Info : Done reading 'test3.geo' Info : Meshing 1D... Info : Meshing curve 1 (Line) Info : Meshing curve 2 (Line) Info : Done meshing 1D (0.004 s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Delaunay) Info : Done meshing 2D (0.00479794 s) Warning : Could not orient normal of surface 1 Info : 110 vertices 216 elements Warning : ------------------------------ Warning : Mesh generation error summary Warning : 1 warning Warning : 0 errors Warning : Check the full log for details Warning : ------------------------------ Info : Writing 'test3.msh'... Info : Done writing 'test3.msh' Info : Stopped on Thu Nov 2 15:49:09 2017 element_size = 500000; // Printing points from 1 to 6. Point(1) = {-0.000000, -0.000000, 0, element_size}; Point(2) = {4275156.460308, 542789.171194, 0, element_size}; Point(3) = {2980231.171904, 4007988.825894, 0, element_size}; Point(4) = {507429.795284, 3586291.304193, 0, element_size}; Point(5) = {-1702226.201816, 1330470.794133, 0, element_size}; Point(6) = {-0.000000, -0.000000, 0, element_size}; // Printing Spline points 1 to 6. Line(1) = {1,2,3}; Line(2) = {3,4,5,6,1}; // Printing line loop from spline 1 to 2 Line Loops(1) = {1,2}; Plane Surface(1) = {1}; ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Homolya, Miklós <m.homolya14@imperial.ac.uk> Sent: Thursday, November 2, 2017 3:42:08 PM To: firedrake Subject: Re: [firedrake] problems with importing a mesh It says that your generated mesh has some vertices, edges, or else that does not belong to any cell. I suggest you have a look at the generated mesh to see if you can spot such an artifact, and then figure out how to drive gmsh to not have such a thing. ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 02 November 2017 19:38:34 To: firedrake Subject: [firedrake] problems with importing a mesh Hello, I have created a mesh of an oceanic basin that I have triangulated with gmsh. Finally! I am very excited to solve for the wind-driven gyre solutions with this new mesh (as I discussed back in March at the workshop), unfortunately, when I try to import the mesh I get the following: mesh = Mesh("test3.msh") ValueError: Provided mesh has some entities not reachable by traversing cells (maybe rogue vertices?) Any ideas as to what I might have done wrong and how I can fix this? Cheers, Francis ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

Thanks for the help Miklos. I am not sure if I understand this completely but I have done a couple of things. One problem was that I had the first and last point as the same but I have since removed the last point. 1) For the case of line this seemed to do the trick and I can now import this into firedrake, which I am happy about. 2) For splines, I have also removed the point (updated code below) and when I look at it using gmsh I see that some points aren't actually on the boundary. This does seem odd but I am not sure how to fix it. Should I use a different kind of spline? Sorry for not better understanding gmsh. element_size = 1000000; // Printing points from 1 to 5. Point(1) = {-0.000000, -0.000000, 0, element_size}; Point(2) = {4275156.460308, 542789.171194, 0, element_size}; Point(3) = {2980231.171904, 4007988.825894, 0, element_size}; Point(4) = {507429.795284, 3586291.304193, 0, element_size}; Point(5) = {-1702226.201816, 1330470.794133, 0, element_size}; // Printing Spline points 1 to 5. Spline(1) = {1,2,3}; Spline(2) = {3,4,5,1}; // Printing line loop from spline 1 to 2 Line Loops(1) = {1,2}; Plane Surface(1) = {1}; ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Homolya, Miklós <m.homolya14@imperial.ac.uk> Sent: Thursday, November 2, 2017 4:00:15 PM To: firedrake Subject: Re: [firedrake] problems with importing a mesh So it looks like some of the points -- e.g., Point 4, 5, and 6 -- are only used as support points to the spline, they don't belong to any cell, yet they are generated into the mesh. ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 02 November 2017 19:50:24 To: firedrake Subject: Re: [firedrake] problems with importing a mesh Thank you. This is the .geo file that I am using is copied below. First I tried splines and that seemed to have a line that goes outside of the domain whereas I have now switched to lines. It's a very, very simple set of coastlines but if this works I can do more complicated coastlines. Unfortunately, when I try to compile it I do get a warning about orienting normals. The .msh file looks okay but clearly there is something not quite right. But it still does not important in firedrake though. fpoulin@domlt32:~/Research/2016-Christine/OceanBasins$ gmsh -2 test3.geo Info : Running 'gmsh -2 test3.geo' [Gmsh 2.10.1, 1 node, max. 1 thread] Info : Started on Thu Nov 2 15:49:09 2017 Info : Reading 'test3.geo'... Info : Done reading 'test3.geo' Info : Meshing 1D... Info : Meshing curve 1 (Line) Info : Meshing curve 2 (Line) Info : Done meshing 1D (0.004 s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Delaunay) Info : Done meshing 2D (0.00479794 s) Warning : Could not orient normal of surface 1 Info : 110 vertices 216 elements Warning : ------------------------------ Warning : Mesh generation error summary Warning : 1 warning Warning : 0 errors Warning : Check the full log for details Warning : ------------------------------ Info : Writing 'test3.msh'... Info : Done writing 'test3.msh' Info : Stopped on Thu Nov 2 15:49:09 2017 element_size = 500000; // Printing points from 1 to 6. Point(1) = {-0.000000, -0.000000, 0, element_size}; Point(2) = {4275156.460308, 542789.171194, 0, element_size}; Point(3) = {2980231.171904, 4007988.825894, 0, element_size}; Point(4) = {507429.795284, 3586291.304193, 0, element_size}; Point(5) = {-1702226.201816, 1330470.794133, 0, element_size}; Point(6) = {-0.000000, -0.000000, 0, element_size}; // Printing Spline points 1 to 6. Line(1) = {1,2,3}; Line(2) = {3,4,5,6,1}; // Printing line loop from spline 1 to 2 Line Loops(1) = {1,2}; Plane Surface(1) = {1}; ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Homolya, Miklós <m.homolya14@imperial.ac.uk> Sent: Thursday, November 2, 2017 3:42:08 PM To: firedrake Subject: Re: [firedrake] problems with importing a mesh It says that your generated mesh has some vertices, edges, or else that does not belong to any cell. I suggest you have a look at the generated mesh to see if you can spot such an artifact, and then figure out how to drive gmsh to not have such a thing. ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 02 November 2017 19:38:34 To: firedrake Subject: [firedrake] problems with importing a mesh Hello, I have created a mesh of an oceanic basin that I have triangulated with gmsh. Finally! I am very excited to solve for the wind-driven gyre solutions with this new mesh (as I discussed back in March at the workshop), unfortunately, when I try to import the mesh I get the following: mesh = Mesh("test3.msh") ValueError: Provided mesh has some entities not reachable by traversing cells (maybe rogue vertices?) Any ideas as to what I might have done wrong and how I can fix this? Cheers, Francis ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

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}; ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of G. D. McBain <gdmcbain@protonmail.com> Sent: Thursday, November 2, 2017 8:07:51 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] problems with importing a mesh 2) For splines, I have also removed the point (updated code below) and when I look at it using gmsh I see that some points aren't actually on the boundary. This does seem odd but I am not sure how to fix it. Should I use a different kind of spline? The problem may be this: the Gmsh concept of Physical versus Elementary entities. http://gmsh.info/doc/texinfo/gmsh.html#Elementary-vs-physical-entities Your Gmsh file has only Elementary entities (Point, Line/Spline, Line Loop, Plane Surface). The trick is to mark the domain (for this two-dimensional problem) as a Physical Surface and each part of the boundary as a Physical Line. In three dimensions, those would be Physical Volume and Physical Surface, respectively. When a Gmsh file has no Physical entities, Gmsh puts everything into the mesh, including construction points (centres of arcs, control-points of splines), which often end up as disconnected, but if there is a least one Physical entity, only those finite elements belonging to Physical entities are included in the output .msh. A quirk of Firedrake, maybe inherited from PETSc which I think it uses to read Gmsh meshes, is that every part of the boundary has to be declared as part of a Physical Line/Surface, even if it's only going to get homogeneous natural boundary conditions in the finite element problem. Hope this helps; I can supply a little example, if the above is unclear.

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 ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: G. D. McBain <gdmcbain@protonmail.com> Sent: Thursday, November 2, 2017 8:53:47 PM To: Francis Poulin Cc: firedrake@imperial.ac.uk Subject: Re: [firedrake] problems with importing a mesh 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}; Great! You're welcome. So you didn't need the Physical Line… I wonder if my warning about those being required all round were wrong or outdated… Let me check.

Thanks for all the helpful information. I am starting to understand but not quiet there yet. I was able to generate the mesh but when I try running your script I get the following problem. Is there a different syntax we should be using? (firedrake) fpoulin@domlt32:~/Dropbox/Research_Christine/Notebooks/TestGeometry$ python membrane.py quadrant.msh File "membrane.py", line 11 def main(meshfile: str) -> None: ^ SyntaxError: invalid syntax What I understand is I should be specifying a label to the boundary and use this label to impose the BC's. I tried importing your geometry into my notebook and that seemed to go okay but when I try plotting I get the Assertion error. Any more help would be greatly appreciated. ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: G. D. McBain <gdmcbain@protonmail.com> Sent: Friday, November 3, 2017 2:09 AM To: Francis Poulin Cc: firedrake@imperial.ac.uk Subject: Re: [firedrake] problems with importing a mesh I am curious to see what you find does the trick. It seems that my warning was obsolete. Whereas, I find from my old notes, that failure label each part of the boundary in Gmsh raised: AssertionError: Every marker has to be contained in unique_markers and in Feb. 2016, L. M. wrote: ‘If you provide a marker for one of the external boundaries of your mesh, you need to do so for all of them.’ https://mailman.ic.ac.uk/mailman/htdig/firedrake/2016-February/001809.html I find today that this isn't the case. I hadn't noticed this being fixed as I'd simply complied by labelling everything. Anyway thank you to whoever fixed that, in Firedrake or upstream. 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. So to satisfy myself of the above, I prepared a simple example with a Dirichlet condition on part of the boundary and natural conditions on the rest. The same technique might work for your problem. Consider the plane Poisson equation with constant unit forcing on a circle with homogeneous Dirichlet conditions, but exploiting two lines of symmetry using just a quadrant with natural conditions on the two radii. The exact answer on the circle is that the integral of the solution divided by the square of the area should be 1/8 pi, or four times on the quarter-model, so 1/2pi. %<---quadrant.geo Point(1) = {1, 0, 0, 1.0}; Point(2) = {0, 0, 0, 1.0}; Point(3) = {0, 1, 0, 1.0}; Circle(1) = {1, 2, 3}; Line(2) = {3, 2}; Line(3) = {2, 1}; Line Loop(1) = {1, 2, 3}; Plane Surface(1) = {1}; Physical Surface("membrane", 2) = {1}; Physical Line("rim", 1) = {1}; --->% Note that the last line specifies the quarter-circle perimeter from "Circle(1)" but not the radii of symmetry "Line(2)" or "Line(3)". The Dirichlet boundary condition is imposed on the perimeter of the circle using "[1]" for the 1 of the Physical Line: bc = DirichletBC(V, 0.0, [1]) Running this with $ gmsh -2 -clscale 6.25e-2 quadrant.geo $ python membrane.py quadrant.msh I get 0.15905285125 which is 1/2 pi to four decimal places.

Sorry, I thought I had upgraded to the latest version of firedrake but I had not. I have since fixed that and am now using the version of firedrake based on python3, which is using 3.5.2. I am happy to say that when I try using this example in my script it seems to work very nicely. It does inf act solve the weak form and imposes the correct BCS. Now for my oceanic basin case, below are the last two lines in my code. From this I gather that the name for the physical surface is [2], and that consists of the entire domain. Plane Surface(1) = {1}; Physical Surface(2) = {1}; If I want to impose Dirichlet BCs everywhere then should I do the following? #bc = DirichletBC(Z.sub(0), 0.0, "on_boundary") bc = DirichletBC(Z.sub(0), 0.0, [2]) I commented out what I had before, which I thought would impose the BCs everywhere but it didn't seem to do that. ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Homolya, Miklós <m.homolya14@imperial.ac.uk> Sent: Friday, November 3, 2017 11:39:23 AM To: G. D. McBain; firedrake Subject: Re: [firedrake] problems with importing a mesh These look like the optional type annotations. I think they were introduced in Python 3.6, while Firedrake requires Python 3.5+ You are probably running Python 3.5.x then. Just remove those colon types, or alternatively try Python 3.6. ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 03 November 2017 15:11:34 To: G. D. McBain Cc: firedrake Subject: Re: [firedrake] problems with importing a mesh Thanks for all the helpful information. I am starting to understand but not quiet there yet. I was able to generate the mesh but when I try running your script I get the following problem. Is there a different syntax we should be using? (firedrake) fpoulin@domlt32:~/Dropbox/Research_Christine/Notebooks/TestGeometry$ python membrane.py quadrant.msh File "membrane.py", line 11 def main(meshfile: str) -> None: ^ SyntaxError: invalid syntax What I understand is I should be specifying a label to the boundary and use this label to impose the BC's. I tried importing your geometry into my notebook and that seemed to go okay but when I try plotting I get the Assertion error. Any more help would be greatly appreciated. ------------------ Francis Poulin Associate Dean, Undergraduate Studies Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: G. D. McBain <gdmcbain@protonmail.com> Sent: Friday, November 3, 2017 2:09 AM To: Francis Poulin Cc: firedrake@imperial.ac.uk Subject: Re: [firedrake] problems with importing a mesh I am curious to see what you find does the trick. It seems that my warning was obsolete. Whereas, I find from my old notes, that failure label each part of the boundary in Gmsh raised: AssertionError: Every marker has to be contained in unique_markers and in Feb. 2016, L. M. wrote: ‘If you provide a marker for one of the external boundaries of your mesh, you need to do so for all of them.’ https://mailman.ic.ac.uk/mailman/htdig/firedrake/2016-February/001809.html I find today that this isn't the case. I hadn't noticed this being fixed as I'd simply complied by labelling everything. Anyway thank you to whoever fixed that, in Firedrake or upstream. 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. So to satisfy myself of the above, I prepared a simple example with a Dirichlet condition on part of the boundary and natural conditions on the rest. The same technique might work for your problem. Consider the plane Poisson equation with constant unit forcing on a circle with homogeneous Dirichlet conditions, but exploiting two lines of symmetry using just a quadrant with natural conditions on the two radii. The exact answer on the circle is that the integral of the solution divided by the square of the area should be 1/8 pi, or four times on the quarter-model, so 1/2pi. %<---quadrant.geo Point(1) = {1, 0, 0, 1.0}; Point(2) = {0, 0, 0, 1.0}; Point(3) = {0, 1, 0, 1.0}; Circle(1) = {1, 2, 3}; Line(2) = {3, 2}; Line(3) = {2, 1}; Line Loop(1) = {1, 2, 3}; Plane Surface(1) = {1}; Physical Surface("membrane", 2) = {1}; Physical Line("rim", 1) = {1}; --->% Note that the last line specifies the quarter-circle perimeter from "Circle(1)" but not the radii of symmetry "Line(2)" or "Line(3)". The Dirichlet boundary condition is imposed on the perimeter of the circle using "[1]" for the 1 of the Physical Line: bc = DirichletBC(V, 0.0, [1]) Running this with $ gmsh -2 -clscale 6.25e-2 quadrant.geo $ python membrane.py quadrant.msh I get 0.15905285125 which is 1/2 pi to four decimal places.

Managed I managed to figure it out. Below is an example where I have a domain with 5 boundaries and I am trying to set the Physical Line like someone suggested. Does this look right?

If yes, then I think I'm in good shape. One problem is that the boundary layer is very narrow for my choice of parameters so it's hard to know whether the boundary value is really zero. I don't suppose there is an easy way to find out the maximum of the boundary value in the numerical solution?

element_size = 0.05; // Printing points from 1 to 6. Point(1) = {-0.000000, -0.000000, 0, element_size}; Point(2) = {4.27, 0.54, 0, element_size}; Point(3) = {2.98, 4.00, 0, element_size}; Point(4) = {0.50, 3.58, 0, element_size}; Point(5) = {-1.70, 1.33, 0, element_size};

// Printing Spline points 1 to 6. Line(1) = {1,2}; Line(2) = {2,3}; Line(3) = {3,4}; Line(4) = {4,5}; Line(5) = {5,1};

// Printing line loop from spline 1 to 2 Line Loops(1) = {1,2,3,4,5};

Plane Surface(6) = {1}; Physical Surface(7) = {6}; Physical Line(8) = {1};

Hello. Yes, sorry about the PEP 3107 function annotations (I fell for these immediately they were released) https://www.python.org/dev/peps/pep-3107/ and yes, it's Physical Line that's wanted for boundary conditions in two-dimensional problems. However, I'm sorry I don't know any answer to your question about how to find extrema on boundaries in Firedrake, but I would like to know that too. I do know a couple of tricks for extracting that information in FreeFem++: * the dirty trick of evaluating a line integral over the boundary patch with the integrand being an impure function that remembers the highest value it has seen and returns an irrelevant value, especially used with a quadrature rule with abscissæ at the nodes; * a technical trick using the machinery for imposing Dirichlet conditions by penalization http://www.um.es/freefem/ff++/pmwiki.php?n=Main.ExtractTheDofsOnAPieceOfBoun... The latter might be a bit specific to the implementation to port to Firedrake, the former I'd rather avoid if there were a proper way to do it in Firedrake, if anyone knows how, but I can't help you myself just yet. If you know where the boundary is geometrically and can independently generate a list of coordinate-points on it, you can use the .at method of point-evaluation. [http://www.firedrakeproject.org/point-evaluation.html](http://www.firedrakeproject.org/point-evaluation.html?highlight=values) I guess the operation of finding extrema of finite element functions on the mesh is not really in the spirit of the finite element method, which is why it's awkward. The method is more about evaluating functionals. The FEniCS Q&A has the same query and an equally awkward procedure. That one might be portable to Firedrake. https://fenicsproject.org/qa/3074/get-function-values-on-boundaries Kind regards, Geordie McBain

On 08/11/17 09:13, Lawrence Mitchell wrote:

On 03/11/17 20:15, Francis Poulin wrote:

...

If yes, then I think I'm in good shape.  One problem is that the boundary layer is very narrow for my choice of parameters so it's hard to know whether the boundary value is really zero.  I don't suppose there is an easy way to find out the maximum of the boundary value in the numerical solution?

If you are imposing strong (Dirichlet) conditions, then the solution on the boundary nodes is, by definition, whatever you set it to.

If you have some other means of enforcing boundary conditions and just want to know that the values are on the boundary, you could do:

V = FunctionSpace(...)

selector = DirichletBC(V, 0, "on_boundary")

f = Function(V)

# solve into f

boundary_values = f.dat.data_ro[selector.nodes, ...]

Note this assumes that inspecting the point values of your space makes sense: this is fine for CG and DG, not so for other spaces. Lawrence