evaluatebasisderivatives.py

finite element library for mathematic majored research

        for j in range(poly_dim):
            code += [(Indent.indent(format_float + format_coeff(i) + format_secondary_index(j)),\
                      format_floating_point(0.0))]
    code += [""]

    code += [Indent.indent(format_comment("Declare new coefficients"))]
    for i in range(num_components):
        for j in range(poly_dim):
            code += [(Indent.indent(format_float + format_new_coeff(i) + format_secondary_index(j)),\
                      format_floating_point(0.0))]
    code += [""]

    code += [Indent.indent(format_comment("Loop possible derivatives"))]
    code += [Indent.indent(format_loop("deriv_num", 0, format_num_derivatives))]
    code += [Indent.indent(format_block_begin)]
    # Increase indentation
    Indent.increase()

    code += [Indent.indent(format_comment("Get values from coefficients array"))]
    for i in range(num_components):
        for j in range(poly_dim):
            code += [(Indent.indent(format_new_coeff(i) + format_secondary_index(j)),\
                      format_coefficients(i) + format_matrix_access(format_dof, j))]
    code += [""]

    code += [Indent.indent(format_comment("Loop derivative order"))]
    code += [Indent.indent(format_loop("j", 0, format_n))]
    code += [Indent.indent(format_block_begin)]
    # Increase indentation
    Indent.increase()

    # Update old coefficients
    code += [Indent.indent(format_comment("Update old coefficients"))]
    for i in range(num_components):
        for j in range(poly_dim):
            code += [(Indent.indent(format_coeff(i) + format_secondary_index(j)),\
                      format_new_coeff(i) + format_secondary_index(j))]
    code += [""]

    # Update new coefficients
    code += multiply_coeffs(element, Indent, format)

    # Decrease indentation
    Indent.decrease()
    code += [Indent.indent(format_block_end)]

    # Compute derivatives on reference element
    code += [Indent.indent(format_comment\
    ("Compute derivatives on reference element as dot product of coefficients and basisvalues"))]
    mapping = pick_first([element.value_mapping(dim) for dim in range(element.value_dimension(0))])
    if mapping == Mapping.PIOLA:
        code += [Indent.indent(format["comment"]("Correct values by the Piola transform"))]
    value_code = []
    for i in range(num_components):
        if (i == 0):
            name = format_derivatives + format_array_access("deriv_num")
        elif (i == 1):
            name = format_derivatives + format_array_access(format_add([format_num_derivatives, "deriv_num"] ))
        else:
            name = format_derivatives + format_array_access(format_add([format_multiply\
                   (["%d" %i, format_num_derivatives]), "deriv_num"] ))

        value = format_add([ format_multiply([format_new_coeff(i) + format_secondary_index(k),\
                             format_basisvalue(k)]) for k in range(poly_dim) ])

        # Use Piola transform to map basisfunctions back to physical element if needed
        if mapping == Mapping.PIOLA:
            value_code.insert(i,(Indent.indent(format_tmp(0, i)), value))
            basis_col = [format_tmp_access(0, j) for j in range(element.cell_dimension())]
            jacobian_row = [format["transform"](Transform.J, j, i, None) for j in range(element.cell_dimension())]
            inner = [format_multiply([jacobian_row[j], basis_col[j]]) for j in range(element.cell_dimension())]
            sum = format_group(format_add(inner))
            value = format_multiply([format_inv(format_det(None)), sum])

        value_code += [(Indent.indent(name), value)]
    code += value_code
    # Decrease indentation
    Indent.decrease()
    code += [Indent.indent(format_block_end)]

    return code + [""]

