Incorrect functional evaluation involving a QuadratureFunction

[Tianyi Li]

This issue is initally posed at FEniCS Q&A.

Consider a trivial 2-d elastic problem defined on a UnitSquareMesh composed of two triangles where the diagonal side appears as "/". All nodes are fixed, except at the right bottom corner where U=[1.0, 0.0] is prescribed. The yield stress is set to an extremely high value. With E=1 and nu=0 The total elastic energy should be M=0.375. The upper left element is stress free due to BC, and in the lower right one we have s=[1, 0, -0.5].

However using the FEniCS Solid Mechanics interface, I found that M=0. I checked that the displacement field is correct, however the total elastic energy functional is incorrectly evaluated. Through an output of the stress in the ConstitutiveUpdate.cpp

_w_stress[num_ip_per_cell*0  + ip] = trial_stress(0);
_w_stress[num_ip_per_cell*1  + ip] = trial_stress(1);
_w_stress[num_ip_per_cell*2  + ip] = trial_stress(3);
std::cout << "stress(0): " << trial_stress(0) << std::endl;
std::cout << "stress(1): " << trial_stress(1) << std::endl;
std::cout << "stress(3): " << trial_stress(3) << std::endl;

I found

step begin: 0
time: 0
stress(0): 0
stress(1): 0
stress(3): 0
stress(0): 0
stress(1): 0
stress(3): 0
Newton iteration 0: r (abs) = 1.000e+00 (tol = 1.000e-15) r (rel) = 1.000e+00 (tol = 1.000e-06)
stress(0): 1
stress(1): 0
stress(3): -0.5
stress(0): 0
stress(1): 0
stress(3): 0
Newton iteration 1: r (abs) = 0.000e+00 (tol = 1.000e-15) r (rel) = 0.000e+00 (tol = 1.000e-06)
Newton solver finished in 1 iterations and 1 linear solver iterations.
M = 0.000e+00

The stress is well calculated in both elements. However the non-zero stress values are overwritten by the zero ones, producing thus a zero elastic energy.