From 2394725a617948857f3dc12dec3a989bb81ff564 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Tue, 2 Oct 2018 15:58:55 -0400 Subject: [PATCH 1/8] Elasticity: add a local buffer struct to preallocate things needed in the assembly loop --- REQUIRE | 1 + src/JuliaFEM.jl | 2 +- src/problems_elasticity.jl | 90 +++++++++++++++++++++++--------------- 3 files changed, 57 insertions(+), 36 deletions(-) diff --git a/REQUIRE b/REQUIRE index 5cbec28..1c5b39d 100644 --- a/REQUIRE +++ b/REQUIRE @@ -14,3 +14,4 @@ HeatTransfer MortarContact2D MortarContact2DAD FEMBeam +Parameters diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index b4a4812..71d5ef5 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -106,7 +106,7 @@ about JuliaFEM, please visit our website at module JuliaFEM using SparseArrays, LinearAlgebra, Statistics -using Reexport, ForwardDiff, LightXML, HDF5 +using Reexport, ForwardDiff, LightXML, HDF5, Parameters @reexport using FEMBase import FEMBase: get_unknown_field_name, get_unknown_field_dimension, diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index faff70b..3717593 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -98,6 +98,55 @@ function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}}) end end +Parameters.@with_kw struct Elasticity3DLocalBuffers{B, T} + ndofs :: Int + dim :: Int + bi :: BasisInfo{B, T} + BL :: Matrix{T} = zeros(6, ndofs) + BNL :: Matrix{T} = zeros(9, ndofs) + Km :: Matrix{T} = zeros(ndofs, ndofs) + Kg :: Matrix{T} = zeros(ndofs, ndofs) + f_int :: Vector{T} = zeros(ndofs) + f_ext :: Vector{T} = zeros(ndofs) + gradu :: Matrix{T} = zeros(dim, dim) + strain :: Matrix{T} = zeros(dim, dim) + strain_vec :: Vector{T} = zeros(6) + stress_vec :: Vector{T} = zeros(6) + F :: Matrix{T} = zeros(dim, dim) + D :: Matrix{T} = zeros(6, 6) + Dtan :: Matrix{T} = zeros(6, 6) + Bt_mul_D :: Matrix{T} = zeros(ndofs, 6) + Bt_mul_D_mul_B :: Matrix{T} = zeros(ndofs, ndofs) + Bt_mul_S :: Vector{T} = zeros(ndofs) +end + +function reset_element!(buf::Elasticity3DLocalBuffers) + fill!(buf.Km, 0.0) + fill!(buf.Kg, 0.0) + fill!(buf.f_int, 0.0) + fill!(buf.f_ext, 0.0) + return +end + +function reset_integration_point!(buf::Elasticity3DLocalBuffers) + fill!(buf.F, 0.0) + fill!(buf.strain, 0.0) + fill!(buf.D, 0.0) + fill!(buf.BL, 0.0) + fill!(buf.BNL, 0.0) + return +end + +function to_voigt!(strain_vec, strain) + strain_vec[1] = strain[1,1] + strain_vec[2] = strain[2,2] + strain_vec[3] = strain[3,3] + strain_vec[4] = 2.0*strain[1,2] + strain_vec[5] = 2.0*strain[2,3] + strain_vec[6] = 2.0*strain[1,3] + return +end + """ Assemble 3d continuum elements in general solid mechanics problem. """ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, @@ -108,34 +157,17 @@ function assemble!(assembly::Assembly, nnodes = length(El) ndofs = dim*nnodes - BL = zeros(6, ndofs) - BNL = zeros(9, ndofs) - Km = zeros(ndofs, ndofs) - Kg = zeros(ndofs, ndofs) - f_int = zeros(ndofs) - f_ext = zeros(ndofs) - bi = BasisInfo(El) - gradu = zeros(dim, dim) - strain = zeros(dim, dim) - strain_vec = zeros(6) - stress_vec = zeros(6) - F = zeros(dim, dim) - D = zeros(6, 6) - Dtan = zeros(6, 6) - Bt_mul_D = zeros(ndofs, 6) - Bt_mul_D_mul_B = zeros(ndofs, ndofs) - Bt_mul_S = zeros(ndofs) + buffer = Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi = BasisInfo(El)) for element in elements - + Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, gradu, strain, + strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = buffer u = element("displacement", time) - fill!(Km, 0.0) - fill!(Kg, 0.0) - fill!(f_int, 0.0) - fill!(f_ext, 0.0) + reset_element!(buffer) for ip in get_integration_points(element) + reset_integration_point!(buffer) X = element("geometry", time) eval_basis!(bi, X, ip) w = ip.weight*bi.detJ @@ -144,8 +176,6 @@ function assemble!(assembly::Assembly, grad!(bi, gradu, u) # displacement gradient ∇u # calculate strain tensor and deformation gradient - fill!(strain, 0.0) - fill!(F, 0.0) F[:,:] += I if props.finite_strain strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) @@ -154,16 +184,10 @@ function assemble!(assembly::Assembly, strain[:,:] = 1/2 * (gradu + gradu') end - strain_vec[1] = strain[1,1] - strain_vec[2] = strain[2,2] - strain_vec[3] = strain[3,3] - strain_vec[4] = 2.0*strain[1,2] - strain_vec[5] = 2.0*strain[2,3] - strain_vec[6] = 2.0*strain[1,3] + to_voigt!(strain_vec, strain) # material stiffness start - fill!(BL, 0.0) if props.finite_strain for i=1:nnodes BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i] @@ -201,7 +225,6 @@ function assemble!(assembly::Assembly, # calculate stress - fill!(D, 0.0) E = element("youngs modulus", ip, time) nu = element("poissons ratio", ip, time) la = E*nu/((1.0+nu)*(1.0-2.0*nu)) @@ -230,7 +253,6 @@ function assemble!(assembly::Assembly, calculate_stress! = plastic_def["type"] yield_surface_ = plastic_def["yield_surface"] params = plastic_def["params"] - initialize_internal_params!(params, ip, Val{:type_3d}) if time == 0.0 @@ -273,8 +295,6 @@ function assemble!(assembly::Assembly, if props.geometric_stiffness # take geometric stiffness into account - fill!(BNL, 0.0) - for i=1:size(dN, 2) BNL[1, 3*(i-1)+1] = dN[1,i] BNL[2, 3*(i-1)+1] = dN[2,i] From cd22d1a5710e690857d06fa0075891206e326ed4 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Fri, 19 Oct 2018 17:35:30 -0400 Subject: [PATCH 2/8] use buffer and call assemble for one elements, also some tweaks for using sparsity pattern in FEMBase.jl --- src/problems_elasticity.jl | 388 +++++++++--------- src/solvers.jl | 10 +- .../test_elasticity_tet10_stiffness_matrix.jl | 4 +- 3 files changed, 209 insertions(+), 193 deletions(-) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 3717593..4993479 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -70,6 +70,15 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, end end +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, + elements::Vector{<:Element}, time, formulation) + local_buffer = allocate_buffer(problem, elements) + assembler = FEMBase.start_assemble(assembly.K) + for element in elements + assemble_element!(assembly, assembler, problem, element, local_buffer, time, formulation) + end +end + include("problems_elasticity_2d.jl") const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9} @@ -120,6 +129,16 @@ Parameters.@with_kw struct Elasticity3DLocalBuffers{B, T} Bt_mul_S :: Vector{T} = zeros(ndofs) end +function allocate_buffer(problem::Problem{Elasticity}, ::Vector{Element{El}}) where El<:Elasticity3DVolumeElements + dim = get_unknown_field_dimension(problem) + + nnodes = length(El) + ndofs = dim*nnodes + + return Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi = BasisInfo(El)) +end + + function reset_element!(buf::Elasticity3DLocalBuffers) fill!(buf.Km, 0.0) fill!(buf.Kg, 0.0) @@ -148,9 +167,11 @@ function to_voigt!(strain_vec, strain) end """ Assemble 3d continuum elements in general solid mechanics problem. """ -function assemble!(assembly::Assembly, +function assemble_element!(assembly::Assembly, + assembler::FEMBase.AssemblerSparsityPattern, problem::Problem{Elasticity}, - elements::Vector{Element{El}}, + element::Element{El}, + local_buffer::Elasticity3DLocalBuffers, time, ::Type{Val{:continuum}}) where El<:Elasticity3DVolumeElements props = problem.properties dim = get_unknown_field_dimension(problem) @@ -158,206 +179,201 @@ function assemble!(assembly::Assembly, nnodes = length(El) ndofs = dim*nnodes - buffer = Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi = BasisInfo(El)) + Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, gradu, strain, + strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = local_buffer + u = element("displacement", time) + reset_element!(local_buffer) - for element in elements - Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, gradu, strain, - strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = buffer - u = element("displacement", time) - reset_element!(buffer) + for ip in get_integration_points(element) + reset_integration_point!(local_buffer) + X = element("geometry", time) + eval_basis!(bi, X, ip) + w = ip.weight*bi.detJ + N = bi.N + dN = bi.grad # deriatives of basis functions w.r.t. X, i.e. ∂N/∂X + grad!(bi, gradu, u) # displacement gradient ∇u - for ip in get_integration_points(element) - reset_integration_point!(buffer) - X = element("geometry", time) - eval_basis!(bi, X, ip) - w = ip.weight*bi.detJ - N = bi.N - dN = bi.grad # deriatives of basis functions w.r.t. X, i.e. ∂N/∂X - grad!(bi, gradu, u) # displacement gradient ∇u + # calculate strain tensor and deformation gradient + F[:,:] += I + if props.finite_strain + strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) + F[:,:] += gradu + else + strain[:,:] = 1/2 * (gradu + gradu') + end - # calculate strain tensor and deformation gradient - F[:,:] += I - if props.finite_strain - strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) - F[:,:] += gradu - else - strain[:,:] = 1/2 * (gradu + gradu') + to_voigt!(strain_vec, strain) + + # material stiffness start + + if props.finite_strain + for i=1:nnodes + BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i] + BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i] + BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i] + BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i] + BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i] + BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i] + BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i] + BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i] + BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i] + BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i] + BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i] + BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i] + BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i] + BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i] + BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i] + BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i] + BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i] + BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i] + end + else + for i=1:nnodes + BL[1, 3*(i-1)+1] = dN[1,i] + BL[2, 3*(i-1)+2] = dN[2,i] + BL[3, 3*(i-1)+3] = dN[3,i] + BL[4, 3*(i-1)+1] = dN[2,i] + BL[4, 3*(i-1)+2] = dN[1,i] + BL[5, 3*(i-1)+2] = dN[3,i] + BL[5, 3*(i-1)+3] = dN[2,i] + BL[6, 3*(i-1)+1] = dN[3,i] + BL[6, 3*(i-1)+3] = dN[1,i] + end + end + + # calculate stress + + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) + la = E*nu/((1.0+nu)*(1.0-2.0*nu)) + mu = E/(2.0*(1.0+nu)) + D[1,1] = D[2,2] = D[3,3] = 2*mu + la + D[4,4] = D[5,5] = D[6,6] = mu + D[1,2] = D[2,1] = D[2,3] = D[3,2] = D[1,3] = D[3,1] = la + + # determine material model + + material_model = :linear_elasticity + if haskey(element, "plasticity") + material_model = :ideal_plasticity + end + + # calculate stress vector based on material model + + if material_model == :linear_elasticity + Dtan[:,:] = D[:,:] + stress_vec[:] = Dtan * strain_vec + end + + if material_model == :ideal_plasticity + plastic_def = element("plasticity")[ip.id] + + calculate_stress! = plastic_def["type"] + yield_surface_ = plastic_def["yield_surface"] + params = plastic_def["params"] + initialize_internal_params!(params, ip, Val{:type_3d}) + + if time == 0.0 + error("Given step time = $(time). Please select time > 0.0") end - to_voigt!(strain_vec, strain) - - # material stiffness start - - if props.finite_strain - for i=1:nnodes - BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i] - BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i] - BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i] - BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i] - BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i] - BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i] - BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i] - BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i] - BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i] - BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i] - BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i] - BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i] - BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i] - BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i] - BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i] - BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i] - BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i] - BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i] - end - else - for i=1:nnodes - BL[1, 3*(i-1)+1] = dN[1,i] - BL[2, 3*(i-1)+2] = dN[2,i] - BL[3, 3*(i-1)+3] = dN[3,i] - BL[4, 3*(i-1)+1] = dN[2,i] - BL[4, 3*(i-1)+2] = dN[1,i] - BL[5, 3*(i-1)+2] = dN[3,i] - BL[5, 3*(i-1)+3] = dN[2,i] - BL[6, 3*(i-1)+1] = dN[3,i] - BL[6, 3*(i-1)+3] = dN[1,i] - end - end - - # calculate stress - - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - la = E*nu/((1.0+nu)*(1.0-2.0*nu)) - mu = E/(2.0*(1.0+nu)) - D[1,1] = D[2,2] = D[3,3] = 2*mu + la - D[4,4] = D[5,5] = D[6,6] = mu - D[1,2] = D[2,1] = D[2,3] = D[3,2] = D[1,3] = D[3,1] = la - - # determine material model - - material_model = :linear_elasticity - if haskey(element, "plasticity") - material_model = :ideal_plasticity - end - - # calculate stress vector based on material model - - if material_model == :linear_elasticity - Dtan[:,:] = D[:,:] - stress_vec[:] = Dtan * strain_vec - end - - if material_model == :ideal_plasticity - plastic_def = element("plasticity")[ip.id] - - calculate_stress! = plastic_def["type"] - yield_surface_ = plastic_def["yield_surface"] - params = plastic_def["params"] - initialize_internal_params!(params, ip, Val{:type_3d}) - - if time == 0.0 - error("Given step time = $(time). Please select time > 0.0") - end - - t_last = ip("prev_time", time) - update!(ip, "prev_time", time => t_last) - dt = time - t_last - stress_last = ip("stress", t_last) - strain_last = ip("strain", t_last) - dstrain_vec = strain_vec - strain_last - fill!(stress_vec, 0.0) - fill!(Dtan, 0.0) - plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) - - end - - :strain in props.store_fields && update!(ip, "strain", time => strain_vec) - :stress in props.store_fields && update!(ip, "stress", time => stress_vec) - :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) - :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) - :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3]) - :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4]) - :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) - :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) - :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain) - - #Km += w*BL'*Dtan*BL - mul!(Bt_mul_D, transpose(BL), Dtan) - mul!(Bt_mul_D_mul_B, Bt_mul_D, BL) - rmul!(Bt_mul_D_mul_B, w) - for i=1:ndofs^2 - @inbounds Km[i] += Bt_mul_D_mul_B[i] - end - - # material stiffness end - - if props.geometric_stiffness - # take geometric stiffness into account - - for i=1:size(dN, 2) - BNL[1, 3*(i-1)+1] = dN[1,i] - BNL[2, 3*(i-1)+1] = dN[2,i] - BNL[3, 3*(i-1)+1] = dN[3,i] - BNL[4, 3*(i-1)+2] = dN[1,i] - BNL[5, 3*(i-1)+2] = dN[2,i] - BNL[6, 3*(i-1)+2] = dN[3,i] - BNL[7, 3*(i-1)+3] = dN[1,i] - BNL[8, 3*(i-1)+3] = dN[2,i] - BNL[9, 3*(i-1)+3] = dN[3,i] - end - - S3 = zeros(3*dim, 3*dim) - S3[1,1] = stress_vec[1] - S3[2,2] = stress_vec[2] - S3[3,3] = stress_vec[3] - S3[1,2] = S3[2,1] = stress_vec[4] - S3[2,3] = S3[3,2] = stress_vec[5] - S3[1,3] = S3[3,1] = stress_vec[6] - S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] - - Kg += w*BNL'*S3*BNL - - end - - # internal load - mul!(Bt_mul_S, transpose(BL), stress_vec) - rmul!(Bt_mul_S, w) - for i=1:ndofs - @inbounds f_int[i] += Bt_mul_S[i] - end - - # external load start - - if haskey(element, "displacement load") - T = element("displacement load", ip, time) - f_ext += w*vec(T*N) - end - - for i=1:dim - if haskey(element, "displacement load $i") - b = element("displacement load $i", ip, time) - f_ext[i:dim:end] += w*vec(b*N) - end - end - - # external load end + t_last = ip("prev_time", time) + update!(ip, "prev_time", time => t_last) + dt = time - t_last + stress_last = ip("stress", t_last) + strain_last = ip("strain", t_last) + dstrain_vec = strain_vec - strain_last + fill!(stress_vec, 0.0) + fill!(Dtan, 0.0) + plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) end - gdofs = get_gdofs(problem, element) + :strain in props.store_fields && update!(ip, "strain", time => strain_vec) + :stress in props.store_fields && update!(ip, "stress", time => stress_vec) + :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) + :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) + :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3]) + :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4]) + :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) + :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) + :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain) - # add contributions to K, Kg, f - add!(assembly.K, gdofs, gdofs, Km) + #Km += w*BL'*Dtan*BL + mul!(Bt_mul_D, transpose(BL), Dtan) + mul!(Bt_mul_D_mul_B, Bt_mul_D, BL) + rmul!(Bt_mul_D_mul_B, w) + for i=1:ndofs^2 + @inbounds Km[i] += Bt_mul_D_mul_B[i] + end + + # material stiffness end if props.geometric_stiffness - add!(assembly.Kg, gdofs, gdofs, Kg) + # take geometric stiffness into account + + for i=1:size(dN, 2) + BNL[1, 3*(i-1)+1] = dN[1,i] + BNL[2, 3*(i-1)+1] = dN[2,i] + BNL[3, 3*(i-1)+1] = dN[3,i] + BNL[4, 3*(i-1)+2] = dN[1,i] + BNL[5, 3*(i-1)+2] = dN[2,i] + BNL[6, 3*(i-1)+2] = dN[3,i] + BNL[7, 3*(i-1)+3] = dN[1,i] + BNL[8, 3*(i-1)+3] = dN[2,i] + BNL[9, 3*(i-1)+3] = dN[3,i] + end + + S3 = zeros(3*dim, 3*dim) + S3[1,1] = stress_vec[1] + S3[2,2] = stress_vec[2] + S3[3,3] = stress_vec[3] + S3[1,2] = S3[2,1] = stress_vec[4] + S3[2,3] = S3[3,2] = stress_vec[5] + S3[1,3] = S3[3,1] = stress_vec[6] + S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] + + Kg += w*BNL'*S3*BNL + end - add!(assembly.f, gdofs, f_ext - f_int) + # internal load + mul!(Bt_mul_S, transpose(BL), stress_vec) + rmul!(Bt_mul_S, w) + for i=1:ndofs + @inbounds f_int[i] += Bt_mul_S[i] + end + + # external load start + + if haskey(element, "displacement load") + T = element("displacement load", ip, time) + f_ext += w*vec(T*N) + end + + for i=1:dim + if haskey(element, "displacement load $i") + b = element("displacement load $i", ip, time) + f_ext[i:dim:end] += w*vec(b*N) + end + end + + # external load end end + gdofs = get_gdofs(problem, element) + + # add contributions to K, Kg, f + FEMBase.assemble_local_matrix!(assembler, gdofs, Km) + + if props.geometric_stiffness + add!(assembly.Kg, gdofs, gdofs, Kg) + end + + add!(assembly.f, gdofs, f_ext - f_int) + return nothing end diff --git a/src/solvers.jl b/src/solvers.jl index 2b58d37..98399ba 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -57,15 +57,16 @@ problems must have unique node ids. function get_field_assembly(solver::Solver) problems = get_field_problems(solver) + problem = problems[1] M = SparseMatrixCOO() - K = SparseMatrixCOO() + K = spzeros(size(problems[1].assembly.K)...) Kg = SparseMatrixCOO() f = SparseMatrixCOO() fg = SparseMatrixCOO() for problem in problems append!(M, problem.assembly.M) - append!(K, problem.assembly.K) + K += problem.assembly.K append!(Kg, problem.assembly.Kg) append!(f, problem.assembly.f) append!(fg, problem.assembly.fg) @@ -74,7 +75,6 @@ function get_field_assembly(solver::Solver) N = size(K, 1) M = sparse(M, N, N) - K = sparse(K, N, N) if nnz(K) == 0 @warn("Field assembly seems to be empty. Check that elements are ", "pushed to problem and formulation is correct.") @@ -142,7 +142,7 @@ function get_boundary_assembly(solver::Solver, N) g = spzeros(N, 1) for problem in get_boundary_problems(solver) assembly = problem.assembly - K_ = sparse(assembly.K, N, N) + # K_ = assembly.K C1_ = sparse(assembly.C1, N, N) C2_ = sparse(assembly.C2, N, N) D_ = sparse(assembly.D, N, N) @@ -167,7 +167,7 @@ function get_boundary_assembly(solver::Solver, N) error("overconstrained dofs, not solving problem.") end - K += K_ + #K += K_ C1 += C1_ C2 += C2_ D += D_ diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl index d85cf3b..a60575f 100644 --- a/test/test_elasticity_tet10_stiffness_matrix.jl +++ b/test/test_elasticity_tet10_stiffness_matrix.jl @@ -25,8 +25,8 @@ update!(element, "geometry", X) update!(element, "displacement", u) problem = Problem(Elasticity, "tet10", 3) add_element!(problem, element) -time = 0.0 -assemble!(problem, time) +ttime = 0.0 +assemble!(problem, ttime) eigs = real(eigvals(Matrix(problem.assembly.K))) eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85, 1632.49, 1264.32, 1212.42, 817.905, From 32ac97f4b1378efc2b56889b0e72b71b041be7c5 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Sun, 28 Oct 2018 20:40:55 -0400 Subject: [PATCH 3/8] keep working on csc assembly --- src/JuliaFEM.jl | 2 ++ src/problems_dirichlet.jl | 6 +++--- src/problems_elasticity.jl | 23 +++++++++-------------- src/solvers.jl | 19 +++++++++---------- 4 files changed, 23 insertions(+), 27 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 71d5ef5..0254349 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -108,6 +108,8 @@ module JuliaFEM using SparseArrays, LinearAlgebra, Statistics using Reexport, ForwardDiff, LightXML, HDF5, Parameters +import FEMSparse + @reexport using FEMBase import FEMBase: get_unknown_field_name, get_unknown_field_dimension, assemble!, update!, initialize! diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl index 441e7d9..183a924 100644 --- a/src/problems_dirichlet.jl +++ b/src/problems_dirichlet.jl @@ -88,9 +88,9 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0; end end for (k, v) in field_vals - add!(problem.assembly.C1, k, k, 1.0) - add!(problem.assembly.C2, k, k, 1.0) - add!(problem.assembly.g, k, 1, v) + FEMBase.add!(problem.assembly.C1, k, k, 1.0) + FEMBase.add!(problem.assembly.C2, k, k, 1.0) + FEMBase.add!(problem.assembly.g, k, 1, v) end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 4993479..9730fc8 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -73,7 +73,7 @@ end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{<:Element}, time, formulation) local_buffer = allocate_buffer(problem, elements) - assembler = FEMBase.start_assemble(assembly.K) + assembler = FEMSparse.start_assemble(assembly.K, assembly.f) for element in elements assemble_element!(assembly, assembler, problem, element, local_buffer, time, formulation) end @@ -168,7 +168,7 @@ end """ Assemble 3d continuum elements in general solid mechanics problem. """ function assemble_element!(assembly::Assembly, - assembler::FEMBase.AssemblerSparsityPattern, + assembler::FEMSparse.AssemblerSparsityPattern, problem::Problem{Elasticity}, element::Element{El}, local_buffer::Elasticity3DLocalBuffers, @@ -341,12 +341,9 @@ function assemble_element!(assembly::Assembly, # internal load mul!(Bt_mul_S, transpose(BL), stress_vec) rmul!(Bt_mul_S, w) - for i=1:ndofs - @inbounds f_int[i] += Bt_mul_S[i] - end + f_int .+= Bt_mul_S # external load start - if haskey(element, "displacement load") T = element("displacement load", ip, time) f_ext += w*vec(T*N) @@ -358,22 +355,21 @@ function assemble_element!(assembly::Assembly, f_ext[i:dim:end] += w*vec(b*N) end end - # external load end - end gdofs = get_gdofs(problem, element) + # Update f_ext in place to be f_ext - f_int + f_ext .-= f_int + # add contributions to K, Kg, f - FEMBase.assemble_local_matrix!(assembler, gdofs, Km) + FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) if props.geometric_stiffness - add!(assembly.Kg, gdofs, gdofs, Kg) + FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) end - add!(assembly.f, gdofs, f_ext - f_int) - return nothing end @@ -424,8 +420,7 @@ function assemble!(assembly::Assembly, end gdofs = get_gdofs(problem, element) - add!(assembly.f, gdofs, f) - + FEMSparse.assemble_local_vector!(assembly.f, gdofs, f) end end diff --git a/src/solvers.jl b/src/solvers.jl index 98399ba..552feb7 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -58,17 +58,17 @@ function get_field_assembly(solver::Solver) problems = get_field_problems(solver) problem = problems[1] - M = SparseMatrixCOO() - K = spzeros(size(problems[1].assembly.K)...) - Kg = SparseMatrixCOO() - f = SparseMatrixCOO() - fg = SparseMatrixCOO() + M = problem.assembly.M + K = problem.assembly.K + f = problem.assembly.f + Kg = problem.assembly.Kg + fg = problem.assembly.fg - for problem in problems + for problem in problems[2:end] append!(M, problem.assembly.M) K += problem.assembly.K append!(Kg, problem.assembly.Kg) - append!(f, problem.assembly.f) + f += problem.assembly.f append!(fg, problem.assembly.fg) end @@ -80,7 +80,6 @@ function get_field_assembly(solver::Solver) "pushed to problem and formulation is correct.") end Kg = sparse(Kg, N, N) - f = sparse(f, N, 1) fg = sparse(fg, N, 1) return M, K, Kg, f, fg @@ -146,7 +145,7 @@ function get_boundary_assembly(solver::Solver, N) C1_ = sparse(assembly.C1, N, N) C2_ = sparse(assembly.C2, N, N) D_ = sparse(assembly.D, N, N) - f_ = sparse(assembly.f, N, 1) + # f_ = assembly.f g_ = sparse(assembly.g, N, 1) for dof in assembly.removed_dofs @info("$(problem.name): removing dof $dof from assembly") @@ -171,7 +170,7 @@ function get_boundary_assembly(solver::Solver, N) C1 += C1_ C2 += C2_ D += D_ - f += f_ + #f += f_ g += g_ end return K, C1, C2, D, f, g From 25866860ab65d25c4a5663071c8c795806ee0c02 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Thu, 8 Nov 2018 14:27:44 -0500 Subject: [PATCH 4/8] wip --- src/problems_elasticity.jl | 21 ++++++++++++++------- test/test_heat_4.jl | 2 +- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 9730fc8..1b32356 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -53,6 +53,8 @@ function get_formulation_type(problem::Problem{Elasticity}) return :incremental end + +using InteractiveUtils """ assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time) @@ -73,9 +75,10 @@ end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{<:Element}, time, formulation) local_buffer = allocate_buffer(problem, elements) - assembler = FEMSparse.start_assemble(assembly.K, assembly.f) - for element in elements - assemble_element!(assembly, assembler, problem, element, local_buffer, time, formulation) + assemblers = [FEMSparse.start_assemble(assembly.K, assembly.f) for i in 1:Threads.nthreads()] + @time @Threads.threads for i in 1:length(elements) + tid = Threads.threadid() + assemble_element!(assembly, assemblers[tid], problem, elements[i], local_buffer, time, formulation) end end @@ -182,11 +185,11 @@ function assemble_element!(assembly::Assembly, Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, gradu, strain, strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = local_buffer u = element("displacement", time) + X = element("geometry", time) reset_element!(local_buffer) for ip in get_integration_points(element) reset_integration_point!(local_buffer) - X = element("geometry", time) eval_basis!(bi, X, ip) w = ip.weight*bi.detJ N = bi.N @@ -243,8 +246,10 @@ function assemble_element!(assembly::Assembly, # calculate stress - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) + #E = element("youngs modulus", ip, time)::Float64 + #nu = element("poissons ratio", ip, time)::Float64 + E = 200e3 + nu = 0.3 la = E*nu/((1.0+nu)*(1.0-2.0*nu)) mu = E/(2.0*(1.0+nu)) D[1,1] = D[2,2] = D[3,3] = 2*mu + la @@ -265,6 +270,7 @@ function assemble_element!(assembly::Assembly, stress_vec[:] = Dtan * strain_vec end + #= if material_model == :ideal_plasticity plastic_def = element("plasticity")[ip.id] @@ -289,6 +295,7 @@ function assemble_element!(assembly::Assembly, calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) end + =# :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -345,7 +352,7 @@ function assemble_element!(assembly::Assembly, # external load start if haskey(element, "displacement load") - T = element("displacement load", ip, time) + T = element("displacement load", ip, time)::Vector{Float64} f_ext += w*vec(T*N) end diff --git a/test/test_heat_4.jl b/test/test_heat_4.jl index 5dd86db..e613541 100644 --- a/test/test_heat_4.jl +++ b/test/test_heat_4.jl @@ -39,7 +39,7 @@ run!(analysis) # two increments, nonlinear solver -delete!(element.fields, "temperature") +delete!(element, "temperature") analysis = Analysis(Nonlinear) add_problems!(analysis, problem, bc) From 3c67c9c2c9689a13d273d5242558efacfac1cbc1 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Thu, 15 Nov 2018 12:26:17 -0500 Subject: [PATCH 5/8] single threaded --- src/problems_elasticity.jl | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 1b32356..65cdc5c 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -75,10 +75,9 @@ end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{<:Element}, time, formulation) local_buffer = allocate_buffer(problem, elements) - assemblers = [FEMSparse.start_assemble(assembly.K, assembly.f) for i in 1:Threads.nthreads()] - @time @Threads.threads for i in 1:length(elements) - tid = Threads.threadid() - assemble_element!(assembly, assemblers[tid], problem, elements[i], local_buffer, time, formulation) + assembler = FEMSparse.start_assemble(assembly.K, assembly.f) + for i in 1:length(elements) + assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation) end end @@ -246,10 +245,8 @@ function assemble_element!(assembly::Assembly, # calculate stress - #E = element("youngs modulus", ip, time)::Float64 - #nu = element("poissons ratio", ip, time)::Float64 - E = 200e3 - nu = 0.3 + E = element("youngs modulus", ip, time)::Float64 + nu = element("poissons ratio", ip, time)::Float64 la = E*nu/((1.0+nu)*(1.0-2.0*nu)) mu = E/(2.0*(1.0+nu)) D[1,1] = D[2,2] = D[3,3] = 2*mu + la From ee7532a749899aaf3f08793bff6eddbc4a7ffaa5 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Mon, 19 Nov 2018 07:44:56 -0500 Subject: [PATCH 6/8] updates --- src/preprocess.jl | 19 ++++++++----------- src/problems_elasticity.jl | 31 ++++++++++++++++++++++++++----- test/test_mesh_coloring.jl | 4 ++-- 3 files changed, 36 insertions(+), 18 deletions(-) diff --git a/src/preprocess.jl b/src/preprocess.jl index 9d04d36..57e1bea 100644 --- a/src/preprocess.jl +++ b/src/preprocess.jl @@ -21,11 +21,10 @@ mutable struct Mesh element_sets :: Dict{Symbol, Set{Int}} surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}} surface_types :: Dict{Symbol, Symbol} - coloring::Union{Nothing, Vector{Vector{Int}}} # Each vector contains a list of elements that do not share nodes end function Mesh() - return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), nothing) + return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict()) end """ @@ -322,15 +321,14 @@ function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_n end """ - create_coloring!(mesh::Mesh) + create_coloring!(mesh::Mesh) -> Dict{Int, Int} Greedy algorithm for coloring a grid such that no two cells with the same node have the same color. -This function sets the `coloring` field in `mesh` to a `Vector{Vector{Int}}` where -each vector contains vectors of elements that do not share any nodes. -It is therefore safe to assemble in parallel each element vector by vector. +The returned value is a mapping between an element id and its color. +It is safe to assemble elements with the same color in parallel """ -function create_coloring!(mesh::Mesh) +function create_coloring(mesh::Mesh) # Contains the elements that each node contain cell_containing_node = Dict{Int, Set{Int}}() for (cellid, nodes) in mesh.elements @@ -359,7 +357,7 @@ function create_coloring!(mesh::Mesh) # cell -> color of cell cell_colors = Dict{Int, Int}() # color -> list of cells - final_colors = Vector{Int}[] + final_colors = Set{Int}[] occupied_colors = Set{Int}() # Zero represents no color set yet for (cellid, _) in mesh.elements @@ -389,13 +387,12 @@ function create_coloring!(mesh::Mesh) if free_color == 0 # no free color found, need to bump max colors total_colors += 1 free_color = total_colors - push!(final_colors, Int[]) + push!(final_colors, Set{Int}()) end cell_colors[cellid] = free_color push!(final_colors[free_color], cellid) end - mesh.coloring = final_colors - return mesh + return cell_colors end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 65cdc5c..e2b773d 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -73,11 +73,30 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - elements::Vector{<:Element}, time, formulation) - local_buffer = allocate_buffer(problem, elements) - assembler = FEMSparse.start_assemble(assembly.K, assembly.f) - for i in 1:length(elements) - assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation) + elements::Vector{T}, time, formulation) where {T <: Element} + + if problem.assemble_parallel + # Threaded assembly + + assemblers = [FEMSparse.start_assemble(assembly.K, assembly.f) for i in 1:Threads.nthreads()] + local_buffers = [allocate_buffer(problem, elements) for i in 1:Threads.nthreads()] + #TODO: We have to be a bit careful here, the index of the element is no longer + # + # should only loop over elements that exist in `elements` here + for (color, elements) in FEMBase.get_color_ranges(elements) + Threads.@threads for i in 1:length(elements) + element = elements[i] + tid = Threads.threadid() + assemble_element!(assembly, assemblers[tid], problem, element, local_buffers[tid], time, formulation) + end + end + else + # Normal assembly + local_buffer = allocate_buffer(problem, elements) + assembler = FEMSparse.start_assemble(assembly.K, assembly.f) + for i in 1:length(elements) + assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation) + end end end @@ -247,6 +266,8 @@ function assemble_element!(assembly::Assembly, E = element("youngs modulus", ip, time)::Float64 nu = element("poissons ratio", ip, time)::Float64 + #E = 200e3 + #nu = 0.3 la = E*nu/((1.0+nu)*(1.0-2.0*nu)) mu = E/(2.0*(1.0+nu)) D[1,1] = D[2,2] = D[3,3] = 2*mu + la diff --git a/test/test_mesh_coloring.jl b/test/test_mesh_coloring.jl index b1be331..5e43966 100644 --- a/test/test_mesh_coloring.jl +++ b/test/test_mesh_coloring.jl @@ -6,8 +6,8 @@ datadir = first(splitext(basename(@__FILE__))) fn = joinpath(datadir, "cube_tet4.inp") mesh = JuliaFEM.Mesh(open(parse_abaqus, fn)) - JuliaFEM.create_coloring!(mesh) - for colors in mesh.coloring + coloring = JuliaFEM.create_coloring(mesh) + for colors in coloring for ele_i in colors for ele_j in colors if ele_i == ele_j From 09c355b9161986ac043f93ff9976a756fe23ed5e Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Thu, 29 Nov 2018 17:40:55 -0500 Subject: [PATCH 7/8] make backwards compatible --- Manifest.toml | 340 +++++++++++++++++++++++++++++++++++++ Project.toml | 33 ++++ src/problems_elasticity.jl | 40 ++--- src/solvers.jl | 29 ++-- test/runtests.jl | 2 +- test/test_heat_4.jl | 2 +- 6 files changed, 414 insertions(+), 32 deletions(-) create mode 100644 Manifest.toml create mode 100644 Project.toml diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 0000000..4286329 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,340 @@ +[[AbaqusReader]] +deps = ["Nullables", "Test"] +git-tree-sha1 = "0c6f5373eeb3f10f8f8405038c071b18d23907f0" +uuid = "bc6b9049-e460-56d6-94b4-a597b2c0390d" +version = "0.2.2" + +[[Arpack]] +deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Random", "SparseArrays", "Test"] +git-tree-sha1 = "1ce1ce9984683f0b6a587d5bdbc688ecb480096f" +uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +version = "0.3.0" + +[[AsterReader]] +deps = ["FEMBase", "HDF5", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "4c43c84e501df1886aa10495e032c5c7da398d0b" +uuid = "cb1a753d-6b7b-52e3-9bfe-57d99d445c39" +version = "0.2.0" + +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[BinDeps]] +deps = ["Compat", "Libdl", "SHA", "URIParser"] +git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9" +uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee" +version = "0.8.10" + +[[BinaryProvider]] +deps = ["Libdl", "Pkg", "SHA", "Test"] +git-tree-sha1 = "055eb2690182ebc31087859c3dd8598371d3ef9e" +uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" +version = "0.5.3" + +[[Blosc]] +deps = ["BinaryProvider", "CMakeWrapper", "Compat", "Libdl"] +git-tree-sha1 = "71fb23581e1f0b0ae7be8ccf0ebfb3600e23ca41" +uuid = "a74b3585-a348-5f62-a45c-50e91977d574" +version = "0.5.1" + +[[BufferedStreams]] +deps = ["Compat", "Test"] +git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f" +uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d" +version = "1.0.0" + +[[CMake]] +deps = ["BinDeps", "Libdl", "Test"] +git-tree-sha1 = "74853a75c26a4a73ac391ee26ee29ebeb5583d9f" +uuid = "631607c0-34d2-5d66-819e-eb0f9aa2061a" +version = "1.1.0" + +[[CMakeWrapper]] +deps = ["BinDeps", "CMake", "Libdl", "Parameters", "Test"] +git-tree-sha1 = "2b43d451639984e3571951cc687b8509b0a86c6d" +uuid = "d5fb7624-851a-54ee-a528-d3f3bac0b4a0" +version = "0.2.2" + +[[Calculus]] +deps = ["Compat"] +git-tree-sha1 = "f60954495a7afcee4136f78d1d60350abd37a409" +uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +version = "0.4.1" + +[[CommonSubexpressions]] +deps = ["Test"] +git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.2.0" + +[[Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "ec61a16eed883ad0cfa002d7489b3ce6d039bb9a" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "1.4.0" + +[[Crayons]] +deps = ["Test"] +git-tree-sha1 = "3017c662a988bcb8a3f43306a793617c6524d476" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "1.0.0" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[DiffResults]] +deps = ["Compat", "StaticArrays"] +git-tree-sha1 = "db8acf46717b13d6c48deb7a12007c7f85a70cf7" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "0.0.3" + +[[DiffRules]] +deps = ["Random", "Test"] +git-tree-sha1 = "c49ec69428ffea0c1d1bbdc63d1a70f5df5860ad" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "0.0.7" + +[[Distributed]] +deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[FEMBase]] +deps = ["FEMBasis", "FEMQuad", "LinearAlgebra", "SparseArrays", "Statistics", "Test", "TimerOutputs"] +path = "../FEMBase" +uuid = "fbcbbc08-f1bf-5204-9233-b69f5d396135" +version = "0.2.1+" + +[[FEMBasis]] +deps = ["Calculus", "LinearAlgebra", "Test"] +git-tree-sha1 = "aa935fc1c3daeca832945d5282c8dbcad8b1c252" +uuid = "353fb843-c566-51e6-ba49-78b3e3d5ebb5" +version = "0.2.0" + +[[FEMBeam]] +deps = ["FEMBase", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "c848335c3b921897b0419dba7d1297eb98df03ec" +uuid = "962f5c4a-ca9e-56d3-a659-14f0d9baaf78" +version = "0.3.0" + +[[FEMQuad]] +deps = ["Test"] +git-tree-sha1 = "7bb06e83f551212fdbd9b1197616159e414ef3bb" +uuid = "be8e8821-3f6f-54c2-987c-d2773c3a52cb" +version = "0.2.0" + +[[FEMSparse]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "76d4c475da886db1d76ef96e368c3b67a37b7910" +repo-rev = "kc/csc2" +repo-url = "https://github.com/JuliaFEM/FEMSparse.jl" +uuid = "55713501-a877-5f50-80b5-148fff7ff4b3" +version = "0.0.0" + +[[ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] +git-tree-sha1 = "d8f3e0f19d0d546aa92eb1cd67cd3e515768d9f7" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.0" + +[[HDF5]] +deps = ["BinDeps", "Blosc", "Distributed", "Homebrew", "Libdl", "LinearAlgebra", "Mmap", "Pkg", "Test", "WinRPM"] +git-tree-sha1 = "8c3bcdb44db436cd20106e2381e1c1ac96aa0ee3" +uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +version = "0.10.2" + +[[HTTPClient]] +deps = ["Compat", "LibCURL"] +git-tree-sha1 = "161d5776ae8e585ac0b8c20fb81f17ab755b3671" +uuid = "0862f596-cf2d-50af-8ef4-f2be67dfa83f" +version = "0.2.1" + +[[HeatTransfer]] +deps = ["FEMBase", "LinearAlgebra", "Pkg", "SparseArrays", "Test"] +git-tree-sha1 = "2a2d51f4eba4ceea7f00111a3cf0f0f91c6b1b78" +uuid = "4030f512-cedb-5907-ac7f-4ab05ad75ee7" +version = "0.3.0" + +[[Homebrew]] +deps = ["BinDeps", "InteractiveUtils", "JSON", "Libdl", "Test", "Unicode"] +git-tree-sha1 = "5582ec74f735cf8d12e562a2e65c47f34063bd51" +uuid = "d9be37ee-ecc9-5288-90f1-b9ca67657a75" +version = "0.7.0" + +[[InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[JSON]] +deps = ["Dates", "Distributed", "Mmap", "Sockets", "Test", "Unicode"] +git-tree-sha1 = "1f7a25b53ec67f5e9422f1f551ee216503f4a0fa" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.20.0" + +[[LibCURL]] +deps = ["BinaryProvider", "Compat", "Libdl", "Printf"] +git-tree-sha1 = "6339c87cb76923a3cf947fcd213cbc364355c9c9" +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.4.1" + +[[LibExpat]] +deps = ["Compat"] +git-tree-sha1 = "fde352ec13479e2f90e57939da2440fb78c5e388" +uuid = "522f3ed2-3f36-55e3-b6df-e94fee9b0c07" +version = "0.5.0" + +[[LibGit2]] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[Libz]] +deps = ["BufferedStreams", "Random", "Test"] +git-tree-sha1 = "d405194ffc0293c3519d4f7251ce51baac9cc871" +uuid = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7" +version = "1.0.0" + +[[LightXML]] +deps = ["BinaryProvider", "Libdl", "Test"] +git-tree-sha1 = "aeec7a341652d47bc773475a42952fa78eccd7cc" +uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179" +version = "0.8.0" + +[[LinearAlgebra]] +deps = ["Libdl"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[MortarContact2D]] +deps = ["FEMBase", "LinearAlgebra", "SparseArrays", "Statistics", "Test"] +git-tree-sha1 = "7ab9145a90cdb71b785a62a60b5b619838580a9f" +uuid = "048d6160-1a0b-53cd-a5b3-316946cc8d80" +version = "0.3.0" + +[[MortarContact2DAD]] +deps = ["FEMBase", "ForwardDiff", "LinearAlgebra", "SparseArrays", "Statistics", "Test"] +git-tree-sha1 = "eb2eabdfb19a3adafe33b66e8d3c2d45b94e9ed0" +uuid = "c1673bdb-6aff-560b-99da-c78ea6da9af3" +version = "0.2.0" + +[[NaNMath]] +deps = ["Compat"] +git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "0.3.2" + +[[Nullables]] +deps = ["Compat"] +git-tree-sha1 = "ae1a63457e14554df2159b0b028f48536125092d" +uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd" +version = "0.0.8" + +[[OrderedCollections]] +deps = ["Random", "Serialization", "Test"] +git-tree-sha1 = "85619a3f3e17bb4761fe1b1fd47f0e979f964d5b" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.0.2" + +[[Parameters]] +deps = ["Markdown", "OrderedCollections", "REPL", "Test"] +git-tree-sha1 = "40f540ec96e50c0b2b9efdb11b5e4d0c63f90923" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.10.1" + +[[Pkg]] +deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + +[[Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[Reexport]] +deps = ["Pkg"] +git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "0.2.0" + +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[SpecialFunctions]] +deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"] +git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "0.7.2" + +[[StaticArrays]] +deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"] +git-tree-sha1 = "97c4bf0f647488dd7ac01ea12be5885f88762938" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "0.10.0" + +[[Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[Test]] +deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[TimerOutputs]] +deps = ["Crayons", "Printf", "Test", "Unicode"] +git-tree-sha1 = "89a9bd610d6bfd62a7c2b85112762b99b979fe5f" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.4.0" + +[[URIParser]] +deps = ["Test", "Unicode"] +git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69" +uuid = "30578b45-9adc-5946-b283-645ec420af67" +version = "0.4.0" + +[[UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[WinRPM]] +deps = ["BinDeps", "Compat", "HTTPClient", "LibExpat", "Libdl", "Libz", "URIParser"] +git-tree-sha1 = "2a889d320f3b77d17c037f295859fe570133cfbf" +uuid = "c17dfb99-b4f7-5aad-8812-456da1ad7187" +version = "0.4.2" diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..d010ff2 --- /dev/null +++ b/Project.toml @@ -0,0 +1,33 @@ +name = "JuliaFEM" +uuid = "f80590ac-b429-510a-8a99-e7c46989f22d" +version = "0.5.0" + +[deps] +AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d" +Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +AsterReader = "cb1a753d-6b7b-52e3-9bfe-57d99d445c39" +FEMBase = "fbcbbc08-f1bf-5204-9233-b69f5d396135" +FEMBasis = "353fb843-c566-51e6-ba49-78b3e3d5ebb5" +FEMBeam = "962f5c4a-ca9e-56d3-a659-14f0d9baaf78" +FEMQuad = "be8e8821-3f6f-54c2-987c-d2773c3a52cb" +FEMSparse = "55713501-a877-5f50-80b5-148fff7ff4b3" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +HeatTransfer = "4030f512-cedb-5907-ac7f-4ab05ad75ee7" +LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MortarContact2D = "048d6160-1a0b-53cd-a5b3-316946cc8d80" +MortarContact2DAD = "c1673bdb-6aff-560b-99da-c78ea6da9af3" +Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" + +[extras] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "Pkg", "Documenter"] diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index e2b773d..b91ea78 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -53,8 +53,6 @@ function get_formulation_type(problem::Problem{Elasticity}) return :incremental end - -using InteractiveUtils """ assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time) @@ -76,26 +74,23 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{T}, time, formulation) where {T <: Element} if problem.assemble_parallel + @assert problem.assemble_csc # Threaded assembly - - assemblers = [FEMSparse.start_assemble(assembly.K, assembly.f) for i in 1:Threads.nthreads()] + assemblers = [FEMSparse.start_assemble(assembly.K_csc, assembly.f_csc) for i in 1:Threads.nthreads()] local_buffers = [allocate_buffer(problem, elements) for i in 1:Threads.nthreads()] - #TODO: We have to be a bit careful here, the index of the element is no longer - # - # should only loop over elements that exist in `elements` here for (color, elements) in FEMBase.get_color_ranges(elements) - Threads.@threads for i in 1:length(elements) + for i in 1:length(elements) element = elements[i] tid = Threads.threadid() - assemble_element!(assembly, assemblers[tid], problem, element, local_buffers[tid], time, formulation) + assemble_element!(assembly, assemblers[tid], problem, element, local_buffers[tid], time, formulation, true) end end else # Normal assembly local_buffer = allocate_buffer(problem, elements) - assembler = FEMSparse.start_assemble(assembly.K, assembly.f) + assembler = FEMSparse.start_assemble(assembly.K_csc, assembly.f_csc) for i in 1:length(elements) - assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation) + assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation, problem.assemble_csc) end end end @@ -193,7 +188,8 @@ function assemble_element!(assembly::Assembly, problem::Problem{Elasticity}, element::Element{El}, local_buffer::Elasticity3DLocalBuffers, - time, ::Type{Val{:continuum}}) where El<:Elasticity3DVolumeElements + time, ::Type{Val{:continuum}}, + use_csc = false) where El<:Elasticity3DVolumeElements props = problem.properties dim = get_unknown_field_dimension(problem) @@ -288,7 +284,6 @@ function assemble_element!(assembly::Assembly, stress_vec[:] = Dtan * strain_vec end - #= if material_model == :ideal_plasticity plastic_def = element("plasticity")[ip.id] @@ -313,7 +308,6 @@ function assemble_element!(assembly::Assembly, calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) end - =# :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -388,11 +382,19 @@ function assemble_element!(assembly::Assembly, # Update f_ext in place to be f_ext - f_int f_ext .-= f_int - # add contributions to K, Kg, f - FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) + if use_csc + # add contributions to K, Kg, f + FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) - if props.geometric_stiffness - FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) + if props.geometric_stiffness + FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) + end + else + add!(assembly.f, gdofs, f_ext) + add!(assembly.K, gdofs, gdofs, Km) + if props.geometric_stiffness + add!(assembly.Kg, gdofs, gdofs, Kg) + end end return nothing @@ -445,7 +447,7 @@ function assemble!(assembly::Assembly, end gdofs = get_gdofs(problem, element) - FEMSparse.assemble_local_vector!(assembly.f, gdofs, f) + add!(assembly.f, gdofs, f) end end diff --git a/src/solvers.jl b/src/solvers.jl index 552feb7..e3992ea 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -61,28 +61,35 @@ function get_field_assembly(solver::Solver) M = problem.assembly.M K = problem.assembly.K f = problem.assembly.f + K_csc = problem.assembly.K_csc + f_csc = problem.assembly.f_csc Kg = problem.assembly.Kg fg = problem.assembly.fg for problem in problems[2:end] append!(M, problem.assembly.M) - K += problem.assembly.K + append!(K, problem.assembly.K) append!(Kg, problem.assembly.Kg) - f += problem.assembly.f + append!(f, problem.assembly.f) + # Use in place addition with .+= ? + K_csc += problem.assembly.K_csc + f_csc += problem.assembly.f_csc append!(fg, problem.assembly.fg) end N = size(K, 1) M = sparse(M, N, N) + K = sparse(K, N, N) if nnz(K) == 0 @warn("Field assembly seems to be empty. Check that elements are ", "pushed to problem and formulation is correct.") end + f = sparse(f, N, 1) Kg = sparse(Kg, N, N) fg = sparse(fg, N, 1) - return M, K, Kg, f, fg + return M, problem.assemble_csc ? K_csc : K, Kg, problem.assemble_csc ? f_csc : f, fg end """ Loop through boundary assemblies and check for possible overconstrain situations. """ @@ -141,11 +148,11 @@ function get_boundary_assembly(solver::Solver, N) g = spzeros(N, 1) for problem in get_boundary_problems(solver) assembly = problem.assembly - # K_ = assembly.K + K_ = sparse(assembly.K, N, N) C1_ = sparse(assembly.C1, N, N) C2_ = sparse(assembly.C2, N, N) D_ = sparse(assembly.D, N, N) - # f_ = assembly.f + f_ = sparse(assembly.f, N, 1) g_ = sparse(assembly.g, N, 1) for dof in assembly.removed_dofs @info("$(problem.name): removing dof $dof from assembly") @@ -166,12 +173,12 @@ function get_boundary_assembly(solver::Solver, N) error("overconstrained dofs, not solving problem.") end - #K += K_ - C1 += C1_ - C2 += C2_ - D += D_ - #f += f_ - g += g_ + K .+= K_ + C1 .+= C1_ + C2 .+= C2_ + D .+= D_ + f .+= f_ + g .+= g_ end return K, C1, C2, D, f, g end diff --git a/test/runtests.jl b/test/runtests.jl index 7d4ff3e..afba444 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using JuliaFEM, Test -include("../docs/make.jl") +# include("../docs/make.jl") @testset "JuliaFEM.jl" begin @testset "test_dirichlet.jl" begin diff --git a/test/test_heat_4.jl b/test/test_heat_4.jl index e613541..5dd86db 100644 --- a/test/test_heat_4.jl +++ b/test/test_heat_4.jl @@ -39,7 +39,7 @@ run!(analysis) # two increments, nonlinear solver -delete!(element, "temperature") +delete!(element.fields, "temperature") analysis = Analysis(Nonlinear) add_problems!(analysis, problem, bc) From 2aca9f0d0c2710cbcb131a57f120ced99884e24c Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Mon, 3 Dec 2018 15:04:16 -0500 Subject: [PATCH 8/8] alloc improvements --- src/problems_elasticity.jl | 59 +++++++++++++++++++++++++------------- 1 file changed, 39 insertions(+), 20 deletions(-) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index b91ea78..d0aab28 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -79,7 +79,7 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, assemblers = [FEMSparse.start_assemble(assembly.K_csc, assembly.f_csc) for i in 1:Threads.nthreads()] local_buffers = [allocate_buffer(problem, elements) for i in 1:Threads.nthreads()] for (color, elements) in FEMBase.get_color_ranges(elements) - for i in 1:length(elements) + Threads.@threads for i in 1:length(elements) element = elements[i] tid = Threads.threadid() assemble_element!(assembly, assemblers[tid], problem, element, local_buffers[tid], time, formulation, true) @@ -133,6 +133,9 @@ Parameters.@with_kw struct Elasticity3DLocalBuffers{B, T} Kg :: Matrix{T} = zeros(ndofs, ndofs) f_int :: Vector{T} = zeros(ndofs) f_ext :: Vector{T} = zeros(ndofs) + f_buffer :: Vector{T} = zeros(ndofs) + f_buffer_dim :: Vector{T} = zeros(div(ndofs, dim)) + gdofs :: Vector{Int} = zeros(Int, ndofs) gradu :: Matrix{T} = zeros(dim, dim) strain :: Matrix{T} = zeros(dim, dim) strain_vec :: Vector{T} = zeros(6) @@ -182,6 +185,9 @@ function to_voigt!(strain_vec, strain) return end +const u = ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]) +const X = ([-93.7197, -93.7197, 150.883], [-91.657, -85.8251, 157.885], [-100.523, -88.8309, 157.883], [-91.6593, -88.8309, 157.883], [-92.6883, -89.7724, 154.384], [-96.0902, -87.328, 157.883], [-97.1216, -91.2753, 154.383], [-92.6895, -91.2753, 154.383], [-91.6581, -87.328, 157.883], [-96.0914, -88.8309, 157.883]) +const displacement_load_string = [string("displacement load ", i) for i in 1:3] """ Assemble 3d continuum elements in general solid mechanics problem. """ function assemble_element!(assembly::Assembly, assembler::FEMSparse.AssemblerSparsityPattern, @@ -190,16 +196,19 @@ function assemble_element!(assembly::Assembly, local_buffer::Elasticity3DLocalBuffers, time, ::Type{Val{:continuum}}, use_csc = false) where El<:Elasticity3DVolumeElements + cheating = false props = problem.properties dim = get_unknown_field_dimension(problem) nnodes = length(El) ndofs = dim*nnodes - Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, gradu, strain, + Parameters.@unpack bi, BL, BNL, Km, Kg, f_int, f_ext, f_buffer, f_buffer_dim, gdofs, gradu, strain, strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = local_buffer - u = element("displacement", time) - X = element("geometry", time) + if !cheating + u = element("displacement", time) + X = element("geometry", time) + end reset_element!(local_buffer) for ip in get_integration_points(element) @@ -211,12 +220,15 @@ function assemble_element!(assembly::Assembly, grad!(bi, gradu, u) # displacement gradient ∇u # calculate strain tensor and deformation gradient - F[:,:] += I + #F[:,:] += I + for i in 1:dim + F[i, i] += 1.0 + end if props.finite_strain strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) F[:,:] += gradu else - strain[:,:] = 1/2 * (gradu + gradu') + strain[:,:] .= 1/2 .* (gradu .+ gradu') end to_voigt!(strain_vec, strain) @@ -260,10 +272,13 @@ function assemble_element!(assembly::Assembly, # calculate stress - E = element("youngs modulus", ip, time)::Float64 - nu = element("poissons ratio", ip, time)::Float64 - #E = 200e3 - #nu = 0.3 + if cheating + E = 200e3 + nu = 0.3 + else + E = element("youngs modulus", ip, time)::Float64 + nu = element("poissons ratio", ip, time)::Float64 + end la = E*nu/((1.0+nu)*(1.0-2.0*nu)) mu = E/(2.0*(1.0+nu)) D[1,1] = D[2,2] = D[3,3] = 2*mu + la @@ -272,6 +287,7 @@ function assemble_element!(assembly::Assembly, # determine material model + material_model = :linear_elasticity if haskey(element, "plasticity") material_model = :ideal_plasticity @@ -280,8 +296,8 @@ function assemble_element!(assembly::Assembly, # calculate stress vector based on material model if material_model == :linear_elasticity - Dtan[:,:] = D[:,:] - stress_vec[:] = Dtan * strain_vec + copyto!(Dtan, D) + mul!(stress_vec, Dtan, strain_vec) end if material_model == :ideal_plasticity @@ -328,7 +344,6 @@ function assemble_element!(assembly::Assembly, end # material stiffness end - if props.geometric_stiffness # take geometric stiffness into account @@ -365,29 +380,33 @@ function assemble_element!(assembly::Assembly, # external load start if haskey(element, "displacement load") T = element("displacement load", ip, time)::Vector{Float64} - f_ext += w*vec(T*N) + mul!(f_buffer, w, vec(T*N)) + f_ext .+= f_buffer end for i=1:dim - if haskey(element, "displacement load $i") - b = element("displacement load $i", ip, time) - f_ext[i:dim:end] += w*vec(b*N) + if haskey(element, displacement_load_string[i]) + b = element(displacement_load_string[i], ip, time)::Float64 + mul!(f_buffer_dim, w, N) + for (i, j) in enumerate(1:dim:length(f_ext)) + f_ext[j] = b * f_buffer_dim[i] + end end end # external load end end - gdofs = get_gdofs(problem, element) + FEMBase.get_gdofs!(gdofs, problem, element) # Update f_ext in place to be f_ext - f_int f_ext .-= f_int if use_csc # add contributions to K, Kg, f - FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) + @inbounds FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) if props.geometric_stiffness - FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) + @inbounds FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg) end else add!(assembly.f, gdofs, f_ext)