def multiply_coeffs(element, Indent, format):
    "Auxilliary function that multiplies coefficients with directional derivatives."

    code = []

    # Prefetch formats to speed up code generation
    format_if               = format["if"]
    format_group            = format["grouping"]
    format_isequal          = format["is equal"]
    format_combinations     = format["derivative combinations"]
    format_block_begin      = format["block begin"]
    format_block_end        = format["block end"]
    format_coeff            = format["coefficient scalar"]
    format_new_coeff        = format["new coefficient scalar"]
    format_secondary_index  = format["secondary index"]
    format_add              = format["add"]
    format_multiply         = format["multiply"]
    format_dmats            = format["dmats table"]
    format_matrix_access    = format["matrix access"]

    # Get number of components, must change for tensor valued elements
    num_components = element.value_dimension(0)

    # Get polynomial dimension of basis
    poly_dim = len(element.basis().base.bs)

    # Get the shape of the element
    cell_shape = element.cell_shape()

    for i in range(cell_shape):

# not language-independent
        code += [Indent.indent(format_if + format_group( format_combinations + \
                 format_matrix_access("deriv_num","j") + format_isequal + "%d" %(i) ) )]
        code += [Indent.indent(format_block_begin)]
        # Increase indentation
        Indent.increase()
        for j in range(num_components):
            for k in range(poly_dim):
                name = format_new_coeff(j) + format_secondary_index(k)
                value = format_add( [format_multiply([format_coeff(j) + format_secondary_index(l),\
                        format_dmats(i) + format_matrix_access(l,k)]) for l in range(poly_dim)])
                code += [(Indent.indent(name), value)]

        # Decrease indentation
        Indent.decrease()
        code += [Indent.indent(format_block_end)]

    return code + [""]

def transform_derivatives(element, sum_value_dim, Indent, format):
    """This function computes the value of the basisfunction as the dot product of the
    coefficients and basisvalues """

    code = []

    # Prefetch formats to speed up code generation
    format_loop             = format["loop"]
    format_num_derivatives  = format["num derivatives"]
    format_derivatives      = format["reference derivatives"]
    format_values           = format["argument values"] 
    format_multiply         = format["multiply"]
    format_add              = format["add"]
    format_array_access     = format["array access"]
    format_matrix_access    = format["matrix access"]
    format_transform        = format["transform matrix"]

    code += [Indent.indent(format["comment"]("Transform derivatives back to physical element"))]
    
    code += [Indent.indent(format_loop("row", 0, format_num_derivatives))]
    code += [Indent.indent(format["block begin"])]
    # Increase indentation
    Indent.increase()
    code += [Indent.indent(format_loop("col", 0, format_num_derivatives))]
    code += [Indent.indent(format["block begin"])]
    # Increase indentation
    Indent.increase()

    # Get number of components, must change for tensor valued elements
    num_components = element.value_dimension(0)

    # Compute offset in array values if any
    for i in range(num_components):
        if (sum_value_dim + i == 0):
            name = format_values + format_array_access("row")
        elif (sum_value_dim + i == 1):
            name = format_values + format_array_access(format_add([format_num_derivatives, "row"]))
        else:
            offset_name = format_multiply(["%d" %(sum_value_dim + i), format_num_derivatives])
            name = format_values + format_array_access(format_add([offset_name, "row"]))
        if (i == 0):
            value = format_multiply([format_transform + format_matrix_access("row","col"),\
                    format_derivatives + format_array_access("col")])
        elif (i == 1):
            value = format_multiply([format_transform + format_matrix_access("row","col"),\
                    format_derivatives + format_array_access(format_add([format_num_derivatives, "col"]))])
        else:
            offset_value = format_multiply(["%d" %i, format_num_derivatives])

            value = format_multiply([format_transform + format_matrix_access("row","col"),\
                    format_derivatives + format_array_access(format_add([offset_value, "col"]))])

        # Compute values
        code += [Indent.indent(format["add equal"](name, value))]

    # Decrease indentation
    Indent.decrease()
    code += [Indent.indent(format["block end"])]

    # Decrease indentation
    Indent.decrease()

    code += [Indent.indent(format["block end"])]

    return code

