[firedrake] Prescribe pressure at a corner for Darcy's equation

[Justin Chang]

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 at 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 at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
>
-------------- next part --------------
HTML attachment scrubbed and removed