Remove component tensors

[pbrubeck]

Consider the bilinear form
a(v, u) = inner(u, v)*dx + inner(div(u), div(v))*dx
where u is a (doubly) Piola mapped tensor in H(div, S) discretized with linear Johnson-Mercier macroelements.

Here we introduce two enhacements targeted to optimize code generation for a(v, u):

1. We eagerly simplify RefGrad(cellwise_constant) to avoid expressions like RefGrad(Jacobian) -> 0 on affine simplex meshes. These zero terms would arise from the product rule on the doubly Piola mapped tensor.
2. We eagerly remove ComponentTensors, to make the UFL easier to traverse by the Form compiler.

Since a Form always maps to a scalar, it is always possible to remove ComponentTensor nodes by pushing Indexed inside expressions and simplifying Indexed(ComponentTensor) with the approapiate index replacements.

This PR makes more agressive simplifications within the Indexed constructor, so it's pushed inside sums, and adds a multifunction to handle index replacements. The removal of ComponentTensors is applied as a final preprocessing step in compute_form_data.

[jorgensd]

In general this looks sensible to me. Some minor comments/suggestions only.

[michalhabera]

You've introduced here that is_cellwise_constant logic is based on degree estimation, which is the only place in the code. How does it work, since the degree estimation for geometric quantities calls is_cellwise_constant itself?

ufl/ufl/algorithms/estimate_degrees.py

Lines 46 to 56 in 1ab30c2

	def geometric_quantity(self, v):
	"""Apply to geometric_quantity.
	Some geometric quantities are cellwise constant. Others are
	nonpolynomial and thus hard to estimate.
	"""
	if is_cellwise_constant(v):
	return 0
	else:
	# As a heuristic, just returning domain degree to bump up degree somewhat
	return extract_unique_domain(v).ufl_coordinate_element().embedded_superdegree

If there is no infinite loop, could this be used in the generic cell-wise constant check? Or is it that only GradRuleset does not cause infinite loop?

[pbrubeck]

I first tried to add this logic to the generic is_cellwise_constant in #334, but simplifying grad(cell_wise_constant) on instantiation turns out to be problematic for nested grad expression simply because grad(f) expects f to have a domain and Zero does not have a domain, so grad(grad(grad(linear))) would simplify to grad(Zero).

[michalhabera]

Is this correct? Looks like one derivative is lost.

[pbrubeck]

Here we are simplifying grad(o) where o = grad(f), if o is constant then we should return zero with the correct shape.