Source code for ffc.tensor.monomialtransformation
# -*- coding: utf-8 -*-
"Transformation of monomial representations of UFL forms."
# Copyright (C) 2009 Anders Logg
#
# This file is part of FFC.
#
# FFC is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# FFC is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with FFC. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Kristian B. Oelgaard, 2009
# Modified by Marie E. Rognes, 2010
#
# First added: 2009-03-06
# Last changed: 2010-02-17
# UFL modules
from ufl.classes import Argument
from ufl.classes import Coefficient
from ufl.classes import FixedIndex
from ufl.permutation import build_component_numbering
# FFC modules
from ffc.log import info, error, ffc_assert
from ffc.fiatinterface import create_element
from ffc.utils import all_equal
from ffc.representationutils import transform_component
# FFC tensor representation modules
from ffc.tensor.monomialextraction import MonomialForm
from ffc.tensor.monomialextraction import MonomialException
[docs]def transform_monomial_form(monomial_form):
"Transform monomial form to reference element."
info("Transforming monomial form to reference element")
# Check that we get a monomial form
ffc_assert(isinstance(monomial_form, MonomialForm),
"Expecting a MonomialForm.")
# Note that we check if each monomial has been transformed before
# and if so we leave it untouched. This is to prevent repeated
# transformation (which fails) which may sometimes happen as a
# result of extracted integrands being cached by the monomial
# extraction.
# Transform each integral
for integrand in monomial_form:
for (i, monomial) in enumerate(integrand.monomials):
if not isinstance(monomial, TransformedMonomial):
integrand.monomials[i] = TransformedMonomial(monomial)
[docs]class MonomialIndex:
"""
This class represents a monomial index. Each index has a type,
a range and a unique id. Valid index types are listed below.
"""
FIXED = "fixed" # Integer index
PRIMARY = "primary" # Argument basis function index
SECONDARY = "secondary" # Index appearing both inside and outside integral
INTERNAL = "internal" # Index appearing only inside integral
EXTERNAL = "external" # Index appearing only outside integral
def __init__(self, index=None, index_type=None, index_range=None,
index_id=None):
"Create index with given type, range and id."
if isinstance(index, MonomialIndex):
self.index_type = index.index_type
self.index_range = [i for i in index.index_range]
self.index_id = index.index_id
else:
self.index_type = index_type
self.index_range = index_range
self.index_id = index_id
def __lt__(self, other):
"Comparison operator."
return self.index_id < other.index_id
def __call__(self, primary=None, secondary=None, internal=None,
external=None):
"Evaluate index at current index list."
if self.index_type == MonomialIndex.FIXED:
return self.index_range[0]
elif self.index_type == MonomialIndex.PRIMARY:
if not primary:
error("Missing index values for primary indices.")
return primary[self.index_id]
elif self.index_type == MonomialIndex.SECONDARY:
if not secondary:
error("Missing index values for secondary indices.")
return secondary[self.index_id]
elif self.index_type == MonomialIndex.INTERNAL:
if not internal:
error("Missing index values for internal auxiliary indices.")
return internal[self.index_id]
elif self.index_type == MonomialIndex.EXTERNAL:
if not external:
error("Missing index values for external auxiliary indices.")
return external[self.index_id]
else:
error("Unknown index type " + str(self.type))
def __add__(self, offset):
"Add offset to index range."
index = MonomialIndex(self)
index.index_range = [offset + i for i in index.index_range]
return index
def __sub__(self, offset):
"Subtract offset from index range."
return self + (-offset)
def __str__(self):
"Return informal string representation (pretty-print)."
if self.index_type == MonomialIndex.FIXED:
return str(self.index_range[0])
elif self.index_type == MonomialIndex.PRIMARY:
return "i_" + str(self.index_id)
elif self.index_type == MonomialIndex.SECONDARY:
return "a_" + str(self.index_id)
elif self.index_type == MonomialIndex.INTERNAL:
return "g_" + str(self.index_id)
elif self.index_type == MonomialIndex.EXTERNAL:
return "b_" + str(self.index_id)
else:
return "?"
[docs]class MonomialDeterminant:
"This class representes a determinant factor in a monomial."
def __init__(self, power=None, restriction=None):
"Create empty monomial determinant."
if power is None:
self.power = 0
else:
self.power = power
self.restriction = restriction
def __str__(self):
"Return informal string representation (pretty-print)."
# FIXME: This pretty-print is plain misleading b/c of the
# implicit relic factor |det J| etc.
if not self.restriction:
r = ""
else:
r = "(%s)" % self.restriction
if self.power == 0:
return "|det F'%s|" % r
elif self.power == 1:
return "|det F'%s| (det F'%s)" % (r, r)
else:
return "|det F'%s| (det F'%s)^%s" % (r, r, str(self.power))
[docs]class MonomialCoefficient:
"This class represents a coefficient in a monomial."
def __init__(self, index, number):
"Create monomial coefficient for given index and number."
self.index = index
self.number = number
def __str__(self):
"Return informal string representation (pretty-print)."
return "c_" + str(self.index)
[docs]class MonomialTransform:
"This class represents a transform (mapping derivative) in a form."
J = "J"
JINV = "JINV"
def __init__(self, index0, index1, transform_type, restriction, offset):
"Create monomial transform."
# Set data
self.index0 = index0
self.index1 = index1
self.transform_type = transform_type
self.restriction = restriction
self.offset = offset
# Subtract offset for fixed indices. Note that the index
# subtraction creates a new index instance. This is ok here
# since a fixed index does not need to match any other index
# (being the same instance) in index summation and index
# extraction.
if index0.index_type is MonomialIndex.FIXED:
self.index0 = index0 - offset
if index1.index_type is MonomialIndex.FIXED:
self.index1 = index1 - offset
def __str__(self):
"Return informal string representation (pretty-print)."
if self.restriction is None:
r = ""
else:
r = "(%s)" % str(self.restriction)
if self.transform_type == "J":
return "dx_%s/dX_%s%s" % (str(self.index0), str(self.index1), r)
else:
return "dX_%s/dx_%s%s" % (str(self.index0), str(self.index1), r)
[docs]class MonomialArgument:
"""
This class represents a monomial argument, that is, a derivative of
a scalar component of a basis function on the reference element.
"""
def __init__(self, element, index, components, derivatives, restriction):
"Create monomial argument."
self.element = element
self.index = index
self.components = components
self.derivatives = derivatives
self.restriction = restriction
def __str__(self):
"Return informal string representation (pretty-print)."
if len(self.components) == 0:
c = ""
else:
c = "[%s]" % ", ".join(str(c) for c in self.components)
if len(self.derivatives) == 0:
d0 = ""
d1 = ""
else:
d0 = "(" + " ".join("d/dX_%s" % str(d) for d in self.derivatives) + " "
d1 = ")"
if self.restriction is None:
r = ""
else:
r = "(%s)" % str(self.restriction)
v = "V_" + str(self.index)
return d0 + v + r + c + d1
[docs]class TransformedMonomial:
"""
This class represents a monomial form after transformation to the
reference element.
"""
def __init__(self, monomial):
"Create transformed monomial from given monomial."
# Reset monomial data
self.float_value = monomial.float_value
self.determinants = []
self.coefficients = []
self.transforms = []
self.arguments = []
# Reset index counters
_reset_indices()
# Initialize index map
index_map = {}
# Iterate over factors
for f in monomial.factors:
# Create FIAT element
ufl_element = f.element()
fiat_element = create_element(ufl_element)
# Note nifty aspect here: when gdim != tdim, it might be
# (in particular, is for H(div)/H(curl), that the value
# dimension is different for the physical and reference
# elements.
# Get number of components
# FIXME: Can't handle tensor-valued elements: vdim = shape[0]
shape = ufl_element.value_shape()
assert(len(shape) <= 1), \
"MonomialTransformation does not handle tensor-valued elements"
if len(shape) == 0:
vdim = 1
else:
vdim = shape[0]
# Extract dimensions
sdim = fiat_element.space_dimension()
cell = ufl_element.cell()
gdim = cell.geometric_dimension()
tdim = cell.topological_dimension()
# Extract basis function index and coefficients
if isinstance(f.function, Argument):
vindex = MonomialIndex(index_type=MonomialIndex.PRIMARY,
index_range=list(range(sdim)),
index_id=f.index())
elif isinstance(f.function, Coefficient):
vindex = MonomialIndex(index_range=list(range(sdim)))
coefficient = MonomialCoefficient(vindex, f.index())
self.coefficients.append(coefficient)
# Extract components
components = self._extract_components(f, index_map, vdim)
if len(components) > 1:
raise MonomialException("Can only handle rank 0 or rank 1 tensors.")
# Handle non-affine mappings (Piola)
if len(components) > 0:
# We can only handle rank 1 elements for now
component = components[0]
# Get mapping (all need to be equal)
mappings = []
for i in component.index_range:
(offset, ufl_sub_element) = ufl_element.extract_component(i)
fiat_sub_element = create_element(ufl_sub_element)
mappings.extend(fiat_sub_element.mapping())
if not all_equal(mappings):
raise MonomialException("Mappings differ: " + str(mappings))
mapping = mappings[0]
# Get component index relative to its sub element and
# its sub element
(component_index, sub_element) = ufl_element.extract_component(component.index_range[0])
# Get offset
if len(component_index) == 0:
offset = 0
else:
offset = component.index_range[0] - component_index[0]
# MER: Need to handle mappings in special ways if gdim
# != tdim and some Piolas are present. This could
# probably be merged with the offset code above, but I
# was not able to wrap my head around the offsets
# always referring to component.index_range[0].
if (gdim != tdim):
assert len(component.index_range) == 1, \
"Component transform not implemented for this case. Please request this feature."
component, offset = transform_component(component.index_range[0], offset, ufl_element)
component = MonomialIndex(index_type=MonomialIndex.FIXED,
index_range=[component],
index_id=None)
components = [component, ]
# Add transforms where appropriate
if mapping == "affine":
pass
elif mapping == "contravariant piola":
# phi(x) = (det J)^{-1} J Phi(X)
index0 = component
index1 = MonomialIndex(index_range=list(range(tdim))) + offset
transform = MonomialTransform(index0, index1,
MonomialTransform.J,
f.restriction, offset)
self.transforms.append(transform)
determinant = MonomialDeterminant(power=-1,
restriction=f.restriction)
self.determinants.append(determinant)
components[0] = index1
elif mapping == "covariant piola":
# phi(x) = J^{-T} Phi(X)
index0 = MonomialIndex(index_range=list(range(tdim))) + offset
index1 = component
transform = MonomialTransform(index0, index1,
MonomialTransform.JINV,
f.restriction, offset)
self.transforms.append(transform)
components[0] = index0
else:
raise MonomialException("Don't know how to handle mapping='%s'." % mapping)
# Extract derivatives / transforms
derivatives = []
for d in f.derivatives:
index0 = MonomialIndex(index_range=list(range(tdim)))
if d in index_map:
index1 = index_map[d]
elif isinstance(d, FixedIndex):
index1 = MonomialIndex(index_type=MonomialIndex.FIXED,
index_range=[int(d)],
index_id=int(d))
else:
index1 = MonomialIndex(index_range=list(range(gdim)))
index_map[d] = index1
transform = MonomialTransform(index0, index1,
MonomialTransform.JINV,
f.restriction, 0)
self.transforms.append(transform)
derivatives.append(index0)
# Extract restriction
restriction = f.restriction
# Create basis function
v = MonomialArgument(ufl_element, vindex, components, derivatives,
restriction)
self.arguments.append(v)
# Figure out secondary and auxiliary indices
internal_indices = self._extract_internal_indices(None)
external_indices = self._extract_external_indices(None)
for i in internal_indices + external_indices:
# Skip already visited indices
if i.index_type is not None:
continue
# Set index type and id
num_internal = len([j for j in internal_indices if j == i])
num_external = len([j for j in external_indices if j == i])
if num_internal == 1 and num_external == 1:
i.index_type = MonomialIndex.SECONDARY
i.index_id = _next_secondary_index()
elif num_internal == 2 and num_external == 0:
i.index_type = MonomialIndex.INTERNAL
i.index_id = _next_internal_index()
elif num_internal == 0 and num_external == 2:
i.index_type = MonomialIndex.EXTERNAL
i.index_id = _next_external_index()
else:
raise Exception("Summation index does not appear exactly twice: %s" % str(i))
[docs] def extract_unique_indices(self, index_type=None):
"Return all unique indices for monomial w.r.t. type and id (not range)."
indices = []
for index in self._extract_indices(index_type):
if index not in indices:
indices.append(index)
return indices
def _extract_components(self, f, index_map, vdim):
"Return list of components."
components = []
for c in f.components:
if c in index_map:
index = index_map[c]
elif isinstance(c, FixedIndex):
# Map component using component map from UFL.
# KBO: Is this the right place to add, and do we only have
# scalar components in the tensor representation at this stage
# in the representation?
ufl_element = f.element()
comp_map, comp_num = build_component_numbering(ufl_element.value_shape(), ufl_element.symmetry())
comp = comp_map[(int(c),)]
index = MonomialIndex(index_type=MonomialIndex.FIXED,
index_range=[comp],
index_id=None)
else:
index = MonomialIndex(index_range=list(range(vdim)))
index_map[c] = index
components.append(index)
return components
def _extract_internal_indices(self, index_type=None):
"Return list of indices appearing inside integral."
indices = []
for v in self.arguments:
indices += [v.index] + v.components + v.derivatives
return [i for i in indices if i.index_type == index_type]
def _extract_external_indices(self, index_type=None):
"Return list of indices appearing outside integral."
indices = [c.index for c in self.coefficients] + \
[t.index0 for t in self.transforms] + \
[t.index1 for t in self.transforms]
return [i for i in indices if i.index_type == index_type]
def _extract_indices(self, index_type=None):
"Return all indices for monomial."
return self._extract_internal_indices(index_type) + \
self._extract_external_indices(index_type)
def __str__(self):
"Return informal string representation (pretty-print)."
factors = []
if self.float_value != 1.0:
factors.append(self.float_value)
factors += self.determinants
factors += self.coefficients
factors += self.transforms
return " * ".join([str(f) for f in factors]) + " | " + " * ".join([str(v) for v in self.arguments])
# Index counters
_current_secondary_index = 0
_current_internal_index = 0
_current_external_index = 0
def _next_secondary_index():
"Return next available secondary index."
global _current_secondary_index
_current_secondary_index += 1
return _current_secondary_index - 1
def _next_internal_index():
"Return next available internal index."
global _current_internal_index
_current_internal_index += 1
return _current_internal_index - 1
def _next_external_index():
"Return next available external index."
global _current_external_index
_current_external_index += 1
return _current_external_index - 1
def _reset_indices():
"Reset all index counters."
global _current_secondary_index
global _current_internal_index
global _current_external_index
_current_secondary_index = 0
_current_internal_index = 0
_current_external_index = 0