def delete_pointers(element, Indent, format):
    "Delete the pointers to arrays."

    code = []

    # Delete pointers
    code += [Indent.indent(format["comment"]("Delete pointer to array of derivatives on FIAT element"))]
    code += [Indent.indent(format["delete"] + format["array access"]("") + " " +\
                           format["reference derivatives"] + format["end line"])] + [""]

    code += [Indent.indent(format["comment"]("Delete pointer to array of combinations of derivatives and transform"))]

    code += [Indent.indent(format["loop"]("row", 0, format["num derivatives"]))]
    code += [Indent.indent(format["block begin"])]
    # Increase indentation
    Indent.increase()

    code += [Indent.indent(format["delete"] + format["array access"]("") + " " +\
                           format["derivative combinations"] + format["array access"]("row") + format["end line"])]

    code += [Indent.indent(format["delete"] + format["array access"]("") + " " +\
                           format["transform matrix"] + format["array access"]("row") + format["end line"])]

    # Decrease indentation
    Indent.decrease()
    code += [Indent.indent(format["block end"])] + [""]

    code += [Indent.indent(format["delete"] + format["array access"]("") + " " +\
                           format["derivative combinations"] + format["end line"])]

    code += [Indent.indent(format["delete"] + format["array access"]("") + " " +\
                           format["transform matrix"] + format["end line"])]

    return code


def debug_combinations(element, Indent, format):

    code = []

    code += [format["comment"]("Debug code")]
    # Debug code
    code += [Indent.indent('std::cout << "%s = " << std::endl;' % format["derivative combinations"])]
    code += [Indent.indent(format["loop"]("j", 0, format["num derivatives"]))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent(format["loop"]("k", 0, format["argument derivative order"]))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent('std::cout << %s << " ";') % (format["derivative combinations"] + "[j][k]")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]
    code += [Indent.indent("std::cout << std::endl;")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]

    return code

def debug_transform(element, Indent, format):

    code = []

    cell_shape = element.cell_shape()
    num_derivatives = format["num derivatives"]
    Jinv = format["transform Jinv"]
    transform = format["transform matrix"]

    # Debug code
    code += [format["comment"]("Debug code")]
    code += [Indent.indent("std::cout.precision(8);")]

    # Jinv
    code += [Indent.indent('std::cout << "%s = " << std::endl;' % Jinv)]
    code += [Indent.indent(format["loop"]("j", 0, cell_shape))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent(format["loop"]("k", 0, cell_shape))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent('std::cout << %s << " ";') % (Jinv + "[j][k]")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]
    code += [Indent.indent("std::cout << std::endl;")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]

    # Transform matrix
    code += [Indent.indent('std::cout << "%s = " << std::endl;' % transform)]
    code += [Indent.indent(format["loop"]("j", 0, num_derivatives))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent(format["loop"]("k", 0, num_derivatives))]
    code += [Indent.indent(format["block begin"])]
    # Increase indent
    Indent.increase()
    code += [Indent.indent('std::cout << %s << " ";') % (transform + "[j][k]")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]
    code += [Indent.indent("std::cout << std::endl;")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent(format["block end"])]

    return code

def debug_reference_derivatives(element, Indent, format):

    # Debug code
    code = [Indent.indent('std::cout << "%s = " << std::endl;' %format["reference derivatives"])]
    code += [Indent.indent(format["loop"]("j", 0, format["num derivatives"]))]
    # Increase indent
    Indent.increase()
    code += [Indent.indent('std::cout << %s << " ";') % (format["reference derivatives"] + "[j]")]

    # Decrease indent
    Indent.decrease()
    code += [Indent.indent("std::cout << std::endl;")]
    return code

def debug_():
    "misc"

    code = []

    # Debug coefficients
#    value = "std::cout << "
#    for i in range(num_components):
#        for j in range(poly_dim):
#            value += 'new_coeff%d_%d << " " << ' % (i,j)
#    value += "std::endl;"
#    code += [value]
#    code += [""]

            # Debug coefficients
#            value = "std::cout << "
#            for i in range(num_components):
#                for j in range(poly_dim):
#                    value += 'new_coeff%d_%d << " " << ' % (i,j)
#            value += "std::endl;"
#            code += [value]

    return code