Fix stress_tensor formulation (#175)

* Fix stress_tensor quadrature

Fixes bug where all of quadrature points were not summed over properly for each basis function during the element-wise integration of the stress tensor. Also adds multiplication by relevant quadrature weights and the determinant of the jacobian.

* Fixing stress_tensor kernel formuation

Updated the stress tensor kernel to allign with Hooke's law for homogenous isotropic materials.

* Add volume correction in stress_tensor.jl

Author: jbecktt

src/Functions/stress_tensor.jl

@@ -62,12 +62,14 @@ function (f::ElementStressTensor)(u::DisplacementResult; element_dofs=false)
     n_basefuncs = getnbasefunctions(st.cellvalues)
     n_quad = getnquadpoints(st.cellvalues)
     dim = TopOptProblems.getdim(st.problem)
-    return sum(
-        map(1:n_basefuncs, 1:n_quad) do a, q_point
+    V = sum(st.cellvalues.detJdV)
+    return sum(map(1:n_quad) do  q_point
+        dΩ = getdetJdV(st.cellvalues, q_point) 
+        sum(map(1:n_basefuncs) do a   
             _u = cellu[dim * (a - 1) .+ (1:dim)]
             return tensor_kernel(f, q_point, a)(DisplacementResult(_u))
-        end,
-    )
+        end) * dΩ
+    end) ./ V
 end
 
 @params struct ElementStressTensorKernel{T} <: AbstractFunction{T}
@@ -82,8 +84,8 @@ function (f::ElementStressTensorKernel)(u::DisplacementResult)
     @unpack E, ν, q_point, a, cellvalues = f
     ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a))
     ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2
-    c1 = E * ν / (1 - ν^2) * sum(diag(ϵ))
-    c2 = E * ν * (1 + ν)
+    c1 = E * ν / ((1 + ν)*(1 - 2*ν)) * sum(diag(ϵ))
+    c2 = E / (1 + ν)
     return c1 * I + c2 * ϵ
 end
 function ChainRulesCore.rrule(f::ElementStressTensorKernel, u::DisplacementResult)