- changed milestone to 1.6
evaluate_dofs for expressions is sub-optimal for vector, tensor or mixed function spaces
Interpolating an Expression
is currently sub-optimal for two reasons:
-
it iterates over the cells and the over the dofs (by restricting it), but for Lagrange-type basis function it would be more convenient to iterate over the dofs;
-
for vector, tensor, mixed function spaces, it recompute the same function (same cell and coordinates) for each component, when the result is clearly the same.
The first point is DOLFIN-related, but the second is due to FFC.
Here an example:
from fenics import *
code = \
"""\
class MyExpression : public Expression
{
public:
MyExpression() : Expression(2) {};
void eval(dolfin::Array<double>&,
const dolfin::Array<double>&) const
{
std::cout << "call" << std::endl;
}
};
"""
mesh = UnitTriangleMesh()
V = VectorFunctionSpace(mesh, "P", 1)
e = Expression(cppcode = code, element = V.ufl_element())
interpolate(e, V)
From the output you can easily count 6 "call"s (1 cell x 6 dofs). The generated code for the vector space overloads the method ufc::evaluate_dofs
as follows:
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const double* vertex_coordinates,
int cell_orientation,
const ufc::cell& c) const
{
// Declare variables for result of evaluation
double vals[2];
// Declare variable for physical coordinates
double y[2];
y[0] = vertex_coordinates[0];
y[1] = vertex_coordinates[1];
f.evaluate(vals, y, c);
values[0] = vals[0];
y[0] = vertex_coordinates[2];
y[1] = vertex_coordinates[3];
f.evaluate(vals, y, c);
values[1] = vals[0];
y[0] = vertex_coordinates[4];
y[1] = vertex_coordinates[5];
f.evaluate(vals, y, c);
values[2] = vals[0];
y[0] = vertex_coordinates[0];
y[1] = vertex_coordinates[1];
f.evaluate(vals, y, c);
values[3] = vals[1];
y[0] = vertex_coordinates[2];
y[1] = vertex_coordinates[3];
f.evaluate(vals, y, c);
values[4] = vals[1];
y[0] = vertex_coordinates[4];
y[1] = vertex_coordinates[5];
f.evaluate(vals, y, c);
values[5] = vals[1];
}
can be rewritten as
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const double* vertex_coordinates,
int cell_orientation,
const ufc::cell& c) const
{
// Declare variables for result of evaluation
double vals[2];
// Declare variable for physical coordinates
double y[2];
y[0] = vertex_coordinates[0];
y[1] = vertex_coordinates[1];
f.evaluate(vals, y, c);
values[0] = vals[0];
values[3] = vals[1];
y[0] = vertex_coordinates[2];
y[1] = vertex_coordinates[3];
f.evaluate(vals, y, c);
values[1] = vals[0];
values[4] = vals[1];
y[0] = vertex_coordinates[4];
y[1] = vertex_coordinates[5];
f.evaluate(vals, y, c);
values[2] = vals[0];
values[5] = vals[1];
}
Here we save half of the cost of calling eval
, but for tensor-valued functions in 3d the cost would be 1/9th.
Comments (22)
-
-
Consider this in relation to
https://bitbucket.org/fenics-project/dolfin/issue/422/allow-arbitrary-ufl-expressions-passed-to
-
- changed milestone to 1.7
-
If
ufc::function::evaluate
was marked as pure (is there the concept in C++?), compiler might be able to do this optimization. -
- removed milestone
Removing milestone: 1.7 (automated comment)
-
- changed milestone to 2016.2
-
assigned issue to
Will fix this in https://bitbucket.org/fenics-project/ffc/branch/jan/fix-issue-52. The fix will rely on the fact that duplicate eval code is literally same so this will need unit test to check for regression.
-
- changed status to resolved
-
- changed status to open
I believe I've seen few weeks ago when debugging something else, that the factorization of common eval code does not work in non-affine cases (manifolds or higher value ranks). There's some geometry transformation code which differs for components.
-
None of the element generation code works for non-affine meshes. I started work on that last year (...) and want to finish it soon. There are a couple of issues touching the subject, see:
https://bitbucket.org/fenics-project/ffc/issues/87/evaluate_basis-error-for-combination-mixed
https://bitbucket.org/fenics-project/ffc/issues/83/add-a-ufc-coordinate_mapping-class-was-ufc
https://bitbucket.org/fenics-project/ffc/issues/93/replace-code-generation-for-ufc-classes
-
By "non-affine cases" I meant
mapping != affine
, i.e. Piola transformed elements. Exampleelement = VectorElement("Raviart-Thomas", triangle, 1)
emitsy[0] = 0.5*coordinate_dofs[2] + 0.5*coordinate_dofs[4]; y[1] = 0.5*coordinate_dofs[3] + 0.5*coordinate_dofs[5]; f.evaluate(vals, y, c); result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1])); values[0] = result; [...] y[0] = 0.5*coordinate_dofs[2] + 0.5*coordinate_dofs[4]; y[1] = 0.5*coordinate_dofs[3] + 0.5*coordinate_dofs[5]; f.evaluate(vals, y, c); result = (detJ*(K[0]*vals[2] + K[1]*vals[3])) + (detJ*(K[2]*vals[2] + K[3]*vals[3])); values[3] = result;
The first three lines of both blocks should be factored out but it does not work because first four lines are considered.
-
This may still be related:
https://bitbucket.org/fenics-project/ffc/issues/87/evaluate_basis-error-for-combination-mixed
-
Fix for non-affine mapped elements in master by 3d57ab0b63fb50cf83fd0d47a80f467f50138074. This fixes cases like
element = VectorElement("Raviart-Thomas", triangle, 1)
.Remaining is to fix the codepath going through
ffc.evaluatedof._generate_multiple_points_body
. This applies typically to elements with quadrature moments. Reproduce byelement = FiniteElement("Brezzi-Douglas-Marini", triangle, 2)
. This needs some refactoring inevaluatedof.py
. -
- changed milestone to 2017.1
-
- changed component to finite-elements
-
- changed milestone to 2017.2
-
This will hopefully fixed by pull request #65.
-
It's not resolved.
-
Generator for
evaluate_dofs
has been rewritten with new UFLACS generation which again suffers with this problem. See https://bitbucket.org/fenics-project/dolfin/issues/939. -
- changed milestone to 2018.1
-
A fix for the simpler case of single point supported dofs is included in pull request #92 .
-
In the meanwhile, is it possible to work around this bug? I got a problem with a code where interpolating a 3x3 tensor expression takes 3 hours due to this issue, almost half of the total run time. I'd like to avoid resorting to an older version of fenics.
-
- removed responsible
- removed milestone
- Log in to comment