Quad and Hex implementation

ART: I'm going to use Morton ordering for all the hex and quad cells, since that's what p4est uses and it seems to work. For now I will just assume all quads are planar since the higher order mesh routines were removed.

MER: How does the Morton ordering match the UFC numbering convention?

ART: Morton ordering is for quads/hexes. UFC numbering only deals with simplices and won't really work for quads (please prove me wrong (for example what side is opposite the lowest vertex)). I haven't found an extension of Morton for simplicies. The biggest question is how to insure an ordering that doesn't require flipping the normals between elements. I'm guessing we will just have to bite the bullet on this but the jury is stil out.

AL: We always knew that quads and hexes would be different when we wrote the UFC numbering but decided to put it in there anyway. It can be changed to whatever is most practical. Is there a numbering for Morton numbering? What is 'prest'?

ART: Okay I'll go re-read the UFC spec, I was only familiar with the simplices ordering. I meant p4est not prest.

ART: Okay so UFC doesn't number edges and uses the right hand rule. I asked Carsten and some others on where to find the different advantages but the only thing I received was that Morton ordering is easier to implement and extend to 3D and beyond. I need to go see how SyFi uses quads and hexes, if they are extensively using the order then I won't change it, otherwise it makes sense to change it.

MSA: I talked to Kent, feel free to ignore what SyFi does about the quad numbering. Do what seems best going forward.

ART: Okay I'll just implement the Morton ordering for now. It seems that the ordering only comes up in a few places so hopefully it will be easy to switch.

ART: It seems that the mesh editor interface will have to change. Right now it take the length of arguments and determines the shape of the cell. I will add a cell type default field with defaulting to current setting. This way we don't break the api. Let me know if there are other opinions on how to change the interface.

ART: Any suggestions for changing
    void MeshEditor::open(Mesh& mesh, uint tdim, uint gdim)
    I don't see a clean way to keep it, but it's used in lots of documentation in the code.
    Also is there an example of a mesh with tdim less than gdim?

GNW: What's wrong with MeshEditor::open? From the tdim and the number of vertices, the cell shape can be inferred. DOLFIN does support tdim < gdim. There may not be a demo, but it has been done before.

[MER: Sidenote: For an example of tdim < gdim, see for instance https://lists.launchpad.net/fenics/msg00715.html. The mesh code actually still runs!]

ART: There's nothing wrong with MeshEditor::open(Mesh& mesh, CellType::Type type, uint tdim, uint gdim); but open(Mesh& mesh, uint tdim, uint gdim) is ambiguous since we will have two CellType::Type's for tdim == 2 or tdim == 3. The open(Mesh& mesh, uint tdim, uint gdim) does a switch on tdim.

ART: @MER This should definitely be supported since that's the common way of representing a surface mesh, although the Mobius strip example is more of a curiosity.

CJC: We are definitely interested in tdim < gdim because we are doing shallow-water simulations on the sphere, so the simulation is performed on a 2D mesh embedded in 3D. This can't be supported in DOLFIN yet (needs Piola transformation from 2D to 3D) but there are plans afoot to make it so.

ART: @CJC do you need Piola because you want to target RT elements? I would think that scalar elements would remain affine. Also if you are doing shallow water are you planning on just doing a 2D mesh with a projection in h? I know that this is what AD CIRC and/or UTBest do, but I haven't seen it supported in FEniCS. I was hoping to avoid any nonlinear transformations of the mesh itself as this adds quite a bit of (algorithmic) complexity. deal.II handles any non-linear mapped quad and it eats up a lot of time in the assembly loop.

CJC: @ART yes, RT elements for 2D shallow water, and we need RT *and* edge elements eventually in 3D. For our shallow-water code we really do have 2D elements in 3D, the Piola transformation guarantees that velocity fields are tangential to the mesh elements. Non-affine is crucial as you can't tile the sphere with affine quadrilaterals. We are planning on doing something with nonlinear transformations of triangles as well. For the equation systems that we are using, some beautiful differential geometry comes in and cancellation of the geometric terms (J and det J) occurs almost all of the terms (except the mass matrix). If you aren't planning to support non-affine quads, I suppose the plan would be for us to look at adding that capability once we have worked out how to add it for triangles.

ART: @CJC. For now I think we have to stick to getting parity with the current triangles and tets, then add the Piola transforms for the mesh. It's not off the table just not the first thing to get done. We can see how things go towards the end of the sprint.

CJC: @ART. Absolutely, we don't even yet have a blueprint for the non-affine transformations on triangles yet!

AL: The MeshEditor interface definitely needs to be changed. MeshEditor::open is one of the places in DOLFIN where we have made the assumption that only simplices are used. Adding a CellType argument would be the natural way to handle it. We could keep the old one for backwards compatibility, possibly with a deprecation warning, and it would call the new MeshEditor::open function with the proper arguments.

Work items:
Add quadrilateral cell type and hook up to mesh.h: TODO
Add Rectangle, UnitQuadrangle, and UnitSquare quadrilateral support: TODO
Create integrated test (poisson) with Syfi generated assembly: TODO
Add hexahedral cell type and hook up to mesh.h: TODO
Add Box, UnitHexahedral, and UnitCube hexahedral support: TODO
Create integrated test (poisson) with Syfi generated assembly: TODO