diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl index 0761fe6..55d3de9 100644 --- a/src/problems_contact_3d.jl +++ b/src/problems_contact_3d.jl @@ -5,14 +5,14 @@ const ContactElements3D = Union{Tri3,Tri6,Quad4,Quad8,Quad9} function create_orthogonal_basis(n) I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] - k = argmax([norm(cross(n,I[:,k])) for k in 1:3]) - t1 = cross(n, I[:,k])/norm(cross(n, I[:,k])) + k = argmax([norm(cross(n, I[:, k])) for k in 1:3]) + t1 = cross(n, I[:, k]) / norm(cross(n, I[:, k])) t2 = cross(n, t1) return t1, t2 end """ Create rotation matrix Q for element nodes rotating quantities to nt coordinaet system. """ -function create_rotation_matrix(element::Element{Tri3}, time::Float64) +function create_rotation_matrix(element::Element{M,Tri3}, time::Float64) where M n = element("normal", time) t11, t21 = create_orthogonal_basis(n[1]) t12, t22 = create_orthogonal_basis(n[2]) @@ -28,7 +28,7 @@ function create_rotation_matrix(element::Element{Tri3}, time::Float64) return Q end -function create_rotation_matrix(element::Element{Quad4}, time::Float64) +function create_rotation_matrix(element::Element{M,Quad4}, time::Float64) where M n = element("normal", time) t11, t21 = create_orthogonal_basis(n[1]) t12, t22 = create_orthogonal_basis(n[2]) @@ -47,7 +47,7 @@ function create_rotation_matrix(element::Element{Quad4}, time::Float64) return Q end -function create_rotation_matrix(element::Element{Tri6}, time::Float64) +function create_rotation_matrix(element::Element{M,Tri6}, time::Float64) where M n = element("normal", time) t11, t21 = create_orthogonal_basis(n[1]) t12, t22 = create_orthogonal_basis(n[2]) @@ -96,7 +96,7 @@ function create_contact_segmentation(slave_element, master_elements, x0, n0, tim length(P) < 3 && continue # no clipping or shared edge (no volume) check_orientation!(P, n0) N_P = length(P) - P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + P_area = sum([norm(1 / 2 * cross(P[i] - P[1], P[mod(i, N_P)+1] - P[1])) for i = 2:N_P]) if isapprox(P_area, 0.0) error("Polygon P has zero area") end @@ -106,7 +106,7 @@ function create_contact_segmentation(slave_element, master_elements, x0, n0, tim return result end -function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time::Float64) +function assemble!(problem::Problem{Contact}, slave_element::Element{FS,Tri3}, time::Float64) where FS props = problem.properties field_dim = get_unknown_field_dimension(problem) @@ -149,19 +149,19 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) detJ = virtual_element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ x_gauss = virtual_element("geometry", ip, time) xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) N1 = slave_element(xi_s, time) - De += w*Matrix(Diagonal(vec(N1))) - Me += w*N1'*N1 + De += w * Matrix(Diagonal(vec(N1))) + Me += w * N1' * N1 end # integration points done end # integration cells done end # master elements done - Ae = De*inv(Me) + Ae = De * inv(Me) end @@ -175,8 +175,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time De = zeros(nsl, nsl) Me = zeros(nsl, nm) - ce = zeros(field_dim*nsl) - ge = zeros(field_dim*nsl) + ce = zeros(field_dim * nsl) + ge = zeros(field_dim * nsl) # loop integration cells for cell in get_cells(P, C0) @@ -191,18 +191,18 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) detJ = virtual_element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ # add contributions N1 = vec(get_basis(slave_element, xi_s, time)) N2 = vec(get_basis(master_element, xi_m, time)) - Phi = Ae*N1 - De += w*Phi*N1' - Me += w*Phi*N2' + Phi = Ae * N1 + De += w * Phi * N1' + Me += w * Phi * N2' - x_s = interpolate(N1, map(+,X1,u1)) - x_m = interpolate(N2, map(+,X2,u2)) - ge += w*vec((x_m-x_s)*Phi') + x_s = interpolate(N1, map(+, X1, u1)) + x_m = interpolate(N2, map(+, X2, u2)) + ge += w * vec((x_m - x_s) * Phi') end # integration points done @@ -215,16 +215,16 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time nmdofs = length(mdofs) D3 = zeros(nsldofs, nsldofs) M3 = zeros(nsldofs, nmdofs) - for i=1:field_dim + for i = 1:field_dim D3[i:field_dim:end, i:field_dim:end] += De M3[i:field_dim:end, i:field_dim:end] += Me end add!(problem.assembly.C1, sdofs, sdofs, D3) add!(problem.assembly.C1, sdofs, mdofs, -M3) - add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3) - add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3) - add!(problem.assembly.g, sdofs, Q3'*ge) + add!(problem.assembly.C2, sdofs, sdofs, Q3' * D3) + add!(problem.assembly.C2, sdofs, mdofs, -Q3' * M3) + add!(problem.assembly.g, sdofs, Q3' * ge) end # master elements done @@ -232,7 +232,7 @@ end """ Assemble quadratic surface element to contact problem. """ -function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time::Float64) +function assemble!(problem::Problem{Contact}, slave_element::Element{FS,Tri6}, time::Float64) where FS props = problem.properties field_dim = get_unknown_field_dimension(problem) @@ -241,15 +241,15 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time if alp != 0.0 T = [ - 1.0 0.0 0.0 0.0 0.0 0.0 - 0.0 1.0 0.0 0.0 0.0 0.0 - 0.0 0.0 1.0 0.0 0.0 0.0 - alp alp 0.0 1.0-2*alp 0.0 0.0 - 0.0 alp alp 0.0 1.0-2*alp 0.0 - alp 0.0 alp 0.0 0.0 1.0-2*alp - ] + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + alp alp 0.0 1.0-2*alp 0.0 0.0 + 0.0 alp alp 0.0 1.0-2*alp 0.0 + alp 0.0 alp 0.0 0.0 1.0-2*alp + ] else - T = Matrix(1.0*I, 6, 6) + T = Matrix(1.0 * I, 6, 6) end nsl = length(slave_element) @@ -258,7 +258,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time Q3 = create_rotation_matrix(slave_element, time) - Ae = Matrix(1.0*I, nsl, nsl) + Ae = Matrix(1.0 * I, nsl, nsl) if problem.properties.dual_basis # construct dual basis @@ -311,7 +311,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time check_orientation!(P, n0) N_P = length(P) - P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + P_area = sum([norm(1 / 2 * cross(P[i] - P[1], P[mod(i, N_P)+1] - P[1])) for i = 2:N_P]) if isapprox(P_area, 0.0) error("Polygon P has zero area") end @@ -324,12 +324,12 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) detJ = virtual_element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ x_gauss = virtual_element("geometry", ip, time) xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time) - N1 = vec(slave_element(xi_s, time)*T) - De += w*Matrix(Diagonal(N1)) - Me += w*N1*N1' + N1 = vec(slave_element(xi_s, time) * T) + De += w * Matrix(Diagonal(N1)) + Me += w * N1 * N1' end # integration points done end # integration cells done @@ -340,7 +340,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time end # sub slave elements done - Ae = De*inv(Me) + Ae = De * inv(Me) end @@ -386,7 +386,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time check_orientation!(P, n0) N_P = length(P) - P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + P_area = sum([norm(1 / 2 * cross(P[i] - P[1], P[mod(i, N_P)+1] - P[1])) for i = 2:N_P]) if isapprox(P_area, 0.0) error("Polygon P has zero area") end @@ -398,7 +398,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time nm = length(master_element) De = zeros(nsl, nsl) Me = zeros(nsl, nm) - ge = zeros(field_dim*nsl) + ge = zeros(field_dim * nsl) # 4. loop integration cells for cell in get_cells(P, C0) @@ -414,21 +414,21 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time) detJ = virtual_element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ # add contributions - N1 = vec(get_basis(slave_element, xi_s, time)*T) + N1 = vec(get_basis(slave_element, xi_s, time) * T) N2 = vec(get_basis(master_element, xi_m, time)) - Phi = Ae*N1 + Phi = Ae * N1 - De += w*Phi*N1' - Me += w*Phi*N2' + De += w * Phi * N1' + Me += w * Phi * N2' us = slave_element("displacement", time) um = master_element("displacement", time) - xs = interpolate(N1, map(+,Xs,us)) - xm = interpolate(N2, map(+,Xs,um)) - ge += w*vec((xm-xs)*Phi') + xs = interpolate(N1, map(+, Xs, us)) + xm = interpolate(N2, map(+, Xs, um)) + ge += w * vec((xm - xs) * Phi') end # integration points done @@ -441,16 +441,16 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time nmdofs = length(mdofs) D3 = zeros(nsldofs, nsldofs) M3 = zeros(nsldofs, nmdofs) - for i=1:field_dim + for i = 1:field_dim D3[i:field_dim:end, i:field_dim:end] += De M3[i:field_dim:end, i:field_dim:end] += Me end add!(problem.assembly.C1, sdofs, sdofs, D3) add!(problem.assembly.C1, sdofs, mdofs, -M3) - add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3) - add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3) - add!(problem.assembly.g, sdofs, Q3'*ge) + add!(problem.assembly.C2, sdofs, sdofs, Q3' * D3) + add!(problem.assembly.C2, sdofs, mdofs, -Q3' * M3) + add!(problem.assembly.g, sdofs, Q3' * ge) end # sub master elements done @@ -480,7 +480,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T # 1. calculate nodal normals and tangents for slave element nodes j ∈ S normals = calculate_normals(slave_elements, time, Val{2}; - rotate_normals=props.rotate_normals) + rotate_normals=props.rotate_normals) update!(slave_elements, "normal", time => normals) # 2. loop all slave elements @@ -489,13 +489,13 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T end # slave elements done, contact virtual work ready S = sort(collect(keys(normals))) # slave element nodes - weighted_gap = Dict{Int64, Vector{Float64}}() - contact_pressure = Dict{Int64, Vector{Float64}}() - complementarity_condition = Dict{Int64, Vector{Float64}}() - is_active = Dict{Int64, Int}() - is_inactive = Dict{Int64, Int}() - is_slip = Dict{Int64, Int}() - is_stick = Dict{Int64, Int}() + weighted_gap = Dict{Int64,Vector{Float64}}() + contact_pressure = Dict{Int64,Vector{Float64}}() + complementarity_condition = Dict{Int64,Vector{Float64}}() + is_active = Dict{Int64,Int}() + is_inactive = Dict{Int64,Int}() + is_slip = Dict{Int64,Int}() + is_stick = Dict{Int64,Int}() la = problem.assembly.la @@ -519,13 +519,13 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T if problem.properties.alpha != 0.0 alp = problem.properties.alpha Te = [ - 1.0 0.0 0.0 0.0 0.0 0.0 - 0.0 1.0 0.0 0.0 0.0 0.0 - 0.0 0.0 1.0 0.0 0.0 0.0 - alp alp 0.0 1.0-2*alp 0.0 0.0 - 0.0 alp alp 0.0 1.0-2*alp 0.0 - alp 0.0 alp 0.0 0.0 1.0-2*alp - ] + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + alp alp 0.0 1.0-2*alp 0.0 0.0 + 0.0 alp alp 0.0 1.0-2*alp 0.0 + alp 0.0 alp 0.0 0.0 1.0-2*alp + ] invTe = [ 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 @@ -539,7 +539,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T invT = SparseMatrixCOO() for element in slave_elements dofs = get_gdofs(problem, element) - for i=1:field_dim + for i = 1:field_dim ldofs = dofs[i:field_dim:end] add!(T, ldofs, ldofs, Te) add!(invT, ldofs, ldofs, invTe) @@ -556,8 +556,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T #@info("invT == invT2? ", invT == invT2) #maxabsdiff = maximum(abs(invT - invT2)) #@info("max diff = $maxabsdiff") - C1 = C1*invT - C2 = C2*invT + C1 = C1 * invT + C2 = C2 * invT end tol = problem.properties.drop_tolerance @@ -565,7 +565,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T SparseArrays.droptol!(C2, tol) for j in S - dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + dofs = [3 * (j - 1) + 1, 3 * (j - 1) + 2, 3 * (j - 1) + 3] weighted_gap[j] = g[dofs] end @@ -588,7 +588,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T # active / inactive node detection for j in S - dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + dofs = [3 * (j - 1) + 1, 3 * (j - 1) + 2, 3 * (j - 1) + 3] weighted_gap[j] = g[dofs] if length(la) != 0 normal = normals[j] @@ -642,17 +642,17 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T for j in S - dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] - tdofs = [3*(j-1)+2, 3*(j-1)+3] + dofs = [3 * (j - 1) + 1, 3 * (j - 1) + 2, 3 * (j - 1) + 3] + tdofs = [3 * (j - 1) + 2, 3 * (j - 1) + 3] if is_inactive[j] == 1 # remove inactive nodes from assembly - C1[dofs,:] .= 0.0 - C2[dofs,:] .= 0.0 - D[dofs,:] .= 0.0 - g[dofs,:] .= 0.0 + C1[dofs, :] .= 0.0 + C2[dofs, :] .= 0.0 + D[dofs, :] .= 0.0 + g[dofs, :] .= 0.0 elseif (is_active[j] == 1) && (is_slip[j] == 1) # constitutive modelling in tangent direction, frictionless contact - C2[tdofs,:] .= 0.0 + C2[tdofs, :] .= 0.0 g[tdofs] .= 0.0 normal = normals[j] tangent1, tangent2 = create_orthogonal_basis(normal) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index c9702d4..80affa6 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -35,10 +35,10 @@ https://en.wikipedia.org/wiki/Hooke's_law """ mutable struct Elasticity <: FieldProblem # these are found from problem.properties for type Problem{Elasticity} - formulation :: Symbol - finite_strain :: Bool - geometric_stiffness :: Bool - store_fields :: Vector{Symbol} + formulation::Symbol + finite_strain::Bool + geometric_stiffness::Bool + store_fields::Vector{Symbol} end function Elasticity() # formulations: plane_stress, plane_strain, continuum @@ -63,7 +63,7 @@ at time. This makes it possible to pre-allocate matrices common to same type of elements. """ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - elements::Vector{Element}, time) + elements::Vector{Element}, time) formulation = Val{problem.properties.formulation} for (element_type, elements_subset) in group_by_element_type(elements) assemble!(assembly, problem, elements_subset, time, formulation) @@ -71,50 +71,35 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - elements::Vector{T}, time, formulation) where {T <: Element} + elements::Vector{T}, time, formulation) where {T<:Element} - if problem.assemble_parallel - @assert problem.assemble_csc - # Threaded assembly - 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) - 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) - end - end - else - # Normal assembly - local_buffer = allocate_buffer(problem, elements) - 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, problem.assemble_csc) - end + # Normal assembly (parallel assembly disabled for now - needs property fields) + local_buffer = allocate_buffer(problem, elements) + for i in 1:length(elements) + assemble_element!(assembly, problem, elements[i], local_buffer, time, formulation) end end include("problems_elasticity_2d.jl") const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9} -const Elasticity3DVolumeElements = Union{Tet4, Pyr5, Wedge6, Wedge15, Hex8, Tet10, Hex20, Hex27} +const Elasticity3DVolumeElements = Union{Tet4,Pyr5,Wedge6,Wedge15,Hex8,Tet10,Hex20,Hex27} function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}}) param_keys = keys(params) all_keys = ip.fields.keys - ip_fields = filter(x->isassigned(all_keys, x), collect(1:length(all_keys))) + ip_fields = filter(x -> isassigned(all_keys, x), collect(1:length(all_keys))) if !("params_initialized" in ip_fields) for key in param_keys update!(ip, key, 0.0 => params[key]) end if type_ == Val{:type_2d} - update!(ip, "stress", 0.0 => [0.0,0.0,0.0]) - update!(ip, "strain", 0.0 => [0.0,0.0,0.0]) + update!(ip, "stress", 0.0 => [0.0, 0.0, 0.0]) + update!(ip, "strain", 0.0 => [0.0, 0.0, 0.0]) elseif type_ == Val{:type_3d} - update!(ip, "stress", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0]) - update!(ip, "strain", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0]) + update!(ip, "stress", 0.0 => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) + update!(ip, "strain", 0.0 => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) else error("daa") end @@ -123,71 +108,40 @@ 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) - 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) - 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) +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) + 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) + 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 allocate_buffer(problem::Problem{Elasticity}, ::Vector{Element{El}}) where El<:Elasticity3DVolumeElements dim = get_unknown_field_dimension(problem) nnodes = length(El) - ndofs = dim*nnodes + ndofs = dim * nnodes -<<<<<<< HEAD - return Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi = BasisInfo(El)) + return Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi=BasisInfo(El)) end -======= - for element in elements - - u = element("displacement", time) - X = element("geometry", time) - - fill!(Km, 0.0) - fill!(Kg, 0.0) - fill!(f_int, 0.0) - fill!(f_ext, 0.0) - - for ip in get_integration_points(element) - 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 - fill!(strain, 0.0) - fill!(F, 0.0) - F[:,:] += I - if props.finite_strain - strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) - F[:,:] += gradu - else - strain[:,:] = 1/2 * (gradu + gradu') - end ->>>>>>> master - function reset_element!(buf::Elasticity3DLocalBuffers) fill!(buf.Km, 0.0) fill!(buf.Kg, 0.0) @@ -205,47 +159,36 @@ function reset_integration_point!(buf::Elasticity3DLocalBuffers) return end -<<<<<<< HEAD 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] + 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 -======= - fill!(D, 0.0) - 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 - 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 ->>>>>>> master 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, - problem::Problem{Elasticity}, - element::Element{El}, - local_buffer::Elasticity3DLocalBuffers, - time, ::Type{Val{:continuum}}, - use_csc = false) where El<:Elasticity3DVolumeElements + assembler::FEMSparse.AssemblerSparsityPattern, + problem::Problem{Elasticity}, + element::Element{El}, + 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 + ndofs = dim * nnodes 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 + strain_vec, stress_vec, F, D, Dtan, Bt_mul_D, Bt_mul_D_mul_B, Bt_mul_S = local_buffer if !cheating u = element("displacement", time) X = element("geometry", time) @@ -255,7 +198,7 @@ function assemble_element!(assembly::Assembly, for ip in get_integration_points(element) reset_integration_point!(local_buffer) eval_basis!(bi, X, ip) - w = ip.weight*bi.detJ + 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 @@ -266,10 +209,10 @@ function assemble_element!(assembly::Assembly, F[i, i] += 1.0 end if props.finite_strain - strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu) - F[:,:] += gradu + 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) @@ -277,37 +220,37 @@ function assemble_element!(assembly::Assembly, # 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] + 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] + 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 @@ -320,11 +263,11 @@ function assemble_element!(assembly::Assembly, 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 - 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 + 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 @@ -380,7 +323,7 @@ function assemble_element!(assembly::Assembly, 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 + for i = 1:ndofs^2 @inbounds Km[i] += Bt_mul_D_mul_B[i] end @@ -388,28 +331,28 @@ function assemble_element!(assembly::Assembly, 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] + 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] + 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 + Kg += w * BNL' * S3 * BNL end @@ -421,11 +364,11 @@ function assemble_element!(assembly::Assembly, # external load start if haskey(element, "displacement load") T = element("displacement load", ip, time)::Vector{Float64} - mul!(f_buffer, w, vec(T*N)) + mul!(f_buffer, w, vec(T * N)) f_ext .+= f_buffer end - for i=1:dim + for i = 1:dim if haskey(element, displacement_load_string[i]) b = element(displacement_load_string[i], ip, time)::Float64 mul!(f_buffer_dim, w, N) @@ -462,44 +405,44 @@ end """ Elasticity equations, surface traction for continuum formulation. """ function assemble!(assembly::Assembly, - problem::Problem{Elasticity}, - elements::Vector{Element{El}}, - time, ::Type{Val{:continuum}}) where El<:Elasticity3DSurfaceElements + problem::Problem{Elasticity}, + elements::Vector{Element{El}}, + time, ::Type{Val{:continuum}}) where El<:Elasticity3DSurfaceElements props = problem.properties dim = get_unknown_field_dimension(problem) for element in elements nnodes = size(element, 2) - f = zeros(dim*nnodes) + f = zeros(dim * nnodes) has_concentrated_forces = false for ip in get_integration_points(element) detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ N = element(ip, time) if haskey(element, "displacement traction force") T = element("displacement traction force", ip, time) - f += w*vec(T*N) + f += w * vec(T * N) end for i in 1:dim if haskey(element, "displacement traction force $i") T = element("displacement traction force $i", ip, time) - f[i:dim:end] += w*vec(T*N) + f[i:dim:end] += w * vec(T * N) end if haskey(element, "concentrated force $i") has_concentrated_forces = true T = element("concentrated force $i", ip, time) - f[i:dim:end] += w*vec(T*N) + f[i:dim:end] += w * vec(T * N) end end if haskey(element, "surface pressure") J = element(ip, time, Val{:Jacobian})' - n = cross(J[:,1], J[:,2]) + n = cross(J[:, 1], J[:, 2]) n /= norm(n) # sign convention, positive pressure is towards surface p = -element("surface pressure", ip, time) - f += w*p*vec(n*N) + f += w * p * vec(n * N) end end if has_concentrated_forces @@ -518,11 +461,11 @@ Assemble all other elements for continuum elasticity problems. Basically, throw an exception telling to filter invalid elements out from the element set. """ function assemble!(assembly::Assembly, - problem::Problem{Elasticity}, - elements::Vector{Element{El}}, - time, ::Type{Val{:continuum}}) where El - @info("It looks that you are trying to assemble elements of type $El to 3d continuum "* - "problem. However, they are not supported yet. To filter out elements from a "* + problem::Problem{Elasticity}, + elements::Vector{Element{El}}, + time, ::Type{Val{:continuum}}) where El + @info("It looks that you are trying to assemble elements of type $El to 3d continuum " * + "problem. However, they are not supported yet. To filter out elements from a " * "element set, try `filter(element->!isa(element, Element{$El}), elements)`") error("Tried to assemble unsupported elements of type $El to 3d continuum problem.") end @@ -530,7 +473,7 @@ end """ Return strain tensor. """ function get_strain_tensor(problem, element, ip, time) gradu = element("displacement", ip, time, Val{:Grad}) - eps = 0.5*(gradu' + gradu) + eps = 0.5 * (gradu' + gradu) return eps end @@ -539,22 +482,22 @@ function get_stress_tensor(problem, element, ip, time) eps = get_strain_tensor(problem, element, ip, time) E = element("youngs modulus", ip, time) nu = element("poissons ratio", ip, time) - mu = E/(2.0*(1.0+nu)) - la = E*nu/((1.0+nu)*(1.0-2.0*nu)) - S = la*tr(eps)*I + 2.0*mu*eps + mu = E / (2.0 * (1.0 + nu)) + la = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) + S = la * tr(eps) * I + 2.0 * mu * eps return S end """ Return stain vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ function get_strain_vector(problem, element, ip, time) eps = get_strain_tensor(problem, element, ip, time) - return [eps[1,1], eps[2,2], eps[3,3], eps[1,2], eps[2,3], eps[1,3]] + return [eps[1, 1], eps[2, 2], eps[3, 3], eps[1, 2], eps[2, 3], eps[1, 3]] end """ Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ function get_stress_vector(problem, element, ip, time) S = get_stress_tensor(problem, element, ip, time) - return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]] + return [S[1, 1], S[2, 2], S[3, 3], S[1, 2], S[2, 3], S[1, 3]] end """ Make least squares fit for some field to nodes. """ @@ -566,22 +509,22 @@ function lsq_fit(problem, elements, field, time) gdofs = get_connectivity(element) for ip in get_integration_points(element) detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ + w = ip.weight * detJ N = element(ip, time) f = field(problem, element, ip, time) - add!(A, gdofs, gdofs, w*kron(N', N)) - for i=1:length(f) - add!(b, gdofs, w*f[i]*N, i) + add!(A, gdofs, gdofs, w * kron(N', N)) + for i = 1:length(f) + add!(b, gdofs, w * f[i] * N, i) end volume += w end end A = sparse(A) b = sparse(b) - A = 1/2*(A + A') + A = 1 / 2 * (A + A') nz = get_nonzero_rows(A) - F = ldlt(A[nz,nz]) + F = ldlt(A[nz, nz]) x = F \ b[nz, :]