quadraturegenerator_utils.py

finite element library for mathematic majored research

        num_dofs = len(vindex.range)

        # Get number of quadrature points
        num_quadrature_points = len(tensor.quadrature.weights)

        # Get lists of secondary indices and auxiliary_0 indices
        aindices = tensor.a.indices
        b0indices = tensor.b0.indices

        names = []
        multi_indices = []

        # Loop secondary and auxiliary indices to generate names and entries in the psi table
        for a in aindices:
            for b in b0indices:

                (name, multi_index) = generate_psi_declaration(tensor_number, indices, vindex,\
                                      a, b, num_quadrature_points, num_dofs, format)

                names += [name]
                multi_indices += [multi_index]

        # Remove redundant names and entries in the psi table
        names, multi_indices = extract_unique(names, multi_indices)

        # Loop names and tabulate psis
        for i in range(len(names)):

            # Get name
            name = names[i]

            # Get values from psi tensor, should have format values[dofs][quad_points]
            vals = values[tuple(multi_indices[i])]

            # Check if the values have the correct dimensions, otherwise quadrature is not correctly
            # implemented for the given form!!
            if numpy.shape(vals) != (num_dofs, num_quadrature_points):
                raise RuntimeError, "Quadrature is not correctly implemented for the given form!"

            # Generate values (FIAT returns [dof, quad_points] transpose to [quad_points, dof])
            value = numpy.transpose(vals)
            tables[name] = value

    return tables

def generate_load_table(tensors):
    "Generate header to load psi tables"

    function = """
void load_table(double A[][%d], const char* name, int m) const
{
  std::ifstream in(name, std::ios::in);
  for (int i = 0; i < m; i++)
    for (int j = 0; j < %d; j++)
      in >> A[i][j];
  in.close();
}\n"""

#    load_int = """
#void load_table(unsigned int& load)
#{
#  std::ifstream in("load.table", std::ios::in);
#  in >> load;
#  in.close();
#}"""
#    save_int = """
#void save()
#{
#  std::ofstream out("load_f%d%d.table", std::ios::out);
#  out << 0;
#  out.close();
#}\n""" %(facet0, facet1)

    group_num_dofs = []

    for tensor in tensors:

        # Get list of psis
        psis = tensor.Psis

        # Loop psis
        for psi_number in range(len(psis)):

            # Get psi
            psi = psis[psi_number]

            # Get the number of dofs
            num_dofs = len(psi[2].range)

            if not num_dofs in group_num_dofs:
                group_num_dofs += [num_dofs]

#    try:
#        file = open("load_table.h", "r")
#        code = file.read()
#        code = code.split("\n")
#        file.close()
#    except IOError:
#        code = False
#    if not code:
        # Create code
#        code = "#include <fstream>"
#        for dof in group_num_dofs:
#            code += function %(dof,dof)
#    else:
#        for dof in group_num_dofs:
#            new_lines = function %(dof,dof)
#            new_lines = new_lines.split("\n")
#            if not new_lines[1] in code:
#                code += new_lines
#        code = "\n".join(code)
    code = []
    for dof in group_num_dofs:
        new_lines = function %(dof,dof)
        new_lines = new_lines.split("\n")
        if not new_lines[1] in code:
            code += new_lines
    return "\n".join(code)

    # Write code
#    file = open("load_table.h", "w")
#    file.write(code)
#    file.close()


def save_psis(tensors, facet0, facet1, Indent, format, tables):
    "Save psis tables instead of tabulating them"

    load_table = "load_table(%s, \"%s\", %d)"
    load_int = "load_table(%s)"
    out = """
std::cout.precision(17);
std::cout << "%s" << std::endl;
for (int j = 0; j < %d; j++)
{
  for (int k = 0; k < %d; k++)
    std::cout << %s[j][k] << " ";
  std::cout << std::endl;
}
std::cout << std::endl;
"""
    code = []

    format_floating_point = format["floating point"]
    format_epsilon        = format["epsilon"]
    format_end_line       = format["end line"]
    format_table          = format["table declaration"]
    format_float          = format["static float declaration"]
#    format_float          = format["const float declaration"]

    # Get list of psis
#    psis = tensor.Psis
    irank = "_i%d" %tensors[0].i.rank
    if (facet0 == None) and (facet1 == None):
        facets = ""
    elif (facet0 != None) and (facet1 == None):
        facets = "_f%d" %(facet0)
    else:
        facets = "_f%d%d" %(facet0, facet1)

    try:
        os.system("mkdir -p tables")

    except IOError:
        dirs = True

    loads = []
