wip (#10)

* cahce stresses at mid surface to determine cell upgrade

* store stresses in qp in own structure

Author: lijas

docs/examples/doubly_curved.jl

@@ -164,7 +164,7 @@ using Test
 
 #For debug plotting
 if false
-    using Plots; pyplot(); PyPlot.pygui(true)
+    using Plots; plotly()
     using JLD2; using FileIO
 
 
@@ -183,11 +183,19 @@ if false
     σxx_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 1, 1)
     σyy_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 2, 2)
     σxy_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 1, 2)
+    σxz_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 1, 3)./ σ₀
+    σyz_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 2, 3)./ σ₀
+    σzz_lumped = getindex.(output.outputdata["Stress at 50%"].data[end][1].stresses, 3, 3)./ σ₀
     zcoords = getindex.(output.outputdata["Stress at 50%"].data[end][1].local_coords, 3)
+    zcoords .-= mean(zcoords)
 
-    plot!(fig[1], σxx_lumped, zcoords, label = "Lumped")
-    plot!(fig[2], σyy_lumped, zcoords, label = "Lumped")
-    plot!(fig[3], σxy_lumped, zcoords, label = "Lumped")
+    plot!(fig[1], σxx_lumped, zcoords, label = "layered")
+    plot!(fig[2], σyy_lumped, zcoords, label = "layered")
+    plot!(fig[3], σxy_lumped, zcoords, label = "layered")
+
+    plot!(fig[6], σxz_lumped, zcoords, label = "layered")
+    plot!(fig[5], σyz_lumped, zcoords, label = "layered")
+    plot!(fig[4], σzz_lumped, zcoords, label = "layered")
 
     #Reference
     #σ_dr_xx = getindex.(data_dr.stresses, 1, 3)./ σ₀
@@ -211,9 +219,9 @@ if false
     #plot!(fig[1], σ_dr_xx./ σ₀, z_dr, label = "Ref")
     #plot!(fig[2], σ_dr_yy, z_dr, label = "Ref")
     #plot!(fig[3], σ_dr_xy, z_dr, label = "Ref")
-    plot!(fig[4], σ_dr_zz, z_dr, label = "Ref")
-    plot!(fig[5], σ_dr_zy, z_dr, label = "Ref")
-    plot!(fig[6], σ_dr_zx, z_dr, label = "Ref")
+    #plot!(fig[4], σ_dr_zz, z_dr, label = "Ref")
+    #plot!(fig[5], σ_dr_zy, z_dr, label = "Ref")
+    #plot!(fig[6], σ_dr_zx, z_dr, label = "Ref")
 
 end
 

src/igashell_integrationdata.jl

@@ -117,7 +117,7 @@ struct IGAShellIntegrationData{dim_p,dim_s,T,ISV<:IGAShellValues,M}
     cell_values_interface::ISV
     cell_values_vertices::ISV
 
-    active_cell_values::Base.RefValue{ISV}
+    cell_value_mid_interfaces::ISV #For getting stress at interface
 
     iqr::QuadratureRule{dim_p,RefCube,T}
     oqr::QuadratureRule{1,RefCube,T}
@@ -128,6 +128,8 @@ struct IGAShellIntegrationData{dim_p,dim_s,T,ISV<:IGAShellValues,M}
     qr_vertices::Array{QuadratureRule{dim_p,RefCube,T}, 1}
 
     cache_values::CachedOOPBasisValues{dim_p,T,M}
+    interfacestresses::Array{SymmetricTensor{2,3,Float64,6},2} #rows: cells,  cols: stress on interface/qp
+    qpstresses::Vector{Vector{SymmetricTensor{2,3,Float64,6}}} #Vector of cells of qps stresses
 
     extraction_operators::Vector{IGA.BezierExtractionOperator{T}}
     oop_order::Int
