(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 (meg@xxxxxxxxx)"
__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 (meg@xxxxxxxxx)"
__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!"
