tensorgenerator.py

finite element library for mathematic majored research

"Code generator for tensor representation"

__author__ = "Anders Logg ([e-mail])"
__date__ = "2004-11-03 -- 2007-06-11"
__copyright__ = "Copyright (C) 2004-2007 Anders Logg"
__license__  = "GNU GPL version 3 or any later version"

# Modified by Kristian B. Oelgaard 2007
# Modified by Marie Rognes ([e-mail]) 2007

# Python modules
from sets import Set

# FFC common modules
from ffc.common.constants import *

# FFC language modules
from ffc.compiler.language.index import *

# FFC code generation common modules
from ffc.compiler.codegeneration.common.codegenerator import *

# FFC format modules
from ffc.compiler.format.removeunused import *

class TensorGenerator(CodeGenerator):
    "Code generator for for tensor representation"

    def __init__(self):
        "Constructor"

        # Initialize common code generator
        CodeGenerator.__init__(self)

    def generate_cell_integral(self, form_data, form_representation, sub_domain, format):
        """Generate dictionary of code for cell integral from the given
        form representation according to the given format"""

        # Extract terms
        terms = form_representation.cell_tensor
        if len(terms) == 0:
            return None

        # Generate element code + set of used geometry terms
        element_code, geo_set = self.__generate_element_tensor(terms, format)

        # Generate geometry code + set of used coefficients + set of jacobi terms
        geo_code, coeff_set, trans_set = self.__generate_geometry_tensors(terms, geo_set, format)

        # Generate code for manipulating coefficients
        coeff_code = self.__generate_coefficients(terms, coeff_set, format) 

        # Get Jacobian snippet
        jacobi_code = [format["generate jacobian"](form_data.cell_dimension, Integral.CELL)]

        # Remove unused declarations
        code = self.__remove_unused(jacobi_code, trans_set, format)

        # Add coefficient and geometry tensor declarations
        code += coeff_code + geo_code

        # Add element code
        code += [""] + [format["comment"]("Compute element tensor")]
        code += element_code

        return {"tabulate_tensor": code, "members":""}

    def generate_exterior_facet_integral(self, form_data, form_representation, sub_domain, format):
        """Generate dictionary of code for exterior facet integral from the given
        form representation according to the given format"""

        # Extract terms
        terms = form_representation.exterior_facet_tensors
        if len(terms) == 0:
            return None

        num_facets = len(terms)
        cases = [None for i in range(num_facets)]

        # Generate element code + set of used geometry terms
        geo_set = Set()
        for i in range(num_facets):
            case, g_set = self.__generate_element_tensor(terms[i], format)
            cases[i] = case
            geo_set = geo_set | g_set

        # Generate code for geometry tensor (should be the same so pick first)
        # Generate set of used coefficients + set of jacobi terms
        geo_code, coeff_set, trans_set = self.__generate_geometry_tensors(terms[0], geo_set, format)

        # Generate code for manipulating coefficients (should be the same so pick first)
        coeff_code = self.__generate_coefficients(terms[0], coeff_set, format)

        # Get Jacobian snippet
        jacobi_code = [format["generate jacobian"](form_data.cell_dimension, Integral.EXTERIOR_FACET)]

        # Remove unused declarations
        code = self.__remove_unused(jacobi_code, trans_set, format)

        # Add coefficient and geometry tensor declarations
        code += coeff_code + geo_code

        # Add element code
        code += [""] + [format["comment"]("Compute element tensor for all facets")]

        return {"tabulate_tensor": (code, cases), "members":""}
    
    def generate_interior_facet_integral(self, form_data, form_representation, sub_domain, format):
        """Generate dictionary of code for interior facet integral from the given
        form representation according to the given format"""

        # Extract terms
        terms = form_representation.interior_facet_tensors
        if len(terms) == 0:
            return None

        num_facets = len(terms)
        cases = [[None for j in range(num_facets)] for i in range(num_facets)]

        # Generate element code + set of used geometry terms
        geo_set = Set()
        for i in range(num_facets):
            for j in range(num_facets):
                case, g_set = self.__generate_element_tensor(terms[i][j], format)
                cases[i][j] = case
                geo_set = geo_set | g_set

        # Generate code for geometry tensor (should be the same so pick first)
        # Generate set of used coefficients + set of jacobi terms
        geo_code, coeff_set, trans_set = self.__generate_geometry_tensors(terms[0][0], geo_set, format)

        # Generate code for manipulating coefficients (should be the same so pick first)
        coeff_code = self.__generate_coefficients(terms[0][0], coeff_set, format)

        # Get Jacobian snippet
        jacobi_code = [format["generate jacobian"](form_data.cell_dimension, Integral.INTERIOR_FACET)]

        # Remove unused declarations
        code = self.__remove_unused(jacobi_code, trans_set, format)

        # Add coefficient and geometry tensor declarations
        code += coeff_code + geo_code

        # Add element code
        code += [""] + [format["comment"]("Compute element tensor for all facet-facet combinations")]

        return {"tabulate_tensor": (code, cases), "members":""}

    def __generate_coefficients(self, terms, coeff_set, format):
        "Generate code for manipulating coefficients"

        # Generate code as a list of declarations
        code = []

        # Add comment
        code += [format["comment"]("Compute coefficients")]

        # A coefficient is identified by 4 numbers:
        #
        #   0 - the number of the function
        #   1 - the position of the (factored) monomial it appears in
        #   2 - the position of the coefficient inside the monomial
        #   3 - the position of the expansion coefficient

        # Iterate over all terms
        j = 0
        for term in terms:
            for G in term.G:
                for k in range(len(G.coefficients)):
                    coefficient = G.coefficients[k]
                    index = coefficient.n0.index
                    if term.monomial.integral.type == Integral.INTERIOR_FACET:
                        space_dimension = 2*len(coefficient.index.range)
                    else:
                        space_dimension = len(coefficient.index.range)
                    for l in range(space_dimension):
                        # If coefficient is not used don't declare it
                        if not format["modified coefficient access"](index, j, k, l) in coeff_set:
                            continue
                        name = format["modified coefficient declaration"](index, j, k, l)
                        value = format["coefficient"](index, l)
                        for l in range(len(coefficient.ops)):
                            op = coefficient.ops[len(coefficient.ops) - 1 - l]
                            if op == Operators.INVERSE:
                                value = format["inverse"](value)
                            elif op == Operators.MODULUS:
                                value = format["absolute value"](value)
                            elif op == Operators.SQRT:
                                value = format["sqrt"](value)
                        code += [(name, value)]
                j += 1

        # Don't add code if there are no coefficients
        if len(code) == 1:
            return []

        # Add newline
        code += [""]
        return code

    def __generate_geometry_tensors(self, terms, geo_set, format):
        "Generate list of declarations for computation of geometry tensors"

        # Generate code as a list of declarations
        code = []    
        
        # Add comment