diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 3090950..3fc2b6e 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -8,91 +8,50 @@ This is JuliaFEM -- Finite Element Package """ module JuliaFEM +using FEMBase +using FEMBase: SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo, + Discrete, Variable, TimeVariant, TimeInvariant, Field, + DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment, + IP, AbstractProblem, IntegrationPoint +using FEMBase: is_field_problem, is_boundary_problem, get_elements, + get_connectivity, assemble_prehook!, assemble_posthook!, + get_parent_field_name, get_reference_coordinates, + get_assembly, get_nonzero_rows, get_nonzero_columns, + eval_basis!, get_basis, get_dbasis, grad!, get_dualbasis, + assemble_mass_matrix!, get_local_coordinates, inside, + get_element_type, filter_by_element_type, get_element_id, + optimize!, resize_sparse, resize_sparsevec +import FEMBase: get_unknown_field_name, get_unknown_field_dimension, + assemble!, update!, initialize! + + +# from other packages TimerOutputs.jl and Logging.jl using TimerOutputs export @timeit, print_timer - -import Base: getindex, setindex!, convert, length, size, isapprox, similar, - start, first, next, done, last, endof, vec, ==, +, -, *, /, haskey, copy, - push!, isempty, empty!, append!, sparse, full, read - -using FEMBasis -using FEMBasis: AbstractBasis -using FEMQuad -using AbaqusReader -using AsterReader - using Logging - -Logging.configure(level=INFO) - -if haskey(ENV, "JULIAFEM_LOGLEVEL") - Logging.configure(level=LogLevel(ENV["JULIAFEM_LOGLEVEL"])) -end - export info, debug +import Base: getindex, setindex!, convert, length, size, isapprox, + similar, start, first, next, done, last, endof, vec, + ==, +, -, *, /, haskey, copy, push!, isempty, empty!, + append!, sparse, full, read + module Testing - - using Base.Test - - export @test, @testset, @test_throws - +using Base.Test +export @test, @testset, @test_throws end -include("fields.jl") -export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment - -include("types.jl") # data types: Point, IntegrationPoint, ... -export AbstractPoint, Point, IntegrationPoint, IP, Node - -### ELEMENTS ### -include("elements.jl") # common element routines -export Node, Element, update!, get_connectivity, get_basis, - get_dbasis, inside, get_local_coordinates, get_element_type, - filter_by_element_type, get_element_id - -include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements -export get_reference_coordinates, - get_interpolation_polynomial, - description -export Poi1, - Seg2, Seg3, - Tri3, Tri6, Tri7, - Quad4, Quad8, Quad9, - Tet4, Tet10, - Pyr5, - Wedge6, - Hex8, Hex20, Hex27 - -include("integrate.jl") # default integration points for elements -export get_integration_points - -include("sparse.jl") -export add!, SparseMatrixCOO, SparseVectorCOO, get_nonzero_rows, get_nonzero_columns, optimize!, resize_sparse, resize_sparsevec - -include("problems.jl") # common problem routines -export Problem, AbstractProblem, FieldProblem, BoundaryProblem, - get_unknown_field_dimension, get_gdofs, Assembly, - get_parent_field_name, get_elements, add_elements! - +using AbaqusReader +using AsterReader include("problems_elasticity.jl") export Elasticity - include("materials_plasticity.jl") export plastic_von_mises - include("problems_dirichlet.jl") export Dirichlet - include("problems_heat.jl") export Heat - export assemble!, postprocess! - -function assemble!(problem::Problem, element::Element, time=0.0) - assemble!(problem.assembly, problem, element, time) -end - ### Mortar methods ### include("problems_mortar.jl") include("problems_mortar_2d.jl") @@ -101,12 +60,8 @@ include("problems_mortar_2d_autodiff.jl") export calculate_normals, calculate_normals!, project_from_slave_to_master, project_from_master_to_slave, Mortar, get_slave_elements, get_polygon_clip - include("io.jl") export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, update_xdmf!, save! - -### ASSEMBLY + SOLVE ### -include("assembly.jl") include("solvers.jl") export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver, get_unknown_field_name, get_formulation_type, get_problems, @@ -116,8 +71,6 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver, is_field_problem, is_boundary_problem include("solvers_modal.jl") export Modal - -### Mortar methods, contact mechanics extension ### include("problems_contact.jl") include("problems_contact_2d.jl") include("problems_contact_3d.jl") @@ -125,22 +78,26 @@ include("problems_contact_2d_autodiff.jl") #include("problems_contact_3d_autodiff.jl") export Contact +# Preprocess module + module Preprocess +using FEMBase include("preprocess.jl") export create_elements, Mesh, add_node!, add_nodes!, add_element!, add_elements!, add_element_to_element_set!, add_node_to_node_set!, find_nearest_nodes, find_nearest_node, reorder_element_connectivity!, create_node_set_from_element_set!, filter_by_element_set - include("preprocess_abaqus_reader.jl") export abaqus_read_mesh, create_surface_elements, create_nodal_elements - include("preprocess_aster_reader.jl") export aster_read_mesh - end +# Postprocess module + module Postprocess +using FEMBase +using FEMBase: get_elements include("postprocess_utils.jl") export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!, calculate_area, calculate_center_of_mass, @@ -149,4 +106,18 @@ end include("deprecations.jl") +export SparseMatrixCOO, SparseVectorCOO, optimize!, resize_sparse +export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment +export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly +export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9, + Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27 +export update!, add_elements!, get_unknown_field_name, add!, + is_field_problem, is_boundary_problem, get_gdofs, + initialize!, get_integration_points, group_by_element_type, + get_unknown_field_dimension, get_connectivity + +export get_nonzero_rows, get_local_coordinates, inside, IP, get_element_type, + get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type, + get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec + end diff --git a/src/assembly.jl b/src/assembly.jl deleted file mode 100644 index 4756283..0000000 --- a/src/assembly.jl +++ /dev/null @@ -1,173 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -function isapprox(a1::Assembly, a2::Assembly) - T = isapprox(a1.K, a2.K) - T &= isapprox(a1.C1, a2.C1) - T &= isapprox(a1.C2, a2.C2) - T &= isapprox(a1.D, a2.D) - T &= isapprox(a1.f, a2.f) - T &= isapprox(a1.g, a2.g) - return T -end - -function assemble_prehook! -end - -function assemble_posthook! -end - -function assemble!(problem::Problem, time=0.0; auto_initialize=true) - if !isempty(problem.assembly) - warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?") - end - if isempty(problem.elements) - warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?") - else - first_element = first(problem.elements) - unknown_field_name = get_unknown_field_name(problem) - if !haskey(first_element, unknown_field_name) - warn("Assemble problem $(problem.name): seems that problem is uninitialized.") - if auto_initialize - info("Initializing problem $(problem.name) at time $time automatically.") - initialize!(problem, time) - end - end - end - if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64}) - assemble_prehook!(problem, time) - end - assemble!(get_assembly(problem), problem, get_elements(problem), time) - if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64}) - assemble_posthook!(problem, time) - end - return true -end - -function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time) - warn("assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!") - for element in elements - assemble!(assembly, problem, element, time) - end - return nothing -end - -function assemble_mass_matrix!(problem::Problem, time) - if !isempty(problem.assembly.M) - info("Mass matrix for $(problem.name) is already assembled, skipping assemble routine") - return - end - elements = get_elements(problem) - for (element_type, elements) in group_by_element_type(get_elements(problem)) - assemble_mass_matrix!(problem::Problem, elements, time) - end - return -end - -function assemble_mass_matrix!{Basis}(problem::Problem, elements::Vector{Element{Basis}}, time) - nnodes = length(Basis) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - for element in elements - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", time) - w = ip.weight*rho*detJ - eval_basis!(Basis, N, ip) - N = element(ip, time) - At_mul_B!(NtN, N, N) - scale!(NtN, w) - for i=1:nnodes^2 - M[i] += NtN[i] - end - end - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j-1)*dim - end - for i=1:dim - add!(problem.assembly.M, ldofs+i, ldofs+i, M) - end - end - return -end - -""" - assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time) - -Assemble Tet10 mass matrices using special method. If Tet10 has constant metric -if can be integrated analytically to gain performance. -""" -function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Tet10}}, time) - nnodes = length(Tet10) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - - M_CM = 1.0/2520.0 * [ - 6 1 1 1 -4 -6 -4 -4 -6 -6 - 1 6 1 1 -4 -4 -6 -6 -4 -6 - 1 1 6 1 -6 -4 -4 -6 -6 -4 - 1 1 1 6 -6 -6 -6 -4 -4 -4 - -4 -4 -6 -6 32 16 16 16 16 8 - -6 -4 -4 -6 16 32 16 8 16 16 - -4 -6 -4 -6 16 16 32 16 8 16 - -4 -6 -6 -4 16 8 16 32 16 16 - -6 -4 -6 -4 16 16 8 16 32 16 - -6 -6 -4 -4 8 16 16 16 16 32] - - function is_CM(element::Element{Tet10}, X; rtol=1.0e-6) - isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false - isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false - isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false - isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false - isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false - isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false - return true - end - - - n_CM = 0 - for element in elements - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j-1)*dim - end - - X = element("geometry", time) - rho = element("density", time) - if is_CM(element, X) && length(rho) == 1 - ip = (1.0/3.0, 1.0/3.0, 1.0/3.0) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - CM_s = detJ*rho - n_CM += 1 - for i=1:dim - add!(problem.assembly.M, ldofs+i, ldofs+i, CM_s * M_CM) - end - else - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - w = ip.weight*rho*detJ - eval_basis!(Tet10, N, ip) - N = element(ip, time) - At_mul_B!(NtN, N, N) - scale!(NtN, w) - for i=1:nnodes^2 - M[i] += NtN[i] - end - end - for i=1:dim - add!(problem.assembly.M, ldofs+i, ldofs+i, M) - end - end - end - info("$n_CM of $(length(elements)) was constant metric.") - return -end diff --git a/src/deprecations.jl b/src/deprecations.jl index 1c3fb7f..29de2b4 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -5,6 +5,10 @@ function assemble!(problem::Problem, time, ::Type{Val{:mass_matrix}}) assemble_mass_matrix!(problem, time) end +function assemble!(problem::Problem, element::Element, time=0.0) + assemble!(problem.assembly, problem, element, time) +end + module Abaqus using JuliaFEM.Preprocess: create_surface_elements end diff --git a/src/elements.jl b/src/elements.jl deleted file mode 100644 index 4b48e67..0000000 --- a/src/elements.jl +++ /dev/null @@ -1,418 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -type Element{E<:AbstractBasis} - id :: Int - connectivity :: Vector{Int} - integration_points :: Vector{IP} - fields :: Dict{String, Field} - properties :: E -end - -""" - Element(element_type, connectivity_vector) - -Construct a new element where element_type is the type of the element -and connectivity_vector is the vector of nodes that the element is connected to. - -Examples --------- -In the example a new element (E in the figure below) of type Tri3 is created. -This spesific element connects to nodes 89, 43, 12 in the finite element mesh. - -```@example -element = Element(Tri3, [89, 43, 12]) -``` -![img](figs/mesh.png) -""" -function Element{E<:AbstractBasis}(::Type{E}, connectivity::Vector{Int}) - return Element{E}(-1, connectivity, [], Dict(), E()) -end - -""" - length(element::Element) - -Return the number of nodes in element. -""" -function length{B}(element::Element{B}) - return length(B) -end - -function size{B}(element::Element{B}) - return size(B) -end - -function getindex(element::Element, field_name::AbstractString) - return element.fields[field_name] -end - -function setindex!(element::Element, data::Field, field_name) - element.fields[field_name] = data -end - -function get_element_type{E}(element::Element{E}) - return E -end - -function get_element_id{E}(element::Element{E}) - return element.id -end - -function is_element_type{E}(element::Element{E}, element_type) - return E === element_type -end - -function filter_by_element_type(element_type, elements) - return filter(element -> is_element_type(element, element_type), elements) -end - -""" - group_by_element_type(elements::Vector{Element}) - -Given a vector of elements, group elements by element type to several vectors. -Returns a dictionary, where key is the element type and value is a vector -containing all elements of type `element_type`. -""" -function group_by_element_type(elements::Vector{Element}) - results = Dict{DataType, Any}() - basis_types = map(element -> typeof(element.properties), elements) - for basis in unique(basis_types) - element_type = Element{basis} - subset = filter(element -> isa(element, element_type), elements) - results[element_type] = convert(Vector{element_type}, subset) - end - return results -end - -function setindex!(element::Element, data::Function, field_name) - if method_exists(data, Tuple{Element, Vector, Float64}) - # create enclosure to pass element as argument - function wrapper_(ip, time) - return data(element, ip, time) - end - field = Field(wrapper_) - else - field = Field(data) - end - element.fields[field_name] = field -end - -function setindex!(element::Element, data, field_name) - element.fields[field_name] = Field(data) -end - -""" Return a Field object from element. - -Examples --------- ->>> element = Element(Seg2, [1, 2]) ->>> data = Dict(1 => 1.0, 2 => 2.0) ->>> update!(element, "my field", data) ->>> element("my field") - -""" -function (element::Element)(field_name::String) - return element[field_name] -end - -""" Return a Field object from element and interpolate in time direction. - -Examples --------- ->>> element = Element(Seg2, [1, 2]) ->>> data1 = Dict(1 => 1.0, 2 => 2.0) ->>> data2 = Dict(1 => 2.0, 2 => 3.0) ->>> update!(element, "my field", 0.0 => data1, 1.0 => data2) ->>> element("my field", 0.5) - -""" -function (element::Element)(field_name::String, time::Float64) - return element[field_name](time) -end - -function last(element::Element, field_name::String) - return last(element[field_name]) -end - -function (element::Element)(ip, time::Float64=0.0) - return get_basis(element, ip, time) -end - -""" -Examples - -julia> el = Element(Quad4, [1, 2, 3, 4]); - -julia> el([0.0, 0.0], 0.0, 1) -1x4 Array{Float64,2}: - 0.25 0.25 0.25 0.25 - -julia> el([0.0, 0.0], 0.0, 2) -2x8 Array{Float64,2}: - 0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0 - 0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25 - -""" -function (element::Element)(ip, time::Float64, dim::Int) - dim == 1 && return get_basis(element, ip, time) - Ni = vec(get_basis(element, ip, time)) - N = zeros(dim, length(element)*dim) - for i=1:dim - N[i,i:dim:end] += Ni - end - return N -end - -function (element::Element)(ip, time::Float64, ::Type{Val{:Jacobian}}) - X = element("geometry", time) - dN = get_dbasis(element, ip, time) - nbasis = length(element) - if isa(X.data, Vector) - J = sum([kron(dN[:,i], X[i]') for i=1:nbasis]) - else - c = get_connectivity(element) - J = sum([kron(dN[:,i], X[c[i]]') for i=1:nbasis]) - end - return J -end - -function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}}) - J = element(ip, time, Val{:Jacobian}) - n, m = size(J) - if n == m # volume element - return det(J) - end - JT = transpose(J) - if size(JT, 2) == 1 # boundary of 2d problem, || ∂X/∂ξ || - return norm(JT) - else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ || - return norm(cross(JT[:,1], JT[:,2])) - end -end - -function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}}) - J = element(ip, time, Val{:Jacobian}) - return inv(J)*get_dbasis(element, ip, time) -end - -function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}}) - return element(ip, time, Val{:Grad})*element[field_name](time) -end - -function (element::Element)(field_name::String, ip, time::Float64) - field = element[field_name] - return element(field, ip, time) -end - -function (element::Element)(field::DCTI, ip, time::Float64) - return field.data -end - -function (element::Element)(field::DCTV, ip, time::Float64) - return field(time).data -end - -function (element::Element)(field::CVTV, ip, time::Float64) - return field(ip, time) -end - -function (element::Element)(field::Field, ip, time::Float64) - field_ = field(time) - basis = element(ip, time) - n = length(element) - if isa(field_.data, Vector) - m = length(field_) - if n != m - error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!") - end - return sum([field_[i]*basis[i] for i=1:n]) - else - c = get_connectivity(element) - return sum([field_[c[i]]*basis[i] for i=1:n]) - end -end - -function size(element::Element, dim) - return size(element)[dim] -end - -""" Update element field based on a dictionary of nodal data and connectivity information. - -Examples --------- -julia> data = Dict(1 => [0.0, 0.0], 2 => [1.0, 2.0]) -julia> element = Seg2([1, 2]) -julia> update!(element, "geometry", data) - -As a result element now have time invariant (variable) vector field "geometry" with data ([0.0, 0.0], [1.0, 2.0]). - -""" -function update!{E}(element::Element{E}, field_name::AbstractString, data::Dict) - #element[field_name] = Field(data) - element_id = element.id - local_connectivity = get_connectivity(element) - for i in local_connectivity - if !haskey(data, i) - ndata = length(data) - critical("Unable to set field data $field_name for element $E with - id $element_id and connectivity $local_connectivity: no data for - node id $i found. Length of data dict = $ndata") - end - end - local_data = [data[i] for i in local_connectivity] - element[field_name] = local_data -end - -function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}}) - time, field_data = data - element_data = V[field_data[i] for i in get_connectivity(element)] - update!(element, field_name, time => element_data) -end - -function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Any}}) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end -end - -function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Int64}}) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end -end - -function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Float64}}) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end -end - -function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Vector{Float64}}}) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end -end - -function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Float64}) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end -end - -function update!(element::Element, field_name::AbstractString, data::Union{Float64, Vector}) - if haskey(element, field_name) - update!(element[field_name], data) - else - if length(data) != length(element) - update!(element, field_name, DCTI(data)) - else - element[field_name] = data - end - end -end - -function update!(element::Element, datas::Pair...) - for (field_name, data) in datas - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = data - end - end -end - -function update!(element::Element, field_name::String, data::Function) - element[field_name] = data -end - -function update!(element::Element, field_name::String, field::Field) - element[field_name] = field -end - -function update!(elements::Vector, field_name::String, data) - for element in elements - update!(element, field_name, data) - end -end - -""" Check existence of field. """ -function haskey(element::Element, field_name::String) - haskey(element.fields, field_name) -end - -function get_connectivity(element::Element) - return element.connectivity -end - -function get_integration_points{E}(element::Element{E}) - # first time initialize default integration points - if length(element.integration_points) == 0 - ips = get_integration_points(element.properties) - if E in (Poi1, Seg2, Seg3, NSeg) - element.integration_points = [IP(i, w, (xi,)) for (i, (w, xi)) in enumerate(ips)] - else - element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)] - end - end - return element.integration_points -end - -""" This is a special case, temporarily change order -of integration scheme mainly for mass matrix. -""" -function get_integration_points{E}(element::Element{E}, change_order::Int) - ips = get_integration_points(element.properties, Val{change_order}) - if E in (Poi1, Seg2, Seg3, NSeg) - return [IP(i, w, (xi,)) for (i, (w, xi)) in enumerate(ips)] - else - return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)] - end -end - -""" Return dual basis transformation matrix Ae. """ -function get_dualbasis(element::Element, time::Float64, order=1) - nnodes = length(element) - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(element, order) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - De += w*diagm(vec(N)) - Me += w*N'*N - end - return De, Me, De*inv(Me) -end - -""" Find inverse isoparametric mapping of element. """ -function get_local_coordinates(element::Element, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6) - haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping") - dim = size(element, 1) - dim == length(X) || error("manifolds not supported.") - xi = zeros(dim) - dX = element("geometry", xi, time) - X - for i=1:max_iterations - J = element(xi, time, Val{:Jacobian})' - xi -= J \ dX - dX = element("geometry", xi, time) - X - norm(dX) < tolerance && return xi - end - info("X = $X, dX = $dX, xi = $xi") - error("Unable to find inverse isoparametric mapping for element $element for X = $X") -end - -""" Test is X inside element. """ -function inside{E}(element::Element{E}, X, time) - xi = get_local_coordinates(element, X, time) - return inside(E, xi) -end diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl deleted file mode 100644 index 121c97e..0000000 --- a/src/elements_lagrange.jl +++ /dev/null @@ -1,70 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using FEMBasis - -# "Poi1" => (0, 1), -"1 node discrete point element", -type Poi1 <: AbstractBasis -end - -function get_basis(element::Element{Poi1}, ip, time) - return [1] -end - -function get_dbasis(element::Element{Poi1}, ip, time) - return [0] -end - -function (element::Element{Poi1})(ip, time::Float64, ::Type{Val{:detJ}}) - return 1.0 -end - -function get_integration_order(element::Poi1) - return 1 -end - -function get_integration_points(element::Poi1, order::Int64) - return [ (1.0, 0.0) ] -end - -function size(::Type{Poi1}) - return (0, 1) -end - -function length(::Type{Poi1}) - return 1 -end - -function FEMBasis.get_reference_element_coordinates(::Type{Poi1}) - Vector{Float64}[[0.0]] -end - -function get_basis{B}(element::Element{B}, ip, time) - T = typeof(first(ip)) - N = zeros(T, 1, length(B)) - eval_basis!(B, N, tuple(ip...)) - return N -end - -function get_dbasis{B}(element::Element{B}, ip, time) - T = typeof(first(ip)) - dN = zeros(T, size(B)...) - eval_dbasis!(B, dN, tuple(ip...)) - return dN -end - -function inside(::Union{Type{Seg2}, Type{Seg3}, Type{Quad4}, Type{Quad8}, - Type{Quad9}, Type{Pyr5}, Type{Hex8}, Type{Hex20}, - Type{Hex27}}, xi) - return all(-1.0 .<= xi .<= 1.0) -end - -function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet10}}, xi) - return all(xi .>= 0.0) && (sum(xi) <= 1.0) -end - -function get_reference_coordinates{B}(element::Element{B}) - return get_reference_element_coordinates(B) -end - diff --git a/src/fields.jl b/src/fields.jl deleted file mode 100644 index 560e0a8..0000000 --- a/src/fields.jl +++ /dev/null @@ -1,457 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -abstract type AbstractField end - -abstract type Discrete<:AbstractField end -abstract type Continuous<:AbstractField end -abstract type Constant<:AbstractField end -abstract type Variable<:AbstractField end -abstract type TimeVariant<:AbstractField end -abstract type TimeInvariant<:AbstractField end - -type Field{A<:Union{Discrete, Continuous}, - B<:Union{Constant, Variable}, - C<:Union{TimeVariant, TimeInvariant}, - T} - data :: T -end - -const FieldSet = Dict{String, Field} - -### Different field combinations and other typealiases - -const DCTI{T} = Field{Discrete, Constant, TimeInvariant, T} -const DVTI{T} = Field{Discrete, Variable, TimeInvariant, T} -const DCTV{T} = Field{Discrete, Constant, TimeVariant, T} -const DVTV{T} = Field{Discrete, Variable, TimeVariant, T} -const CCTI{T} = Field{Continuous, Constant, TimeInvariant, T} -const CVTI{T} = Field{Continuous, Variable, TimeInvariant, T} -const CCTV{T} = Field{Continuous, Constant, TimeVariant, T} -const CVTV{T} = Field{Continuous, Variable, TimeVariant, T} - -# Discrete fields - - -""" Discrete, constant, time-invariant field. This is constant in both spatial -direction and time direction, i.e. df/dX = 0 and df/dt = 0. - -This is the most basic type of field having no anything special functionality. - -Examples --------- - -julia> f = DCTI() -julia> update!(f, 1.0) - -Multiplying by constant works: - -julia> 2*f -2.0 - -Interpolation in time direction gives the same constant: - -julia> f(1.0) -1.0 - -By default, when calling Field with scalar, DCTI is assumed, i.e. - -julia> Field(0.0) == DCTI(0.0) -true - -""" -function DCTI() - return DCTI(nothing) -end - -function DCTI{T}(a::T) - return DCTI{T}(a) -end - -function DVTI{T}(a::T) - return DVTI{T}(a) -end - -function DCTV{T}(a::T) - return DCTV{T}(a) -end - -function DVTV{T}(a::T) - return DVTV{T}(a) -end - -function CCTI{T}(a::T) - return CCTI{T}(a) -end - -function CVTI{T}(a::T) - return CVTI{T}(a) -end - -function CCTV{T}(a::T) - return CCTV{T}(a) -end - -function CVTV{T}(a::T) - return CVTV{T}(a) -end - -function Field() - return DCTI() -end - -function Field{T}(data::T) - return DCTI{T}(data) -end - -function ==(x::DCTI, y::DCTI) - return ==(x.data, y.data) -end - -function ==(x::DCTI, y) - return ==(x.data, y) -end - -function isapprox(x::DCTI, y::DCTI) - isapprox(x.data, y.data) -end - -function isapprox(x::DCTI, y) - isapprox(x.data, y) -end - -function length(f::DCTI) - return 1 -end - -function *(c::Number, f::DCTI) - return c*f.data -end - -""" Kind of spatial interpolation of DCTI. """ -function *(N::Matrix, f::DCTI) - @assert length(N) == 1 - return N[1]*f.data -end - -function update!(field::DCTI, data) - field.data = data -end - -""" Interpolate time-invariant field in time direction. """ -function (field::DCTI)(time::Float64) - return field.data -end - -""" Discrete, variable, time-invariant field. This is constant in time direction, -but not in spatial direction, i.e. df/dt = 0 but df/dX != 0. The basic structure -of data is Vector, and it is implicitly assumed that length of field matches to -the number of shape functions, so that interpolation in spatial direction works. - -Examples --------- -""" -function DVTI() - return DVTI([]) -end - -""" For vector data, DVTI is automatically created. - -julia> DVTI([1.0, 2.0]) == Field([1.0, 2.0]) -true - -""" -function Field(data::Vector) - return DVTI(data) -end - -""" For dictionary data, DVTI is automatically created. - -Define e.g. nodal coordinates in dictionary -julia> X = Dict(1 => [1.0, 2.0], 2 => [3.0, 4.0]) -julia> Field(X) == DVTI(X) - -""" -function Field(data::Dict) - return DVTI(data) -end - -function ==(x::DVTI, y::DVTI) - return ==(x.data, y.data) -end - -function isapprox(x::DVTI, y) - return isapprox(x.data, y) -end - -""" Default slicing of field. - -julia> f = DVTI([1.0, 2.0]) -julia> f[1] -1.0 - -""" -function getindex(field::DVTI, i::Int64) - return field.data[i] -end - -""" Multi-slicing of field. - -julia> f = DVTI([1.0, 2.0, 3.0]) -julia> f[[1, 3]] -[1.0, 3.0] - -""" -function getindex(field::DVTI, I::Array{Int64, 1}) - return [field.data[i] for i in I] -end - -function length(field::DVTI) - return length(field.data) -end - -function start(field::DVTI) - return 1 -end - -function +(f1::DVTI, f2::DVTI) - return DVTI(f1.data + f2.data) -end - -function -(f1::DVTI, f2::DVTI) - return DVTI(f1.data - f2.data) -end - -function update!(field::DVTI, data::Union{Vector, Dict}) - field.data = data -end - -""" Take scalar product of DVTI and constant T. """ -function *(T::Number, field::DVTI) - return DVTI(T*field.data) -end - -""" Take dot product of DVTI field and vector T. Vector length must match to the -field length and this can be used mainly for interpolation purposes, i.e., u = ∑ Nᵢuᵢ. -""" -function *(T::Union{Vector, RowVector}, f::DVTI) - @assert length(T) <= length(f) - return sum([T[i]*f[i] for i=1:length(T)]) -end - -""" Take outer product of DVTI field and matrix T. """ -function *(T::Matrix, f::DVTI) - n, m = size(T) - return sum([kron(T[:,i], f[i]') for i=1:m])' -end - -function vec(field::DVTI) - return [field.data...;] -end - -""" Interpolate time-invariant field in time direction. """ -function (field::DVTI)(time::Float64) - return field -end - -""" Create a similar DVTI field from vector data. - -julia> f1 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]]) -julia> f2 = similar(f1, [2.0, 3.0, 4.0, 5.0]) -julia> f2 == DVTI(Vector[[2.0, 3.0], [4.0, 5.0]]) -true - -""" -function similar(field::DVTI, data::Vector) - n = length(field) - m = length(data) - dim = round(Int, m/n) - @assert dim*n == m - new_data = reshape(data, dim, n) - new_field = DVTI() - new_field.data = [new_data[:,i] for i=1:n] - return new_field -end - -function next(f::DVTI, state) - return f.data[state], state+1 -end - -function done(f::DVTI, s) - return s > length(f.data) -end - -""" Simple time frame / increment to contain both time and data. """ -type Increment{T} - time :: Float64 - data :: T -end - -""" Discrete, constant, time variant field. This is constant in spatial -direction but non-constant in time direction, i.e. df/dX = 0 but df/dt != 0. - -Examples --------- -julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0 -julia> f = DCTV(t0 => y0, t1 => y1) - -""" -function DCTV(data::Pair...) - return DCTV([Increment(d[1],d[2]) for d in data]) -end - -function Field{T}(data::Pair{Float64, T}...) - return DCTV([Increment{T}(d[1], d[2]) for d in data]) -end - -function getindex(field::DCTV, i::Int64) - return field.data[i] -end - -function length(field::DCTV) - return length(field.data) -end - -function first(field::DCTV) - return field[1] -end - -""" Interpolate constant time-variant field in time direction. """ -function (field::DCTV)(time::Number) - time < first(field).time && return DCTI(first(field).data) - time > last(field).time && return DCTI(last(field).data) - for i=reverse(1:length(field)) - isapprox(field[i].time, time) && return DCTI(field[i].data) - end - for i=reverse(2:length(field)) - t0 = field[i-1].time - t1 = field[i].time - if t0 < time < t1 - y0 = field[i-1].data - y1 = field[i].data - dt = t1-t0 - new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) - return DCTI(new_data) - end - end -end - -function endof(field::DCTV) - return endof(field.data) -end - -""" Discrete, variable, time variant fields. """ -function DVTV() - return DVTV(Increment[]) -end - -function DVTV{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...) - return DVTV([Increment{T}(d[1], d[2]) for d in data]) -end - -function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...) - return DVTV([Increment{T}(d[1], d[2]) for d in data]) -end - -function length(field::DVTV) - return length(field.data) -end - -function getindex(field::DVTV, i::Int64) - return field.data[i] -end - -function first(field::DVTV) - return field[1] -end - -function endof(field::DVTV) - return endof(field.data) -end - -""" Interpolate discrete, variable, time-variant field in time direction. """ -function (field::DVTV)(time::Float64) - time < first(field).time && return DVTI(first(field).data) - time > last(field).time && return DVTI(last(field).data) - for i=reverse(1:length(field)) - isapprox(field[i].time, time) && return DVTI(field[i].data) - end - for i=reverse(2:length(field)) - t0 = field[i-1].time - t1 = field[i].time - if t0 < time < t1 - y0 = field[i-1].data - y1 = field[i].data - dt = t1-t0 - new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) - return DVTI(new_data) - end - end -end - -""" Update time-dependent fields with new values. - -Examples --------- - -julia> f = Field(0.0 => 1.0) -julia> update!(f, 1.0 => 2.0) - -Now field has two (time, value) pairs: (0.0, 1.0) and (1.0, 2.0) - -Notes ------ -Time vector is assumed to be ordered t_i-1 < t_i < t_i+1. If updating -field with already existing time the old value is replaced with new one. - -""" -function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T}) - time, data = val - if isapprox(last(field).time, time) - last(field).data = data - else - push!(field.data, Increment(val...)) - end -end - -### Basic data structure for continuous field - -type Basis - basis :: Function - dbasis :: Function -end - -### Convenient functions to create fields - -function Field(func::Function) - if method_exists(func, Tuple{}) - return CCTI(func) - elseif method_exists(func, Tuple{Float64}) - return CCTV(func) - elseif method_exists(func, Tuple{Vector}) - return CVTI(func) - elseif method_exists(func, Tuple{Vector, Float64}) - return CVTV(func) - else - error("no proper definition found for function: check methods.") - end -end - -### Accessing continuous fields - -function (field::CCTI)(xi::Vector, time::Number) - return field.data() -end - -function (field::CVTI)(xi::Vector, time::Number) - return field.data(xi) -end - -function (field::CCTV)(xi::Vector, time::Number) - return field.data(time) -end - -function (field::CVTV)(xi, time) - return field.data(xi, time) -end - diff --git a/src/integrate.jl b/src/integrate.jl deleted file mode 100644 index 2dda5bc..0000000 --- a/src/integrate.jl +++ /dev/null @@ -1,58 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using FEMQuad: get_quadrature_points - -# Default number of integration points for each element. First rule is the -# default integration rule returned by `get_integration_points(element)`. -# Sometimes we want to increase integration order, e.g. when integrating mass -# matrix or boundary conditions. For that reason, additional rules are provied -# in list, so e.g. `get_integration_points(element, 1)` returns the second rule, -# `get_integration_points(element, 2)` third rule and so on. Rules should be -# ordered so that picking next one integrates more accurately. -integration_rule_mapping = ( - :Seg2 => (:GLSEG1, :GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5), - :Seg3 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5), - :NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5), - :Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25), - :Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :Hex8 => (:GLHEX8, :GLHEX27, :GLHEX81, :GLHEX243), - :Hex20 => (:GLHEX27, :GLHEX81, :GLHEX243), - :Hex27 => (:GLHEX27, :GLHEX81, :GLHEX243), - :NSolid => (:GLHEX27, :GLHEX81, :GLHEX243), - :Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15), - :Tet10 => (:GLTET4, :GLTET5, :GLTET15), - :Pyr5 => (:GLPYR5, ), - :Wedge6 => (:GLWED6, :GLWED21), - :Wedge15 => (:GLWED21, )) - -for (E, R) in integration_rule_mapping - for i in 1:length(R) - P = Val{R[i]} - order = Val{i-1} - if i == 1 - code = quote - function get_integration_points(element::$E) - return get_quadrature_points($P) - end - end - else - code = quote - function get_integration_points(element::$E, ::Type{$order}) - return get_quadrature_points($P) - end - end - end - eval(code) - end -end - -# All good codes needs a special case. Here we have it: Poi1 -function get_integration_points(element::Poi1) - [ (1.0, 0.0) ] -end diff --git a/src/problems.jl b/src/problems.jl deleted file mode 100644 index cd26991..0000000 --- a/src/problems.jl +++ /dev/null @@ -1,442 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -abstract type AbstractProblem end -abstract type FieldProblem<:AbstractProblem end -abstract type BoundaryProblem<:AbstractProblem end -abstract type MixedProblem<:AbstractProblem end - -""" -General linearized problem to solve - (K₁+K₂)Δu + C1'*Δλ = f₁+f₂ - C2Δu + D*Δλ = g -""" -type Assembly - - M :: SparseMatrixCOO # mass matrix - - # for field assembly - K :: SparseMatrixCOO # stiffness matrix - Kg :: SparseMatrixCOO # geometric stiffness matrix - f :: SparseMatrixCOO # force vector - fg :: SparseMatrixCOO # - - # for boundary assembly - C1 :: SparseMatrixCOO - C2 :: SparseMatrixCOO - D :: SparseMatrixCOO - g :: SparseMatrixCOO - c :: SparseMatrixCOO - - u :: Vector{Float64} # solution vector u - u_prev :: Vector{Float64} # previous solution vector u - u_norm_change :: Real # change of norm in u - - la :: Vector{Float64} # solution vector la - la_prev :: Vector{Float64} # previous solution vector u - la_norm_change :: Real # change of norm in la - - removed_dofs :: Vector{Int64} # manually remove dofs from assembly -end - -function Assembly() - return Assembly( - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - [], [], Inf, - [], [], Inf, - []) -end - -function empty!(assembly::Assembly) - empty!(assembly.K) - empty!(assembly.Kg) - empty!(assembly.f) - empty!(assembly.fg) - empty!(assembly.C1) - empty!(assembly.C2) - empty!(assembly.D) - empty!(assembly.g) - empty!(assembly.c) -end - -function isempty(assembly::Assembly) - T = isempty(assembly.K) - T &= isempty(assembly.Kg) - T &= isempty(assembly.f) - T &= isempty(assembly.fg) - T &= isempty(assembly.C1) - T &= isempty(assembly.C2) - T &= isempty(assembly.D) - T &= isempty(assembly.g) - T &= isempty(assembly.c) - return T -end - -""" -Defines types for Problem variables. - -# Examples - -The type of 'elements' is Vector{Element} - -Add elements into the Problem element list. -```@example -a = [1, 2, 3] -Problem.elements = a -``` - -""" -type Problem{P<:AbstractProblem} - name :: AbstractString # descriptive name for the problem - dimension :: Int # degrees of freedom per node - parent_field_name :: AbstractString # (optional) name of the parent field e.g. "displacement" - elements :: Vector{Element} - dofmap :: Dict{Element, Vector{Int64}} # connects the element local dofs to the global dofs - assembly :: Assembly - fields :: Dict{AbstractString, Field} - postprocess_fields :: Vector{String} - properties :: P -end - -""" - Problem(problem_type, problem_name::String, problem_dimension) - -Construct a new field problem where `problem_type` is the type of the problem -(Elasticity, Dirichlet, etc.), `problem_name` is the name of the problem and -`problem_dimension` is the number of DOF:s in one node (2 in a 2D problem, 3 -in an elastic 3D problem, 6 in a 3D beam problem, etc.). - -# Examples - -Create a vector-valued (dim=3) elasticity problem: - -```@example -prob1 = Problem(Elasticity, "this is my problem", 3) -``` - -""" -function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64) - return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P()) -end - -""" -Construct a new boundary problem. - -Examples --------- -Create a Dirichlet boundary problem for a vector-valued (dim=3) elasticity problem. - -julia> bc1 = Problem(Dirichlet, "support", 3, "displacement") -solver. -""" -function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name) - return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P()) -end - -function get_formulation_type(problem::Problem) - return :incremental -end - -function get_unknown_field_name{P<:BoundaryProblem}(::Type{P}) - return "lambda" -end - -function get_assembly(problem) - return problem.assembly -end - -""" - update!(problem.properties, attr...) - -Update properties for a problem. - -# Example - -```julia -update!(body.properties, "finite_strain" => "false") -``` -""" -function update!{P<:AbstractProblem}(problem::P, attr::Pair{String, String}...) - for (name, value) in attr - debug("$P: set $name to $value") - setfield!(problem, parse(name), parse(value)) - end -end - -""" - function initialize!(problem_type, element_name, time) - -Initialize the element ready for calculation, where `problem_type` is the type -of the problem (Elasticity, Dirichlet, etc.), `element_name` is the name of a -constructed element (see Element(element_type, connectivity_vector)) and `time` -is the starting time of the initializing process. -""" -function initialize!(problem::Problem, element::Element, time::Float64) - field_name = get_unknown_field_name(problem) - field_dim = get_unknown_field_dimension(problem) - nnodes = length(element) - - # initialize primary field - if !haskey(element, field_name) - if field_dim == 1 - update!(element, field_name, time => zeros(nnodes)) - else - update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes]) - end - end - - # if a boundary problem, initialize also a field for the main problem - is_boundary_problem(problem) || return - field_name = get_parent_field_name(problem) - if !haskey(element, field_name) - if field_dim == 1 - update!(element, field_name, time => zeros(nnodes)) - else - update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes]) - end - end -end - -function initialize!(problem::Problem, time::Float64=0.0) - for element in get_elements(problem) - initialize!(problem, element, time) - end -end - -""" - update!(problem, assembly, u, la) - -Update the problem solution vector for assembly. -""" -function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector) - - # resize & fill with zeros vectors if length mismatch with current solution - - if length(u) != length(assembly.u) - info("resizing solution vector u") - resize!(assembly.u, length(u)) - fill!(assembly.u, 0.0) - end - - if length(la) != length(assembly.la) - info("resizing lagrange multiplier vector la") - resize!(assembly.la, length(la)) - fill!(assembly.la, 0.0) - end - - # copy current solutions to previous ones and add/replace new solution - # TODO: here we have couple of options and they need to be clarified - # for total formulation we are solving total quantity Ku = f while in - # incremental formulation we solve KΔu = f and u = u + Δu - assembly.u_prev = copy(assembly.u) - assembly.la_prev = copy(assembly.la) - - if get_formulation_type(problem) == :total - assembly.u = u - assembly.la = la - elseif get_formulation_type(problem) == :incremental - assembly.u += u - assembly.la = la - elseif get_formulation_type(problem) == :forwarddiff - assembly.u += u - assembly.la += la - else - info("$(problem.name): unknown formulation type, don't know what to do with results") - error("serious failure with problem formulation: $(get_formulation_type(problem))") - end - - # calculate change of norm - assembly.u_norm_change = norm(assembly.u - assembly.u_prev) - assembly.la_norm_change = norm(assembly.la - assembly.la_prev) - return assembly.u, assembly.la -end - -""" - get_global_solution(problem, assembly) - -Return a global solution (u, la) for a problem. - -Notes ------ -If the length of solution vector != number of nodes, i.e. the field dimension is -something else than 1, reshape vectors so that their length matches to the -number of nodes. This helps to get nodal results easily. -""" -function get_global_solution(problem::Problem, assembly::Assembly) - u = assembly.u - la = assembly.la - field_dim = get_unknown_field_dimension(problem) - if field_dim == 1 - return u, la - else - nnodes = round(Int, length(u)/field_dim) - u = reshape(u, field_dim, nnodes) - u = Vector{Float64}[u[:,i] for i in 1:nnodes] - la = reshape(la, field_dim, nnodes) - la = Vector{Float64}[la[:,i] for i in 1:nnodes] - return u, la - end -end - -""" - update!(problem, assembly, elements, time) - -Update a solution from the assebly to elements. -""" -function update!{P<:FieldProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) - u, la = get_global_solution(problem, assembly) - field_name = get_unknown_field_name(problem) - # update solution u for elements - for element in elements - connectivity = get_connectivity(element) - update!(element, field_name, time => u[connectivity]) - end -end - -function update!{P<:BoundaryProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) - u, la = get_global_solution(problem, assembly) - parent_field_name = get_parent_field_name(problem) # displacement - field_name = get_unknown_field_name(problem) # lambda - # update solution and lagrange multipliers for boundary elements - for element in elements - connectivity = get_connectivity(element) - update!(element, parent_field_name, time => u[connectivity]) - update!(element, field_name, time => la[connectivity]) - end -end - -""" - add_elements!(problem::Problem, elements) - -Add new elements into the problem. -""" -function add_elements!(problem::Problem, elements) - for element in elements - push!(problem.elements, element) - end -end - -function get_elements(problem::Problem) - return problem.elements -end - -function get_assembly(problem::Problem) - return problem.assembly -end - -function length(problem::Problem) - return length(problem.elements) -end - -function update!(problem::Problem, field_name::AbstractString, data) - #if haskey(problem.fields, field_name) - # update!(problem.fields[field_name], field_name::AbstractString, data) - #else - # problem.fields[field_name] = Field(data) - #end - update!(problem.elements, field_name::AbstractString, data) -end - -function haskey(problem::Problem, field_name::AbstractString) - return haskey(problem.fields, field_name) -end - -function getindex(problem::Problem, field_name::AbstractString) - return problem.fields[field_name] -end - -""" Return field calculated to nodal points for elements in problem p. """ -function (problem::Problem)(field_name::AbstractString, time::AbstractFloat=0.0) - #if haskey(problem, field_name) - # return problem[field_name](time) - #end - f = nothing - for element in get_elements(problem) - haskey(element, field_name) || continue - for (c, v) in zip(get_connectivity(element), element(field_name, time)) - if f == nothing - f = Dict(c => v) - continue - end - if haskey(f, c) - if !isapprox(f[c], v) - info("several values for single node when returning field $field_name") - info("already have: $(f[c]), and trying to set $v") - end - else - f[c] = v - end - end - end - #f == nothing && return f - #update!(problem, field_name, time => f) - return f -end - -""" Return the dimension of the unknown field of this problem. """ -function get_unknown_field_dimension(problem::Problem) - return problem.dimension -end - -""" Return the name of the unknown field of this problem. """ -function get_unknown_field_name{P}(problem::Problem{P}) - return get_unknown_field_name(P) -end - -""" Return the name of the parent field of this (boundary) problem. """ -function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P}) - return problem.parent_field_name -end - -function push!(problem::Problem, elements...) - push!(problem.elements, elements...) -end - -function push!(problem::Problem, elements::Vector) - push!(problem.elements, elements...) -end - -function push!(problem::Problem, elements_::Vector...) - for elements in elements_ - push!(problem.elements, elements...) - end -end - -function get_gdofs(element::Element, dim::Int) - conn = get_connectivity(element) - if length(conn) == 0 - error("element connectivity not defined, cannot determine global dofs for element: $element") - end - gdofs = vec([dim*(i-1)+j for j=1:dim, i in conn]) - return gdofs -end - -function empty!(problem::Problem) - empty!(problem.assembly) -end - -""" Return global degrees of freedom for element. - -Notes ------ -First look dofs from problem.dofmap, it not found, update dofmap from -element.element connectivity using formula gdofs = [dim*(nid-1)+j for j=1:dim] -1. look element dofs from problem.dofmap -2. if not found, use element.connectivity to update dofmap and 1. -""" -function get_gdofs(problem::Problem, element::Element) - if !haskey(problem.dofmap, element) - dim = get_unknown_field_dimension(problem) - problem.dofmap[element] = get_gdofs(element, dim) - end - return problem.dofmap[element] -end diff --git a/src/solvers.jl b/src/solvers.jl index f9bf1c7..c0f6172 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -52,12 +52,6 @@ function haskey(solver::Solver, field_name::String) return haskey(solver.fields, field_name) end -# one-liner helpers to identify problem types - -is_field_problem(problem) = false -is_field_problem{P<:FieldProblem}(problem::Problem{P}) = true -is_boundary_problem(problem) = false -is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) = true get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver)) get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver)) diff --git a/src/sparse.jl b/src/sparse.jl deleted file mode 100644 index 3cf5eaa..0000000 --- a/src/sparse.jl +++ /dev/null @@ -1,192 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# Sparse utils to make assembly of local and global matrices easier. -# Unoptimized but should do all necessary stuff for at start. - -type SparseMatrixCOO{T<:Real} - I :: Vector{Int} - J :: Vector{Int} - V :: Vector{T} -end - -const SparseVectorCOO = SparseMatrixCOO - -function SparseMatrixCOO() - return SparseMatrixCOO{Float64}([], [], []) -end - -function SparseVectorCOO(I::Vector, V::Vector) - return SparseVectorCOO(I, ones(I), V) -end - -function convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC) - return SparseMatrixCOO(findnz(A)...) -end - -function convert(::Type{SparseVectorCOO}, A::SparseVector) - return SparseVectorCOO(findnz(A)...) -end - -function convert(::Type{SparseMatrixCOO}, A::Matrix) - return SparseMatrixCOO(findnz(A)...) -end - -""" Convert from COO format to CSC. - -Parameters ----------- -tol - used to drop near zero values less than tol. -""" -function sparse(A::SparseMatrixCOO; tol=1.0e-12) - B = sparse(A.I, A.J, A.V) - SparseArrays.droptol!(B, tol) - return B -end - -function sparse(A::SparseMatrixCOO, n::Int, m::Int; tol=1.0e-12) - B = sparse(A.I, A.J, A.V, n, m) - SparseArrays.droptol!(B, tol) - return B -end - -function sparse(A::SparseMatrixCOO, n::Int, m::Int, f::Function; tol=1.0e-12) - B = sparse(A.I, A.J, A.V, n, m, f) - SparseArrays.droptol!(B, tol) - return B -end - -function push!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64) - push!(A.I, I) - push!(A.J, J) - push!(A.V, V) -end - -function empty!(A::SparseMatrixCOO) - empty!(A.I) - empty!(A.J) - empty!(A.V) -end - -function append!(A::SparseMatrixCOO, B::SparseMatrixCOO) - append!(A.I, B.I) - append!(A.J, B.J) - append!(A.V, B.V) -end - -function isempty(A::SparseMatrixCOO) - return isempty(A.I) && isempty(A.J) && isempty(A.V) -end - -function full(A::SparseMatrixCOO, args...) - return full(sparse(A.I, A.J, A.V, args...)) -end - -""" Add local element matrix to sparse matrix. This basically does: - ->>> A[dofs1, dofs2] = A[dofs1, dofs2] + data - -Example -------- - ->>> S = [3, 4] ->>> M = [6, 7, 8] ->>> data = Float64[5 6 7; 8 9 10] ->>> A = SparseMatrixCOO() ->>> add!(A, S, M, data) ->>> full(A) -4x8 Array{Float64,2}: - 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 5.0 6.0 7.0 - 0.0 0.0 0.0 0.0 0.0 8.0 9.0 10.0 - -""" -function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix) - n, m = size(data) - for j=1:m - for i=1:n - push!(A.I, dofs1[i]) - push!(A.J, dofs2[j]) - end - end - append!(A.V, vec(data)) -end - -""" Add sparse matrix of CSC to COO. """ -function add!(A::SparseMatrixCOO, B::SparseMatrixCSC) - I, J, V = findnz(B) - C = SparseMatrixCOO(I, J, V) - append!(A, C) -end - -""" Add new data to COO Sparse vector. """ -function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1) - if length(dofs) != length(data) - info("dofs = $dofs") - info("data = $(vec(data))") - error("when adding to sparse vector dimension mismatch!") - end - append!(A.I, dofs) - append!(A.J, dim*ones(Int, length(dofs))) - append!(A.V, vec(data)) -end - -""" Add SparseVector to SparseVectorCOO. """ -function add!(a::SparseVectorCOO, b::SparseVector) - I, V = findnz(b) - c = SparseVectorCOO(I, V) - append!(a, c) - return -end - -""" Combine (I,J,V) values if possible to reduce memory usage. """ -function optimize!(A::SparseMatrixCOO) - I, J, V = findnz(sparse(A)) - A.I = I - A.J = J - A.V = V - return -end - -""" Find all nonzero rows from sparse matrix. - -Returns -------- - -Ordered list of row indices. -""" -function get_nonzero_rows(A::SparseMatrixCSC) - return sort(unique(rowvals(A))) -end - -function get_nonzero_columns(A::SparseMatrixCSC) - return get_nonzero_rows(transpose(A)) -end - -function size(A::SparseMatrixCOO) - isempty(A) && return (0, 0) - return maximum(A.I), maximum(A.J) -end - -function size(A::SparseMatrixCOO, idx::Int) - return size(A)[idx] -end - -""" Resize sparse matrix A to (higher) dimension n x m. """ -function resize_sparse(A, n, m) - return sparse(findnz(A)..., n, m) -end - -""" Resize sparse vector b to (higher) dimension n. """ -function resize_sparsevec(b, n) - return sparsevec(findnz(b)..., n) -end - -""" Approximative comparison of two matricse A and B. """ -function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO) - A2 = sparse(A) - B2 = sparse(B, size(A2)...) - return isapprox(A2, B2) -end diff --git a/src/types.jl b/src/types.jl deleted file mode 100644 index e4dcec4..0000000 --- a/src/types.jl +++ /dev/null @@ -1,76 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -const Node = Vector{Float64} - -abstract type AbstractPoint end - -type Point{P<:AbstractPoint} - id :: Int - weight :: Float64 - coords :: Tuple{Vararg{Float64}} - fields :: Dict{AbstractString, Field} - properties :: P -end - -function setindex!{T}(point::Point, val::Pair{Float64, T}, field_name) - point.fields[field_name] = Field(val) -end - -function getindex(point::Point, field_name) - return point.fields[field_name] -end - -function getindex(point::Point, idx::Int) - return point.coords[idx] -end - -function haskey(point::Point, field_name) - return haskey(point.fields, field_name) -end - -function (point::Point)(field_name, time=0.0) - point.fields[field_name](time).data -end - -function start(point::Point) - return start(point.coords) -end - -function done(point::Point, i) - return done(point.coords, i) -end - -function next(point::Point, i) - return next(point.coords, i) -end - -function update!{T}(point::Point, field_name, val::Pair{Float64, T}) - if haskey(point, field_name) - update!(point[field_name], val) - else - point[field_name] = val - end -end - -#= TODO: in future -type Node <: AbstractPoint -end - -type MaterialPoint <: AbstractPoint -end -=# - -type IntegrationPoint <: AbstractPoint -end - -const IP = Point{IntegrationPoint} - -function IP(id, weight, coords::Tuple) - return IP(id, weight, coords, Dict(), IntegrationPoint()) -end - -function IP(id, weight, coords::Vector) - warn("Consider giving coords as Tuple.") - return IP(id, weight, [c for c in coords], Dict(), IntegrationPoint()) -end