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.