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 index a9c3cd3..1b9a7f6 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ 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" @@ -18,12 +19,18 @@ 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" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" 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" + [compat] HDF5 = "≥ 0.7.0" LightXML = "≥ 0.4.0" diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index ddaecc0..82c9676 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -106,7 +106,9 @@ about JuliaFEM, please visit our website at module JuliaFEM using SparseArrays, LinearAlgebra, Statistics -using Reexport, ForwardDiff, LightXML, HDF5 +using Reexport, ForwardDiff, LightXML, HDF5, Parameters + +import FEMSparse @reexport using FEMBase import FEMBase: get_unknown_field_name, get_unknown_field_dimension, 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_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 13e680f..b63154b 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -70,6 +70,31 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, end end +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_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 + end +end + include("problems_elasticity_2d.jl") const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9} @@ -98,35 +123,42 @@ function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}}) end end -""" Assemble 3d continuum elements in general solid mechanics problem. """ -function assemble!(assembly::Assembly, - problem::Problem{Elasticity}, - elements::Vector{Element{El}}, - time, ::Type{Val{:continuum}}) where El<:Elasticity3DVolumeElements - props = problem.properties +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 - 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) +<<<<<<< HEAD + return Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi = BasisInfo(El)) +end +======= for element in elements u = element("displacement", time) @@ -154,54 +186,36 @@ function assemble!(assembly::Assembly, else strain[:,:] = 1/2 * (gradu + gradu') end +>>>>>>> master - 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] +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 - # 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] - 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 +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 +<<<<<<< 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] + return +end +======= fill!(D, 0.0) E = element("youngs modulus", ip, time)::Float64 nu = element("poissons ratio", ip, time)::Float64 @@ -210,133 +224,237 @@ function assemble!(assembly::Assembly, 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 - # determine material model +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 + cheating = false + props = problem.properties + dim = get_unknown_field_dimension(problem) - material_model = :linear_elasticity - if haskey(element, "plasticity") - material_model = :ideal_plasticity + nnodes = length(El) + 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 + if !cheating + u = element("displacement", time) + X = element("geometry", time) + end + reset_element!(local_buffer) + + for ip in get_integration_points(element) + reset_integration_point!(local_buffer) + 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 + 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') + 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 + + 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 + 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 + copyto!(Dtan, D) + mul!(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 - # 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 - - 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] - 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) + #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) + f_int .+= Bt_mul_S + + # external load start + if haskey(element, "displacement load") + T = element("displacement load", ip, time)::Vector{Float64} + mul!(f_buffer, w, vec(T*N)) + f_ext .+= f_buffer + end + + 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) + 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 + + 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 - add!(assembly.K, gdofs, gdofs, Km) + @inbounds FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext) + if props.geometric_stiffness + @inbounds 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 - - add!(assembly.f, gdofs, f_ext - f_int) - end return nothing @@ -390,7 +508,6 @@ function assemble!(assembly::Assembly, gdofs = get_gdofs(problem, element) add!(assembly.f, gdofs, f) - end end diff --git a/src/solvers.jl b/src/solvers.jl index 2b58d37..e3992ea 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -57,17 +57,23 @@ problems must have unique node ids. function get_field_assembly(solver::Solver) problems = get_field_problems(solver) - M = SparseMatrixCOO() - K = SparseMatrixCOO() - Kg = SparseMatrixCOO() - f = SparseMatrixCOO() - fg = SparseMatrixCOO() + problem = problems[1] + 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 + for problem in problems[2:end] append!(M, problem.assembly.M) append!(K, problem.assembly.K) append!(Kg, problem.assembly.Kg) 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 @@ -79,11 +85,11 @@ function get_field_assembly(solver::Solver) @warn("Field assembly seems to be empty. Check that elements are ", "pushed to problem and formulation is correct.") end - Kg = sparse(Kg, N, N) 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. """ @@ -167,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 eb7a173..01a490d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,11 @@ using JuliaFEM, Test +<<<<<<< HEAD +# include("../docs/make.jl") + +======= +>>>>>>> master @testset "JuliaFEM.jl" begin @testset "test_dirichlet.jl" begin include("test_dirichlet.jl") diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl index e0de522..b9b7ee1 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, 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