@@ -0,0 +1,1896 @@
+%22%22%22Solve a mixed Helmholtz problem%0A%0Asigma + grad(u) = 0,%0Au + div(sigma) = f,%0A%0Ausing hybridisation with SLATE performing the forward elimination and%0Abackwards reconstructions. The corresponding finite element variational%0Aproblem is:%0A%0Adot(sigma, tau)*dx - u*div(tau)*dx + lambdar*dot(tau, n)*dS = 0%0Adiv(sigma)*v*dx + u*v*dx = f*v*dx%0Agammar*dot(sigma, n)*dS = 0%0A%0Afor all tau, v, and gammar.%0A%0AThis is solved using broken Raviart-Thomas elements of degree k for%0A(sigma, tau), discontinuous Galerkin elements of degree k - 1%0Afor (u, v), and HDiv-Trace elements of degree k - 1 for (lambdar, gammar).%0A%0AThe forcing function is chosen as:%0A%0A(1+8*pi*pi)*sin(x%5B0%5D*pi*2)*sin(x%5B1%5D*pi*2),%0A%0Awhich reproduces the known analytical solution:%0A%0Asin(x%5B0%5D*pi*2)*sin(x%5B1%5D*pi*2)%0A%22%22%22%0A%0Afrom __future__ import absolute_import, print_function, division%0A%0Afrom firedrake import *%0A%0A# Create a mesh%0Amesh = UnitSquareMesh(120, 120)%0A%0A# Define relevant function spaces%0Adegree = 2%0ART = FiniteElement(%22RT%22, triangle, degree)%0ABRT = FunctionSpace(mesh, RT)%0ADG = FunctionSpace(mesh, %22DG%22, degree - 1)%0A%0AW = BRT * DG%0A%0A# Define the trial and test functions%0Asigma, u = TrialFunctions(W)%0Atau, v = TestFunctions(W)%0A%0A# Define the source function%0Af = Function(DG)%0Ax, y = SpatialCoordinate(mesh)%0Af.interpolate((1+8*pi*pi)*sin(x*pi*2)*sin(y*pi*2))%0A%0A# Define finite element variational forms%0Aa = (dot(sigma, tau) + u * v + div(sigma) * v - div(tau) * u) * dx%0AL = f * v * dx%0Aw = Function(W)%0Aparameters = %7B%22mat_type%22: %22matfree%22,%0A              %22pc_type%22: %22python%22,%0A              %22pc_python_type%22: %22firedrake.HybridizationPC%22,%0A              %22trace_pc_type%22: %22lu%22,%0A              %22trace_ksp_type%22: %22preonly%22,%0A              %22trace_ksp_monitor_true_residual%22: True,%0A              %22ksp_monitor_true_residual%22: True%7D%0Asolve(a == L, w, solver_parameters=parameters)%0A%0Au, p = w.split()%0Au.rename(%22Velocity%22)%0Ap.rename(%22Pressure%22)%0A%0AFile(%22pchybrid.pvd%22).write(u, p)%0A

Add complete hybridization test with python PC