Thanks Lawrence! I will give this (as well as the other mail about projecting data) a try shortly. On Wed, Feb 1, 2017 at 4:02 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 01/02/17 09:53, Justin Chang wrote:
Hi all,
I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners.
Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition?
What you need to do is, before doing anything with your mesh, add additional boundary marker labels that correspond to the appropriate facets.
You do this half through plex:
from firedrake import *
mesh = UnitSquareMesh(...)
plex = mesh._plex
coords = plex.getCoordinates() coord_sec = plex.getCoordinateSection()
if plex.getStratumSize("boundary_faces", 1) > 0: faces = plex.getStratumIS("boundary_faces", 1).getIndices() for face in faces: face_coords = plex.vecGetClosure(coord_sec, coords, face) # now face_coords is a numpy array of the vertices of # the face (XY XY) if face_coords are in corner face: plex.setLabelValue("boundary_ids", face, ID)
V = FunctionSpace(mesh, "RT", 1) ...
now you can use ID just as you do the built in mesh ids.
(untested, but this is general idea).
If you look in utility_meshes.py at (say) RectangleMesh, you can see how we mark the full faces of the mesh and copy that logic.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake