Hi Anna,

On 7 Jan 2016, at 10:58, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote:

Dear all,

I have the mesh defined in the file attached. I can successfully import it and use it with Firedrake, but how can I define a surface integral? For example, I want to somehow define a "ds" which will denote the left boundary of my domain, corresponding to y=0. The nodes that will define this "ds" are conveniently the first 20 nodes in the mesh file attached, with y=0 (the second column).

The way you have to do this is by marking the faces in the plex object before you pass it to the Mesh constructor.

I think what you're doing at the moment is something like:

plex = _from_cell_list(...)
mesh = Mesh(plex)


Before creating the mesh, you can do something like:

plex.createLabel("boundary_ids")
plex.markBoundaryFaces("boundary_faces")
coords = plex.getCoordinates()
coord_sec = plex.getCoordinateSection()
if plex.getStratumSize("boundary_faces", 1) > 0:
   faces = plex.getStratumIS("boundary_faces", 1).getIndices()
   eps = 1e-5
   for face in faces:
       face_coords = plex.vecGetClosure(coord_sec, coords, face)
       # face_coords has 4 entries x0, y0, x1, y1.
       if abs(face_coords[1]) < eps and abs(face_coords[3]) < eps:
           plex.setLabelValue("boundary_ids", face, 1)

mesh = Mesh(plex)

Now you can use ds(1) to integrate over the boundary with y == 0.

Lawrence
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake