ufcformat.py

finite element library for mathematic majored research

"Code generation for the UFC 1.0 format"

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

# Modified by Kristian B. Oelgaard 2007

# UFC code templates
from ufc import *

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

# FFC language modules
from ffc.compiler.language.restriction import *
from ffc.compiler.language.tokens import *
from ffc.compiler.language.integral import *

# FFC format modules
from codesnippets import *
from removeunused import *

# Choose map from restriction
choose_map = {Restriction.PLUS: "0", Restriction.MINUS: "1", Restriction.CONSTANT: "0", None: ""}
transform_options = {Transform.JINV: lambda m, j, k: "Jinv%s_%d%d" % (m, j, k),
                     Transform.J: lambda m, j, k: "J%s_%d%d" % (m, k, j)}
# Options for the printing q or 1.0/(q) for q string:
power_options = {True: lambda q: q, False: lambda q: "1.0/(%s)" % q}

# Specify formatting for code generation
format = { "add": lambda v: " + ".join(v),
           "subtract": lambda v: " - ".join(v),
           "multiply": lambda v: "*".join(v),
           "grouping": lambda v: "(%s)" % v,
           "block": lambda v: "{%s}" % v,
           "block begin": "{",
           "block end": "}",
           "separator": ", ",
           "return": lambda v: "return %s;" % v,
           "bool": lambda v: {True: "true", False: "false"}[v],
           "floating point": lambda v: "<not defined>",
           "epsilon": "<not defined>",
           "tmp declaration": lambda j, k: "const double " + format["tmp access"](j, k),
           "tmp access": lambda j, k: "tmp%d_%d" % (j, k),
           "comment": lambda v: "// %s" % v,
           "exception": lambda v: "throw std::runtime_error(\"%s\");" % v,
           "determinant": lambda r: "detJ%s" % choose_map[r],
           "scale factor": "det",
           "power": lambda base, exp: power_options[exp >= 0](format["multiply"]([str(base)]*abs(exp))),
           "constant": lambda j: "c%d" % j,
           "coefficient table": lambda j, k: "w[%d][%d]" % (j, k),
           "coefficient": lambda j, k: "w[%d][%d]" % (j, k),
           "coeff": "w",
           "modified coefficient declaration": lambda i, j, k, l: "const double c%d_%d_%d_%d" % (i, j, k, l),
           "modified coefficient access": lambda i, j, k, l: "c%d_%d_%d_%d" % (i, j, k, l),
           "transform": lambda type, j, k, r: "%s" % (transform_options[type](choose_map[r], j, k)),
           "reference tensor" : lambda j, i, a: None,
           "geometry tensor declaration": lambda j, a: "const double " + format["geometry tensor access"](j, a),
           "geometry tensor access": lambda j, a: "G%d_%s" % (j, "_".join(["%d" % index for index in a])),
           "element tensor": lambda i, k: "A[%d]" % k,
           "sign tensor": lambda type, i, k: "S%s%s_%d" % (type, i, k),
           "sign tensor declaration": lambda s: "const int " + s,
           "signs": "S",
           "vertex values": lambda i: "vertex_values[%d]" % i,
           "dof values": lambda i: "dof_values[%d]" % i,
           "dofs": lambda i: "dofs[%d]" % i,
           "entity index": lambda d, i: "c.entity_indices[%d][%d]" % (d, i),
           "num entities": lambda dim : "m.num_entities[%d]" % dim,
           "offset declaration": "unsigned int offset",
           "offset access": "offset",
           "cell shape": lambda i: {1: "ufc::line", 2: "ufc::triangle", 3: "ufc::tetrahedron"}[i],
# quadrature, evalutate_basis(), evaluate_basis_derivatives()
# declarations
           "float declaration": "double ",
           "const float declaration": "const double ",
           "static float declaration": "static double ",
           "uint declaration": "unsigned int ",
           "const uint declaration": "const unsigned int ",
           "static uint declaration": "static unsigned int ",
           "table declaration": "const static double ",
# access
           "array access": lambda i: "[%s]" %(i),
           "matrix access": lambda i,j: "[%s][%s]" %(i,j),
           "secondary index": lambda i: "_%s" %(i),
# program flow
           "dof map if": lambda i,j: "if (%d <= %s && %s <= %d)" %(i,\
                         format["argument basis num"], format["argument basis num"], j),
           "loop": lambda i,j,k: "for (unsigned int %s = %s; %s < %s; %s++)"% (i, j, i, k, i),
           "if": "if",
# operators
           "times equal": lambda i,j: "%s *= %s;" %(i,j),
           "add equal": lambda i,j: "%s += %s;" % (i,j),
           "is equal": " == ",
           "inverse": lambda v: "(1.0/%s)" % v,
           "absolute value": lambda v: "std::abs(%s)" % v,
           "sqrt": lambda v: "std::sqrt(%s)" % v,
# variable names
           "element tensor quad": "A",
           "integration points": "ip",
           "first free index": "j",
           "second free index": "k",
           "free secondary indices":["r","s","t","u"],
           "derivatives": lambda i,j,k,l: "dNdx%d_%d[%s][%s]" % (i,j,k,l),
           "element coordinates": lambda i,j: "x[%s][%s]" % (i,j),
           "weights": lambda i,j: "W%d[%s]" % (i,j),
           "psis": "P",
           "argument coordinates": "coordinates",
           "argument values": "values",
           "argument basis num": "i",
           "argument derivative order": "n",
           "local dof": "dof",
           "x coordinate": "x",
           "y coordinate": "y",
           "z coordinate": "z",
           "scalings": lambda i,j: "scalings_%s_%d" %(i,j),
           "coefficients table": lambda i: "coefficients%d" %(i),
           "dmats table": lambda i: "dmats%d" %(i),
           "coefficient scalar": lambda i: "coeff%d" %(i),
           "new coefficient scalar": lambda i: "new_coeff%d" %(i),
           "psitilde_a": "psitilde_a",
           "psitilde_bs": lambda i: "psitilde_bs_%d" %(i),
           "psitilde_cs": lambda i,j: "psitilde_cs_%d%d" %(i,j),
           "basisvalue": lambda i: "basisvalue%d" %(i),
           "num derivatives": "num_derivatives",
           "reference derivatives": "derivatives",
           "derivative combinations": "combinations",
           "transform matrix": "transform",
           "transform Jinv": "Jinv",
# snippets
           "coordinate map": lambda i: {2:map_coordinates_2D, 3:map_coordinates_3D}[i],
           "facet sign": lambda e: "sign_facet%d" % e,
           "snippet facet signs": lambda d: eval("facet_sign_snippet_%dD" % d),
           "snippet dof map": evaluate_basis_dof_map,
           "snippet eta_triangle": eta_triangle_snippet,
           "snippet eta_tetrahedron": eta_tetrahedron_snippet,
           "snippet jacobian": lambda d: eval("jacobian_%dD" % d),
           "snippet combinations": combinations_snippet,
           "snippet transform2D": transform2D_snippet,
           "snippet transform3D": transform3D_snippet,
#           "snippet inverse 2D": inverse_jacobian_2D,
#           "snippet inverse 3D": inverse_jacobian_3D,
           "snippet evaluate_dof": lambda d : {2: evaluate_dof_2D, 3: evaluate_dof_3D}[d],
           "get cell vertices" : "const double * const * x = c.coordinates;",
           "generate jacobian": lambda d,i: __generate_jacobian(d,i),
           "generate body": lambda d: __generate_body(d),
# misc
           "block separator": ",\n",
#           "block separator": ",",
           "new line": "\\\n",
           "end line": ";",
           "space": " ",
           "pointer": "*",
           "new": "new ",
           "delete": "delete "}

