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). Thanks, Anna.
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
Hi Lawrence, I had already tried this based on examples from firedrake/utility_meshes.py, but I get the following error: Traceback (most recent call last): File "mesh.py", line 39, in <module> V = FunctionSpace(mesh, "CG", 1) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/functionspace.py", line 636, in __new__ mesh.init() File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/mesh.py", line 750, in init self._callback(self) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/mesh.py", line 1040, in callback dim=geometric_dim) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/functionspace.py", line 695, in __new__ self = super(VectorFunctionSpace, cls).__new__(cls, mesh_t, element, name=name, shape=(dim,)) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/functionspace.py", line 140, in __new__ dmplex.get_facet_nodes(mesh.exterior_facets.facet_cell, File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", line 64, in __get__ obj.__dict__[self.__name__] = result = self.fget(obj) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/mesh.py", line 471, in exterior_facets boundary_ids, unique_markers=unique_ids) File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/mesh.py", line 61, in __init__ "Every marker has to be contained in unique_markers" AssertionError: Every marker has to be contained in unique_markers A sample code can be found attached. Best, Anna.
On 8 Jan 2016, at 12:20, Anna Kalogirou <A.Kalogirou@leeds.ac.uk> wrote:
Hi Lawrence,
I had already tried this based on examples from firedrake/utility_meshes.py, but I get the following error:
Ah, the problem is that we expect that all boundary facets are marked with something! You could add an else branch to you code that marks with (say) 0. Cheers, Lawrence
participants (3)
-
Anna Kalogirou
-
Anna Kalogirou
-
Lawrence Mitchell