@@ -260,6 +262,15 @@ function IGAShellIntegrationData(data::IGAShellData{dim_p,dim_s,T}, C::Vector{IG
     cell_values_sr = IGAShellValues(data.thickness, QuadratureRule{dim_p,RefCube}(1), oqr, mid_ip, getnbasefunctions(cache.basis_values_discont), oop_ip)
     set_oop_basefunctions!(cell_values_sr, initial_basis_values)
 
+    #
+    oop_ip = [ip_layered for i in 1:getnbasefunctions(mid_ip)]
+    cell_values_mid_interface = IGAShellValues(data.thickness, QuadratureRule{dim_p,RefCube}(1), oqr_cohesive[1], mid_ip, getnbasefunctions(ip_layered), oop_ip)
+    set_oop_basefunctions!(cell_values_mid_interface, cache.basis_cohesive_layered[1])
+
+    #Store stresses (cauchy) in gauss points and at interfaces
+    interfacestresses = zeros(SymmetricTensor{2,3,T,6}, ninterfaces(data), getncells(data)) 
+    qpstresses        = [zeros(SymmetricTensor{2,3,T,6}, getnquadpoints(cell_values)) for _ in 1:getncells(data)]
+
     ###
     # Cohesive data
     cell_values_cohesive_top = IGAShellValues(data.thickness, iqr_inp_cohesive, oqr_cohesive[2], mid_ip, getnbasefunctions(ip_discont))
@@ -272,10 +283,9 @@ function IGAShellIntegrationData(data::IGAShellData{dim_p,dim_s,T}, C::Vector{IG
         cell_values_sr,
         cell_values_cohesive_top, cell_values_cohesive_bot,
         face_values, side_values, interface_values, vertices_values,
-        #cell_values_face_top, cell_values_face_bot,
-        Ref(cell_values), #Ref(face_values),
+        cell_values_mid_interface,
         iqr, oqr, iqr_inp_cohesive, oqr_face, oqr_cohesive, iqr_sides, iqr_vertices,
-        cache,
+        cache, interfacestresses, qpstresses,
         C, order)
 end
 

src/igashell_main.jl

@@ -401,6 +401,11 @@ function Five.assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDof
     _assemble_stiffnessmatrix_and_forcevector!(dh, igashell, state, IGASHELL_STIFFMAT)
 end
 
+
+function Five.assemble_massmatrix!(dh::Ferrite.AbstractDofHandler, part::IGAShell, state::StateVariables)
+    return nothing
+end
+
 function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHandler, 
                                                      igashell::IGAShell{dim_p,dim_s,T},  
                                                      state::StateVariables, 
@@ -457,14 +462,15 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
             ke = zeros(T, ndofs_layer, ndofs_layer)
 
             states = @view materialstates[:, ilay]
+            stress_state = igashell.integration_data.qpstresses[ic]
 
             resize_cache2!(igashell.cache.cache2, ndofs_layer)
 
             if assemtype == IGASHELL_STIFFMAT
                 @timeit "integrate shell" _get_layer_forcevector_and_stiffnessmatrix!(
                                                     cv, 
                                                     ke, fe, 
-                                                    getmaterial(layerdata(igashell)), states, 
+                                                    getmaterial(layerdata(igashell)), states, stress_state,
                                                     ue_layer, ilay, nlayers(igashell), active_dofs, 
                                                     is_small_deformation_theory(layerdata(igashell)), IGASHELL_STIFFMAT, getwidth(layerdata(igashell)), igashell.cache.cache2)
 
@@ -475,7 +481,7 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
                 @timeit "integrate shell" _get_layer_forcevector_and_stiffnessmatrix!(
                                         cv, 
                                         ke, fe, 
-                                        getmaterial(layerdata(igashell)), ⁿstates, 
+                                        getmaterial(layerdata(igashell)), ⁿstates, stress_state,
                                         ue_layer, ilay, nlayers(igashell), active_dofs, 
                                         is_small_deformation_theory(layerdata(igashell)), IGASHELL_FSTAR, getwidth(layerdata(igashell)), igashell.cache.cache2)
 
@@ -538,7 +544,7 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
 
             states = @view interfacestates[:, iint]
 
-            active_dofs = active_interface_dofs[iint] #1:Ferrite.ndofs_per_cell(dh,ic)#
+            active_dofs = 1:Ferrite.ndofs_per_cell(dh,ic)#active_interface_dofs[iint] #
             ndofs_interface = length(active_dofs)
 
             resize!(ue_interface, ndofs_interface)
@@ -599,10 +605,6 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
 
 end
 
-function Five.assemble_massmatrix!( dh::Ferrite.AbstractDofHandler, igashell::IGAShell{dim_p,dim_s,T}, system_arrays::StateVariables) where {dim_p,dim_s,T}
-
-end
-
 function Five.post_part!(dh, igashell::IGAShell{dim_p,dim_s,T}, states) where {dim_s, dim_p, T}
     #if dim_s == 2
     #    return
@@ -612,50 +614,88 @@ function Five.post_part!(dh, igashell::IGAShell{dim_p,dim_s,T}, states) where {d
 
         cellstate = getcellstate(adapdata(igashell), ic)
 
-        if !is_lumped(cellstate) && !is_layered(cellstate)
+        if is_mixed(cellstate) || is_fully_discontiniuos(cellstate)
             continue
         end
 
         #Get cellvalues for cell
         Ce = get_extraction_operator(intdata(igashell), ic)
-        
-        #Extract stresses from states
-        σ_states = states.partstates[ic].materialstates[:]
-        σ_states = getproperty.(σ_states, :σ)
+
         #Data for cell
         _celldofs = celldofs(dh, cellid)
         ue = states.d[_celldofs]
 
         nnodes = Ferrite.nnodes_per_cell(igashell)
         X = zeros(Vec{dim_s,T}, nnodes)
         Ferrite.cellcoords!(X, dh, cellid)
-        Xᵇ= IGA.compute_bezier_points(Ce, X)
-        celldata = (celldofs = _celldofs, 
-                    Xᵇ=Xᵇ, X=X, ue=ue, 
-                    nlayers=nlayers(igashell), ninterfaces=ninterfaces(igashell), 
-                    cellid=cellid, ic=ic)
+        Xᵇ = IGA.compute_bezier_points(Ce, X)
 
-        #Build basis_values for cell
-        cv = build_cellvalue!(igashell, cellstate)
-        IGA.set_bezier_operator!(cv, Ce)
-        reinit!(cv, Xᵇ)
+        if is_lumped(cellstate)
+            _post_lumped(igashell, Xᵇ, X, ue, Ce, cellstate, ic, cellid)
+        elseif is_layered(cellstate)
+            _post_layered(igashell, Xᵇ, X, ue, Ce, cellstate, ic, cellid)
+        else
+            continue
+        end
 
-        #Build basis_values for stress_recovory
-        cv_sr = intdata(igashell).cell_values_sr
-        oop_values = _build_oop_basisvalue!(igashell, cellstate)
-        set_oop_basefunctions!(cv_sr, oop_values)
-        IGA.set_bezier_operator!(cv_sr, Ce)
-        reinit!(cv_sr, Xᵇ)
+    end
 
-        recover_cell_stresses(srdata(igashell), σ_states, celldata, cv_sr, cv)
+end
+
+function _post_layered(igashell, Xᵇ, X, ue, Ce, cellstate, ic::Int, cellid::Int)
+
+    #Shape values for evaluating stresses at center of cell
+    cv_mid_interface = igashell.integration_data.cell_value_mid_interfaces
+    set_bezier_operator!(cv_mid_interface, Ce)
+    
+    #oop_values = _build_oop_basisvalue!(igashell, cellstate)
+    #set_oop_basefunctions!(cv_mid_interface, oop_values)
+    
+    reinit!(cv_mid_interface, Xᵇ)
+    active_layer_dofs = build_active_layer_dofs(igashell, cellstate)
+
+    iqp = 0
+    for ilay in 1:nlayers(igashell)-1
+        iqp += 1
+        active_dofs = active_layer_dofs[ilay]
+        ue_layer = ue[active_dofs]
+        
+        #Only one quad points per layer 
+        σ, _, _ = _eval_stress_center(cv_mid_interface, igashell.layerdata.layer_materials[ilay], iqp, Xᵇ, ue_layer, active_dofs, is_small_deformation_theory(igashell.layerdata))
+
+        igashell.integration_data.interfacestresses[ilay, ic] = σ
     end
+end
+
+function _post_lumped(igashell, Xᵇ, X, ue, Ce, cellstate, ic::Int, cellid::Int)
 
+    #Extract stresses from states
+    σ_states = igashell.integration_data.qpstresses[ic]
+
+    celldata = (Xᵇ=Xᵇ, X=X, ue=ue, 
+                nlayers=nlayers(igashell), ninterfaces=ninterfaces(igashell), 
+                cellid=cellid, ic=ic)
+
+    #Build basis_values for cell
+    cv = build_cellvalue!(igashell, cellstate)
+    IGA.set_bezier_operator!(cv, Ce)
+    reinit!(cv, Xᵇ)
+
+    #Build basis_values for stress_recovory
+    cv_sr = intdata(igashell).cell_values_sr
+    oop_values = _build_oop_basisvalue!(igashell, cellstate)
+    set_oop_basefunctions!(cv_sr, oop_values)
+    IGA.set_bezier_operator!(cv_sr, Ce)
+    reinit!(cv_sr, Xᵇ)
+
+    recover_cell_stresses(srdata(igashell), σ_states, celldata, cv_sr, cv)
+    
 end
 
 function _get_layer_forcevector_and_stiffnessmatrix!(
                                 cv::IGAShellValues{dim_s,dim_p,T}, 
                                 ke::AbstractMatrix, fe::AbstractVector,
-                                material, materialstate, 
+                                material, materialstate, stress_state, 
                                 ue_layer::AbstractVector{T}, ilay::Int, nlayers::Int, active_dofs::Vector{Int}, 
                                 is_small_deformation_theory::Bool, calculate_what::IGASHELL_ASSEMBLETYPE, width::T, cache::IGAShellCacheSolid{dim_s,T}) where {dim_s,dim_p,T}
 
@@ -704,13 +744,13 @@ function _get_layer_forcevector_and_stiffnessmatrix!(
                     _calculate_linear_forces!(fe, ke, cv, 
                                                 ilay, qpᴸ, qp, width,
                                                 F, R, δF, δɛ,
-                                                material, materialstate, 
+                                                material, materialstate, stress_state,
                                                 ndofs_layer)
                 else
                     _calculate_nonlinear_forces!(fe, ke, cv, 
                                                 ilay, qpᴸ, qp, width,
                                                 F, R, δF, δɛ,
-                                                material, materialstate, 
+                                                material, materialstate, stress_state,
                                                 ndofs_layer)
                 end
             elseif calculate_what === IGASHELL_FSTAR
@@ -727,7 +767,7 @@ function _get_layer_forcevector_and_stiffnessmatrix!(
 
 end
 
-function _calculate_linear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F::Tensor{2,dim_s}, R, δF, δɛ, material, materialstates, ndofs_layer) where {dim_s}
+function _calculate_linear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F::Tensor{2,dim_s}, R, δF, δɛ, material, materialstates, stress_state, ndofs_layer) where {dim_s}
     ɛ = symmetric(F) - one(SymmetricTensor{2,dim_s})
 
     δɛ .= symmetric.(δF)
@@ -740,8 +780,7 @@ function _calculate_linear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F::Ten
     ∂σ∂ɛ = otimesu(R,R) ⊡ ∂̂σ∂ɛ ⊡ otimesu(R',R')
     σ = R⋅_̂σ⋅R'
 
-    #σ, ∂σ∂ɛ, new_matstate = constitutive_driver(material[ilay], ɛ, ⁿmaterialstates[layer_qp])
-    #materialstates[layer_qp] = new_matstate
+    stress_state[qp] = _to3d(_̂σ)
 
     for i in 1:ndofs_layer
 
@@ -762,7 +801,7 @@ function _calculate_linear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F::Ten
 
 end
 
-function _calculate_nonlinear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F, R, δF, δE, material, materialstates, ndofs_layer)
+function _calculate_nonlinear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F, R, δF, δE, material, materialstates, stress_state, ndofs_layer)
     dΩ = getdetJdV(cv,qp)*width
 
     E = symmetric(1/2 * (F' ⋅ F - one(F)))
@@ -774,6 +813,10 @@ function _calculate_nonlinear_forces!(fe, ke, cv, ilay, layer_qp, qp, width, F,
     ∂S∂E = otimesu(R,R) ⊡ _∂S∂E ⊡ otimesu(R',R')
     S = R⋅_S⋅R'
 
+    σ = inv(det(F)) * symmetric(F ⋅ S ⋅ F')
+    _̂σ = symmetric(R'⋅σ⋅R)
+    stress_state[qp] = _to3d(_̂σ)
+
     #σ = inv(det(F)) * F ⋅ S ⋅ F'
 
     # Hoist computations of δE