def init(options):
    "Initialize code generation for given options"

    # Set number of digits for floating point and machine precision
    format_string = "%%.%dg" % eval(options["precision="])
    format["floating point"] = lambda v : format_string % v
    format["epsilon"] = 10.0*eval("1e-%s" % options["precision="])

def write(generated_forms, prefix, options):
    "Generate UFC 1.0 code for a given list of pregenerated forms"
    debug("Generating code for UFC 1.0")

    # Generate code for header
    output = ""
    output += generate_header(prefix, options)
    output += "\n"

    # Generate UFC code
    output += generate_ufc(generated_forms, prefix, options)

    # Generate code for footer
    output += generate_footer(prefix, options)

    # Write file
    filename = "%s.h" % prefix
    file = open(filename, "w")
    file.write(output)
    file.close()
    debug("Output written to " + filename)

def generate_header(prefix, options):
    "Generate code for header"

    # Check if BLAS is required
    if options["blas"]:
        blas_include = "\n#include <cblas.h>"
        blas_warning = "\n// Warning: This code was generated with '-f blas' and requires cblas.h."
    else:
        blas_include = ""
        blas_warning = ""
        
    return """\
// This code conforms with the UFC specification version 1.0
// and was automatically generated by FFC version %s.%s

#ifndef __%s_H
#define __%s_H

#include <cmath>
#include <stdexcept>
#include <ufc.h>%s
""" % (FFC_VERSION, blas_warning, prefix.upper(), prefix.upper(), blas_include)

