#!/usr/bin/env python

from os.path import basename, splitext

from firedrake import (Mesh, FunctionSpace, TestFunction, DirichletBC,
                       interpolate, Constant, assemble, solve, File)

from ufl import dx, inner, grad


def main(meshfile: str) -> None:

    mesh = Mesh(meshfile)

    V = FunctionSpace(mesh, 'CG', 1)
    v = TestFunction(V)
    bc = DirichletBC(V, 0.0, [1])

    w = interpolate(Constant(1.0, mesh), V)
    area = assemble(w * dx)
    print('Area = ', area)

    F = (inner(grad(w), grad(v)) - v) * dx
    solve(F == 0, w, bcs=bc, solver_parameters={'ksp_type': 'cg'})

    volume = assemble(w * dx)
    print('Boussinesq k-factor = ', volume / area**2)

    File(splitext(basename(meshfile))[0] + '.pvd').write(w)


if __name__ == '__main__':

    from sys import argv

    main(argv[1])