Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-18 17:47:29 +00:00
Code Issues Packages Projects Releases Wiki Activity

code cleanup. inverse isoparametric mapping. function field wrapper.

Browse Source
This commit is contained in:
Jukka Aho
2016-07-07 18:00:36 +03:00
parent 6417b76f8c
commit 6c9b76ae54
1 changed files with 33 additions and 358 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elements.jl
+33 -358
View File
@@ -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
=#