# Source code for ufl.algorithms.change_to_reference

# -*- coding: utf-8 -*-
"""Algorithm for replacing gradients in an expression with reference gradients and coordinate mappings."""

# Copyright (C) 2008-2016 Martin Sandve Alnæs
#
# This file is part of UFL (https://www.fenicsproject.org)
#

from ufl.log import error

from ufl.core.multiindex import indices
from ufl.corealg.multifunction import MultiFunction
from ufl.corealg.map_dag import map_expr_dag

from ufl.classes import (FormArgument, GeometricQuantity,
Jacobian, JacobianInverse, JacobianDeterminant,
Indexed, MultiIndex, FixedIndex)

from ufl.constantvalue import as_ufl
from ufl.tensors import as_tensor
from ufl.permutation import compute_indices

from ufl.finiteelement import MixedElement

from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks
from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering
from ufl.checks import is_cellwise_constant

"""
# Some notes:
# Below, let v_i mean physical coordinate of vertex i and V_i mean the reference cell coordinate of the same vertex.

# Add a type for CellVertices? Note that vertices must be computed in linear cell cases!
triangle_vertices[i,j] = component j of vertex i, following ufc numbering conventions

# DONE Add a type for CellEdgeLengths? Note that these are only easy to define in the linear cell case!
triangle_edge_lengths    = [v1v2, v0v2, v0v1] # shape (3,)
tetrahedron_edge_lengths = [v0v1, v0v2, v0v3, v1v2, v1v3, v2v3] # shape (6,)

# DONE Here's how to compute edge lengths from the Jacobian:
J =[ [dx0/dX0, dx0/dX1],
[dx1/dX0, dx1/dX1] ]
# First compute the edge vector, which is constant for each edge: the vector from one vertex to the other
reference_edge_vector_0 = V2 - V1 # Example! Add a type ReferenceEdgeVectors?
# Then apply J to it and take the length of the resulting vector, this is generic for affine cells
edge_length_i = || dot(J, reference_edge_vector_i) ||

e2 = || J[:,0] . < 1, 0> || = || J[:,0] || = || dx/dX0 || = edge length of edge 2 (v0-v1)
e1 = || J[:,1] . < 0, 1> || = || J[:,1] || = || dx/dX1 || = edge length of edge 1 (v0-v2)
e0 = || J[:,:] . <-1, 1> || = || < J[0,1]-J[0,0], J[1,1]-J[1,0] > || = || dx/dX <-1,1> || = edge length of edge 0 (v1-v2)

trev = triangle_reference_edge_vector
evec0 = J00 * trev[edge][0] + J01 * trev[edge][1]  =  J*trev[edge]
evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
elen[edge] = sqrt(evec0*evec0 + evec1*evec1)  = sqrt((J*trev[edge])**2)

trev = triangle_reference_edge_vector
evec0 = J00 * trev[edge][0] + J01 * trev[edge][1]  =  J*trev
evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
evec2 = J20 * trev[edge][0] + J21 * trev[edge][1] # Manifold: triangle in 3D
elen[edge] = sqrt(evec0*evec0 + evec1*evec1 + evec2*evec2)  = sqrt((J*trev[edge])**2)

trev = tetrahedron_reference_edge_vector
evec0 = sum(J[0,k] * trev[edge][k] for k in range(3))
evec1 = sum(J[1,k] * trev[edge][k] for k in range(3))
evec2 = sum(J[2,k] * trev[edge][k] for k in range(3))
elen[edge] = sqrt(evec0*evec0 + evec1*evec1 + evec2*evec2)  = sqrt((J*trev[edge])**2)

# DONE Here's how to compute min/max facet edge length:
triangle:
r = facetarea
tetrahedron:
min(elen[edge] for edge in range(6))
or
min( min(elen[0], min(elen[1], elen[2])), min(elen[3], min(elen[4], elen[5])) )
or
min1 = min_value(elen[0], min_value(elen[1], elen[2]))
min2 = min_value(elen[3], min_value(elen[4], elen[5]))
r = min_value(min1, min2)
(want proper Min/Max types for this!)

# DONE Here's how to compute circumradius for an interval:

# DONE Here's how to compute circumradius for a triangle:
e0 = elen[0]
e1 = elen[1]
e2 = elen[2]

# DONE Here's how to compute circumradius for a tetrahedron:
# v1v2 = edge length between vertex 1 and 2
# la,lb,lc = lengths of the sides of an intermediate triangle
la = v1v2 * v0v3
lb = v0v2 * v1v3
lc = v0v1 * v2v3
# p = perimeter
p = (la + lb + lc)
# s = semiperimeter
s = p / 2
# area of intermediate triangle with Herons formula
tmp_area = sqrt(s * (s - la) * (s - lb) * (s - lc))

"""

