# Steady-state continuity equation on an extruded mesh
# ====================================================
#
# This demo showcases the use of extruded meshes, including the new regions of
# integration and the construction of sophisticated finite element spaces.
#
# We now consider the equation
#
# .. math::
#
# \nabla\cdot(\vec{u}q) &= 0 \\
# q &= q_\mathrm{in} \quad \text{on} \ \Gamma_\mathrm{inflow},
#
# in a domain :math:`\Omega`, where :math:`\vec{u}` is a prescribed vector field,
# and :math:`q` is an unknown scalar field. The value of :math:`q` is known on the
# 'inflow' part of the boundary :math:`\Gamma`, where :math:`\vec{u}` is directed
# towards the interior of the domain. :math:`q` can be interpreted as the
# steady-state distribution of a passive tracer carried by a fluid with velocity
# field :math:`\vec{u}`.
#
# We apply an upwind DG method, as we saw in the
# :doc:`previous example `. Denoting the
# upwind value of :math:`q` on interior facets by :math:`\widetilde{q}`, the full
# set of equations are then
#
# .. math::
#
# -\int_\Omega \! q \vec{u_0} \cdot \nabla \phi \, \mathrm{d} x
# + \int_{\Gamma_{\mathrlap{\mathrm{ext, outflow}}}} \! \phi q \vec{u} \cdot \vec{n}
# \, \mathrm{d} s
# + \int_{\Gamma_\mathrm{int}} \! (\phi_+ \vec{u} \cdot \vec{n}_+ +
# \phi_- \vec{u} \cdot \vec{n}_-) \widetilde{q} \, \mathrm{d} S
# \quad = \quad
# -\int_{\Gamma_{\mathrlap{\mathrm{ext, inflow}}}} \phi q_\mathrm{in} \vec{u} \cdot
# \vec{n} \, \mathrm{d} s \quad \forall \ \phi \in V,
#
# We will take the domain :math:`\Omega` to be the cuboid
# :math:`\Omega = [0,1] \times [0,1] \times [0,0.2]`. We will use the uniform
# velocity field :math:`\vec{u} = (0, 0, 1)`. :math:`\Gamma_\mathrm{inflow}`
# is therefore the base of the cuboid, while :math:`\Gamma_\mathrm{outflow}`
# is the top. The four vertical sides can be ignored, since
# :math:`\vec{u} \cdot \vec{n} = 0` on these faces.
#
# We use an *extruded* mesh, where the base mesh is a 20 by 20 unit square,
# divided into triangles, with 10 evenly-spaced vertical layers. This gives
# prism-shaped cells. ::
from firedrake import *
m = UnitSquareMesh(20, 20)
mesh = ExtrudedMesh(m, layers=10, layer_height=0.02)
# We will use a simple piecewise-constant function space for the unknown scalar
# :math:`q`: ::
V = FunctionSpace(mesh, "DG", 0)
# Our velocity will live in a low-order Raviart-Thomas space. The construction of
# this is more complicated than element spaces that have appeared previously. The
# horizontal and vertical components of the field are specified separately. They
# are combined into a single element which is used to build a FunctionSpace. ::
# RT1 element on a prism
W0_h = FiniteElement("RT", "triangle", 1)
W0_v = FiniteElement("DG", "interval", 0)
W0 = HDivElement(TensorProductElement(W0_h, W0_v))
W1_h = FiniteElement("DG", "triangle", 0)
W1_v = FiniteElement("CG", "interval", 1)
W1 = HDivElement(TensorProductElement(W1_h, W1_v))
W_elt = W0 + W1
W = FunctionSpace(mesh, W_elt)
# As an aside, since our prescibed velocity is purely in the vertical direction, a
# simpler space would have sufficed: ::
# Vertical part of RT1 element
# W_h = FiniteElement("DG", "triangle", 0)
# W_v = FiniteElement("CG", "interval", 1)
# W_elt = HDivElement(TensorProductElement(W_h, W_v))
# W = FunctionSpace(mesh, W_elt)
# Or even: ::
# Why can't everything in life be this easy?
# W = VectorFunctionSpace(mesh, "CG", 1)
# Next, we set the prescribed velocity field: ::
velocity = as_vector((0.0, 0.0, 1.0))
u = project(velocity, W)
# if we had used W = VectorFunctionSpace(mesh, "CG", 1), we could have done
# u = Function(W)
# u.interpolate(velocity)
# Next, we will set the boundary value on our scalar to be a simple indicator
# function over part of the bottom of the domain: ::
x, y, z = SpatialCoordinate(mesh)
inflow = conditional(And(z < 0.02, x > 0.5), 1.0, -1.0)
q_in = Function(V)
q_in.interpolate(inflow)
# Now we will define our forms. We use the same trick as in the
# :doc:`previous example ` of defining ``un`` to aid
# with the upwind terms: ::
n = FacetNormal(mesh)
un = 0.5*(dot(u, n) + abs(dot(u, n)))
# We define our trial and test functions in the usual way: ::
q = TrialFunction(V)
phi = TestFunction(V)
# Since we are on an extruded mesh, we have several new integral types at our
# disposal. An integral over the cells of the domain is still denoted by ``dx``.
# Boundary integrals now come in several varieties: ``ds_b`` denotes an integral
# over the base of the mesh, while ``ds_t`` denotes an integral over the top of
# the mesh. ``ds_v`` denotes an integral over the sides of a mesh, though we will
# not use that here.
#
# Similiarly, interior facet integrals are split into ``dS_h`` and ``dS_v``, over
# *horizontal* interior facets and *vertical* interior facets respectively. Since
# our velocity field is purely in the vertical direction, we will omit the
# integral over vertical interior facets, since we know
# :math:`\vec{u} \cdot \vec{n}` is zero for these. ::
a1 = -q*dot(u, grad(phi))*dx
a2 = dot(jump(phi), un('+')*q('+') - un('-')*q('-'))*dS_h
a3 = dot(phi, un*q)*ds_t # outflow at top wall
a = a1 + a2 + a3
L = -q_in*phi*dot(u, n)*ds_b # inflow at bottom wall
# Finally, we will compute the solution: ::
out = Function(V)
solve(a == L, out)
# By construction, the exact solution is quite simple: ::
exact = Function(V)
exact.interpolate(conditional(x > 0.5, 1.0, -1.0))
# We finally compare our solution to the expected solution: ::
assert max(abs(out.dat.data - exact.dat.data)) < 1e-10
# This demo can be found as a script in
# :demo:`extruded_continuity.py `.