def generate_footer(prefix, options):
    "Generate code for footer"
    return """\
#endif
"""

def generate_ufc(generated_forms, prefix, options):
    "Generate code for body"

    output = ""
    
    # Iterate over forms
    for i in range(len(generated_forms)):

        # Get pregenerated code, form data and prefix
        (form_code, form_data) = generated_forms[i]
        form_prefix = compute_prefix(prefix, generated_forms, i)

        # Generate code for ufc::finite_element(s)
        for (label, sub_element) in form_code["finite_elements"]:
            output += __generate_finite_element(sub_element, form_data, options, form_prefix, label)
            output += "\n"

        # Generate code for ufc::dof_map(s)
        for (label, sub_dof_map) in form_code["dof_maps"]:
            output += __generate_dof_map(sub_dof_map, form_data, options, form_prefix, label)
            output += "\n"

        # Generate code for ufc::cell_integral
        for j in range(form_data.num_cell_integrals):
            output += __generate_cell_integral(form_code[("cell_integral", j)], form_data, options, form_prefix, j)
            output += "\n"

        # Generate code for ufc::exterior_facet_integral
        for j in range(form_data.num_exterior_facet_integrals):
            output += __generate_exterior_facet_integral(form_code[("exterior_facet_integral", j)], form_data, options, form_prefix, j)
            output += "\n"
    
        # Generate code for ufc::interior_facet_integral
        for j in range(form_data.num_interior_facet_integrals):
            output += __generate_interior_facet_integral(form_code[("interior_facet_integral", j)], form_data, options, form_prefix, j)
            output += "\n"

        # Generate code for ufc::form
        if "form" in form_code:
            output += __generate_form(form_code["form"], form_data, options, form_prefix)
            output += "\n"

    return output

def compute_prefix(prefix, generated_forms, i):
    "Compute prefix for form i"

    # Get form ranks
    ranks = [form_data.rank for (form_code, form_data) in generated_forms]

    # Return prefixFunctional, prefixLinearForm or prefixBilinearForm
    # when we have exactly one form of ranks 0, 1 or 2
    count = [ranks.count(0), ranks.count(1), ranks.count(2)]
    if len(ranks) <= 3 and sum(count) > 0 and min(count) >= 0 and max(count) <= 1:
        postfixes = ["Functional", "LinearForm", "BilinearForm"]
        return "%s%s" % (prefix, postfixes[ranks[i]])

    # Return prefix_i if we have more than one rank
    if len(ranks) > 1:
        return "%s_%d" % (prefix, i)

    # Else, just return prefix
    return prefix

def __generate_finite_element(code, form_data, options, prefix, label):
    "Generate code for ufc::finite_element"

    ufc_code = {}

    # Set class name
    ufc_code["classname"] = "%s_finite_element_%s" % (prefix, "_".join([str(i) for i in label]))

    # Generate code for members
    ufc_code["members"] = ""

    # Generate code for constructor
    ufc_code["constructor"] = "// Do nothing"

    # Generate code for destructor
    ufc_code["destructor"] = "// Do nothing"

    # Generate code for signature
    ufc_code["signature"] = "return \"%s\";" % code["signature"]

    # Generate code for cell_shape
    ufc_code["cell_shape"] = "return %s;" % code["cell_shape"]
    
    # Generate code for space_dimension
    ufc_code["space_dimension"] = "return %s;" % code["space_dimension"]

    # Generate code for value_rank
    ufc_code["value_rank"] = "return %s;" % code["value_rank"]

    # Generate code for value_dimension
    cases = ["return %s;" % case for case in code["value_dimension"]]