~chaffra/ufl/main-old

"Basic algebra operations."

 

#

# This file is part of UFL.

#

 

class Sum(AlgebraOperator):

    __slots__ = ("_operands",)

    def __new__(cls, *operands): # TODO: This whole thing seems a bit complicated... Can it be simplified? Maybe we can merge some loops for efficiency?

        ufl_assert(operands, "Can't take sum of nothing.")

        #if not operands:

        #    return Zero() # Allowing this leads to zeros with invalid type information in other places, need indices and shape

        # make sure everything is an Expr

        operands = [as_ufl(o) for o in operands]

        # Got one operand only? Do nothing then.

        if len(operands) == 1:

            return operands[0]

        # assert consistent tensor properties

        sh = operands[0].shape()

        fi = operands[0].free_indices()

            error("Shape mismatch in Sum.")

        if any((set(fi) ^ set(o.free_indices())) for o in operands[1:]):

            error("Can't add expressions with different free indices.")

        # sort operands in a canonical order

        operands = sorted(operands, cmp=cmp_expr)

        # purge zeros

        operands = [o for o in operands if not isinstance(o, Zero)]

        # sort scalars to beginning and merge them

        scalars = [o for o in operands if isinstance(o, ScalarValue)]

        if scalars:

                operands = nonscalars

            else:

                operands = [f] + nonscalars

        if not operands:

            return Zero(sh, fi, fid)

        # left with one operand only?

        if len(operands) == 1:

            return operands[0]

        # Replace n-repeated operands foo with n*foo

        newoperands = []

        op = operands[0]

                op = o

                n = 1

        operands = newoperands

        # left with one operand only?

        if len(operands) == 1:

            return operands[0]

        # construct and initialize a new Sum object

        self = AlgebraOperator.__new__(cls)

        self._init(*operands)

 

    def _init(self, *operands):

        self._operands = operands

    def __init__(self, *operands):

        AlgebraOperator.__init__(self)

    def operands(self):

        return self._operands

    def free_indices(self):

        return self._operands[0].free_indices()

    def index_dimensions(self):

        return self._operands[0].index_dimensions()

    def shape(self):

        return self._operands[0].shape()

    def evaluate(self, x, mapping, component, index_values):

        return sum(o.evaluate(x, mapping, component, index_values) for o in self.operands())

    def __str__(self):

        ops = [parstr(o, self) for o in self._operands]

        if False:

class Product(AlgebraOperator):

    """The product of two or more UFL objects."""

    __slots__ = ("_operands", "_free_indices", "_index_dimensions",)

    def __new__(cls, *operands):

        # Make sure everything is an Expr

        operands = [as_ufl(o) for o in operands]

        # Make sure everything is scalar

        #ufl_assert(not any(o.shape() for o in operands),

        #    "Product can only represent products of scalars.")

        if any(o.shape() for o in operands):

            error("Product can only represent products of scalars.")

        # No operands? Return one.

        if not operands:

            return IntValue(1)

        # Got one operand only? Just return it.

        if len(operands) == 1:

            return operands[0]

        # Got any zeros? Return zero.

        if any(isinstance(o, Zero) for o in operands):

            free_indices     = unique_indices(tuple(chain(*(o.free_indices() for o in operands))))

            index_dimensions = subdict(mergedicts([o.index_dimensions() for o in operands]), free_indices)

            return Zero((), free_indices, index_dimensions)

        # Merge scalars, but keep nonscalars sorted

        scalars = []

        nonscalars = []

                    return nonscalars[0]

            else:

                scalars = [p]

        # Sort operands in a canonical order (NB! This is fragile! Small changes here can have large effects.)

        operands = scalars + sorted(nonscalars, cmp=cmp_expr)

        # Replace n-repeated operands foo with foo**n

        newoperands = []

        op, nop = operands[0], 1

                # Reset op as o

                op, nop = o, 1

        operands = newoperands

        # Left with one operand only after simplifications?

        if len(operands) == 1:

            return operands[0]

        # Construct and initialize a new Product object

        self = AlgebraOperator.__new__(cls)

        self._init(*operands)

        return self

    def _init(self, *operands):

        "Constructor, called by __new__ with already checked arguments."

        # Store basic properties

        self._operands = operands

        # Extract indices

        self._free_indices     = unique_indices(tuple(chain(*(o.free_indices() for o in operands))))

        self._index_dimensions = mergedicts([o.index_dimensions() for o in operands]) or EmptyDict

    def __init__(self, *operands):

        AlgebraOperator.__init__(self)

    def operands(self):

        return self._operands

    def free_indices(self):

        return self._free_indices

    def index_dimensions(self):

        return self._index_dimensions

    def shape(self):

        return ()

    def evaluate(self, x, mapping, component, index_values):

        ops = self.operands()

        sh = self.shape()

        for o in ops:

            tmp *= o.evaluate(x, mapping, (), index_values)

        return tmp

    def __str__(self):

        ops = [parstr(o, self) for o in self._operands]

        if False:

            return s

        # Implementation with no line splitting:

        return "%s" % " * ".join(ops)

    def __repr__(self):

        return "Product(%s)" % ", ".join(repr(o) for o in self._operands)

 

