bad indexing of basisFunctionsGrad

I have tested your code on a single hexahedron of size 1x1x1.

• The C matrix seems correct (in this case).
• The K matrix is wrong. After some tests I found that this is due to the fact that the indexing of your shape function gradient vector (basisFunctionsGrad) is wrong

bFGrad1(k,0) = basisFunctionsGrad[(numNodes*m+l)*3 + k];

=> read again the documentation of gmsh::model::mesh::getBasisFunctions

=> Modify you code so that you are able to test it on a single element.