# Von mises stress

**URL:** <https://discourse.julialang.org/t/von-mises-stress/127393>\
**Category:** New to Julia\
**Tags:** gridap\
**Created:** [March 26, 2025, 9:34pm UTC](https://discourse.julialang.org/t/von-mises-stress/127393 "2025-03-26T21:34:55Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mary](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mary/32/53139_2.png) [@mary](https://discourse.julialang.org/u/mary)\
**Post date:** [March 26, 2025, 9:34pm UTC](https://discourse.julialang.org/t/von-mises-stress/127393/1 "2025-03-26T21:34:55Z")

</div>

Hi, I am using gridap in julia to find the vonMises stress. As it is for topology optimization I penalized the E for material distribution. It is working correctly when I have unifor material distribution but it is not showing correct results when material distribution is not uniform

```julia
const E_mat = 1.0
const ν_mat = 0.3 
const ηe = 1e-9
const penal = 3

function ElasFourthOrderConstTensor(E,ν,PlanarState)
    # 1 for Plane Stress and 2 Plane Strain Condition
    if PlanarState == 1
        C1111 =E/(1-ν*ν)
        C1122 = (ν*E)/(1-ν*ν)
        C1112 = 0.0
        C2222 =E/(1-ν*ν)
        C2212 = 0.0
        C1212 =E/(2*(1+ν))
    elseif PlanarState == 2
        C1111 = (E*(1-ν*ν))/((1+ν)*(1-ν-2*ν*ν))
        C1122 = (ν*E)/(1-ν-2*ν*ν)
        C1112 = 0.0
        C2222 = (E*(1-ν))/(1-ν-2*ν*ν)
        C2212 = 0.0
        C1212 =E/(2*(1+ν))
    end
    C_ten = SymFourthOrderTensorValue(C1111 ,C1112 ,C1122 ,C1112 ,C1212 ,C2212 ,C1122 ,C2212 ,C2222)
    return C_ten
end 
const C_mat = ElasFourthOrderConstTensor(E_mat ,ν_mat ,1)
function σfun(ε_in,pth)
        σ = (C_mat⊙ε_in)
    return σ
    
end

function Em(p)
    Em = (ηe + (1-ηe)*(p ^ penal))
    return Em
end 
f = VectorValue(0,-0.125) 
function stepDisp(fem_params,pth)
    A_Disp(u,v,pth) = ((p->Em(p))∘pth) * ε(v) ⊙ (σfun∘(ε(u),pth))
    a_Disp(u,v) = ∫(A_Disp(u,v,pth))fem_params.dΩ
    b_Disp(v) = ∫(v ⋅ f)fem_params.dΓ_Load
    op_Disp = AffineFEOperator(a_Disp ,b_Disp ,U_Disp ,fem_params.V0_Disp)
    uh_out = solve(op_Disp)
    return uh_out
end 
I2 = SymTensorValue{2,Float64}(1.0 ,0.0 ,1.0)
I4 = I2⊗I2
I4_sym = one(SymFourthOrderTensorValue{2,Float64})
P_vol = (1.0/2)*I4
P_dev = I4_sym - P_vol
function σvon(ε_in,pth)
    
    σ = (C_mat⊙ε_in)
    von_dev = P_dev ⊙ σ
    j2 = 0.5*(von_dev ⊙ von_dev)
    von_stress = sqrt(3*j2)
    
    return von_stress
end
p0 = rand(fem_params.np)
pf_vec = pf_p0(p0; r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
uh = stepDisp(fem_params,pth)
von_stress = (σvon∘(ε(uh),pth))

```

this is the result

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/7/07ece89aa05851bbc789a8cec64e40bd08a7eeef.png)  
while there is almost no material in high stress region. graph below shows material distribution  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/5/25b16f67f72c2c1755fea74e8a1217f0e24d381d.png)

I did my best to include most of the code. I appreciate your help.