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