Subject: Re: [Ffc] [Branch ~ffc-core/ffc/dev] Rev 1494: Somewhat evil and dead-slow (but very useful!) testing of generated

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 20 Jan 2010 22:52:43 +0100 (CET)

2010/1/20  <[email-masked]>:
> ------------------------------------------------------------
> revno: 1494
> committer: Marie E. Rognes <[email-masked]>
> branch nick: ffc-unstable
> timestamp: Wed 2010-01-20 21:43:48 +0100
> message:
>  Somewhat evil and dead-slow (but very useful!) testing of generated
>  element and dofmap code. A lot of stuff fails, including:
>
>  (a) element_map fails for mixed elements (-> AL)
>  (b) generated evaluate_basis_derivatives code does not compile (-> KBO)

I just pushed a bug fix that should? fix this.

When I try to run the tests in verify_element_code I get:
Traceback (most recent call last):
  File "test.py", line 16, in <module>
    from REFERENCE import references
ImportError: No module named REFERENCE

Did you forget to add the file?

Kristian

>  (c) BDM/RT on tetrahedra needs inspection (-> MER)
> modified:
>  test/verify_element_code/cppcode.py
>  test/verify_element_code/generate_reference.py
>  test/verify_element_code/test.py
>
>
> --
> lp:~ffc-core/ffc/dev
> https://code.launchpad.net/~ffc-core/ffc/dev
>
> You are subscribed to branch lp:~ffc-core/ffc/dev.
> To unsubscribe from this branch go to https://code.launchpad.net/~ffc-core/ffc/dev/+edit-subscription.
>
> === modified file 'test/verify_element_code/cppcode.py'
> --- test/verify_element_code/cppcode.py 2010-01-18 13:37:16 +0000
> +++ test/verify_element_code/cppcode.py 2010-01-20 20:43:48 +0000
> @@ -5,39 +5,104 @@
>  __author__ = "Marie E. Rognes ([email-masked])"
>  __license__  = "GNU GPL version 3 or any later version"
>
> +import numpy
>
>  # Last modified: 15.01.2010
>
> -def element_code(element):
> -
> -    ir = {}
> -
> -    ir["space_dimension"] = _space_dimension()
> -    ir["value_dimension"] = _value_dimension()
> -    #ir["evaluate_basis"] = _evaluate_basis(element)
> -    #ir["evaluate_basis_derivatives"] = _evaluate_basis_derivatives(element)
> -    #ir["evaluate_dof"] = _evaluate_dof(element)
> -
> -    return ir
> -
> -def dofmap_code(element):
> -
> -    ir = {}
> -    #ir["tabulate_dofs"] = _tabulate_dofs(element)
> -    return ir
> +def test_code(element):
> +
> +    test = {}
> +
> +    # Element stuff
> +    test["space_dimension"] = _space_dimension()
> +    test["value_dimension"] = _value_dimension()
> +    #test["evaluate_basis"] = _evaluate_basis(element)
> +    #test["evaluate_basis_derivatives"] = _evaluate_basis_derivatives(element)
> +    test["evaluate_dof"] = _evaluate_dof(element)
> +    test["interpolate_vertex_values"] = _interpolate_vertex_values(element)
> +    #test["num_sub_elements"] = _num_sub_elements(element)
> +
> +    # Dofmap stuff
> +    test["num_facet_dofs"] = _num_facet_dofs(element)
> +    #test["num_entity_dofs"] = _num_entity_dofs(element)
> +    test["tabulate_dofs"] = _tabulate_dofs(element)
> +    test["tabulate_facet_dofs"] = _tabulate_facet_dofs(element)
> +    test["tabulate_coordinates"] = _tabulate_coordinates(element)
> +    return test
> +
> +
> +# --------- general utilities ---------------
> +
> +def extract_value_dimension(element):
> +    try:
> +        N = element.value_shape()[0]
> +    except:
> +        N = 1
> +    return N
> +
> +def extract_cell_dimension(element):
> +    try:
> +        return element.cell().geometric_dimension()
> +    except:
> +        return element.ref_el.get_spatial_dimension()
> +
> +
> +def _mesh(d):
> +    if d == 1:
> +        return "dolfin::UnitInterval m(5);\nm.init();\ndolfin::UFCMesh mesh(m);"
> +    elif d == 2:
> +        return "dolfin::UnitSquare m(5, 7);\nm.init();\ndolfin::UFCMesh mesh(m);"
> +    elif d == 3:
> +        return "dolfin::UnitCube m(5, 6, 3);\nm.init();\ndolfin::UFCMesh mesh(m);"
> +
> +def _cell(i=None):
> +    i = i or 3
> +    return "dolfin::UFCCell cell(dolfin::Cell(m, %s));" % i
> +
> +def _function(n):
> +
> +    values = ["values[%d] = sin(%d*x[0]);" % (i, i+1) for i in range(n)]
> +    if n == 1:
> +        value_declaration = ""
> +    else:
> +        value_declaration = "public:\n  Source() : Expression(%d) {}\n\n" % n
> +
> +    code = """
> +class Source : public dolfin::Expression
> +{
> +  %s
> +
> +  void eval(dolfin::Array<double>& values, const dolfin::Array<const double>& x) const
> +  {
> +    %s
> +   }
> +};
> +
> +Source f;
> +""" % (value_declaration, "\n".join(values))
> +    return code
> +
> +dolfin_includes = """
> +#include <dolfin/common/Array.h>
> +#include <dolfin/mesh/Cell.h>
> +#include <dolfin/mesh/UnitInterval.h>
> +#include <dolfin/mesh/UnitSquare.h>
> +#include <dolfin/mesh/UnitCube.h>
> +#include <dolfin/function/Expression.h>
> +#include <dolfin/fem/UFC.h>
> +"""
>
>  # -------- element ------------------
>
> -
>  template = """\
>  #include <iostream>
> -#include <ufc.h>
> +%s
>  #include "element.h"
>
>  int main()
>  {
>
> -element_finite_element_0 element;
> +element_0_finite_element_0 element;
>
>  %s
>
> @@ -45,70 +110,98 @@
>  }
>  """
>
> -cell_template = """
> -// Initialize ufc cell
> -ufc::cell c;
> -"""
> -
> -function_template = """
> -// Initialize ufc function
> -ufc::function f;
> -"""
> -
>  def _space_dimension():
> -    return template % """\
> -    int dim = element.space_dimension();
> -    std::cout << dim << std::endl;
> -    """
> +    "Test code for space_dimension."
> +    code = """\
> +int dim = element.space_dimension();
> +std::cout << dim << std::endl;"""
> +    return template % ("", code)
>
>  def _value_dimension():
> -    return template % """\
> -    int dim = element.value_dimension(0);
> -    std::cout << dim << std::endl;
> -    """
> +    "Test code for value_dimension."
> +
> +    code = """\
> +int dim = element.value_dimension(0);
> +std::cout << dim << std::endl;"""
> +    return template % ("", code)
>
>  def _evaluate_basis(element):
> -    return template % ""
> +    "Test code for evaluate_basis."
> +
> +    return template % ("", code)
>
>  def _evaluate_basis_derivatives(element):
> -    return template % ""
> +    "Test code for evaluate_basis_derivatives."
> +
> +    return template % ("", code)
>
>  def _evaluate_dof(element):
> -    " c++ code for testing evaluate_dof"
> -
> -    # Declarations
> -    code = [cell_template, function_template]
> -    code += ["""
> -    // Initialize result:
> -    double result;
> -
> -    // Initialize iterator
> -    unsigned int i;
> -    """]
> +    "Test code for evaluate_dof."
> +
> +    # Initialize cell
> +    d = extract_cell_dimension(element)
> +    code = [_mesh(d), _cell()]
> +
> +    # Initialize function f
> +    n = extract_value_dimension(element)
> +    code += [_function(n)]
>
>     # Evaluate_dof code
> +    code += ["double result;"]
>     dofs_code = ["""\
> -    // Evaluate dof number %(i)d
> -    i = %(i)d;
> -    result = element.evaluate_dof(i, f, c);
> -    std::cout << result << std::endl;
> -    """ % {"i": i} for i in range(element.space_dimension())]
> +result = element.evaluate_dof(%d, f, cell);
> +std::cout << result << std::endl;""" % i
> +                 for i in range(element.space_dimension())]
>
>     code += dofs_code
> -    return template % "".join(code)
> -
> +    return template % (dolfin_includes, "\n".join(code))
> +
> +
> +def _interpolate_vertex_values(element):
> +    "Test code for interpolate_vertex_values."
> +
> +    N = extract_value_dimension(element)
> +    d = extract_cell_dimension(element)
> +    n = element.space_dimension()
> +
> +    # Initialize default mesh and cell
> +    code = [_mesh(d), _cell()]
> +
> +    # Initialize dof values
> +    dof_values = ["%d" % i for i in range(n)]
> +    code += ["double dof_values[%d] = {%s};" % (n, ",".join(dof_values))]
> +
> +    # Declare vertex_values:
> +    code += ["double vertex_values[%d];" % (N*(d+1))]
> +
> +    # Call interpolate
> +    func = "element.interpolate_vertex_values"
> +    code += [func + "(vertex_values, dof_values, cell);"]
> +
> +    # Print resulting vertex values
> +    for i in range(N*(d+1)):
> +        code += ["std::cout << vertex_values[%d] << std::endl;" % i]
> +
> +    return template % (dolfin_includes, "\n".join(code))
> +
> +def _num_sub_elements(element):
> +    "Test code for num_sub_elements."
> +
> +    code = ["int num_elements = element.num_sub_elements();"]
> +    code += ["std::cout << num_elements << std::endl;"]
> +    return template % ("", "\n".join(code))
>
>  # -------- dof map ------------------
>
>  dof_template = """\
>  #include <iostream>
> -#include <ufc.h>
> +%s
>  #include "element.h"
>
>  int main()
>  {
>
> -element_dof_map_0 dofmap;
> +element_0_dof_map_0 dofmap;
>
>  %s
>
> @@ -116,14 +209,89 @@
>  }
>  """
>
> +
> +def _num_facet_dofs(element):
> +    "Test code for num_facet_dofs."
> +
> +    code = "unsigned int result = dofmap.num_facet_dofs();"
> +    return dof_template % ("", code)
> +
> +def _num_entity_dofs(element):
> +    "Test code for num_entity_dofs."
> +
> +    d = extract_cell_dimension(element)
> +
> +    code = ["unsigned int result;"]
> +    for i in range(d):
> +        code += ["result = dofmap.num_entity_dofs(%d);" % i]
> +        code += ["std::cout << result << std::endl;"]
> +    return dof_template % ("", "\n".join(code))
> +
> +
>  def _tabulate_dofs(element):
> -    return dof_template % """\
> -    //unsigned int* dofs;
> -    //const ufc::mesh m;
> -
> -    ufc::cell cell;
> -    cell.coordinates = new double * [3]
> -    """
> -
> -
> -
> +    "Test code for tabulate_dofs."
> +
> +    d = extract_cell_dimension(element)
> +    code = []
> +
> +    # Initialize default mesh and cell
> +    code += [_mesh(d), _cell()]
> +
> +    # Declare dof array
> +    n = element.space_dimension()
> +    code += ["unsigned int dofs[%d];" % n]
> +
> +    # Tabulate dofs
> +    code += ["std::cout << 1 << std::endl;"]
> +    code += ["dofmap.tabulate_dofs(dofs, mesh, cell);"]
> +
> +    # Print dofs
> +    code += ["std::cout << dofs[%d] << std::endl;" % i for i in range(n)]
> +
> +    return dof_template % (dolfin_includes, "\n".join(code))
> +
> +def _tabulate_coordinates(element):
> +    "Test code for tabulate_coordinates."
> +
> +    n = element.space_dimension()
> +    d = extract_cell_dimension(element)
> +    code = []
> +
> +    # Initialize default mesh and cell
> +    code += [_mesh(d), _cell()]
> +
> +    # Initialize coords
> +    code += ["double** coords = new double*[%d];" % n]
> +    for i in range(n):
> +        code += ["coords[%d] = new double[%d];" % (i, d)]
> +
> +    # Tabulate coordinates and print
> +    code += ["dofmap.tabulate_coordinates(coords, cell);"]
> +    for i in range(n):
> +        for j in range(d):
> +            code += ["std::cout << coords[%d][%d] << std::endl;" % (i, j)]
> +
> +    return dof_template % (dolfin_includes, "\n".join(code))
> +
> +def _tabulate_facet_dofs(element):
> +    "Test code for tabulate_facet_dofs."
> +
> +    code = []
> +
> +    # Map (spatial dimension - 1) to number of facets
> +    facets = [2, 3, 6]
> +
> +    # Declare dofs
> +    code += ["int n = dofmap.num_facet_dofs();"]
> +    code += ["unsigned int* dofs = new unsigned int[n];"]
> +
> +    D = extract_cell_dimension(element) - 1
> +    for d in range(facets[D]):
> +        code += ["""
> +        dofmap.tabulate_facet_dofs(dofs, %d);
> +        for (unsigned int i=0; i < n; i++) {
> +          std::cout << dofs[i] << std::endl;
> +        }
> +        """ % d]
> +
> +    return dof_template % ("", "\n".join(code))
>
> === modified file 'test/verify_element_code/generate_reference.py'
> --- test/verify_element_code/generate_reference.py      2010-01-18 13:37:16 +0000
> +++ test/verify_element_code/generate_reference.py      2010-01-20 20:43:48 +0000
> @@ -6,9 +6,8 @@
>  __license__  = "GNU GPL version 3 or any later version"
>
>  import commands
> -from ufl import FiniteElement
> -from ffc.fiatinterface import create_element
> -from cppcode import element_code, dofmap_code
> +from ufl import FiniteElement, MixedElement, VectorElement
> +from cppcode import test_code
>
>  def run_command(command):
>     "Run system command and collect any errors."
> @@ -16,53 +15,97 @@
>     return (status == 0, output)
>
>  def generate_reference(element):
> -
> +    print "Generating tests for %s" % element
>     reference = {}
>
>     # Write ufl code to file:
>     open("element.ufl", "w").write("element = %s" % element)
>
>     # Compile element with ffc
> -    c = "ffc element.ufl"
> +    c = "ffc -r tensor element.ufl"
>     (ok, output) = run_command(c)
>     if not ok:
>         print "FFC compilation failed for %s" % element
>
> -    # Fetch codes for this element
> -    fiat_element = create_element(eval(element))
> -    codes = element_code(fiat_element)
> +    # Fetch test code for this element
> +    try:
> +        from ffc.fiatinterface import create_element
> +    except:
> +        from ffc.createelement import create_element
> +    codes = test_code(create_element(eval(element)))
>
>     for (function, code) in codes.iteritems():
> -
>         # Create test (c++)
> +        print "\tGenerating, compiling and testing %s" % function
>         open("test.cpp", "w").write(code)
>
>         # Compile c++ code
> -        c = "g++ `pkg-config --cflags ufc-1` -Werror -o test test.cpp"
> +        pkg_config = "`pkg-config --cflags ufc-1` `pkg-config --libs dolfin`"
> +        c = "g++ %s -Werror -o test test.cpp" % pkg_config
>         (ok, output) = run_command(c)
>         if not ok:
> -            raise Exception, "GCC compilation failed:\n %s" % str(output)
> +            raise Exception, "GCC compilation for %s failed:\n %s" \
> +                  % (function, str(output))
>
>         # Run test
>         (ok, result) = run_command("./test")
> +        if not ok:
> +            raise Exception, "test for %s did not run!\n" % function
>
> -        # Add result to reference dictionary
> +        # Add test code and result to reference dictionary
>         reference[function] = (code, result)
>
>     return reference
>
> -elements = ["FiniteElement('CG', 'triangle', 1)"]
> -references = {}
> -
> -for element in elements:
> -
> -    # Generate reference
> -    references[element] = generate_reference(element)
> -    print "Reference: ", references[element]
> -
> -
> -filename = "REFERENCE.py"
> -file = open(filename, 'w')
> -file.write("references = %s" % str(references))
> +def look_at_reference():
> +    from REFERENCE import references
> +
> +    for (element, reference) in references.iteritems():
> +        print "-------------------------"
> +        print element
> +        for function in reference:
> +            print "\t %s: %s" % (function, reference[function][1])
> +
> +
> +
> +if __name__ == "__main__":
> +
> +    cells = ['triangle', 'tetrahedron']
> +    families = ["DG", "CG", "RT", "BDM", "N1curl"]
> +    max_degree = 0
> +    elements = []
> +
> +    # List of basic elements to be tested
> +    elements += ["FiniteElement('%s', '%s', %d)" % (family, cell, k)
> +                 for family in families
> +                 for k in range(1, max_degree+1)
> +                 for cell in cells]
> +    elements += ["FiniteElement('DG', '%s', 0)" % cell for cell in cells]
> +    elements += ["FiniteElement('%s', 'interval', %d)" % (family, k)
> +                 for family in ['DG', 'CG'] for k in range(1, max_degree+1)]
> +
> +    # Some mixed elements:
> +    mixed = ["MixedElement([VectorElement('CG', 'triangle', 3), FiniteElement('DG', 'triangle', 2)])",
> +             "MixedElement([FiniteElement('CG', 'tetrahedron', 2), FiniteElement('DG', 'tetrahedron', 1)])",
> +             "MixedElement([FiniteElement('DG', 'tetrahedron', 2), FiniteElement('CG', 'tetrahedron', 1)])",
> +             "MixedElement([FiniteElement('RT', 'tetrahedron', 1), FiniteElement('DG', 'tetrahedron', 0)])",
> +             "MixedElement([FiniteElement('DG', 'tetrahedron', 1), FiniteElement('BDM', 'tetrahedron', 2)])",
> +             "MixedElement([FiniteElement('CG', 'triangle', 2), FiniteElement('N1curl', 'triangle', 2)])",
> +             "MixedElement([FiniteElement('CG', 'tetrahedron', 2), FiniteElement('N1curl', 'tetrahedron', 3)])",
> +             "MixedElement([VectorElement('CG', 'tetrahedron', 1), FiniteElement('N1curl', 'tetrahedron', 2)])"]
> +
> +    last = ["MixedElement([VectorElement('CG', 'tetrahedron', 2), FiniteElement('N1curl', 'tetrahedron', 3), FiniteElement('BDM', 'tetrahedron', 3), FiniteElement('DG', 'tetrahedron', 0)])"]
> +
> +    #elements += mixed
> +
> +    # Generate reference for each element and store in dictionary
> +    references = {}
> +    for element in elements:
> +        references[element] = generate_reference(element)
> +
> +    # Dump reference data:
> +    filename = "REFERENCE.py"
> +    file = open(filename, 'w')
> +    file.write("references = %s" % str(references))
>
>
>
> === modified file 'test/verify_element_code/test.py'
> --- test/verify_element_code/test.py    2010-01-18 13:37:16 +0000
> +++ test/verify_element_code/test.py    2010-01-20 20:43:48 +0000
> @@ -5,7 +5,6 @@
>  __author__ = "Marie E. Rognes ([email-masked])"
>  __license__  = "GNU GPL version 3 or any later version"
>
> -
>  # Last modified: 15.01.2010
>  import commands
>
> @@ -24,10 +23,12 @@
>     open("element.ufl", "w").write("element = %s" % element)
>
>     # Compile element with ffc
> -    c = "ffc element.ufl"
> +    c = "ffc -r tensor element.ufl"
>     (ok, output) = run_command(c)
>     if not ok:
>         print "FFC compilation failed for %s" % element
> +        failures += [(element, "FFC compilation", output, "")]
> +        continue
>
>     reference = references[element]
>
> @@ -36,27 +37,50 @@
>
>         (code, reference_result) = reference[function]
>
> -        # Write test code to file
> +        # TEMPORARY: Changes in generated code...
> +        code = code.replace("element_0", "element", 1)
> +
> +        # Write c++ test code to file
>         open("test.cpp", "w").write(code)
>
> -        # Compile file with ufc
> -        c = "g++ `pkg-config --cflags ufc-1` -Werror -o test test.cpp"
> +        # Compile c++ test
> +        pkg_config = "`pkg-config --cflags ufc-1` `pkg-config --libs dolfin`"
> +        c = "g++ %s -Werror -o test test.cpp" % pkg_config
>         (ok, output) = run_command(c)
>         if not ok:
> -            raise Exception, "GCC compilation failed"
> +            print "GCC compilation for %s failed:\n %s" \
> +                   % (function, str(output))
> +            failures += [(element, "GCC compilation for %s" % function, output, "")]
> +            continue
>
>         # Run code
>         (ok, result) = run_command("./test")
> +        if not ok:
> +            print "Running test for for %s failed:\n %s" \
> +                   % (function, str(result))
> +            failures += [(element, "Running test for %s" % function, output, "")]
> +            continue
> +
> +        # Check result against reference
>         if result != reference_result:
> -            failures += [(element, function, output)]
> +            failures += [(element, function, result, reference_result)]
>             print
>         else:
> -            print "OK: %s: %s" % (function, result)
> +            print "\tResult matches reference for %s." % (function)
> +
>
>  if failures:
> -    print "\n%d failure:" % len(failures)
> -    for (element, function, output) in failures:
> -        print "%s for %s failed." % (function, element)
> +    print "\n%d failure(s) [out of possibly %d]:" % (len(failures), len(references.keys())*len(reference.keys()))
> +    file = open("failures.txt", 'w')
> +    for (element, function, result, reference_result) in failures:
> +        msg = "Result does not match reference for %s for %s" % (function, element)
> +        print "Result does not match reference for %s for %s" % (function, element)
> +        file.write(msg + ".\n")
> +        file.write("refere = %s\n" % str(reference_result))
> +        file.write("result = %s\n\n" % str(result))
> +    print "\n%d failure(s) [out of possibly %d]:" % (len(failures), len(references.keys())*len(reference.keys()))
> +    print "Check failures.txt for more info. Good luck and may the force be with you."
> +    file.close()
>
>  else:
>     print "All ok!"
>
>
>

---

From: Marie Rognes <[email-masked]>
Date: Wed, 20 Jan 2010 22:58:40 +0100

> Did you forget to add the file?

Yes (on purpose ;) ) The testing is not quite done yet.
But, running generate_reference.py should generate the REFERENCE.py.

--
Marie

---

From: Marie Rognes <[email-masked]>
Date: Wed, 20 Jan 2010 23:21:12 +0100

> I just pushed a bug fix that should? fix this.

Still fails for BDM/RT. Note that for the piola mapped elements it might
be easier to use format["J"] and format ["inv(J)"] directly than the
format["transform"]. (The latter constantly confuses me with regard to
indices for dx/dX and dX/dx). Also detJ has been given an extra_
(by someone?)

There is a doc-string in _change_variables in evaluatedof.py
on how the Piolas should work in case reminding is useful.

--
Marie

---

From: Marie Rognes <[email-masked]>
Date: Wed, 20 Jan 2010 23:25:14 +0100

> on how the Piolas should work in case reminding is useful.

Actually, the one in interpolatevertex values is probably the relevant one.

--
Marie

---

From: Anders Logg <[email-masked]>
Date: Wed, 20 Jan 2010 23:43:29 +0100

> indices for dx/dX and dX/dx). Also detJ has been given an extra_

Maybe the extra _ is not necessary. Feel free to remove it.

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 21 Jan 2010 09:05:04 +0100 (CET)

> indices for dx/dX and dX/dx). Also detJ has been given an extra_

I'm still using format_old, will fix quadrature code generation first and then change everything to use the new format (evaluate_basis*, quadraturegenerators) extending the dictionary as I go along.

> on how the Piolas should work in case reminding is useful.

OK, thanks.

Kristian