From 6c9b76ae545be22b47d238f40028d194776e65e9 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 7 Jul 2016 18:00:36 +0300 Subject: [PATCH] code cleanup. inverse isoparametric mapping. function field wrapper. --- src/elements.jl | 391 ++++-------------------------------------------- 1 file changed, 33 insertions(+), 358 deletions(-) diff --git a/src/elements.jl b/src/elements.jl index c5ff893..154ae72 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -25,6 +25,19 @@ function setindex!(element::Element, data::Field, field_name::String) element.fields[field_name] = data end +function setindex!(element::Element, data::Function, field_name::String) + 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::String) element.fields[field_name] = Field(data) end @@ -229,16 +242,6 @@ function update!(elements::Vector, field_name::String, data) end end -#= -dbasis_cache = ForwardDiff.jacobian -""" Evaluate partial derivatives of basis functions using ForwardDiff. """ -function get_dbasis(element::Element, ip, time) - xi = isa(ip, IP) ? ip.coords : ip - basis(xi) = vec(get_basis(element, xi, time)) - return ForwardDiff.jacobian(basis, xi)' -end -=# - """ Check existence of field. """ function haskey(element::Element, field_name) haskey(element.fields, field_name) @@ -272,11 +275,11 @@ function get_gdofs(element::Element) end """ Return dual basis transformation matrix Ae. """ -function get_dualbasis(element::Element, time) +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) + for ip in get_integration_points(element, order) detJ = element(ip, time, Val{:detJ}) w = ip.weight*detJ N = element(ip, time) @@ -286,354 +289,26 @@ function get_dualbasis(element::Element, time) return De, Me, De*inv(Me) end -#= - -type Element{E} - connectivity :: Vector{Int} - fields :: Dict{String, Field} - # matrices to construct dual basis - D :: Matrix{Float64} - M :: Matrix{Float64} - A :: Matrix{Float64} -end - -function Base.size{E}(::Element{E}) - return size(E) -end - -function Base.size{E}(::Element{E}, i::Int64) - return size(E)[i] -end - -function convert{E}(::Type{Element{E}}, connectivity::Vector{Int}) - return Element{E}(connectivity, Dict(), Matrix(), Matrix(), Matrix()) -end - -function get_integration_points{E}(element::Element{E}, args...) - return get_integration_points(E, args...) -end - -function update_gauss_fields!(element::Element, data::Vector{IntegrationPoint}, time::Real) - if haskey(element, "integration points") - # push or update - if !isapprox(last(element["integration points"]).time, time) - push!(element["integration points"], time => data) - else - last(element["integration points"]).data = data - end - else - # create - element["integration points"] = Field(time => data) +""" 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 -""" Get FieldSet from element. """ -function Base.getindex(element::Element, field_name) - return element.fields[field_name] +""" 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 -function Base.length{E}(element::Element{E}) - size(E)[2] -end - -"""Add new Field to element. - -Examples --------- ->>> element["temperature"] = [1, 2, 3, 4] ->>> element["temperature"] = (0.0, [0, 0, 0, 0]), (1.0, [1, 2, 3, 4]) ->>> element["temperature"] = (0.0 => [0, 0, 0, 0], 1.0 => [1, 2, 3, 4]) -""" -function Base.setindex!(element::Element, data, name::String) - element.fields[name] = Field(data) -end -function Base.setindex!(element::Element, field::Field, name::String) - element.fields[name] = field -end -function Base.setindex!(element::Element, data::Tuple, name::String) - element.fields[name] = Field(data...) -end - - -typealias VecOrIP Union{Vector, IntegrationPoint} - -function call(element::Element, field_name::String, time::Real, variation=nothing) - return isa(variation, Void) ? element[field_name](time) : variation -end - -function call(element::Element, field_name::String, xi::VecOrIP, time::Number, variation=nothing) - field = element(field_name, time, variation) -# field = isa(variation, Void) ? element[field_name](time) : variation - basis = get_basis(element) - return basis(field, xi) -end - -function call(element::Element, field_name::String, xi::VecOrIP, time::Number, ::Type{Val{:grad}}, variation=nothing) -# field = isa(variation, Void) ? element[field_name](time) : variation - field = element(field_name, time, variation) - basis = get_basis(element) - geom = element["geometry"](time) - return basis(geom, field, xi, Val{:grad}) -end - -function call(element::Element, field_name::String, xi::VecOrIP) - field = element[field_name] - basis = get_basis(element) - return basis(element[field_name], xi) -end - -function call(element::Element, field_name::String, xi::VecOrIP, ::Type{Val{:grad}}) - field = element[field_name] - geom = element["geometry"] - basis = get_basis(element) - return basis(geom, field, xi, Val{:grad}) -end - -function call(element::Element, field_name::String, time::Number) - return element[field_name](time) -end - - -function get_basis{E}(element::Element{E}, ip::IntegrationPoint) - return get_basis(E, ip.xi) -end - -function get_basis{E}(::Type{Element{E}}, xi::Vector{Float64}) - return get_basis(E, xi) -end - -function get_basis{E}(element::Element{E}, xi::Vector) - return get_basis(E, xi) -end - -function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0) - return get_basis(element, xi) -end - -""" Given a list of elementa and nodes, find a subset of elements -containing nodes. -""" -function find_elements(elements, nodes) - s = Set{Element}() - for element in elements - conn = get_connectivity(element) - for j in nodes - if j in conn - push!(s, element) - break - end - end - end - return collect(s) -end - - -function get_dbasis{E}(element::Element{E}, ip::IntegrationPoint) - return get_dbasis(E, ip.xi) -end - -function get_basis{E, T<:Real}(element::Element{E}, xi::T) - return get_basis(E, xi) -end - - -function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}}) - De, Me, Ae = get_dualbasis(element, time) - N = get_basis(element, xi) - Phi = Ae*N' - return Phi' -end - -function get_basis{E}(element::Element{E}) - basis = CVTI( - (xi::Vector) -> get_basis(E, xi), - (xi::Vector) -> get_dbasis(E, xi)) - return basis -end - -function call{E}(element::Element{E}, xi::VecOrIP, ::Type{Val{:grad}}) - basis = get_basis(element) - geom = element["geometry"] - return basis(geom, xi, Val{:grad}) -end - -function call{E}(element::Element{E}, xi::VecOrIP, time::Float64, ::Type{Val{:grad}}) - basis = get_basis(element) - return basis(element["geometry"](time), xi, Val{:grad}) -end - -function call(element::Element, field_name::String) - return element[field_name] -end - - -""" Return the jacobian of element. """ -function get_jacobian{E}(element::Element{E}, xi::Vector{Float64}, time::Real) - X = element("geometry", time) - dN = get_dbasis(E, xi) - J = sum([kron(dN[:,i], X[i]') for i=1:length(X)]) - return J -end -function get_jacobian{E}(element::Element{E}, ip::IntegrationPoint, time::Real) - return get_jacobian(element, ip.xi, time) -end - -""" Return Jacobian of element in deformed state. """ -function get_jacobian{E}(element::Element{E}, xi::Vector{Float64}, time::Real, ::Type{Val{:deformed}}) - x = element("geometry", time) - if haskey(element, "displacement") - x += element("displacement", time) - end - dN = get_dbasis(E, xi) - j = sum([kron(dN[:,i], x[i]') for i=1:length(x)]) - return j -end -function get_jacobian{E}(element::Element{E}, ip::IntegrationPoint, time::Real, ::Type{Val{:deformed}}) - return get_jacobian(element, ip.xi, time, Val{:deformed}) -end - - - -""" Calculate local normal-tangential coordinates for element. """ -function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real) - ntcoords = Matrix[] - normals = Vector{Float64}[] - refcoords = get_reference_element_coordinates(E) - x = element("geometry", time) - for xi in refcoords - dN = get_dbasis(E, xi)*x - n, m = size(dN) - @assert n != m # if n == m -> this is not manifold - if m == 1 # plane case - tangent = dN / norm(dN) - normal = [-tangent[2] tangent[1]]' - push!(normals, vec(normal)) - push!(ntcoords, [normal tangent]) - elseif m == 2 - normal = cross(dN[:,1], dN[:,2]) - normal /= norm(normal) - u1 = normal - j = indmax(abs(u1)) - v2 = zeros(3) - v2[mod(j,3)+1] = 1.0 - u2 = v2 - dot(u1, v2) / dot(v2, v2) * v2 - u3 = cross(u1, u2) - tangent1 = u2/norm(u2) - tangent2 = u3/norm(u3) - push!(ntcoords, [normal tangent1 tangent2]) - push!(normals, vec(normal)) - else - error("calculate_normal_tangential_coordinates!(): n=$n, m=$m") - end - end - element["normal-tangential coordinates"] = ntcoords - element["normals"] = normals -end - -""" Return list of nodes / connectivity points from a set of elements. -""" -function get_nodes(elements::Vector) - nodes = Set{Int64}() - for element in elements - push!(nodes, get_connectivity(element)...) - end - nodes = sort(collect(nodes)) - return nodes -end - -""" Calculate normal-tangential coordinates for a set of elements. - -Notes ------ -Average normals so that normals are unique in nodes. -""" - -function calculate_normal_tangential_coordinates!(elements::Vector, time::Real, configuration::Symbol=:deformed) - if size(elements[1], 1) == 1 - return calculate_normal_tangential_coordinates!(elements, time, Val{2}, configuration) - else - return calculate_normal_tangential_coordinates!(elements, time, Val{3}, configuration) - end -end - -""" Calculate normal-tangential coordinates for 2d case. - -Notes ------ -n = (e₃×∂X/∂ξ) / || e₃×∂X/∂ξ || and e₃ = [0 0 1] -""" -function calculate_normal_tangential_coordinates!(elements::Vector, time::Real, ::Type{Val{2}}, configuration::Symbol) - nodes = get_nodes(elements) - n = zeros(2, maximum(nodes)) - Q = [0 -1; 1 0] - for element in elements - gdofs = get_gdofs(element, 1) - for ip in get_integration_points(element, Val{3}) - if configuration == :deformed - J = get_jacobian(element, ip, time, Val{:deformed}) - else - J = get_jacobian(element, ip, time) - end - N = element(ip, time) - n[:, gdofs] += ip.weight*Q*J'*N - end - end - t = zeros(n) - for i=1:size(n,2) - n[:,i] = n[:,i] / norm(n[:,i]) - t[:,i] = [-n[2,i], n[1,i]] - end - for element in elements - node_ids = get_connectivity(element) - Q = Matrix{Float64}[ [n[:,i] t[:,i]] for i in node_ids] - element["normal-tangential coordinates"] = (time => Q) - element["normals"] = (time => Vector{Float64}[n[:,i] for i in node_ids]) - end -end - -""" Calculate normal-tangential coordinates for 3d case. -""" -function calculate_normal_tangential_coordinates!(elements::Vector, time::Real, ::Type{Val{3}}) - nodes = get_nodes(elements) - n = zeros(3, maximum(nodes)) - for element in elements - gdofs = get_gdofs(element, 1) - for ip in get_integration_points(element, Val{3}) - J = transpose(get_jacobian(element, ip, time, Val{:deformed})) - N = element(ip, time) - c = reshape(cross(J[:,1], J[:,2]), 3, 1) - n[:, gdofs] += ip.weight*c*N - end - end - t1 = zeros(n) - t2 = zeros(n) - for i=1:size(n,2) - i in nodes || continue - n[:,i] = n[:,i] / norm(n[:,i]) - u1 = n[:,i] - j = indmax(abs(n[:,i])) - v2 = zeros(3) - v2[mod(j,3)+1] = 1.0 - u2 = v2 - dot(u1, v2) / dot(v2, v2) * v2 - u3 = cross(u1, u2) - t1[:,i] = u2/norm(u2) - t2[:,i] = u3/norm(u3) - end - for element in elements - node_ids = get_connectivity(element) - Q = Matrix{Float64}[ [n[:,i] t1[:,i] t2[:,i]] for i in node_ids] - element["normal-tangential coordinates"] = (time => Q) - element["normals"] = (time => Vector{Float64}[n[:,i] for i in node_ids]) - end -end - - -""" Update values for several elements at once. """ -# FIXME: with or without {T} ? -function update!{T}(elements::Vector{Element{T}}, field_name::String, data...) - for element in elements - update!(element, field_name, data...) - end -end - -=#