#    outputs = []

    if tables:
#        print "save tables"

        for names in tables:
#            print "names: ", names
            # Get value and save as an array
            value = tables[names]
            shape = numpy.shape(value)
            name = names.split(" ")[-1].split("[")[0]
            a = []
            for j in range(shape[0]):
                for k in range(shape[1]):
                    val = value[j][k]
                    if abs(val) < format_epsilon:
                        val = 0.0
                    a += [format_floating_point(val)]
            a = " ".join(a)
            file = open("tables/" + name + facets + irank + ".table", "w")
            file.write(a)
            file.close

            # Generate code to load table
            loads += [load_table %(name, "tables/" + name + facets + irank + ".table", shape[0])]
#            outputs += [out %(name, num_quadrature_points, num_dofs, name)]
            code += [Indent.indent(names.replace(format_table,format_float) + format_end_line)]
    else:
        for tensor_number in range(len(tensors)):
            tensor = tensors[tensor_number]
            tables = get_names_tables(tensor, tensor_number, format)

            for names in tables:
                # Get value and save as an array
                value = tables[names]
                shape = numpy.shape(value)
                name = names.split(" ")[-1].split("[")[0]
                a = []
                for j in range(shape[0]):
                    for k in range(shape[1]):
                        val = value[j][k]
                        if abs(val) < format_epsilon:
                            val = 0.0
                        a += [format_floating_point(val)]
                a = " ".join(a)
                file = open("tables/" + name + facets + irank + ".table", "w")
                file.write(a)
                file.close

                # Generate code to load table
                loads += [load_table %(name, "tables/" + name + facets + irank + ".table", shape[0])]
#                outputs += [out %(name, num_quadrature_points, num_dofs, name)]
                code += [Indent.indent(names.replace(format_table,format_float) + format_end_line)]

    code += [(Indent.indent(format["static uint declaration"] + "load"),"1")]
#    code += [Indent.indent(load_int %("load") + format_end_line)]
    code += [Indent.indent(format["if"] + format["grouping"]("load"))]
    code += [Indent.indent(format["block begin"])]

    load_code = [Indent.indent("std::cout << \"Loading tables\" << std::endl;")]
    for load in loads:
        code += [Indent.indent(load + format_end_line)]
#        load_code += [Indent.indent(load + format_end_line)]

    code += [(Indent.indent("load"),"0")]
#    code += [Indent.indent("save();")]
    code += [Indent.indent(format["block end"])]

#    for output in outputs:
#        code += [Indent.indent(output)]

#    file = open("load.table", "w")
#    file.write("0")
#    file.close

    return code

def unique_tables(tensors, format):
    # Determine if some tensors have the same tables (and same names)

    name_map = {}
    tables = {}

    for tensor_number in range(len(tensors)):

        tensor = tensors[tensor_number]
        table = get_names_tables(tensor, tensor_number, format)
        for name in table:
            tables[name] = table[name]

#    print "tables 2: ", tables
    # Initialise name map
    for name in tables:
        name_map[name] = []

    # Loop all tables to see if some are redundant
    for name0 in tables:
        val0 = numpy.array(tables[name0])
        for name1 in tables:
            # Don't compare values with self
            if not name0 == name1:
                val1 = numpy.array(tables[name1])
                # Check if dimensions match
                if numpy.shape(val0) == numpy.shape(val1):
#                    print "numpy.shape(val0): ", numpy.shape(val0)
#                    print "numpy.shape(val1): ", numpy.shape(val1)
#                    print "val0: ", val0
#                    print "val1: ", val1
                    diff = val1 - val0
                    # Check if values are the same
                    if not diff.any():
                        if name0 in name_map:
                            name_map[name0] += [name1]
                            if name1 in name_map:
                                del name_map[name1]

    inverse_name_map = {}
    for name in name_map:
#        print "name: ", name
        # Get name of psi
        name_strip = name.split()[-1].split("[")[0]
#        print "name_strip: ", name_strip

        maps = name_map[name]
        for m in maps:
#            print "m: ", m
            # Get name of psi
            m_strip = m.split()[-1].split("[")[0]
#            print "m_strip: ", m_strip
            inverse_name_map[m_strip] = name_strip
            del tables[m]
#    print "inverse_name_map: ", inverse_name_map

    return (inverse_name_map, tables)