class Division(AlgebraOperator):

    __slots__ = ("_a", "_b",)

    def __new__(cls, a, b):

        a = as_ufl(a)

        b = as_ufl(b)

        self = AlgebraOperator.__new__(cls)

        self._init(a, b)

        return self

    def _init(self, a, b):

        #ufl_assert(isinstance(a, Expr) and isinstance(b, Expr), "Expecting Expr instances.")

        if not (isinstance(a, Expr) and isinstance(b, Expr)):

 

    def __init__(self, a, b):

        AlgebraOperator.__init__(self)

    def operands(self):

        return (self._a, self._b)

    def free_indices(self):

        return self._a.free_indices()

    def index_dimensions(self):

        return self._a.index_dimensions()

    def shape(self):

        return () # self._a.shape()

        a, b = self.operands()

        a = a.evaluate(x, mapping, component, index_values)

        b = b.evaluate(x, mapping, component, index_values)

 

class Power(AlgebraOperator):

    __slots__ = ("_a", "_b",)

    def __new__(cls, a, b):

        a = as_ufl(a)

        b = as_ufl(b)

        if not is_true_ufl_scalar(a): error("Cannot take the power of a non-scalar expression.")

        if not is_true_ufl_scalar(b): error("Cannot raise an expression to a non-scalar power.")

        if isinstance(a, ScalarValue) and isinstance(b, ScalarValue):

            return as_ufl(a._value ** b._value)

        if b == 1:

            return a

        if b == 0:

            return IntValue(1)

        # construct and initialize a new Power object

        self = AlgebraOperator.__new__(cls)

        self._init(a, b)

        return self

    def _init(self, a, b):

        #ufl_assert(isinstance(a, Expr) and isinstance(b, Expr), "Expecting Expr instances.")

        if not (isinstance(a, Expr) and isinstance(b, Expr)):

 

    def __init__(self, a, b):

        AlgebraOperator.__init__(self)

    def operands(self):

        return (self._a, self._b)

    def free_indices(self):

        return self._a.free_indices()

    def index_dimensions(self):

        return self._a.index_dimensions()

    def shape(self):

        return ()

        a, b = self.operands()

        a = a.evaluate(x, mapping, component, index_values)

        b = b.evaluate(x, mapping, component, index_values)

        return a**b

    def __str__(self):

        return "%s ** %s" % (parstr(self._a, self), parstr(self._b, self))

 

 

class Abs(AlgebraOperator):

    __slots__ = ("_a",)

    def __init__(self, a):

        AlgebraOperator.__init__(self)

        ufl_assert(isinstance(a, Expr), "Expecting Expr instance.")

        if not isinstance(a, Expr): error("Expecting Expr instances.")

        self._a = a

    def operands(self):

        return (self._a, )

    def free_indices(self):

        return self._a.free_indices()

    def index_dimensions(self):

        return self._a.index_dimensions()

    def shape(self):

        return self._a.shape()

        a = self._a.evaluate(x, mapping, component, index_values)

        return abs(a)

    def __str__(self):

        return "| %s |" % parstr(self._a, self)