quadraturegenerator.py

finite element library for mathematic majored research

    def __tabulate_weights(self, weights, tensor_number, Indent, format):
        "Generate table of quadrature weights"

        code = []
        
        # Prefetch formats to speed up code generation
        format_floating_point = format["floating point"]

        code += [Indent.indent(format["comment"]\
                ("Array of quadrature weights (tensor/monomial term %d)" %(tensor_number,) ))]

        # Create variable name
        name = format["table declaration"] + format["weights"](tensor_number, str(len(weights)))
        value = format["block"](format["separator"].join([format_floating_point(w)\
                 for w in weights]))

        code += [(Indent.indent(name), value)]

        return code + [""]

    def __tabulate_psis(self, tensors, Indent, format, tables):
        "Tabulate values of basis functions and their derivatives at quadrature points"

        code = []

        format_psis = format["psis"]
        format_floating_point = format["floating point"]
        format_block_begin = format["block begin"]
        format_block_end = format["block end"]

        if tables:
            for names in tables:
                # Get value and save as an array
                vals = tables[names]
                # Generate array of values (FIAT returns [dof, quad_points] transpose to [quad_points, dof])
                value = tabulate_matrix(vals, format)
                code += [(Indent.indent(names), Indent.indent(value))]# + [""]
        else:
            for tensor_number in range(len(tensors)):
                tensor = tensors[tensor_number]
                tables = get_names_tables(tensor, tensor_number, format)

                code += [Indent.indent(format["comment"]\
                        ("Values of shapefunctions and their derivatives at quadrature points"))]
                code += [Indent.indent(format["comment"]\
                        ("Format: [quadrature points][dofs] (tensor/monomial term %d)" % (tensor_number,)))]
                for names in tables:
                    # Get value and save as an array
                    vals = tables[names]
                    # Generate array of values (FIAT returns [dof, quad_points] transpose to [quad_points, dof])
                    value = tabulate_matrix(vals, format)
                    code += [(Indent.indent(names), Indent.indent(value))]# + [""]
        return code

    def __reset_element_tensor(self, tensor, Indent, format):
        "Reset the entries of the local element tensor"

        code = []

        # Get rank and dims of primary indices
        irank = tensor.i.rank
        idims = tensor.i.dims

        # Get monomial and compute macro dimensions in case of restricted basisfunctions
        monomial = tensor.monomial
        macro_idims = compute_macro_idims(monomial, idims, irank)

        # Generate value
        value = format["floating point"](0.0)

        # FIXME: quadrature only support Functionals and Linear and Bilinear forms
        if (irank == 0):
            code += [Indent.indent(format["comment"]("Reset value"))]

            # Generate entry and name
            entry = "0"
            name =  format["element tensor quad"] + format["array access"](entry)
            code += [(Indent.indent(name),value)]

        elif (irank == 1):
            code += [Indent.indent(format["comment"]\
                    ("Reset values of the element tensor block"))]

            # Generate entry and name
            entry = format["first free index"]
            name =  format["element tensor quad"] + format["array access"](entry)

            # Create boundaries for loop
            boundaries = [0, macro_idims[0]]
            loop_vars = [[format["first free index"]] + boundaries]
            code += generate_loop(name, value, loop_vars, Indent, format)

        elif (irank == 2):
            code += [Indent.indent(format["comment"]\
                    ("Reset values of the element tensor block"))]

            # Generate entry and name
            entry = format["add"]([format["multiply"]([format["first free index"], "%d" %macro_idims[1]]),\
                                                       format["second free index"]])
            name =  format["element tensor quad"] + format["array access"](entry)

            # Create boundaries for loop
            boundaries = [[0, macro_idims[0]], [0, macro_idims[1]]]
            loop_vars = [[format["first free index"]] + boundaries[0],\
                         [format["second free index"]] + boundaries[1]]

            code += generate_loop(name, value, loop_vars, Indent, format)
        else:
            raise RuntimeError, "Quadrature only support Functionals and Linear and Bilinear forms"

        return code + [""]

    def __element_tensor(self, tensor, tensor_number, Indent, format, name_map):
        "Generate loop over primary indices"

        code = []
        s_set = Set()

        # Prefetch formats to speed up code generation
        format_add      = format["add"]
        format_multiply = format["multiply"]
        format_group    = format["grouping"]
        format_new_line = format["new line"]

        dic_indices = {0:format["first free index"], 1:format["second free index"]}

        # Get rank and dims of primary indices
        irank, idims = tensor.i.rank, tensor.i.dims

        # Get monomial and compute macro dimensions in case of restricted basisfunctions
        monomial = tensor.monomial
        macro_idims = compute_macro_idims(monomial, idims, irank)

        # Get Psi indices, list of primary and secondary indices e.g. [[i0, a0], [i1, a1]]
        indices = [psi[1] for psi in tensor.Psis]
        vindices = [psi[2] for psi in tensor.Psis]

        # Get secondary and auxiliary multiindices
        aindices, b0indices, bgindices = tensor.a.indices, tensor.b0.indices, tensor.bg.indices

        # Compute scaling
        weight = [format["weights"](tensor_number, format["integration points"])]

        # Choose level of optimisation
        if self.optimise_level == 0:
            values, secondary_loop, t_set = values_level_0(indices, vindices, aindices, b0indices,\
                                            bgindices, tensor, tensor_number, weight, format, name_map)
        elif self.optimise_level == 1:
            values, secondary_loop, t_set = values_level_1(indices, vindices, aindices, b0indices,\
                                            bgindices, tensor, tensor_number, weight, format, name_map)
        elif self.optimise_level == 2:
            values, secondary_loop, t_set = values_level_2(indices, vindices, aindices, b0indices,\
                                            bgindices, tensor, tensor_number, weight, format, name_map)
        elif self.optimise_level == 3:
            values, secondary_loop, t_set = values_level_3(indices, vindices, aindices, b0indices,\
                                            bgindices, tensor, tensor_number, weight, format, name_map)
        else:
            raise RuntimeError, "Optimisation level not implemented!"

        value = format_add(values)

        # FIXME: quadrature only support Functionals and Linear and Bilinear forms
        if (irank == 0):

            # Entry is zero because functional is a scalar value
            entry = "0"
            # Generate name
            name =  format["element tensor quad"] + format["array access"](entry)
            code += [Indent.indent(format["comment"]\
                    ("Compute value (tensor/monomial term %d)" % (tensor_number,)))]

            # Create boundaries for loop
            loop_vars = secondary_loop
            if secondary_loop:
                code += generate_loop(name, value, loop_vars, Indent, format, format["add equal"])
            else:
                code += [format["add equal"](Indent.indent(name), value)]

        elif (irank == 1):
            # Generate entry
            for i in range(irank):
                for v in monomial.basisfunctions:
                    if v.index.type == Index.PRIMARY and v.index.index == i:
                        if v.restriction == Restriction.MINUS:
                            entry = format_add([dic_indices[i], str(idims[0])])
                        else:
                            entry = dic_indices[i]
                        break

            # Generate name
            name =  format["element tensor quad"] + format["array access"](entry)
            code += [Indent.indent(format["comment"]\
                    ("Compute block entries (tensor/monomial term %d)" % (tensor_number,)))]

            # Create boundaries for loop
            loop_vars = secondary_loop
            if secondary_loop:
                code += generate_loop(name, value, loop_vars, Indent, format, format["add equal"])
            else:
                code += [format["add equal"](Indent.indent(name), value)]

        elif (irank == 2):
            entry = []
            # Generate entry
            for i in range(irank):
                for v in monomial.basisfunctions:
                    if v.index.type == Index.PRIMARY and v.index.index == i:
                        if v.restriction == Restriction.MINUS:
                            entry += [format["grouping"](format_add([dic_indices[i], str(idims[i])]))]
                        else:
                            entry += [dic_indices[i]]
                        break

            entry[0] = format_multiply([entry[0], str(macro_idims[1])])
            name =  format["element tensor quad"] + format["array access"](format_add(entry))

            code += [Indent.indent(format["comment"]\
                    ("Compute block entries (tensor/monomial term %d)" % (tensor_number,)))]

            # Create boundaries for loop
            loop_vars = secondary_loop
            if secondary_loop:
                code += generate_loop(name, value, loop_vars, Indent, format, format["add equal"])
            else:
                code += [format["add equal"](Indent.indent(name), value)]

        else:
            raise RuntimeError, "Quadrature only support Functionals and Linear and Bilinear forms"

        return (code, t_set)

    def __remove_unused(self, code, set, format):

        if code:
            # Generate body of code, using the format
            lines = format["generate body"](code)

            # Generate auxiliary code line that uses all members of the set (to trick remove_unused)
            line_set = format["add equal"]("A", format["multiply"](set))
            lines += "\n" + line_set

            # Remove unused Jacobi declarations
            code = remove_unused(lines)

            # Delete auxiliary line
            code = code.replace("\n" + line_set, "")

            return [code]
        else:
            return code