Implementation of CutFEM in Firedrake

[miguelcoolchips]

Starting this discussion thread to have a common place for everything related to the implementation of CutFEM in Firedrake. The critical component required is the custom quadrature with runtime tabulation of points. This routine would need to accept the points in physical space. In my discussion with @ReubenHill on slack, he pointed out that he implemented the infrastructure, but the API is not there yet. Here is his response:

Reuben Nixon-Hill:
If you take a look at what happens when you interpolate onto a vertex only mesh you see what we have to do to get runtime points

firedrake/firedrake/interpolation.py

Lines 255 to 275 in 2ec31f9

	if target_mesh is not source_mesh:
	if not isinstance(target_mesh.topology, firedrake.mesh.VertexOnlyMeshTopology):
	raise NotImplementedError("Can only interpolate onto a Vertex Only Mesh")
	if target_mesh.geometric_dimension() != source_mesh.geometric_dimension():
	raise ValueError("Cannot interpolate onto a mesh of a different geometric dimension")
	if not hasattr(target_mesh, "_parent_mesh") or target_mesh._parent_mesh is not source_mesh:
	raise ValueError("Can only interpolate across meshes where the source mesh is the parent of the target")
	# For trans-mesh interpolation we use a FInAT QuadratureElement as the
	# (base) target element with runtime point set expressions as their
	# quadrature rule point set and weights from their dual basis.
	# NOTE: This setup is useful for thinking about future design - in the
	# future this `rebuild` function can be absorbed into FInAT as a
	# transformer that eats an element and gives you an equivalent (which
	# may or may not be a QuadratureElement) that lets you do run time
	# tabulation. Alternatively (and this all depends on future design
	# decision about FInAT how dual evaluation should work) the
	# to_element's dual basis (which look rather like quadrature rules) can
	# have their pointset(s) directly replaced with run-time tabulated
	# equivalent(s) (i.e. finat.point_set.UnknownPointSet(s))
	rt_var_name = 'rt_X'
	to_element = rebuild(to_element, expr, rt_var_name)

Specifically this comment

firedrake/firedrake/interpolation.py

Lines 262 to 264 in 2ec31f9

	# For trans-mesh interpolation we use a FInAT QuadratureElement as the
	# (base) target element with runtime point set expressions as their
	# quadrature rule point set and weights from their dual basis.

Quadrature elements in FInAT (https://github.com/FInAT/FInAT/blob/master/finat/quadrature_element.py) require an AbstractQuadratureRule (see https://github.com/FInAT/FInAT/blob/4112fb88b56d6f740c6de0ff105927ca1492a4c7/finat/quadrature.py) which have point_set and weight_expression properties. An option for a point set is an UnknownPointSet (see https://github.com/FInAT/FInAT/blob/master/finat/point_set.py#L86)
In TSFC we create a kernel for interpolation (dual evaluation) using the compile_expression_dual_evaluation function (see https://github.com/firedrakeproject/tsfc/blob/master/tsfc/driver.py#L161)

I will expand this thread as we make progress on the implementation.