simplify.py

finite element library for mathematic majored research

_author__ = "Marie Rognes ([e-mail])"
__date__ = "2006-10-23 -- 2007-08-30"
__copyright__ = "Copyright (C) 2006"
__license__  = "GNU GPL version 3 or any later version"

# Modified by Anders Logg 2007
# Modified by Kristian Oelgaard 2007

# Python modules
import sys
from sets import Set

# FFC common modules
from ffc.common.debug import *
from ffc.common.exceptions import *

# FFC compiler.language modules
from ffc.compiler.language.index import *
from ffc.compiler.language.algebra import *
from ffc.compiler.language.tokens import *
from ffc.compiler.language.integral import *

def simplify(form):
    """ Simplification of a form"""

    if not isinstance(form, Form):
        raise FormError, "simplify assumes a Form as input."

    debug("Simplifying form...")
    # Handle restrictions on exterior facets before simplifying
    restriction_exterior(form)

    # Change constant restrictions on interior facets before simplifying
    change_constant_restrictions(form)

    previous = str(form)
    simplified = ""

    while(previous != simplified):
        reassign_indices(form)
        previous = str(form)
        simplify_form(form)
        simplified = str(form)
    reassign_indices(form)
    debug("done")

def simplify_form(f):
    # First: contract indices and factorize
    f.monomials = contract_list(contract_monomials, f.monomials)
    reassign_indices(f)
    #f.monomials = contract_list(factorize_monomials, f.monomials)

    for monomial in f.monomials:
        # Second, remove monomials with numeric = 0.0:
        if monomial.numeric == 0.0:
            f.monomials.remove(monomial)
            continue
        # Third, contract determinants:
        # FIX: Should move this.
        if monomial.determinants:
            monomial.determinants = contract_list(contract_determinants,
                                                  monomial.determinants)
        
        # Fourth, simplify each monomial with regard to derivatives.
        monomial = simplify_monomial(monomial)

        
def contract_list(contraction, monomials):
    """ Given a list of ..., run contraction on (all) pairs of these,
    and return the new list. contraction should return a tuple
    (result, contracted) where result is a list containing either the
    contracted result or the original input.""" 
    # meg: Again, please replace this if there is an easier/prettier
    # way of doing this.
    if len(monomials) < 2:
        return monomials 
    current = listcopy(monomials)
    i = 0
    while i < len(current):
        j = i+1
        while j < len(current):
            (q, contracted) = contraction(current[i], current[j])
            if contracted:
                current.remove(current[j])
                current[i] = q[0]      # (Note that j > i.)
            else:
                j +=1
        i += 1
    return current

def factorize_monomials(m, n):
    """ Given a two monomials, factorize any common factors and return
    the new monomials."""
    # meg: Not quite finished yet.
    if contraction_likely(m,n):
        differences = diff(m, n)
        if is_empty(differences):
            q = Monomial(m)
            q.numeric += n.numeric
            return ([q], True)
    return ([m, n], False)
    
def contract_monomials(m, n):
    if not (isinstance(m, Monomial) and isinstance(n, Monomial)):
        raise FormError, "contract_monomials can only contract monomials"
    # First, check to see if it is at all likely that these monomials
    # are contractable
    if not contraction_likely(m, n):
        return ([m, n], False)
    q = contract_indices(m, n)
    return (q, len(q) == 1)

def contract_determinants(d0, d1):
    if not (isinstance(d0, Determinant) and isinstance(d1, Determinant)):
        raise FormError,"contract_determinants can only contract determinants!"
    # Determinants are contracted by trying to multiply them:
    d = d0*d1
    if d: return ([d], True)
    else: return ([d0, d1], False)

def contraction_likely(m, n):
    """ Given two monomials/basisfunctions, check if they have the
    same numbers of basisfunctions, transforms, derivatives,
    components etc. (This is just intended as a preliminary check.)"""
    # Comparing monomials:
    if isinstance(m, Monomial) and isinstance(n, Monomial):
        if m.integral != n.integral:
            return False

        if m.numeric != n.numeric:
            return False

        if len(m.determinants) != len(n.determinants):
            return False
        for i in range (len(m.determinants)):
            if not m.determinants[i] == n.determinants[i]:
                return False

        if len(m.coefficients) != len(n.coefficients):
            return False
        for i in range(len(m.coefficients)):
            if not contraction_likely(m.coefficients[i], n.coefficients[i]):
                return False

        if len(m.transforms) != len(n.transforms):
            return False
        for i in range(len(m.transforms)):
            if not contraction_likely(m.transforms[i], n.transforms[i]):
                return False

        if len(m.basisfunctions) != len(n.basisfunctions):
            return False
        for i in range(len(m.basisfunctions)):
            if not contraction_likely(m.basisfunctions[i],
                                      n.basisfunctions[i]):
                return False
        return True

    # Comparing basis functions:
    elif isinstance(m, BasisFunction) and isinstance(n, BasisFunction):
        if m.element != n.element:
            return False
        if m.index != n.index:
            return False
        if m.restriction != n.restriction:
            return False
        if len(m.component) != len(n.component):
            return False
        if len(m.derivatives) != len(n.derivatives):
            return False
        return True

    # Comparing transforms:
    elif isinstance(m, Transform) and isinstance(n, Transform):
        if m.type != n.type or m.restriction != n.restriction:
            return False
        return True

    # Comparing coefficients:
    elif isinstance(m, Coefficient) and isinstance(n, Coefficient):
        if m.n0 != n.n0:
            return False
        if m.n1 != n.n1:
            return False
        return True

    # Others, not implemented
    else:
        return True
            
def contract_indices(m, n):
    """ Given two monomials, contract indices in the following way: If
    m and n only differ by one index, contract and replace this index."""
    indices = [{}, {}]

    # Extract the differences between the monomials
    differences = abbreviate(diff(m, n))

    # Extract the different index values:
    for attribute in differences:
        for [mv, nv] in differences[attribute]:
            # Each key/label in mv is also in nv by construction.
            for label in mv:  
                if indices[0].has_key(str(mv[label])):
                    indices[0][str(mv[label])] += [(label, mv[label])]
                else:
                    indices[0][str(mv[label])] = [(label, mv[label])]