import firedrake as fd
import numpy as np

# Define mesh
Delta_x = 0.1
Delta_y = 0.1
Delta_z = 0.1

Nx = int(10)
Ny = int(20)

Nz = int(10)

Lx = Nx*Delta_x
Ly = Ny*Delta_y
Lz = Nz*Delta_z

mesh = fd.utility_meshes.RectangleMesh(nx=Nx, ny=Ny, Lx=Lx, Ly=Ly,
                                       quadrilateral=True)
mesh = fd.ExtrudedMesh(mesh, layers=Nz, layer_height=Delta_z)

# Create extruded mesh element
horiz_elt = fd.FiniteElement("DG", fd.quadrilateral, 0)
vert_elt =  fd.FiniteElement("DG", fd.interval, 0)
elt = fd.TensorProductElement(horiz_elt, vert_elt)

# Create function space
V = fd.FunctionSpace(mesh, elt)

# Create function f and point x
f = fd.interpolate(fd.Constant(0.25), V)
x = np.array([Delta_x/2., Delta_y/2., Delta_z/2.])

# Evaluate f at point x
fd.File("test.pvd").write(f)

print(f.at(x, tolerance=1e-8))