# FIXME: This implementation semeed to work last year but lead to performance problems. Look through and test again now.
[docs]class NEWChangeToReferenceGrad(MultiFunction): def __init__(self): MultiFunction.__init__(self) self._ngrads = 0 self._restricted = '' self._avg = ''
[docs] def expr(self, o, *ops): return o._ufl_expr_reconstruct_(*ops)
[docs] def terminal(self, o): return o
[docs] def coefficient_derivative(self, o, *dummy_ops): error("Coefficient derivatives should be expanded before applying change to reference grad.")
[docs] def restricted(self, o, *dummy_ops): "Store modifier state." if self._restricted != '': error("Not expecting nested restrictions.") self._restricted = o.side() f, = o.ufl_operands r = self(f) self._restricted = '' return r
[docs] def grad(self, o, *dummy_ops): "Store modifier state." self._ngrads += 1 f, = o.ufl_operands r = self(f) self._ngrads -= 1 return r
[docs] def facet_avg(self, o, *dummy_ops): if self._avg != '': error("Not expecting nested averages.") self._avg = "facet" f, = o.ufl_operands r = self(f) self._avg = "" return r
[docs] def cell_avg(self, o, *dummy_ops): if self._avg != '': error("Not expecting nested averages.") self._avg = "cell" f, = o.ufl_operands r = self(f) self._avg = "" return r
[docs] def form_argument(self, t): return self._mapped(t)
[docs] def geometric_quantity(self, t): if self._restricted or self._ngrads or self._avg: return self._mapped(t) else: return t
def _mapped(self, t): # Check that we have a valid input object if not isinstance(t, Terminal): error("Expecting a Terminal.") # Get modifiers accumulated by previous handler calls ngrads = self._ngrads restricted = self._restricted avg = self._avg if avg != "": error("Averaging not implemented.") # FIXME # These are the global (g) and reference (r) values if isinstance(t, FormArgument): g = t r = ReferenceValue(g) elif isinstance(t, GeometricQuantity): g = t r = g else: error("Unexpected type {0}.".format(type(t).__name__)) # Some geometry mapping objects we may need multiple times below domain = t.ufl_domain() J = Jacobian(domain) detJ = JacobianDeterminant(domain) K = JacobianInverse(domain) # Restrict geometry objects if applicable if restricted: J = J(restricted) detJ = detJ(restricted) K = K(restricted) # Create Hdiv mapping from possibly restricted geometry objects Mdiv = (1.0 / detJ) * J # Get component indices of global and reference terminal objects gtsh = g.ufl_shape # rtsh = r.ufl_shape gtcomponents = compute_indices(gtsh) # rtcomponents = compute_indices(rtsh) # Create core modified terminal, with eventual # layers of grad applied directly to the terminal, # then eventual restriction applied last for i in range(ngrads): g = Grad(g) r = ReferenceGrad(r) if restricted: g = g(restricted) r = r(restricted) # Get component indices of global and reference objects with # grads applied gsh = g.ufl_shape # rsh = r.ufl_shape # gcomponents = compute_indices(gsh) # rcomponents = compute_indices(rsh) # Get derivative component indices dsh = gsh[len(gtsh):] dcomponents = compute_indices(dsh) # Create nested array to hold expressions for global # components mapped from reference values def ndarray(shape): if len(shape) == 0: return [None] elif len(shape) == 1: return [None] * shape[-1] else: return [ndarray(shape[1:]) for i in range(shape[0])] global_components = ndarray(gsh) # Compute mapping from reference values for each global component for gtc in gtcomponents: if isinstance(t, FormArgument): # Find basic subelement and element-local component # ec, element, eoffset = t.ufl_element().extract_component2(gtc) # FIXME: Translate this correctly eoffset = 0 ec, element = t.ufl_element().extract_reference_component(gtc) # Select mapping M from element, pick row emapping = # M[ec,:], or emapping = [] if no mapping if isinstance(element, MixedElement): error("Expecting a basic element here.") mapping = element.mapping() if mapping == "contravariant Piola": # S == HDiv: # Handle HDiv elements with contravariant piola # mapping contravariant_hdiv_mapping = (1/det J) * # J * PullbackOf(o) ec, = ec emapping = Mdiv[ec, :] elif mapping == "covariant Piola": # S == HCurl: # Handle HCurl elements with covariant piola mapping # covariant_hcurl_mapping = JinvT * PullbackOf(o) ec, = ec emapping = K[:, ec] # Column of K is row of K.T elif mapping == "identity": emapping = None else: error("Unknown mapping {0}".format(mapping)) elif isinstance(t, GeometricQuantity): eoffset = 0 emapping = None else: error("Unexpected type {0}.".format(type(t).__name__)) # Create indices # if rtsh: # i = Index() if len(dsh) != ngrads: error("Mismatch between derivative shape and ngrads.") if ngrads: ii = indices(ngrads) else: ii = () # Apply mapping row to reference object if emapping: # Mapped, always nonscalar terminal Not # using IndexSum for the mapping row dot product to # keep it simple, because we don't have a slice type emapped_ops = [emapping[s] * Indexed(r, MultiIndex((FixedIndex(eoffset + s),) + ii)) for s in range(len(emapping))] emapped = sum(emapped_ops[1:], emapped_ops[0]) elif gtc: # Nonscalar terminal, unmapped emapped = Indexed(r, MultiIndex((FixedIndex(eoffset),) + ii)) elif ngrads: # Scalar terminal, unmapped, with derivatives emapped = Indexed(r, MultiIndex(ii)) else: # Scalar terminal, unmapped, no derivatives emapped = r for di in dcomponents: # Multiply derivative mapping rows, parameterized by # free column indices dmapping = as_ufl(1) for j in range(ngrads): dmapping *= K[ii[j], di[j]] # Row of K is column of JinvT # Compute mapping from reference values for this # particular global component global_value = dmapping * emapped # Apply index sums # if rtsh: # global_value = IndexSum(global_value, MultiIndex((i,))) # for j in range(ngrads): # Applied implicitly in the dmapping * emapped above # global_value = IndexSum(global_value, MultiIndex((ii[j],))) # This is the component index into the full object # with grads applied gc = gtc + di # Insert in nested list comp = global_components for i in gc[:-1]: comp = comp[i] comp[0 if gc == () else gc[-1]] = global_value # Wrap nested list in as_tensor unless we have a scalar # expression if gsh: tensor = as_tensor(global_components) else: tensor, = global_components return tensor
[docs]class OLDChangeToReferenceGrad(MultiFunction): def __init__(self): MultiFunction.__init__(self) expr = MultiFunction.reuse_if_untouched
[docs] def terminal(self, o): return o