[spec] Fluidity Boundary Conditions

Branch

Fluidity mesh objects contain face objects
which in turn contain a list of boundary_ids
boundary_ids are (in FEniCS terms) a MeshFunction over facets, associating each facet with exactly one boundary ID.
For interior facets and for exterior facets to which no other ID has been associated, this ID is 0
(Boundary ID is always an integer)
Now at Field level, there are boundary condition objects
these associate a boundary condition with one or more boundary ids
among other things, the boundary condition objects contain a list of the surface facets to which they apply.
so that's your list of surface elements

does that mean that the easiest thing is to associate any boundary condition with a boundary, then retrieve the list of surface elements through the mechanism you just described?
since I imagine we want to be able to write in the "UFL" box in diamond something like: DirichletBC(mesh, bdry_id, condition)
and not actually use the fluidity infrastructure for computing boundary conditions
I think that's correct. And actually we can eventually unify the two by that route
is there any danger that i'll end up overwriting some other boundary condition specified in the flml that the user wanted by inserting my own "dummy" BC?
If we have a FEniCS-compatible way of doing it, then we could eventually call our routines driven by the flml to cause the fluidity interface.
I think for the moment, the rule is that UFL equations do not honour flml boundary conditions so there is no conflict
ah, that would be the "unified" case, and there would be no conflict

The weak BCs way of doing this is described at: http://fenicsproject.org/documentation/dolfin/1.0.0/python/demo/pde/subdomains-poisson/python/documentation.html
if i understand correctly, fluidity has already created boundaries and assigned them IDs based on the input mesh - is that right?
however, that doesn't preclude us from supporting the fenics way of doing it, since we could always compute the subdomain specified anyway

mesh%faces%face_list: face (i,j) is the face in element i adjacent to element j. the i-th row of that sparsity is the list of facets of that element in local facet number order. j will be negative for a boundary facet, but you probably don't need to care about that.
mesh%faces%face_element_list gives you the element corresponding to a face.
You just need to map from the global face number to the element local facet id which you need to get by finding the index of that facet in the ith row of face list.
Given that this is the way that we will always do integrals over ds, it may be worth cacheing that lookup.

We need:

An OP2 Set of elements for each boundary.
A Dat on this set containing the local face number of the face that is on the boundary, which agrees with the UFC specification on facet numbering.
A mapping from the local element numbers of each boundary element to the global node numbers.

To build the boundary data structures:

First uniquify mesh.faces.boundary_ids to get a list of the boundaries we need to build.
Scan through the column indices for face_list. If an entry is negative, then that element/face is on the boundary.
For each element that is on the boundary, we add it to the correct boundary element set. The correct set is looked up by looking up face_list[i,j] where i is the element and j is the negative index with the boundary, and using this to index into the boundary_ids array. The mappings to the global node numbers are simply looked up through ndglno for that element. The facet number is the index of j into the i-th row of face_list - this is add to the facet number dat.