simplify.py

finite element library for mathematic majored research

                if indices[1].has_key(str(nv[label])):
                    indices[1][str(nv[label])] += [(label, nv[label])]
                else:
                    indices[1][str(nv[label])] = [(label, nv[label])]

    # We can contract the indices and thus the monomials, if there is
    # only a difference of one index and if there are at least two
    # occurances of this index. (Summation over _repeated_ indices.)
    if len(indices[0].keys()) == len(indices[1].keys()) == 1 \
           and len(indices[0].values()[0]) == 2:
        # Constructing the new index:
        i0 = indices[0].values()[0][0][1]
        i1 = indices[1].values()[0][0][1]

        # We only want to contract if each index value occurs once.
        common_indices = Set(i0.range) & Set(i1.range)
        if not common_indices:
            # Constucting the new monomial based on the old m:
            index = Index(i0) + Index(i1)
            q = Monomial(m)
            for (label, i) in indices[0][str(i0)]:
                s = "q.%s = index" % str(label)
                exec(s)
            return [q]

    return [m, n]

def simplify_monomial(monomial):
    """ Simpliy monomials with construction of the form:
    (dx_j/dX_i)(dX_l/dx_j) | (d/dX_l) => (d/dX_i)"""
    
    for basis in monomial.basisfunctions:
        success = 0
        first = None
        second = None
        
        # The derivatives of this basis function is a good
        # starting point so we run through these. We use index
        # notation since we may have to replace some of them.
        for i in range(len(basis.derivatives)):
            derivative = basis.derivatives[i]
            therestriction = basis.restriction
            success = 0
            theindex = derivative.index
            # Now, lets run through the transforms and see whether
            # there are two matching:
            for transform in monomial.transforms:
                if transform.type == Transform.JINV:
                    if (not cmp(transform.index0, theindex)
                        and transform.restriction == therestriction):
                        first = transform
                        break
            # If there are no matching JINV-transforms, no hope of success.
            if not first: break

            for transform in monomial.transforms:
                if transform.type == Transform.J:
                    if (not cmp(transform.index1, first.index1)
                        and transform.restriction == therestriction):
                        second = transform
                        success = 1
                        break
                        
            if success == 1:
                # Now, we should first remove the transforms from
                # the transform list. Second: replace the old
                # derivative with the new one.
                basis.derivatives[i] = Derivative(derivative.element,
                                                  second.index0)
                monomial.transforms.remove(first) 
                monomial.transforms.remove(second)
    return monomial

def diff(m, n, key = None):
    """ Take two elements and returns the difference between these in
    an appropriate manner.""" 

    # Dictionaries containing each version when m and n are different
    mversion = {}
    nversion = {}

    if isinstance(m, list) and isinstance(n, list):
        #The difference of two lists (of equal length!)
        diffs = []
        if len(m) == len(n):
            for i in range(len(m)):
                difference = diff(m[i], n[i], key[i])
                if not is_empty(difference):
                    diffs += [difference]
            return diffs
        else:
            raise FormError("Only know how to diff lists of equal length.")
        
    elif isinstance(m, Monomial) and isinstance(n, Monomial):
        # The difference between two monomials
        coeffids = ["coefficients[%d]" % i for i in range(len(m.coefficients))]
        coeffdiff = diff(m.coefficients, n.coefficients, coeffids)
        tids = ["transforms[%d]" % i for i in range(len(m.transforms))]
        tdiff = diff(m.transforms, n.transforms, tids)
        dictionary = {'coefficients': coeffdiff, 'transforms': tdiff}

        # Treat the basis functions items separately:
        bids = ["basisfunctions[%d]." % i
                for i in range(len(m.basisfunctions))]
        bdiff = diff(m.basisfunctions, n.basisfunctions, bids)
        for dict in bdiff:
            for key in dict:
                if dictionary.has_key(key):
                    dictionary[key] += dict[key]
                else:
                    dictionary[key] = dict[key]
        return abbreviate(dictionary)
        
    elif isinstance(m, BasisFunction) and isinstance(n, BasisFunction):
        # The difference between two basis functions
        idiff = diff(m.index, n.index, key + "index")
        cids = [key + "component[%d]" % i for i in range(len(m.component))]
        cdiff = diff(m.component, n.component, cids)
        dids = [key + "derivatives[%d]" % i for i in range(len(m.derivatives))]
        ddiff = diff(m.derivatives, n.derivatives, dids)
        return abbreviate({'index': idiff, 'component': cdiff,
                           'derivatives': ddiff})

    elif isinstance(m, Index) and isinstance(n, Index):
        # The difference between two indices
        if not m == n:
            mversion[key] = m
            nversion[key] = n

    elif isinstance(m, Transform) and isinstance(n, Transform):
        # The index difference between two transforms
        if not m.index0 == n.index0:
            mversion[key + ".index0"] = m.index0
            nversion[key + ".index0"] = n.index0
        if not m.index1 == n.index1:
            mversion[key + ".index1"] = m.index1
            nversion[key + ".index1"] = n.index1
        
    elif isinstance(m, Derivative) and isinstance(n, Derivative):
        # The index difference between two derivatives:
        if not m.index == n.index:
            mversion[key + ".index"] = m.index
            nversion[key + ".index"] = n.index

    elif isinstance(m, Coefficient) and isinstance(n, Coefficient):
        # The index difference between two coefficients
        if not m.index == n.index:
            mversion[key + ".index"] = m.index
            nversion[key + ".index"] = n.index
        
    else:
        raise FormError, "Diff is not implemented between such elements"

    if mversion:
        return [mversion, nversion]
    else:
        return []

def restriction_exterior(form):
    """Removing restrictions on exterior facets.
       If this makes monomial terms equal, all but one monomial term will be deleted."""

    # List of numbers of monomials with removed restrictions on exterior facets
    removed_restrictions = []
    for i in range(len(form.monomials)):
        p = form.monomials[i]
        for v in p.basisfunctions:
            type = p.integral.type
            if type == Integral.EXTERIOR_FACET:
                if not (v.restriction == None or v.restriction == Restriction.CONSTANT):
                    # Remove restriction and keep track of the monomial number
                    v.restriction = None
                    removed_restrictions += [i]

    # Create a set of the removed restrictions, (remove duplicate numbers) and get monomials
    removed_restrictions = tuple(set(removed_restrictions))
    monomials = [form.monomials[r] for r in removed_restrictions]

    # If any restrictions were moved
    if monomials:
        # The first monomial is always unique
        unique = [monomials[0]]
        for i in range(1,len(monomials)):
            p = monomials[i]
            equals = False
            # Check if monomial already exists
            for p0 in unique:
                if contraction_likely(p,p0):
                    if not diff(p,p0):
                        # If there are no differences the monmials are equal
                        equals = True

            # If monomial is redundant, remove it. Otherwise add it to the list of uniuqe monomials
            if equals:
                form.monomials.remove(p)
            else:
                unique += [p]

def change_constant_restrictions(form):
    """Change Restriction.CONSTANT to Restriction.PLUS on interior facets"""

    for m in form.monomials:
        if m.integral.type == Integral.INTERIOR_FACET:
            for v in m.basisfunctions:
                if v.restriction == Restriction.CONSTANT:
                    v.restriction = Restriction.PLUS