Second derivatives of coefficients create invalid C++ code

[Martin Sandve Alnæs]

cell = tetrahedron
V = FiniteElement("CG", cell, 2)
f = Coefficient(V)
M = f.dx(0,0)*dx

[Martin Sandve Alnæs]

The good solution involves improving the value numbering such that d^2f/dxdy and d^2f/dydx gets the same symbols. Ties in with handling element symmetries at the same time, it's basically the same thing.

[Martin Sandve Alnæs]

Symmetries are now handled in value numbering, both for tensor elements and second derivatives and higher.

The next failure is that somewhere in the pipeline second derivatives of piecewise linear coefficients should dissapear but they don't. Thus this code fails because the element tables are not generated, and opportunities for constant propagation of zero second derivative values are lost:

degree = 1 # fails
#degree = 2 # ok

cell = triangle
V = VectorElement("CG", cell, degree)
U = Coefficient(V)
D = Coefficient(V)
M = div(dot(grad(D), U))*dx

[Jan Blechta]

Is this correct fix

diff --git a/uflacs/generation/integralgenerator.py b/uflacs/generation/integralgenerator.py
index 23c206b..cfafcb8 100644
--- a/uflacs/generation/integralgenerator.py
+++ b/uflacs/generation/integralgenerator.py
@@ -215,6 +215,10 @@ class IntegralGenerator(object):

         # Build loops for each dofrange
         for dofrange in dofranges:
+            # Empty loop
+            if dofrange[1] == dofrange[0]:
+                continue
+
             dofblock = outer_dofblock + (dofrange,)
             body = []

?

[Martin Sandve Alnæs]

This patch is correct. Although probably better to fix earlier somehow.

[Jan Blechta]

Actually there is another occurrence of zero basis tables in IntegralGenerator.generate_partition() so this seems to work for the original example of #35

diff --git a/uflacs/backends/ffc/definitions.py b/uflacs/backends/ffc/definitions.py
index a11a362..67170c1 100644
--- a/uflacs/backends/ffc/definitions.py
+++ b/uflacs/backends/ffc/definitions.py
@@ -64,7 +64,13 @@ class FFCDefinitionsBackend(MultiFunction):
             pass
         else:
             # No need to store basis function value in its own variable, just get table value directly
+            code += [VariableDecl("double", access, "0.0")]
             uname, begin, end = tabledata
+
+            # Empty loop needs to be skipped as zero tables may not be generated
+            if begin >= end:
+                return code
+
             entity = format_entity_name(self.ir["entitytype"], mt.restriction)

             iq = names.iq
@@ -79,7 +85,6 @@ class FFCDefinitionsBackend(MultiFunction):
             body = [AssignAdd(access, prod)]

             # Loop to accumulate linear combination of dofs and tables
-            code += [VariableDecl("double", access, "0.0")]
             code += [ForRange(idof, begin, end, body=body)]

         return code
diff --git a/uflacs/generation/integralgenerator.py b/uflacs/generation/integralgenerator.py
index 23c206b..3032bc8 100644
--- a/uflacs/generation/integralgenerator.py
+++ b/uflacs/generation/integralgenerator.py
@@ -215,6 +215,10 @@ class IntegralGenerator(object):

         # Build loops for each dofrange
         for dofrange in dofranges:
+            # Empty loop needs to be skipped as zero tables may not be generated
+            if dofrange[1] == dofrange[0]:
+                continue
+
             dofblock = outer_dofblock + (dofrange,)
             body = []

Although I'm not sure whether there should not be further elimination of zero terms. The patch just removes assignments with undefined derivatives but keeps zero intermediate results which are used further, like this

        double w2_d11 = 0.0;
        double w2_d01 = 0.0;
        double w2_d00 = 0.0;
...
        sv361[203] = w2_d11 * sp361[4];
        sv361[204] = w2_d01 * sp361[6];
        sv361[207] = w2_d01 * sp361[4];
        sv361[208] = w2_d00 * sp361[6];

instead of

        double w2_d11 = 0.0;
        for (int ic = 3; ic < 3; ++ic)
        {
            w2_d11 += w[2][ic] * FE0_C0_D02[0][iq][ic - 3];
        }
        double w2_d01 = 0.0;
        for (int ic = 3; ic < 3; ++ic)
        {
            w2_d01 += w[2][ic] * FE0_C0_D02[0][iq][ic - 3];
        }
        double w2_d00 = 0.0;
        for (int ic = 3; ic < 3; ++ic)
        {
           w2_d00 += w[2][ic] * FE0_C0_D02[0][iq][ic - 3];
        }
...
        sv361[203] = w2_d11 * sp361[4];
        sv361[204] = w2_d01 * sp361[6];
        sv361[207] = w2_d01 * sp361[4];
        sv361[208] = w2_d00 * sp361[6];

Earlier patch seems to be needed in FFC.

[Jan Blechta]

Sorry, this was not with example from #35 but with another complicated form.

[Jan Blechta]

@martinal, do you like this fix or would you like to get it fixed more upstream? In the second case, do you have some hint where to look?

[Martin Sandve Alnæs]

I want to fix it earlier, the empty range here shouldn't happen. I have a hunch but need to look into the code. Let me get back to you on Monday. Busy with proposal deadline.

[Jan Blechta]

Ok, thanks.

[Martin Sandve Alnæs]

Sorry for the delay, I'm working on this now. I'll use your workarounds in the first round and do a proper fix in a later iteration over the algorithms.

[Jan Blechta]

Issue #35 was marked as a duplicate of this issue.

[Martin Sandve Alnæs]

Original issue resolved. A more optimal solution depends on a later review of core algorithms.

[Martin Sandve Alnæs]

Removing milestone: 1.6 (automated comment)