@@ -0,0 +1,2780 @@
+from firedrake import *%0Afrom firedrake.utils import cached_property%0A%0Aimport base%0A%0A%0Aclass HDGProblem(base.Problem):%0A%0A    name = %22HDG Helmholtz%22%0A%0A    @cached_property%0A    def tau(self):%0A        # Stability parameter for the HDG method%0A        return Constant(100.)%0A%0A    @cached_property%0A    def function_space(self):%0A        U = VectorFunctionSpace(self.mesh, %22DG%22, self.degree)%0A        V = FunctionSpace(self.mesh, %22DG%22, self.degree)%0A        T = FunctionSpace(self.mesh, %22DGT%22, self.degree)%0A        return U * V * T%0A%0A    @cached_property%0A    def forcing(self):%0A        f = Function(self.function_space, name=%22forcing%22)%0A        u = self.analytic_solution%0A        f.interpolate(-div(grad(u)) + u)%0A        return f%0A%0A    @cached_property%0A    def a(self):%0A        W = self.function_space%0A        sigma, u, lambdar = TrialFunctions(W)%0A        tau, v, gamma = TestFunctions(W)%0A        n = FacetNormal(self.mesh)%0A%0A        sigmahat = sigma + self.tau*(u - lambdar)*n%0A%0A        def both(arg):%0A            return 2*avg(arg)%0A%0A        a = (dot(sigma, tau)*dx - div(tau)*u*dx +%0A             both(sigmahat*dot(tau, n))*dS +%0A             sigmahat*dot(tau, n)*ds%0A             - dot(grad(v), sigma)*dx +%0A             both(dot(sigmahat, n)*v)*dS +%0A             dot(sigmahat, n)*v*ds%0A             + both(gamma*dot(sigmahat, n))*dS%0A             + gamma*lambdar*ds)%0A        return a%0A%0A    @cached_property%0A    def L(self):%0A        W = self.function_space%0A        _, v, gamma = TestFunctions(W)%0A        f = self.forcing%0A        u0 = Constant(1.0)%0A        return inner(f, v)*dx + inner(u0, gamma)*ds%0A%0A    @cached_property%0A    def analytic_flux(self):%0A        u = self.analytic_solution%0A        return -grad(u)%0A%0A    @cached_property%0A    def bcs(self):%0A        # All bcs are weakly enforced%0A        return None%0A%0A    @cached_property%0A    def output(self):%0A        sigma, u, _ = self.u.split()%0A        return (sigma, u)%0A%0A    @cached_property%0A    def err(self):%0A        u_a = Function(self.function_space%5B1%5D)%0A        u_a.interpolate(self.analytic_solution)%0A        sigma_a = FunctionSpace(self.function_space%5B0%5D)%0A        sigma_a.project(self.analytic_flux)%0A        u_err = errornorm(self.u.split()%5B1%5D, u_a)%0A        sigma_err = errornorm(self.u.split()%5B0%5D, sigma_a)%0A        return (sigma_err, u_err)%0A%0A    @cached_property%0A    def true_err(self):%0A        u_err = errornorm(self.analytic_solution, self.u.split()%5B1%5D)%0A        sigma_err = errornorm(self.analytic_flux, self.u.split()%5B0%5D)%0A        return (sigma_err, u_err)%0A%0A    @cached_property%0A    def sol(self):%0A        sigma = Function(self.function_space%5B0%5D, name=%22Analytic flux%22)%0A        u = Function(self.function_space%5B1%5D, name=%22Analytic scalar%22)%0A        u.interpolate(self.analytic_solution)%0A        sigma.project(self.analytic_flux)%0A        return (sigma, u)%0A

Set up HDG problem