mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Testing/code coverage (#83)
Change the code coverage to green. * removed duplicate code * Removed unused code * removed unmaintained code * DCTI + DVTI refactored * discrete fields refactored and tested * fields are now tested quite well. * Removed obsolete code not used anywhere * Element descriptions to common dictionary * size in global const dictionary also * Added coverage to sparse tools and removed couple unused functions * get nonzero rows from SparseMatrixCSC * bugfix: extending element basis now working and tested * Removed two unused functions from elements.jl * removed useless function * Useless conversion * remove elasticity assembly using ForwardDiff because it's not used anywhere' * Added basic testing for NURBS. Fixed bug in NSolid interpolation. * removed unused functions * Removed some debug stuff * renamed file * removed field assembly posthook, i think not good idea at all * test for nnz(K) == 0 and automatic determination of dofs * Testing that solver is throwing error if having problems with boundary assembly * Removed some unused options. Refactoring. * Moved solver non-related code to elements.jl * Removed custom exception (no need) * unneeded postprocess code * More tests for NURBS elements. * Removed unfinished .mail parser * proper use of Logging package * also read results * renamed test file * create_surface_elements accepts surface name in String now * bugfix: remove zero rows from constraint matrix after manually removing dofs from some boundary assemblies. * New test, displacement 3d patch test * skip displacement field in surface element splitting if not defined * test element splitting and linear surface elements, fails. * Bugfix: Xdmf, not XDMF * removed nonworking tests, requires bugfix * abaqus_read_results is not working -> bug
This commit is contained in:
Jukka Aho
committed by
Tero Frondelius
Tero Frondelius
parent
ddabc9d82b
commit
c307c1482c
+4
-2
@@ -44,7 +44,9 @@ export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
|
||||
get_dbasis, inside, get_local_coordinates
|
||||
|
||||
include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
|
||||
export get_reference_coordinates, get_interpolation_polynomial
|
||||
export get_reference_coordinates,
|
||||
get_interpolation_polynomial,
|
||||
description
|
||||
export Poi1,
|
||||
Seg2, Seg3,
|
||||
Tri3, Tri6, Tri7,
|
||||
@@ -62,7 +64,7 @@ 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
|
||||
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,
|
||||
|
||||
+5
-1
@@ -703,7 +703,7 @@ end
|
||||
function abaqus_open_results(name)
|
||||
path = abaqus_input_file_path(name)
|
||||
result_file = "$path/$name.xmf"
|
||||
return XDMF(result_file)
|
||||
return Xdmf(result_file)
|
||||
end
|
||||
|
||||
### JuliaFEM-ABAQUS interface entry point
|
||||
@@ -759,3 +759,7 @@ function create_surface_elements(mesh::Mesh, surface_name::Symbol)
|
||||
return elements
|
||||
end
|
||||
|
||||
function create_surface_elements(mesh::Mesh, surface_name::String)
|
||||
return create_surface_elements(mesh, Symbol(surface_name))
|
||||
end
|
||||
|
||||
|
||||
@@ -1,46 +1,6 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
function optimize!(assembly::Assembly)
|
||||
optimize!(assembly.K)
|
||||
optimize!(assembly.Kg)
|
||||
optimize!(assembly.f)
|
||||
optimize!(assembly.fg)
|
||||
optimize!(assembly.C1)
|
||||
optimize!(assembly.C2)
|
||||
optimize!(assembly.D)
|
||||
optimize!(assembly.g)
|
||||
optimize!(assembly.c)
|
||||
end
|
||||
|
||||
function append!(assembly::Assembly, sub_assembly::Assembly)
|
||||
append!(assembly.M, sub_assembly.M)
|
||||
append!(assembly.K, sub_assembly.K)
|
||||
append!(assembly.Kg, sub_assembly.Kg)
|
||||
append!(assembly.f, sub_assembly.f)
|
||||
append!(assembly.fg, sub_assembly.fg)
|
||||
append!(assembly.C1, sub_assembly.C1)
|
||||
append!(assembly.C2, sub_assembly.C2)
|
||||
append!(assembly.D, sub_assembly.D)
|
||||
append!(assembly.g, sub_assembly.g)
|
||||
append!(assembly.c, sub_assembly.c)
|
||||
end
|
||||
|
||||
""" Calculate norm of assembly, i.e., norm of each block of matrix. """
|
||||
function norm(assembly::Assembly, p=2)
|
||||
N1 = norm(assembly.M, p)
|
||||
N2 = norm(assembly.K, p)
|
||||
N3 = norm(assembly.Kg, p)
|
||||
N4 = norm(assembly.f, p)
|
||||
N5 = norm(assembly.fg, p)
|
||||
N6 = norm(assembly.C1, p)
|
||||
N7 = norm(assembly.C2, p)
|
||||
N8 = norm(assembly.D, p)
|
||||
N9 = norm(assembly.g, p)
|
||||
N10 = norm(assembly.c, p)
|
||||
return [N1, N2, N3, N4, N5, N6, N7, N8, N9, N10]
|
||||
end
|
||||
|
||||
function isapprox(a1::Assembly, a2::Assembly)
|
||||
T = isapprox(a1.K, a2.K)
|
||||
T &= isapprox(a1.C1, a2.C1)
|
||||
|
||||
+17
-40
@@ -29,6 +29,22 @@ 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 is(E, element_type)
|
||||
end
|
||||
|
||||
function filter_by_element_type(element_type, elements)
|
||||
return filter(element -> is_element_type(element, element_type), elements)
|
||||
end
|
||||
|
||||
function setindex!(element::Element, data::Function, field_name)
|
||||
if method_exists(data, Tuple{Element, Vector, Float64})
|
||||
# create enclosure to pass element as argument
|
||||
@@ -100,7 +116,7 @@ julia> el([0.0, 0.0], 0.0, 2)
|
||||
"""
|
||||
function (element::Element)(ip, time::Float64, dim::Int)
|
||||
dim == 1 && return get_basis(element, ip, time)
|
||||
Ni = 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
|
||||
@@ -144,14 +160,6 @@ function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:G
|
||||
return element(ip, time, Val{:Grad})*element[field_name](time)
|
||||
end
|
||||
|
||||
function (element::Element)(field::Field, time::Float64)
|
||||
return field(time)
|
||||
end
|
||||
|
||||
function (element::Element)(field::DCTI, time::Float64)
|
||||
return field.data
|
||||
end
|
||||
|
||||
function (element::Element)(field_name::String, ip, time::Float64)
|
||||
field = element[field_name]
|
||||
return element(field, ip, time)
|
||||
@@ -220,35 +228,8 @@ function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K,
|
||||
time, field_data = data
|
||||
element_data = V[field_data[i] for i in get_connectivity(element)]
|
||||
update!(element, field_name, time => element_data)
|
||||
#if haskey(element, field_name)
|
||||
# update!(element[field_name], data)
|
||||
#else
|
||||
# element[field_name] = Field(data)
|
||||
#end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
|
||||
for data in datas
|
||||
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
|
||||
end
|
||||
|
||||
#=
|
||||
function update!(element::Element, field_name, data::Pair...)
|
||||
for data in datas
|
||||
update!(element, field_name, data)
|
||||
end
|
||||
end
|
||||
=#
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Any}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
@@ -353,10 +334,6 @@ function get_integration_points(element::Element, change_order::Int)
|
||||
return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
|
||||
end
|
||||
|
||||
function get_gdofs(element::Element)
|
||||
return get_gdofs(element, 1)
|
||||
end
|
||||
|
||||
""" Return dual basis transformation matrix Ae. """
|
||||
function get_dualbasis(element::Element, time::Float64, order=1)
|
||||
nnodes = length(element)
|
||||
|
||||
+55
-192
@@ -1,23 +1,66 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
global const ELEMENT_DESCRIPTIONS = Dict(
|
||||
"Poi1" => "1 node discrete point element",
|
||||
"Seg2" => "2 node linear segment/line element",
|
||||
"Seg3" => "3 node quadratic segment/line element",
|
||||
"Tri3" => "3 node linear triangle element",
|
||||
"Tri6" => "6 node quadratic triangle element",
|
||||
"Tri7" => "7 node quadratic triangle element (has middle node)",
|
||||
"Quad4" => "4 node linear quadrangle element",
|
||||
"Quad8" => "8 node quadratic quadrangle element (Serendip)",
|
||||
"Quad9" => "9 node quadratic quadrangle element",
|
||||
"Tet4" => "4 node linear tetrahedral element",
|
||||
"Tet10" => "10 node quadratic tetrahedral element",
|
||||
"Wedge6" => "6 node linear prismatic element (wedge)",
|
||||
"Wedge15" => "15 node quadratic prismatic element (wedge)",
|
||||
"Hex8" => "8 node linear hexahedral element",
|
||||
"Hex20" => "20 node biquadratic hexahedral element",
|
||||
"Hex27" => "27 node quadratic hexahedral element")
|
||||
|
||||
global const ELEMENT_SIZES = Dict(
|
||||
"Poi1" => (0, 1),
|
||||
"Seg2" => (1, 2),
|
||||
"Seg3" => (1, 3),
|
||||
"Tri3" => (2, 3),
|
||||
"Tri6" => (2, 6),
|
||||
"Tri7" => (2, 7),
|
||||
"Quad4" => (2, 4),
|
||||
"Quad8" => (2, 8),
|
||||
"Quad9" => (2, 9),
|
||||
"Tet4" => (3, 4),
|
||||
"Tet10" => (3, 10),
|
||||
"Wedge6" => (3, 6),
|
||||
"Wedge15" => (3, 15),
|
||||
"Hex8" => (3, 8),
|
||||
"Hex20" => (3, 20),
|
||||
"Hex27" => (3, 27))
|
||||
|
||||
""" Return description line of element. """
|
||||
function description{T}(element::Element{T})
|
||||
element_type = last(split("$T", '.'))
|
||||
return get(ELEMENT_DESCRIPTIONS, element_type, "Unknown element description")
|
||||
end
|
||||
|
||||
""" Return size of element, i.e. tuple (n, m) where n is dimension of element
|
||||
(0, 1, 2, 3) and m is number of nodes. """
|
||||
function size{T}(element::Element{T})
|
||||
element_type = last(split("$T", '.'))
|
||||
return ELEMENT_SIZES[element_type]
|
||||
end
|
||||
|
||||
""" Return length of element, i.e. number of nodes. """
|
||||
function length{T}(element::Element{T})
|
||||
return size(element)[end]
|
||||
end
|
||||
|
||||
### 0d element
|
||||
|
||||
|
||||
type Poi1 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Poi1})
|
||||
"1 node point"
|
||||
end
|
||||
|
||||
function size(element::Element{Poi1})
|
||||
return (0, 1)
|
||||
end
|
||||
|
||||
function length(element::Element{Poi1})
|
||||
return 1
|
||||
end
|
||||
|
||||
function get_basis(element::Element{Poi1}, ip, time)
|
||||
return [1]
|
||||
end
|
||||
@@ -47,18 +90,6 @@ end
|
||||
type Seg2 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Seg2})
|
||||
"2 node segment"
|
||||
end
|
||||
|
||||
function size(element::Element{Seg2})
|
||||
return (1, 2)
|
||||
end
|
||||
|
||||
function length(element::Element{Seg2})
|
||||
return 2
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Seg2})
|
||||
Vector{Float64}[
|
||||
[-1.0], # N1
|
||||
@@ -78,18 +109,6 @@ end
|
||||
type Seg3 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Seg3})
|
||||
"3 node segment"
|
||||
end
|
||||
|
||||
function size(element::Element{Seg3})
|
||||
return (1, 3)
|
||||
end
|
||||
|
||||
function length(element::Element{Seg3})
|
||||
return 3
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Seg3})
|
||||
Vector{Float64}[
|
||||
[-1.0], # N1
|
||||
@@ -110,18 +129,6 @@ end
|
||||
type Tri3 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Tri3})
|
||||
"3 node triangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Tri3})
|
||||
return (2, 3)
|
||||
end
|
||||
|
||||
function length(element::Element{Tri3})
|
||||
return 3
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Tri3})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0], # N1
|
||||
@@ -147,18 +154,6 @@ end
|
||||
type Tri6 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Tri6})
|
||||
"6 node triangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Tri6})
|
||||
return (2, 6)
|
||||
end
|
||||
|
||||
function length(element::Element{Tri6})
|
||||
return 6
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Tri6})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0], # N1
|
||||
@@ -187,18 +182,6 @@ end
|
||||
type Tri7 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Tri7})
|
||||
"7 node triangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Tri7})
|
||||
return (2, 7)
|
||||
end
|
||||
|
||||
function length(element::Element{Tri7})
|
||||
return 7
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Tri7})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0], # N1
|
||||
@@ -228,18 +211,6 @@ end
|
||||
type Quad4 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Quad4})
|
||||
"4 node quadrangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Quad4})
|
||||
return (2, 4)
|
||||
end
|
||||
|
||||
function length(element::Element{Quad4})
|
||||
return 4
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Quad4})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0], # N1
|
||||
@@ -266,18 +237,6 @@ end
|
||||
type Quad8 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Quad8})
|
||||
"8 node Serendip quadrangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Quad8})
|
||||
return (2, 8)
|
||||
end
|
||||
|
||||
function length(element::Element{Quad8})
|
||||
return 8
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Quad8})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0], # N1
|
||||
@@ -308,18 +267,6 @@ end
|
||||
type Quad9 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Quad9})
|
||||
"9 node quadrangle"
|
||||
end
|
||||
|
||||
function size(element::Element{Quad9})
|
||||
return (2, 9)
|
||||
end
|
||||
|
||||
function length(element::Element{Quad9})
|
||||
return 9
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Quad9})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0], # N1
|
||||
@@ -351,18 +298,6 @@ end
|
||||
type Tet4 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Tet4})
|
||||
"4 node tetrahedral element"
|
||||
end
|
||||
|
||||
function size(element::Element{Tet4})
|
||||
return (3, 4)
|
||||
end
|
||||
|
||||
function length(element::Element{Tet4})
|
||||
return 4
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Tet4})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0, 0.0], # N1
|
||||
@@ -390,18 +325,6 @@ end
|
||||
type Tet10 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Tet10})
|
||||
"10 node tetrahedral element"
|
||||
end
|
||||
|
||||
function size(element::Element{Tet10})
|
||||
return (3, 10)
|
||||
end
|
||||
|
||||
function length(element::Element{Tet10})
|
||||
return 10
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Tet10})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0, 0.0], # N1
|
||||
@@ -435,18 +358,6 @@ end
|
||||
type Wedge6 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Wedge6})
|
||||
"6 node prismatic element (wedge)"
|
||||
end
|
||||
|
||||
function size(element::Element{Wedge6})
|
||||
return (3, 6)
|
||||
end
|
||||
|
||||
function length(element::Element{Wedge6})
|
||||
return 6
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Wedge6})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0, -1.0], # N1
|
||||
@@ -476,18 +387,6 @@ end
|
||||
type Wedge15 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Wedge15})
|
||||
"15 node prismatic element (wedge)"
|
||||
end
|
||||
|
||||
function size(element::Element{Wedge15})
|
||||
return (3, 15)
|
||||
end
|
||||
|
||||
function length(element::Element{Wedge15})
|
||||
return 15
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Wedge15})
|
||||
Vector{Float64}[
|
||||
[0.0, 0.0, -1.0], # N1
|
||||
@@ -526,18 +425,6 @@ end
|
||||
type Hex8 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Hex8})
|
||||
"8 node hexahedral element"
|
||||
end
|
||||
|
||||
function size(element::Element{Hex8})
|
||||
return (3, 8)
|
||||
end
|
||||
|
||||
function length(element::Element{Hex8})
|
||||
return 8
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Hex8})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0, -1.0], # N1
|
||||
@@ -569,18 +456,6 @@ end
|
||||
type Hex20 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Hex20})
|
||||
"20 node hexahedral element"
|
||||
end
|
||||
|
||||
function size(element::Element{Hex20})
|
||||
return (3, 20)
|
||||
end
|
||||
|
||||
function length(element::Element{Hex20})
|
||||
return 20
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Hex20})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0, -1.0], # N1
|
||||
@@ -624,18 +499,6 @@ end
|
||||
type Hex27 <: AbstractElement
|
||||
end
|
||||
|
||||
function description(::Type{Hex27})
|
||||
"27 node hexahedral element"
|
||||
end
|
||||
|
||||
function size(element::Element{Hex27})
|
||||
return (3, 27)
|
||||
end
|
||||
|
||||
function length(element::Element{Hex27})
|
||||
return 27
|
||||
end
|
||||
|
||||
function get_reference_coordinates(::Type{Hex27})
|
||||
Vector{Float64}[
|
||||
[-1.0, -1.0, -1.0], # N1
|
||||
|
||||
@@ -92,12 +92,12 @@ function get_basis(element::Element{NSolid}, xi::Vector, time)
|
||||
tu = element.properties.knots_u
|
||||
tv = element.properties.knots_v
|
||||
tw = element.properties.knots_w
|
||||
w = element.properties.weights
|
||||
nu = length(tu)
|
||||
nv = length(tv)
|
||||
nw = length(tw)
|
||||
weights = element.properties.weights
|
||||
nu = length(tu)-pu-1
|
||||
nv = length(tv)-pv-1
|
||||
nw = length(tw)-pw-1
|
||||
u, v, w = xi
|
||||
N = [w[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw]
|
||||
N = vec([weights[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw])'
|
||||
return N / sum(N)
|
||||
end
|
||||
|
||||
|
||||
+284
-307
@@ -10,48 +10,11 @@ abstract Variable <: AbstractField
|
||||
abstract TimeVariant <: AbstractField
|
||||
abstract TimeInvariant <: AbstractField
|
||||
|
||||
|
||||
type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
|
||||
data
|
||||
end
|
||||
|
||||
typealias FieldSet Dict{AbstractString, Field}
|
||||
|
||||
### Basic data structure for discrete field
|
||||
|
||||
type Increment{T}
|
||||
time :: Float64
|
||||
data :: T
|
||||
end
|
||||
|
||||
function convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
|
||||
return Increment{T}(data[1], data[2])
|
||||
end
|
||||
|
||||
function convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}})
|
||||
time = data[1]
|
||||
content = data[2]
|
||||
return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)])
|
||||
end
|
||||
|
||||
function getindex{T}(increment::Increment{Vector{T}}, i::Int64)
|
||||
return increment.data[i]
|
||||
end
|
||||
|
||||
### Basic data structure for continuous field
|
||||
|
||||
type Basis
|
||||
basis :: Function
|
||||
dbasis :: Function
|
||||
end
|
||||
|
||||
function (basis::Basis)(xi::Vector)
|
||||
basis.basis(xi)
|
||||
end
|
||||
|
||||
function (basis::Basis)(xi::Vector, ::Type{Val{:grad}})
|
||||
basis.dbasis(xi)
|
||||
end
|
||||
typealias FieldSet Dict{String, Field}
|
||||
|
||||
### Different field combinations and other typealiases
|
||||
|
||||
@@ -64,156 +27,158 @@ typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to inter
|
||||
typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time
|
||||
typealias CVTV Field{Continuous, Variable, TimeVariant}
|
||||
|
||||
typealias ScalarIncrement{T} Increment{T}
|
||||
typealias VectorIncrement{T} Increment{Vector{T}}
|
||||
typealias TensorIncrement{T} Increment{Matrix{T}}
|
||||
|
||||
typealias DiscreteField Union{DCTI, DVTI, DCTV, DVTV}
|
||||
typealias ContinuousField Union{CCTI, CVTI, CCTV, CVTV}
|
||||
typealias ConstantField Union{DCTI, DCTV, CCTI, CCTV}
|
||||
typealias VariableField Union{DVTI, DVTV, CVTI, CVTV}
|
||||
typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI}
|
||||
typealias TimeVariantField Union{DCTV, DVTV, CCTV, CVTV}
|
||||
# Discrete fields
|
||||
|
||||
|
||||
### Convenient functions to create 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.
|
||||
|
||||
#function Base.convert(::Type{Field}, data)
|
||||
# return Field(data)
|
||||
#end
|
||||
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 Field()
|
||||
return DCTI()
|
||||
end
|
||||
|
||||
function Field(data)
|
||||
return DCTI(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 Base.:*(c::Number, f::DCTI)
|
||||
return c*f.data
|
||||
end
|
||||
|
||||
""" Kind of spatial interpolation of DCTI. """
|
||||
function Base.:*(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
|
||||
|
||||
function Field{T}(data::Pair{Float64, T}...)
|
||||
return DCTV([Increment{T}(d[1], d[2]) for d in data])
|
||||
end
|
||||
#=
|
||||
function Field{T}(data::Pair{Float64, Vector{T}}...)
|
||||
return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
|
||||
end
|
||||
""" For dictionary data, DVTI is automatically created.
|
||||
|
||||
function Field{T}(data::Pair{Float64, Dict{Int64, T}}...)
|
||||
return DVTV([Increment{Dict{Int64, 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
|
||||
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 convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
|
||||
return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
|
||||
function ==(x::DVTI, y::DVTI)
|
||||
return ==(x.data, y.data)
|
||||
end
|
||||
|
||||
""" Create new discrete, constant, time variant field.
|
||||
function isapprox(x::DVTI, y)
|
||||
return isapprox(x.data, y)
|
||||
end
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0
|
||||
julia> f = DCTV(t0 => y0, t1 => y1)
|
||||
""" Default slicing of field.
|
||||
|
||||
julia> f = DVTI([1.0, 2.0])
|
||||
julia> f[1]
|
||||
1.0
|
||||
|
||||
"""
|
||||
#function convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...)
|
||||
# return DCTV([Increment(d[1],d[2]) for d in data])
|
||||
#end
|
||||
function DCTV(data::Pair...)
|
||||
return DCTV([Increment(d[1],d[2]) for d in data])
|
||||
end
|
||||
|
||||
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, Number})
|
||||
return CVTV(func)
|
||||
else
|
||||
error("no proper definition found for function: check methods.")
|
||||
end
|
||||
end
|
||||
|
||||
function CVTI(basis::Function, dbasis::Function)
|
||||
return CVTI(Basis(basis, dbasis))
|
||||
end
|
||||
|
||||
function Field(basis::Function, dbasis::Function)
|
||||
return CVTI(basis, dbasis)
|
||||
end
|
||||
|
||||
### Accessing and manipulating discrete fields
|
||||
|
||||
function getindex(field::DVTV, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
function push!(field::DCTV, data::Pair)
|
||||
push!(field.data, data)
|
||||
end
|
||||
|
||||
function push!(field::DVTV, data::Pair)
|
||||
push!(field.data, data)
|
||||
end
|
||||
|
||||
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 getindex(field::DCTV, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
function getindex(field::Field, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
function length(field::DVTI)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function length(field::DCTI)
|
||||
function start(field::DVTI)
|
||||
return 1
|
||||
end
|
||||
|
||||
function length(field::DVTV)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function length(field::DCTV)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function first(field::Union{DCTV, DVTV})
|
||||
return field[1]
|
||||
end
|
||||
|
||||
function isapprox(f1::DCTI, f2::DCTI)
|
||||
isapprox(f1.data, f2.data)
|
||||
end
|
||||
|
||||
for op = (:+, :*, :/, :-)
|
||||
@eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data)
|
||||
@eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data)
|
||||
@eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data)
|
||||
@eval ($op)(field::DCTI, k::Number) = ($op)(field.data, k)
|
||||
@eval ($op)(k::Number, field::DCTI) = ($op)(field.data, k)
|
||||
end
|
||||
|
||||
function Base.:+(f1::DVTI, f2::DVTI)
|
||||
return DVTI(f1.data + f2.data)
|
||||
end
|
||||
@@ -222,49 +187,55 @@ function Base.:-(f1::DVTI, f2::DVTI)
|
||||
return DVTI(f1.data - f2.data)
|
||||
end
|
||||
|
||||
function Base.:*{T<:Real}(c::T, field::DVTI)
|
||||
return DVTI(c*field.data)
|
||||
function update!(field::DVTI, data::Union{Vector, Dict})
|
||||
field.data = data
|
||||
end
|
||||
|
||||
function Base.:*(N::Matrix, f::DCTI)
|
||||
return f.data*N'
|
||||
""" Take scalar product of DVTI and constant T. """
|
||||
function Base.:*(T::Number, field::DVTI)
|
||||
return DVTI(T*field.data)
|
||||
end
|
||||
|
||||
|
||||
|
||||
# Multiply DVTI field with another vector T. Vector length
|
||||
# must match to the field length and this can be used mainly
|
||||
# for interpolation purposes, i.e., u = ∑ Nᵢuᵢ
|
||||
""" 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 Base.:*(T::Vector, 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 Base.:*(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
|
||||
|
||||
function vec(field::DCTV)
|
||||
error("trying to vectorize $field does not make sense")
|
||||
""" Interpolate time-invariant field in time direction. """
|
||||
function (field::DVTI)(time::Float64)
|
||||
return field
|
||||
end
|
||||
|
||||
function endof(field::Field)
|
||||
return endof(field.data)
|
||||
end
|
||||
""" Create a similar DVTI field from vector data.
|
||||
|
||||
#function Base.similar{T}(field::DVTI, data::Vector{T})
|
||||
# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
|
||||
#end
|
||||
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{T}(field::DVTI, data::Vector{T})
|
||||
n = length(field.data)
|
||||
data = reshape(data, round(Int, length(data)/n), n)
|
||||
newdata = Vector[data[:,i] for i=1:n]
|
||||
return typeof(field)(newdata)
|
||||
end
|
||||
|
||||
function start(::DVTI)
|
||||
return 1
|
||||
"""
|
||||
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)
|
||||
@@ -275,6 +246,114 @@ 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
|
||||
@@ -300,146 +379,44 @@ function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T})
|
||||
end
|
||||
end
|
||||
|
||||
function update!{T}(field::Union{DCTI, DVTI}, val::T)
|
||||
field.data = val
|
||||
### 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::CVTI)(xi::Vector)
|
||||
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::CVTV)(xi, time::Float64)
|
||||
return field.data(xi, time)
|
||||
end
|
||||
|
||||
function (field::CVTI)(xi::Vector, ::Type{Val{:Grad}})
|
||||
return field.data(xi, Val{:Grad})
|
||||
end
|
||||
|
||||
function (field::CCTV)(time::Float64)
|
||||
function (field::CCTV)(xi::Vector, time::Number)
|
||||
return field.data(time)
|
||||
end
|
||||
|
||||
function convert(::Type{Basis}, field::CVTI)
|
||||
return field.data
|
||||
function (field::CVTV)(xi::Vector, time::Number)
|
||||
return field.data(xi, time)
|
||||
end
|
||||
|
||||
### Interpolation
|
||||
|
||||
""" Interpolate time-invariant field in time direction. """
|
||||
function (field::DVTI)(time::Float64)
|
||||
return field
|
||||
end
|
||||
function (field::DCTI)(time::Float64)
|
||||
return field.data
|
||||
end
|
||||
function (field::CVTI)(time::Float64)
|
||||
return field.data()
|
||||
end
|
||||
function (field::CCTI)(time::Float64)
|
||||
return field.data()
|
||||
end
|
||||
|
||||
""" Interpolate constant time-variant field in time direction. """
|
||||
function (field::DCTV)(time::Real)
|
||||
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
|
||||
error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time")
|
||||
end
|
||||
|
||||
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
|
||||
error("interpolate DVTV: unknown failure when interpolating $(field.data) for time $time")
|
||||
end
|
||||
|
||||
""" Interpolate constant field in spatial dimension. """
|
||||
function (basis::CVTI)(field::DCTI, xi::Vector)
|
||||
return field.data
|
||||
end
|
||||
|
||||
""" Interpolate variable field in spatial dimension. """
|
||||
function (basis::CVTI)(values::DVTI, xi::Vector)
|
||||
N = basis(xi)
|
||||
return sum([N[i]*values[i] for i=1:length(N)])
|
||||
end
|
||||
|
||||
function (basis::CVTI)(geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
|
||||
dbasis = basis(xi, Val{:grad})
|
||||
# J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
|
||||
J = sum([kron(dbasis[:,i], geometry[i]') for i=1:length(geometry)])
|
||||
invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
|
||||
grad = invJ * dbasis
|
||||
return grad
|
||||
end
|
||||
|
||||
function (basis::CVTI)(geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
|
||||
grad = basis(geometry, xi, Val{:grad})
|
||||
# gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
|
||||
gradf = sum([kron(grad[:,i], values[i]') for i=1:length(values)])'
|
||||
return length(gradf) == 1 ? gradf[1] : gradf
|
||||
end
|
||||
|
||||
function (basis::CVTI)(xi::Vector, time::Number)
|
||||
basis(xi)
|
||||
end
|
||||
|
||||
function Base.:*(grad::Matrix, field::DVTI)
|
||||
n, m = size(grad)
|
||||
return sum([kron(grad[:,i], field[i]') for i=1:m])'
|
||||
end
|
||||
|
||||
function DVTV(data::Pair{Float64, Vector}...)
|
||||
return DVTV([Increment(d[1], d[2]) for d in data])
|
||||
end
|
||||
|
||||
function start(f::DVTV)
|
||||
return start(f.data)
|
||||
end
|
||||
|
||||
function next(f::DVTV, state)
|
||||
return next(f.data, state)
|
||||
end
|
||||
|
||||
function done(f::DVTV, state)
|
||||
return done(f.data, state)
|
||||
end
|
||||
|
||||
""" Return time vector from time variable field. """
|
||||
function keys(field::DVTV)
|
||||
return Float64[increment.time for increment in field]
|
||||
end
|
||||
|
||||
function setindex!(field::Field, val, idx::Int64)
|
||||
field.data[idx] = val
|
||||
end
|
||||
|
||||
@@ -26,8 +26,8 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
|
||||
add!(A, gdofs, gdofs, w*kron(N', N))
|
||||
end
|
||||
end
|
||||
nz = get_nonzero_rows(A)
|
||||
A = sparse(A)
|
||||
nz = get_nonzero_rows(A)
|
||||
A = 1/2*(A + A')
|
||||
F = ldltfact(A[nz,nz])
|
||||
end
|
||||
|
||||
@@ -33,27 +33,6 @@ function aster_parse_nodes(section; strip_characters=true)
|
||||
return nodes
|
||||
end
|
||||
|
||||
function parse(mesh, ::Type{Val{:CODE_ASTER_MAIL}})
|
||||
model = Dict()
|
||||
header = nothing
|
||||
data = []
|
||||
for line in split(mesh, '\n')
|
||||
length(line) != 0 || continue
|
||||
info("line: $line")
|
||||
if is_aster_mail_keyword(strip(line))
|
||||
header = parse_aster_header(line)
|
||||
empty!(data)
|
||||
continue
|
||||
end
|
||||
if line == "FINSF"
|
||||
info(data)
|
||||
header = nothing
|
||||
process_aster_section!(model, join(data, ""), header, Val{header[1]})
|
||||
end
|
||||
end
|
||||
return model
|
||||
end
|
||||
|
||||
|
||||
""" Code Aster binary file (.med). """
|
||||
type MEDFile
|
||||
@@ -161,7 +140,7 @@ Returns
|
||||
Dict containing fields "nodes" and "connectivity".
|
||||
|
||||
"""
|
||||
function parse_aster_med_file(fn, mesh_name=nothing; debug=false)
|
||||
function parse_aster_med_file(fn, mesh_name=nothing)
|
||||
med = MEDFile(fn)
|
||||
mesh_names = get_mesh_names(med::MEDFile)
|
||||
all_meshes = join(mesh_names, ", ")
|
||||
@@ -173,12 +152,10 @@ function parse_aster_med_file(fn, mesh_name=nothing; debug=false)
|
||||
end
|
||||
nsets = get_node_sets(med, mesh_name)
|
||||
elsets = get_element_sets(med, mesh_name)
|
||||
if debug
|
||||
elset_names = join(values(elsets), ", ")
|
||||
info("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names")
|
||||
nset_names = join(values(nsets), ", ")
|
||||
info("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names")
|
||||
end
|
||||
elset_names = join(values(elsets), ", ")
|
||||
debug("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names")
|
||||
nset_names = join(values(nsets), ", ")
|
||||
debug("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names")
|
||||
nodes = get_nodes(med, nsets, mesh_name)
|
||||
conn = get_connectivity(med, elsets, mesh_name)
|
||||
result = Dict("nodes" => nodes, "connectivity" => conn)
|
||||
|
||||
@@ -81,10 +81,6 @@ function isempty(assembly::Assembly)
|
||||
return T
|
||||
end
|
||||
|
||||
function get_dofs(assembly::Assembly)
|
||||
return sort(unique(assembly.K.J))
|
||||
end
|
||||
|
||||
type Problem{P<:AbstractProblem}
|
||||
name :: AbstractString # descriptive name for problem
|
||||
dimension :: Int # degrees of freedom per node
|
||||
@@ -359,10 +355,6 @@ function push!(problem::Problem, elements_::Vector...)
|
||||
end
|
||||
end
|
||||
|
||||
function get_connectivity(problem::Problem)
|
||||
return union([get_connectivity(element) for element in get_elements(problem)]...)
|
||||
end
|
||||
|
||||
function get_gdofs(element::Element, dim::Int)
|
||||
conn = get_connectivity(element)
|
||||
if length(conn) == 0
|
||||
@@ -372,10 +364,6 @@ function get_gdofs(element::Element, dim::Int)
|
||||
return gdofs
|
||||
end
|
||||
|
||||
function get_dofs(problem::Problem)
|
||||
return get_dofs(problem.assembly)
|
||||
end
|
||||
|
||||
function empty!(problem::Problem)
|
||||
empty!(problem.assembly)
|
||||
end
|
||||
@@ -396,33 +384,3 @@ function get_gdofs(problem::Problem, element::Element)
|
||||
end
|
||||
return problem.dofmap[element]
|
||||
end
|
||||
|
||||
""" Find dofs corresponding to nodes. """
|
||||
function find_dofs_by_nodes(problem::Problem, nodes)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
return find_dofs_by_nodes(dim, nodes)
|
||||
end
|
||||
function find_dofs_by_nodes(dim::Int, nodes)
|
||||
dofs = Int64[]
|
||||
for node in nodes
|
||||
for j=1:dim
|
||||
push!(dofs, dim*(node-1)+j)
|
||||
end
|
||||
end
|
||||
return dofs
|
||||
end
|
||||
|
||||
""" Find nodes corresponding to dofs. """
|
||||
function find_nodes_by_dofs(problem::Problem, dofs)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
return find_nodes_by_dofs(dim, dofs)
|
||||
end
|
||||
function find_nodes_by_dofs(dim, dofs)
|
||||
nodes = Int64[]
|
||||
for dof in dofs
|
||||
j = Int(ceil(dof/dim))
|
||||
j in nodes && continue
|
||||
push!(nodes, j)
|
||||
end
|
||||
return nodes
|
||||
end
|
||||
|
||||
@@ -287,273 +287,3 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
problem.assembly.g = g
|
||||
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Frictionless 2d small sliding contact without forwarddiff.
|
||||
|
||||
true/false flags: finite_sliding, friction, use_forwarddiff
|
||||
"""
|
||||
function _assemble!(problem::Problem{Contact}, time::Float64,
|
||||
::Type{Val{1}}, ::Type{Val{false}},
|
||||
::Type{Val{false}}, ::Type{Val{false}}; debug=false)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_elements = get_slave_elements(problem)
|
||||
|
||||
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
|
||||
normals, tangents = calculate_normals(slave_elements, time, Val{1};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
update!(slave_elements, "tangent", time => tangents)
|
||||
|
||||
Rn = 0.0
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = slave_element("displacement", time)
|
||||
la1 = slave_element("reaction force", time)
|
||||
n1 = slave_element("normal", time)
|
||||
t1 = slave_element("tangent", time)
|
||||
x1 = X1 + u1
|
||||
Q1_ = [n1[1] t1[1]]
|
||||
Q2_ = [n1[2] t1[2]]
|
||||
Z = zeros(2, 2)
|
||||
Q2 = [Q1_ Z; Z Q2_]
|
||||
contact_area = 0.0
|
||||
contact_error = 0.0
|
||||
|
||||
if "element area" in props.store_fields
|
||||
element_area = 0.0
|
||||
for ip in get_integration_points(slave_element)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
element_area += w
|
||||
end
|
||||
update!(slave_element, "element area", time => element_area)
|
||||
end
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
nm = length(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval
|
||||
continue
|
||||
end
|
||||
|
||||
if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval
|
||||
continue
|
||||
end
|
||||
|
||||
# 3.1 calculate segmentation
|
||||
xi1a = project_from_master_to_slave(slave_element, X2[1], time)
|
||||
xi1b = project_from_master_to_slave(slave_element, X2[2], time)
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
# 3.2. bi-orthogonal basis
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
Ae = zeros(nsl, nsl)
|
||||
if props.dual_basis
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
De += w*diagm(N1)
|
||||
Me += w*N1*N1'
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
else
|
||||
Ae = eye(nsl)
|
||||
end
|
||||
|
||||
# 3.3. loop integration points of one integration segment and calculate
|
||||
# local mortar matrices
|
||||
fill!(De, 0.0)
|
||||
fill!(Me, 0.0)
|
||||
ge = zeros(field_dim*nsl)
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
Phi = Ae*N1
|
||||
|
||||
# project gauss point from slave element to master element in direction n_s
|
||||
X_s = N1*X1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = N1*t1 # tangent condition in gauss point
|
||||
n_s /= norm(n_s)
|
||||
t_s /= norm(t_s)
|
||||
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
X_m = N2*X2
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
x_s = X_s + u_s
|
||||
x_m = X_m + u_m
|
||||
la_s = Phi*la1
|
||||
ge += w*vec((x_m-x_s)*Phi')
|
||||
|
||||
# virtual work
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
contact_area += w
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end
|
||||
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
|
||||
# add contribution to contact virtual work
|
||||
D2 = zeros(field_dim*nsl, field_dim*nsl)
|
||||
M2 = zeros(field_dim*nsl, field_dim*nsl)
|
||||
for i=1:field_dim
|
||||
D2[i:field_dim:end, i:field_dim:end] += De
|
||||
M2[i:field_dim:end, i:field_dim:end] += Me
|
||||
end
|
||||
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D2)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M2)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2)
|
||||
ge = -D2*vec(x1)+M2*vec(x2)
|
||||
add!(problem.assembly.g, sdofs, Q2'*ge)
|
||||
ce = vec(la1) + ge
|
||||
add!(problem.assembly.c, sdofs, Q2'*ce)
|
||||
|
||||
end # master elements done
|
||||
|
||||
if "contact area" in props.store_fields
|
||||
update!(slave_element, "contact area", time => contact_area)
|
||||
end
|
||||
|
||||
if "contact error" in props.store_fields
|
||||
update!(slave_element, "contact error", time => contact_error)
|
||||
end
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
S = sort(collect(keys(normals))) # slave element nodes
|
||||
weighted_gap = Dict{Int64, Vector{Float64}}()
|
||||
contact_pressure = Dict{Int64, Vector{Float64}}()
|
||||
complementarity_condition = Dict{Int64, Vector{Float64}}()
|
||||
is_active = Dict{Int64, Int}()
|
||||
is_inactive = Dict{Int64, Int}()
|
||||
is_slip = Dict{Int64, Int}()
|
||||
is_stick = Dict{Int64, Int}()
|
||||
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C2 = sparse(problem.assembly.C2, ndofs, ndofs)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
c = full(problem.assembly.c, ndofs, 1)
|
||||
g = full(problem.assembly.g, ndofs, 1)
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
weighted_gap[j] = g[dofs]
|
||||
|
||||
if length(la) != 0
|
||||
p = dot(normals[j], la[dofs])
|
||||
t = dot(tangents[j], la[dofs])
|
||||
contact_pressure[j] = [p, t]
|
||||
else
|
||||
contact_pressure[j] = [0.0, 0.0]
|
||||
end
|
||||
|
||||
#complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
complementarity_condition[j] = c[dofs]
|
||||
if complementarity_condition[j][1] < 0
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
is_stick[j] = 0
|
||||
else
|
||||
is_inactive[j] = 0
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
if "weighted gap" in props.store_fields
|
||||
update!(slave_elements, "weighted gap", time => weighted_gap)
|
||||
end
|
||||
if "contact pressure" in props.store_fields
|
||||
update!(slave_elements, "contact pressure", time => contact_pressure)
|
||||
end
|
||||
if "complementarity condition" in props.store_fields
|
||||
update!(slave_elements, "complementarity condition", time => complementarity_condition)
|
||||
end
|
||||
if "active nodes" in props.store_fields
|
||||
update!(slave_elements, "active nodes", time => is_active)
|
||||
end
|
||||
if "inactive nodes" in props.store_fields
|
||||
update!(slave_elements, "inactive nodes", time => is_inactive)
|
||||
end
|
||||
if "stick nodes" in props.store_fields
|
||||
update!(slave_elements, "stick nodes", time => is_stick)
|
||||
end
|
||||
if "slip nodes" in props.store_fields
|
||||
update!(slave_elements, "slip nodes", time => is_slip)
|
||||
end
|
||||
|
||||
# info("# | active | inactive | stick | slip | gap | pres | comp")
|
||||
# for j in S
|
||||
# str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
|
||||
# str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
|
||||
# info(str1 * str2)
|
||||
# end
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
# constitutive modelling in tangent direction, frictionless contact
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if (is_active[j] == 1) && (is_slip[j] == 1)
|
||||
# info("$j is in active/slip, removing tangential constraint $(dofs[2])")
|
||||
C2[dofs[2],:] = 0.0
|
||||
g[dofs[2]] = 0.0
|
||||
D[dofs[2], dofs] = tangents[j]
|
||||
end
|
||||
end
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if is_inactive[j] == 1
|
||||
# info("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
g[dofs,:] = 0.0
|
||||
end
|
||||
end
|
||||
|
||||
problem.assembly.C1 = C1
|
||||
problem.assembly.C2 = C2
|
||||
problem.assembly.D = D
|
||||
problem.assembly.g = g
|
||||
|
||||
end
|
||||
|
||||
|
||||
@@ -304,151 +304,6 @@ function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity},
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
""" Elasticity equations, 3d, linear. """
|
||||
function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
|
||||
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = length(element)
|
||||
ndofs = dim*nnodes
|
||||
BL = zeros(6, ndofs)
|
||||
Km = zeros(ndofs, ndofs)
|
||||
Kg = zeros(ndofs, ndofs)
|
||||
f = zeros(ndofs)
|
||||
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
|
||||
fill!(BL, 0.0)
|
||||
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
|
||||
|
||||
E = element("youngs modulus", ip, time)
|
||||
nu = element("poissons ratio", ip, time)
|
||||
|
||||
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
|
||||
1.0-nu nu nu 0.0 0.0 0.0
|
||||
nu 1.0-nu nu 0.0 0.0 0.0
|
||||
nu nu 1.0-nu 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.5-nu 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.5-nu 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.5-nu]
|
||||
|
||||
Km += w*BL'*D*BL
|
||||
|
||||
if haskey(element, "displacement load")
|
||||
T = element("displacement load", ip, time)
|
||||
f += w*vec(T*N)
|
||||
end
|
||||
for i=1:dim
|
||||
if haskey(element, "displacement load $i")
|
||||
b = element("displacement load $i", ip, time)
|
||||
f[i:dim:end] += w*vec(b*N)
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if get_formulation_type(problem) == :incremental
|
||||
if haskey(element, "displacement")
|
||||
u = vec(element["displacement"](time))
|
||||
f -= Kt*u
|
||||
end
|
||||
end
|
||||
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
""" Material and geometric stiffness for linear buckling analysis. """
|
||||
function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_buckling}})
|
||||
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = length(element)
|
||||
ndofs = dim*nnodes
|
||||
BL = zeros(6, ndofs)
|
||||
BNL = zeros(9, ndofs)
|
||||
Km = zeros(ndofs, ndofs)
|
||||
Kg = zeros(ndofs, ndofs)
|
||||
f = zeros(ndofs)
|
||||
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
|
||||
gradu = element("displacement", ip, time, Val{:Grad})
|
||||
strain = 1/2*(gradu' + gradu)
|
||||
|
||||
fill!(BL, 0.0)
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
E = element("youngs modulus", ip, time)
|
||||
nu = element("poissons ratio", ip, time)
|
||||
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
|
||||
1.0-nu nu nu 0.0 0.0 0.0
|
||||
nu 1.0-nu nu 0.0 0.0 0.0
|
||||
nu nu 1.0-nu 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.5-nu 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.5-nu 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.5-nu]
|
||||
|
||||
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
|
||||
stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
|
||||
|
||||
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]
|
||||
|
||||
Km += w*BL'*D*BL
|
||||
Kg += w*BNL'*S3*BNL
|
||||
|
||||
end
|
||||
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
""" Elasticity equations, 3d nonlinear. """
|
||||
function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
|
||||
props = problem.properties
|
||||
@@ -679,185 +534,3 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity},
|
||||
end
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
|
||||
return assemble(problem, element, time, Val{:continuum})
|
||||
end
|
||||
|
||||
""" Elasticity equations using ForwardDiff
|
||||
"""
|
||||
function assemble(problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:forwarddiff}})
|
||||
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = size(element, 2)
|
||||
|
||||
function get_residual_vector(u::Vector)
|
||||
u = reshape(u, dim, nnodes)
|
||||
u = Field([u[:,i] for i=1:nnodes])
|
||||
r = zeros(dim, nnodes)
|
||||
|
||||
for ip in get_integration_points(element)
|
||||
|
||||
JT = transpose(get_jacobian(element, ip, time))
|
||||
n, m = size(JT)
|
||||
if n == m
|
||||
w = ip.weight*det(JT)
|
||||
elseif m == 1
|
||||
w = ip.weight*norm(JT)
|
||||
elseif m == 2
|
||||
w = ip.weight*norm(cross(JT[:,1], JT[:,2]))
|
||||
else
|
||||
error("jacobian $JT")
|
||||
end
|
||||
|
||||
# calculate internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# kinematics
|
||||
F = I + gradu
|
||||
E = 1/2*(F'*F - I)
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if problem.properties.formulation == :plane_stress
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for plane stress
|
||||
end
|
||||
|
||||
# stress
|
||||
S = lambda*trace(E)*I + 2*mu*E
|
||||
|
||||
r += w*F*S*grad
|
||||
end
|
||||
|
||||
# calculate external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= w*b*basis
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
r -= w*T*basis
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
|
||||
end
|
||||
|
||||
field = element("displacement", time)
|
||||
Km, allresults = ForwardDiff.jacobian(get_residual_vector, vec(field),
|
||||
AllResults, cache=autodiffcache)
|
||||
Kg = zeros(Km)
|
||||
f = -ForwardDiff.value(allresults)
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
|
||||
###############################
|
||||
# Plastic material #
|
||||
###############################
|
||||
#=
|
||||
abstract PlaneStressLinearElasticPlasticProblem <: LinearElasticityProblem
|
||||
|
||||
function PlaneStressLinearElasticPlasticProblem(name="plane stress linear elasticity", dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressLinearElasticPlasticProblem}(name, dim, elements)
|
||||
end
|
||||
|
||||
""" Elasticity equations, plane stress. """
|
||||
function assemble!{E<:CG, P<:PlaneStressLinearElasticPlasticProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real)
|
||||
|
||||
gdofs = get_gdofs(element, problem.dim)
|
||||
ndim, nnodes = size(E)
|
||||
B = zeros(3, 2*nnodes)
|
||||
for ip in get_integration_points(element)
|
||||
w = ip.weight
|
||||
J = get_jacobian(element, ip, time)
|
||||
N = element(ip, time)
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
nu = element("poissons ratio", ip, time)
|
||||
E_ = element("youngs modulus", ip, time)
|
||||
C = E_/(1.0 - nu^2) .* [
|
||||
1.0 nu 0.0
|
||||
nu 1.0 0.0
|
||||
0.0 0.0 (1.0-nu)/2.0]
|
||||
dN = element(ip, time, Val{:grad})
|
||||
fill!(B, 0.0)
|
||||
for i=1:size(dN, 2)
|
||||
B[1, 2*(i-1)+1] = dN[1,i]
|
||||
B[2, 2*(i-1)+2] = dN[2,i]
|
||||
B[3, 2*(i-1)+1] = dN[2,i]
|
||||
B[3, 2*(i-1)+2] = dN[1,i]
|
||||
end
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J))
|
||||
end
|
||||
if haskey(element, "displacement load")
|
||||
b = element("displacement load", ip, time)
|
||||
add!(assembly.force_vector, gdofs, w*N'*b*det(J))
|
||||
end
|
||||
if haskey(element, "displacement traction force")
|
||||
T = element("displacement traction force", ip, time)
|
||||
L = w*T*N*norm(J)
|
||||
add!(assembly.force_vector, gdofs, vec(L))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
include("elasticplastic.jl")
|
||||
|
||||
# Elasticity problems
|
||||
abstract ElasticityProblem <: AbstractProblem
|
||||
abstract PlaneStressElasticityProblem <: ElasticityProblem
|
||||
|
||||
function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
|
||||
return "displacement"
|
||||
end
|
||||
|
||||
function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
|
||||
return Vector{Float64}
|
||||
end
|
||||
|
||||
|
||||
|
||||
=#
|
||||
|
||||
function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:E}})
|
||||
haskey(element, "displacement") || return nothing
|
||||
gradu = element("displacement", ip, time, Val{:Grad})
|
||||
eps = 0.5*(gradu + gradu')
|
||||
return eps
|
||||
end
|
||||
|
||||
function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:S}})
|
||||
haskey(element, "displacement") || return nothing
|
||||
props = problem.properties
|
||||
eps = problem(element, ip, time, Val{:E})
|
||||
eps == nothing && return nothing
|
||||
E = element("youngs modulus", ip, time)
|
||||
nu = element("poissons ratio", ip, time)
|
||||
mu = E/(2.0*(1.0+nu))
|
||||
la = E*nu/((1.0+nu)*(1.0-2.0*nu))
|
||||
if props.formulation in [:plane_stress, :plane_strain]
|
||||
la = 2.0*la*mu/(la+2.0*mu)
|
||||
end
|
||||
S = la*trace(eps)*I + 2.0*mu*eps
|
||||
return S
|
||||
end
|
||||
|
||||
function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:COORD}})
|
||||
haskey(element, "geometry") || return nothing
|
||||
return element("geometry", ip, time)
|
||||
end
|
||||
|
||||
@@ -1,237 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Elasticity problems
|
||||
abstract ElasticPlasticProblem <: AbstractProblem
|
||||
abstract PlaneStressElasticPlasticProblem <: ElasticPlasticProblem
|
||||
|
||||
function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return "displacement"
|
||||
end
|
||||
|
||||
function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return Vector{Float64}
|
||||
end
|
||||
|
||||
# 3D Elasticity problems
|
||||
function ElasticPlasticProblem(dim::Int=3, elements=[])
|
||||
return Problem{ElasticPlasticProblem}("elasticplastic problem", dim, elements)
|
||||
end
|
||||
|
||||
# 2D Plane stress elasticity problems
|
||||
function PlaneStressElasticPlasticProblem(dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressElasticPlasticProblem}("plane stress elasticplastic problem", dim, elements)
|
||||
end
|
||||
|
||||
|
||||
function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
# internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
if !haskey(element, "integration points")
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
last_stress = zeros(2,2)
|
||||
last_strain = zeros(2,2)
|
||||
else
|
||||
last_stress = zeros(3,3)
|
||||
last_strain = zeros(3,3)
|
||||
end
|
||||
else
|
||||
for each_ip in element("integration points", time)
|
||||
if isapprox(each_ip.xi, ip.xi)
|
||||
last_stress = ip("stress", time)
|
||||
last_strain = ip("stress", time)
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# deformation gradient
|
||||
F = I + gradu
|
||||
E = 1/2*(F'*F - I)
|
||||
|
||||
#E = 1/2*(gradu + gradu') # finite strain (total)
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
stress_y = element("yield stress", time).data
|
||||
dstrain = E - last_strain
|
||||
de_v = [dstrain[1,1], dstrain[2,2], dstrain[1,2]]
|
||||
material_model = element("material model", time)
|
||||
s = last_stress
|
||||
de = copy(ForwardDiff.get_value(dstrain))
|
||||
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
C = stiffnessTensorPlaneStress(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[1,2]]
|
||||
problem_stress_type = :PlaneStressElasticPlasticProblem
|
||||
else
|
||||
C = stiffnessTensor(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
|
||||
problem_stress_type = :ElasticPlasticProblem
|
||||
end
|
||||
dep = zeros(3)
|
||||
stress_inc, dep = calculate_stress(de_,
|
||||
s_v,
|
||||
C,
|
||||
stress_y,
|
||||
Val{:vonMises},
|
||||
Val{problem_stress_type})
|
||||
info("%% ", dep)
|
||||
s_v += C * (de_v - dep)
|
||||
info("--: ", ForwardDiff.get_value(s_v))
|
||||
# stress
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
S = [s_v[1] s_v[3];
|
||||
s_v[3] s_v[2]]
|
||||
else
|
||||
S = [s_v[1] s_v[6] s_v[5];
|
||||
s_v[6] s_v[2] s_v[4];
|
||||
s_v[5] s_v[4] s_v[3]]
|
||||
end
|
||||
r += F*S*grad*det(J)
|
||||
end
|
||||
|
||||
# external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= b*basis*det(J)
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
JT = transpose(J)
|
||||
s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
|
||||
r -= T*basis*norm(s)
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
end
|
||||
|
||||
|
||||
|
||||
#=
|
||||
function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
info("_____________________")
|
||||
# internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
|
||||
if !haskey(element, "integration points")
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
last_stress = zeros(2,2)
|
||||
last_strain = zeros(2,2)
|
||||
else
|
||||
last_stress = zeros(3,3)
|
||||
last_strain = zeros(3,3)
|
||||
end
|
||||
else
|
||||
for each_ip in element("integration points", time)
|
||||
if isapprox(each_ip.xi, ip.xi)
|
||||
last_stress = ip("stress", time)
|
||||
last_strain = ip("stress", time)
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# deformation gradient
|
||||
F = I + gradu
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if P == PlaneStressElasticityProblem
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems
|
||||
end
|
||||
|
||||
# strain
|
||||
E = 1/2*(F'*F - I)
|
||||
#E = 1/2*(gradu + gradu') # finite strain (total)
|
||||
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
stress_y = element("yield stress", time).data
|
||||
dstrain = E - last_strain
|
||||
material_model = element("material model", time)
|
||||
s = last_stress
|
||||
de = ForwardDiff.get_value(dstrain)
|
||||
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
C = stiffnessTensorPlaneStress(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[1,2]]
|
||||
problem_stress_type = :PlaneStressElasticPlasticProblem
|
||||
else
|
||||
C = stiffnessTensor(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
|
||||
problem_stress_type = :ElasticPlasticProblem
|
||||
end
|
||||
|
||||
stress_inc, lambda = plastic_multiplier = calculate_stress(de_,
|
||||
s_v,
|
||||
C,
|
||||
stress_y,
|
||||
Val{:vonMises},
|
||||
Val{problem_stress_type})
|
||||
|
||||
# dep = lambda * dfds(s)
|
||||
# upate_material_parameters!(...)
|
||||
s_new = s_v + stress_inc
|
||||
#S = [s_v[1] s_v[6] s_v[5];
|
||||
# s_v[6] s_v[2] s_v[4];
|
||||
# s_v[5] s_v[4] s_v[3]]
|
||||
S = [s_new[1] s_new[3];
|
||||
s_new[3] s_new[2]]
|
||||
# S = C * (E - dep)
|
||||
|
||||
|
||||
info("Stress: ", vec(ForwardDiff.get_value(S)))
|
||||
# stress
|
||||
#S = lambda*trace(E)*I + 2*mu*E
|
||||
|
||||
r += F*S*grad*det(J)
|
||||
|
||||
end
|
||||
|
||||
|
||||
# external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= b*basis*det(J)
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
JT = transpose(J)
|
||||
s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
|
||||
r -= T*basis*norm(s)
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
end
|
||||
=# #fff
|
||||
+47
-55
@@ -34,15 +34,15 @@ function vertex_inside_polygon(q, P; atol=1.0e-3)
|
||||
cosa = dot(A,B)/c
|
||||
isapprox(cosa, 1.0; atol=atol) && return false
|
||||
isapprox(cosa, -1.0; atol=atol) && return true
|
||||
try
|
||||
angle += acos(cosa)
|
||||
catch
|
||||
info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
|
||||
info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
|
||||
info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
|
||||
info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
|
||||
rethrow()
|
||||
end
|
||||
#try
|
||||
angle += acos(cosa)
|
||||
#catch
|
||||
# info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
|
||||
# info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
|
||||
# info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
|
||||
# info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
|
||||
# rethrow()
|
||||
#end
|
||||
end
|
||||
return isapprox(angle, 2*pi; atol=atol)
|
||||
end
|
||||
@@ -68,26 +68,8 @@ function get_cells(P, C; allow_quads=false)
|
||||
if N == 4 && allow_quads
|
||||
return Vector[P]
|
||||
end
|
||||
#V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N])
|
||||
#A = 1/2*abs(dot(n, V))
|
||||
#info("A = $A")
|
||||
cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N]
|
||||
return cells
|
||||
|
||||
maxa = 0.0
|
||||
maxj = 0
|
||||
for i=1:N
|
||||
A = P[i] - C
|
||||
B = P[mod(i,N)+1] - C
|
||||
theta = acos(dot(A,B)/(norm(A)*norm(B)))
|
||||
if theta > maxa
|
||||
maxa = theta
|
||||
maxj = i
|
||||
end
|
||||
end
|
||||
info("max angle $(maxa/pi*180) at index $maxj, N=$N")
|
||||
indices = mod(collect(maxj:maxj+N), N)
|
||||
info("indices = $indices")
|
||||
end
|
||||
|
||||
""" Test does P contain q. """
|
||||
@@ -175,7 +157,7 @@ function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
|
||||
return theta[1:2], theta[3]
|
||||
end
|
||||
end
|
||||
|
||||
#=
|
||||
info("failed to project vertex from auxiliary plane back to surface")
|
||||
info("element type: $E")
|
||||
info("element connectivity: $(get_connectivity(element))")
|
||||
@@ -199,7 +181,7 @@ function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
|
||||
info("dtheta = $(dtheta)")
|
||||
theta -= dtheta
|
||||
end
|
||||
|
||||
=#
|
||||
throw(error("project_point_to_surface: did not converge in $max_iterations iterations!"))
|
||||
end
|
||||
|
||||
@@ -272,8 +254,10 @@ function split_quadratic_element(element::Element{Tri6}, time::Float64)
|
||||
new_element = Element(Tri3, connectivity[elmap])
|
||||
X = element("geometry", time)
|
||||
update!(new_element, "geometry", time => X[elmap])
|
||||
u = element("displacement", time)
|
||||
update!(new_element, "displacement", time => u[elmap])
|
||||
if haskey(element, "displacement")
|
||||
u = element("displacement", time)
|
||||
update!(new_element, "displacement", time => u[elmap])
|
||||
end
|
||||
#n = element("normal", time)
|
||||
#update!(new_element, "normal", time => n[elmap])
|
||||
if haskey(element, "master elements")
|
||||
@@ -372,6 +356,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
check_orientation!(P, n0)
|
||||
N_P = length(P)
|
||||
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
|
||||
|
||||
if first_slave_element
|
||||
debug("Polygon clip info for first slave element:")
|
||||
debug("S = $S")
|
||||
@@ -380,11 +365,15 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
debug("N_P = $N_P")
|
||||
debug("P_area = $P_area")
|
||||
end
|
||||
|
||||
if isapprox(P_area, 0.0)
|
||||
info("Polygon P has zero area: $P_area")
|
||||
continue
|
||||
end
|
||||
|
||||
C0 = calculate_centroid(P)
|
||||
|
||||
#=
|
||||
if isnan(C0[1])
|
||||
info("C0 = $C0")
|
||||
info("P = $P")
|
||||
@@ -393,6 +382,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
info("n0 = $n0")
|
||||
error("Calculation of centroid of polygon clip P failed.")
|
||||
end
|
||||
=#
|
||||
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
@@ -412,7 +402,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
x_gauss = nothing
|
||||
try
|
||||
#try
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
@@ -420,17 +410,17 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1*N1'
|
||||
catch
|
||||
info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.")
|
||||
info("x_gauss = $x_gauss")
|
||||
info("cell = $cell")
|
||||
info("C0 = $C0")
|
||||
info("P = $P")
|
||||
info("S = $S")
|
||||
info("M = $M")
|
||||
info("n0 = $n0")
|
||||
rethrow()
|
||||
end
|
||||
#catch
|
||||
# info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.")
|
||||
# info("x_gauss = $x_gauss")
|
||||
# info("cell = $cell")
|
||||
# info("C0 = $C0")
|
||||
# info("P = $P")
|
||||
# info("S = $S")
|
||||
# info("M = $M")
|
||||
# info("n0 = $n0")
|
||||
# rethrow()
|
||||
#end
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
else
|
||||
@@ -449,6 +439,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
# project gauss point from auxiliary plane to master and slave element
|
||||
#x_gauss = N*x_cell
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
#=
|
||||
if isnan(x_gauss[1])
|
||||
info("is nan")
|
||||
info("x_gauss = $x_gauss")
|
||||
@@ -460,25 +451,26 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
|
||||
info("n0 = $n0")
|
||||
error("nan, unable to continue")
|
||||
end
|
||||
=#
|
||||
|
||||
xi_s = nothing
|
||||
xi_m = nothing
|
||||
alpha = nothing
|
||||
|
||||
try
|
||||
#try
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
|
||||
catch
|
||||
info("projecting vertex back to surface has failed.")
|
||||
info("x_gauss = $x_gauss")
|
||||
info("cell = $cell")
|
||||
info("C0 = $C0")
|
||||
info("P = $P")
|
||||
info("S = $S")
|
||||
info("M = $M")
|
||||
info("n0 = $n0")
|
||||
rethrow()
|
||||
end
|
||||
#catch
|
||||
# info("projecting vertex back to surface has failed.")
|
||||
# info("x_gauss = $x_gauss")
|
||||
# info("cell = $cell")
|
||||
# info("C0 = $C0")
|
||||
# info("P = $P")
|
||||
# info("S = $S")
|
||||
# info("M = $M")
|
||||
# info("n0 = $n0")
|
||||
# rethrow()
|
||||
#end
|
||||
|
||||
# add contributions
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
|
||||
+41
-132
@@ -55,24 +55,6 @@ 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))
|
||||
|
||||
"""
|
||||
Posthook for field assembly. By default, do nothing.
|
||||
This can be used to make some modifications for assembly
|
||||
after all elements are assembled.
|
||||
|
||||
Examples
|
||||
--------
|
||||
function field_assembly_posthook!(solver::Solver,
|
||||
K::SparseMatrixCSC,
|
||||
Kg::SparseMatrixCSC,
|
||||
f::SparseMatrixCSC,
|
||||
fg::SpareMatrixCSC)
|
||||
info("doing stuff, size(K) = ", size(K))
|
||||
end
|
||||
"""
|
||||
function field_assembly_posthook!
|
||||
end
|
||||
|
||||
"""Return one combined field assembly for a set of field problems.
|
||||
|
||||
Parameters
|
||||
@@ -120,12 +102,6 @@ function get_field_assembly(solver::Solver; show_info=true)
|
||||
f = sparse(f, solver.ndofs, 1)
|
||||
fg = sparse(fg, solver.ndofs, 1)
|
||||
|
||||
# run any posthook for assembly if defined
|
||||
args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC}
|
||||
if method_exists(field_assembly_posthook!, args)
|
||||
field_assembly_posthook!(solver, K, Kg, fg, fg)
|
||||
end
|
||||
|
||||
return M, K, Kg, f, fg
|
||||
end
|
||||
|
||||
@@ -191,9 +167,11 @@ function get_boundary_assembly(solver::Solver)
|
||||
g_ = sparse(assembly.g, ndofs, 1)
|
||||
for dof in assembly.removed_dofs
|
||||
info("$(problem.name): removing dof $dof from assembly")
|
||||
C1_[:,dof] = 0.0
|
||||
C1_[dof,:] = 0.0
|
||||
C2_[dof,:] = 0.0
|
||||
end
|
||||
SparseArrays.dropzeros!(C1_)
|
||||
SparseArrays.dropzeros!(C2_)
|
||||
|
||||
already_constrained = get_nonzero_rows(C2)
|
||||
new_constraints = get_nonzero_rows(C2_)
|
||||
@@ -203,8 +181,6 @@ function get_boundary_assembly(solver::Solver)
|
||||
warn("already constrained = $already_constrained")
|
||||
warn("new constraints = $new_constraints")
|
||||
overconstrained_dofs = sort(overconstrained_dofs)
|
||||
overconstrained_nodes = find_nodes_by_dofs(problem, overconstrained_dofs)
|
||||
warn("in overconstrained nodes $overconstrained_nodes")
|
||||
error("overconstrained dofs, not solving problem.")
|
||||
end
|
||||
|
||||
@@ -224,10 +200,9 @@ Solve linear system using LDLt factorization (SuiteSparse). This version
|
||||
requires that final system is symmetric and positive definite, so boundary
|
||||
conditions are first eliminated before solution.
|
||||
"""
|
||||
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
|
||||
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}})
|
||||
|
||||
nnz(D) == 0 || return false
|
||||
C1 == C2 || return false
|
||||
|
||||
A = get_nonzero_rows(K)
|
||||
B = get_nonzero_rows(C2)
|
||||
@@ -235,28 +210,14 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug
|
||||
B == B2 || return false
|
||||
I = setdiff(A, B)
|
||||
|
||||
if debug
|
||||
info("# A = $(length(A))")
|
||||
info("# B = $(length(B))")
|
||||
info("# I = $(length(I))")
|
||||
end
|
||||
debug("# A = $(length(A))")
|
||||
debug("# B = $(length(B))")
|
||||
debug("# I = $(length(I))")
|
||||
|
||||
if length(B) == 0
|
||||
warn("No rows in C2, forget to set Dirichlet boundary conditions to model?")
|
||||
else
|
||||
# solver boundary dofs (usually a trivial solution Iu = g
|
||||
try
|
||||
u[B] = lufact(C2[B,B2]) \ full(g[B])
|
||||
catch
|
||||
info("solver #1 failed to solve boundary dofs (you should not see this message).")
|
||||
info("# A = $(length(A))")
|
||||
info("# B = $(length(B))")
|
||||
info("# B2 = $(length(B2))")
|
||||
info("# I = $(length(I))")
|
||||
info("B = $B")
|
||||
info("B2 = $B2")
|
||||
rethrow()
|
||||
end
|
||||
u[B] = lufact(C2[B,B2]) \ full(g[B])
|
||||
end
|
||||
|
||||
# solve interior domain using LDLt factorization
|
||||
@@ -287,12 +248,9 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
|
||||
end
|
||||
|
||||
""" Default linear system solver for solver. """
|
||||
function solve!(solver::Solver; empty_assemblies_before_solution=true,
|
||||
show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false)
|
||||
function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric=true)
|
||||
|
||||
if show_info
|
||||
info("Solving problems ...")
|
||||
end
|
||||
info("Solving problems ...")
|
||||
t0 = Base.time()
|
||||
|
||||
# assemble field & boundary problems
|
||||
@@ -310,21 +268,10 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false)
|
||||
M = 1/2*(M + M')
|
||||
end
|
||||
|
||||
if fill_D_diagonal
|
||||
nz = ones(solver.ndofs)
|
||||
nz[get_nonzero_rows(C2)] = 0.0
|
||||
nz[get_nonzero_rows(D)] = 0.0
|
||||
D += spdiagm(nz)
|
||||
end
|
||||
|
||||
# free up some memory before solution by either emptying field assemblies
|
||||
# or combining values with same indices in sparse COO matrices. Small
|
||||
# boundary problems are untouched.
|
||||
for problem in get_field_problems(solver)
|
||||
if empty_assemblies_before_solution
|
||||
if empty_assemblies_before_solution
|
||||
# free up some memory before solution by emptying field assemblies from problems
|
||||
for problem in get_field_problems(solver)
|
||||
empty!(problem.assembly)
|
||||
elseif optimize
|
||||
optimize!(problem.assembly)
|
||||
end
|
||||
gc()
|
||||
end
|
||||
@@ -332,13 +279,17 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false)
|
||||
ndofs = solver.ndofs
|
||||
u = zeros(ndofs)
|
||||
la = zeros(ndofs)
|
||||
status = false
|
||||
is_solved = false
|
||||
i = 0
|
||||
for i in [1, 2]
|
||||
status = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
|
||||
status && break
|
||||
is_solved = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
|
||||
if is_solved
|
||||
break
|
||||
end
|
||||
end
|
||||
if !is_solved
|
||||
error("Failed to solve linear system!")
|
||||
end
|
||||
status || error("Failed to solve linear system!")
|
||||
t1 = round(Base.time()-t0, 2)
|
||||
norms = (norm(u), norm(la))
|
||||
push!(solver.norms, norms)
|
||||
@@ -346,10 +297,8 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false)
|
||||
solver.u = u
|
||||
solver.la = la
|
||||
|
||||
if show_info
|
||||
info("Solved problems in $t1 seconds using solver $i.")
|
||||
info("Solution norms = $norms.")
|
||||
end
|
||||
info("Solved problems in $t1 seconds using solver $i.")
|
||||
info("Solution norms = $norms.")
|
||||
|
||||
return
|
||||
end
|
||||
@@ -459,22 +408,6 @@ function get_all_elements(solver::Solver)
|
||||
return [elements...;]
|
||||
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 is(E, element_type)
|
||||
end
|
||||
|
||||
function filter_by_element_type(element_type, elements)
|
||||
return filter(element -> is_element_type(element, element_type), elements)
|
||||
end
|
||||
|
||||
function (solver::Solver)(field_name::AbstractString, time::Float64)
|
||||
fields = []
|
||||
for problem in get_problems(solver)
|
||||
@@ -644,47 +577,27 @@ Notes
|
||||
-----
|
||||
Default convergence criteria is obtained by checking each sub-problem convergence.
|
||||
"""
|
||||
function has_converged(solver::Solver{Nonlinear}; show_info=false,
|
||||
check_convergence_for_boundary_problems=false)
|
||||
function has_converged(solver::Solver{Nonlinear})
|
||||
properties = solver.properties
|
||||
converged = true
|
||||
eps = properties.convergence_tolerance
|
||||
for problem in solver.problems
|
||||
has_converged = true
|
||||
if is_field_problem(problem)
|
||||
has_converged = problem.assembly.u_norm_change < eps
|
||||
if isapprox(norm(problem.assembly.u), 0.0)
|
||||
# trivial solution
|
||||
has_converged = true
|
||||
end
|
||||
show_info && info("Details for problem $(problem.name)")
|
||||
show_info && info("Norm: $(norm(problem.assembly.u))")
|
||||
show_info && info("Norm change: $(problem.assembly.u_norm_change)")
|
||||
show_info && info("Has converged? $(has_converged)")
|
||||
end
|
||||
if is_boundary_problem(problem) && check_convergence_for_boundary_problems
|
||||
has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps
|
||||
show_info && info("Details for problem $(problem.name)")
|
||||
show_info && info("Norm: $(norm(problem.assembly.la))")
|
||||
show_info && info("Norm change: $(problem.assembly.la_norm_change)")
|
||||
show_info && info("Has converged? $(has_converged)")
|
||||
for problem in get_field_problems(solver)
|
||||
has_converged = problem.assembly.u_norm_change < eps
|
||||
if isapprox(norm(problem.assembly.u), 0.0)
|
||||
# trivial solution
|
||||
has_converged = true
|
||||
end
|
||||
debug("Details for problem $(problem.name)")
|
||||
debug("Norm: $(norm(problem.assembly.u))")
|
||||
debug("Norm change: $(problem.assembly.u_norm_change)")
|
||||
debug("Has converged? $(has_converged)")
|
||||
converged &= has_converged
|
||||
end
|
||||
return converged
|
||||
end
|
||||
|
||||
type NonlinearConvergenceError <: Exception
|
||||
solver :: Solver
|
||||
end
|
||||
|
||||
function Base.showerror(io::IO, exception::NonlinearConvergenceError)
|
||||
max_iters = exception.solver.properties.max_iterations
|
||||
print(io, "nonlinear iteration did not converge in $max_iters iterations!")
|
||||
end
|
||||
|
||||
""" Default solver for quasistatic nonlinear problems. """
|
||||
function (solver::Solver{Nonlinear})(; show_info=true)
|
||||
function (solver::Solver{Nonlinear})()
|
||||
|
||||
properties = solver.properties
|
||||
|
||||
@@ -693,10 +606,10 @@ function (solver::Solver{Nonlinear})(; show_info=true)
|
||||
|
||||
# 2. start non-linear iterations
|
||||
for properties.iteration=1:properties.max_iterations
|
||||
show_info && info(repeat("-", 80))
|
||||
show_info && info("Starting nonlinear iteration #$(properties.iteration)")
|
||||
show_info && info("Increment time t=$(round(solver.time, 3))")
|
||||
show_info && info(repeat("-", 80))
|
||||
info(repeat("-", 80))
|
||||
info("Starting nonlinear iteration #$(properties.iteration)")
|
||||
info("Increment time t=$(round(solver.time, 3))")
|
||||
info(repeat("-", 80))
|
||||
|
||||
# 2.1 update linearized assemblies
|
||||
assemble!(solver)
|
||||
@@ -714,7 +627,9 @@ function (solver::Solver{Nonlinear})(; show_info=true)
|
||||
end
|
||||
|
||||
# 3. did not converge
|
||||
properties.error_if_no_convergence && throw(NonlinearConvergenceError(solver))
|
||||
if properties.error_if_no_convergence
|
||||
error("nonlinear iteration did not converge in $(properties.iteration) iterations!")
|
||||
end
|
||||
end
|
||||
|
||||
""" Convenience function to call nonlinear solver. """
|
||||
@@ -851,9 +766,3 @@ function Postprocessor(problems::Problem...)
|
||||
end
|
||||
return solver
|
||||
end
|
||||
|
||||
function Postprocessor(name::AbstractString, problems::Problem...)
|
||||
solver = Postprocessor(problems...)
|
||||
solver.name = name
|
||||
return solver
|
||||
end
|
||||
|
||||
+1
-49
@@ -32,10 +32,6 @@ function convert(::Type{SparseMatrixCOO}, A::Matrix)
|
||||
return SparseMatrixCOO(findnz(A)...)
|
||||
end
|
||||
|
||||
function convert(::Type{SparseMatrixCOO}, A::Vector)
|
||||
return SparseMatrixCOO(findnz(sparse(A))...)
|
||||
end
|
||||
|
||||
""" Convert from COO format to CSC.
|
||||
|
||||
Parameters
|
||||
@@ -61,12 +57,6 @@ function empty!(A::SparseMatrixCOO)
|
||||
empty!(A.V)
|
||||
end
|
||||
|
||||
function append!(A::SparseMatrixCOO, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})
|
||||
append!(A.I, I)
|
||||
append!(A.J, J)
|
||||
append!(A.V, V)
|
||||
end
|
||||
|
||||
function append!(A::SparseMatrixCOO, B::SparseMatrixCOO)
|
||||
append!(A.I, B.I)
|
||||
append!(A.J, B.J)
|
||||
@@ -77,17 +67,6 @@ function isempty(A::SparseMatrixCOO)
|
||||
return isempty(A.I) && isempty(A.J) && isempty(A.V)
|
||||
end
|
||||
|
||||
function Base.:+(A::SparseMatrixCOO, B::SparseMatrixCOO)
|
||||
if isempty(A)
|
||||
return B
|
||||
end
|
||||
if isempty(B)
|
||||
return A
|
||||
end
|
||||
C = SparseMatrixCOO([A.I;B.I], [A.J;B.J], [A.V;B.V])
|
||||
return C
|
||||
end
|
||||
|
||||
function full(A::SparseMatrixCOO, args...)
|
||||
return full(sparse(A.I, A.J, A.V, args...))
|
||||
end
|
||||
@@ -156,6 +135,7 @@ function optimize!(A::SparseMatrixCOO)
|
||||
A.I = I
|
||||
A.J = J
|
||||
A.V = V
|
||||
return
|
||||
end
|
||||
|
||||
""" Find all nonzero rows from sparse matrix.
|
||||
@@ -173,20 +153,6 @@ function get_nonzero_columns(A::SparseMatrixCSC)
|
||||
return get_nonzero_rows(transpose(A))
|
||||
end
|
||||
|
||||
function get_nonzero_rows(A::Union{SparseMatrixCOO, Matrix})
|
||||
return get_nonzero_rows(sparse(A))
|
||||
end
|
||||
|
||||
function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix})
|
||||
return get_nonzero_columns(sparse(A))
|
||||
end
|
||||
|
||||
function get_nonzeros(C::Union{SparseMatrixCSC, Matrix})
|
||||
nz1 = get_nonzero_rows(C)
|
||||
nz2 = get_nonzero_columns(C)
|
||||
return (nz1, nz2)
|
||||
end
|
||||
|
||||
function size(A::SparseMatrixCOO)
|
||||
isempty(A) && return (0, 0)
|
||||
return maximum(A.I), maximum(A.J)
|
||||
@@ -206,20 +172,6 @@ function resize_sparsevec(b, n)
|
||||
return sparsevec(findnz(b)..., n)
|
||||
end
|
||||
|
||||
""" Matrix norm. Automatically convert to dense when asking for 2-norm for small matrices. """
|
||||
function norm(A::SparseMatrixCOO, p=Inf; maxdim=1000)
|
||||
dim = size(A, 1)
|
||||
if p == 2 && dim > maxdim
|
||||
warn("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.")
|
||||
return 0.0
|
||||
end
|
||||
if p == 2
|
||||
return norm(full(A), p)
|
||||
else
|
||||
return norm(sparse(A), p)
|
||||
end
|
||||
end
|
||||
|
||||
""" Approximative comparison of two matricse A and B. """
|
||||
function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
|
||||
A2 = sparse(A)
|
||||
|
||||
@@ -57,9 +57,3 @@ typealias IP Point{IntegrationPoint}
|
||||
function IP(id, weight, coords)
|
||||
return IP(id, weight, coords, Dict(), IntegrationPoint())
|
||||
end
|
||||
|
||||
function convert(::Type{IP}, data::Tuple{Float64, Vector{Float64}})
|
||||
weight, coords = data
|
||||
return IP(-1, weight, coords)
|
||||
end
|
||||
|
||||
|
||||
@@ -124,3 +124,10 @@ end
|
||||
@test isa(lst, Vector)
|
||||
end
|
||||
|
||||
@testset "extend basis" begin
|
||||
el = Element(Quad4, [1, 2, 3, 4])
|
||||
expected = [
|
||||
0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00
|
||||
0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25]
|
||||
@test isapprox(el([0.0, 0.0], 0.0, 2), expected)
|
||||
end
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
|
||||
@testset "NSeg interpolate" begin
|
||||
element = Element(NSeg, [1, 2])
|
||||
@test element([0.0], 0.0) == [0.5 0.5]
|
||||
@test size(element) == (1, 2)
|
||||
@test is_nurbs(element)
|
||||
element2 = Element(Seg2, [1, 2])
|
||||
@test !is_nurbs(element2)
|
||||
end
|
||||
|
||||
@testset "NSurf interpolate" begin
|
||||
element = Element(NSurf, [1, 2, 3, 4])
|
||||
@test element([0.0, 0.0], 0.0) == [0.25 0.25 0.25 0.25]
|
||||
@test size(element) == (2, 4)
|
||||
@test is_nurbs(element)
|
||||
end
|
||||
|
||||
@testset "NSolid interpolate" begin
|
||||
element = Element(NSolid, [1, 2, 3, 4, 5, 6, 7, 8])
|
||||
@test element([0.0, 0.0, 0.0], 0.0) == [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125]
|
||||
@test size(element) == (3, 8)
|
||||
@test is_nurbs(element)
|
||||
end
|
||||
+139
-17
@@ -3,13 +3,83 @@
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
using Logging
|
||||
Logging.configure(level=DEBUG)
|
||||
|
||||
@testset "create and manipulate fields" begin
|
||||
@testset "discrete, constant, time invariant field" begin
|
||||
@test isa(DCTI(), DCTI)
|
||||
@test DCTI(0.0).data == 0.0
|
||||
@test isa(Field(0.0), DCTI)
|
||||
@test isa(Field(), DCTI)
|
||||
f = DCTI()
|
||||
update!(f, 1.0)
|
||||
@test f.data == 1.0
|
||||
@test DCTI(1) == 1
|
||||
@test length(DCTI(1)) == 1
|
||||
@test f == DCTI(1.0)
|
||||
@test isapprox(f, DCTI(1.0))
|
||||
@test isapprox(f, 1.0)
|
||||
@test 2*f == 2.0 # multiply by constant
|
||||
@test f(1.0) == 1.0 # time interpolation
|
||||
@test isapprox([2.0]''*f, 2.0) # wanted behavior?
|
||||
end
|
||||
|
||||
@testset "updating time dependent fields" begin
|
||||
@testset "discrete, variable, time invariant field" begin
|
||||
@test isa(DVTI(), DVTI)
|
||||
@test DVTI([1.0, 2.0]).data == [1.0, 2.0]
|
||||
@test isa(Field([1.0, 2.0]), DVTI)
|
||||
|
||||
f = DVTI()
|
||||
update!(f, [2.0, 3.0])
|
||||
@test isapprox(f.data, [2.0, 3.0])
|
||||
@test length(f) == 2
|
||||
|
||||
# slicing
|
||||
@test isapprox(f[1], 2.0)
|
||||
@test isapprox(f[[1, 2]], [2.0, 3.0])
|
||||
|
||||
# boolean comparison and multiplying by a constant
|
||||
@test f == DVTI([2.0, 3.0])
|
||||
@test isapprox(2*f, [4.0, 6.0])
|
||||
|
||||
f3 = 2*f
|
||||
@test isa(f3, DVTI)
|
||||
@test f3+f == 3*f
|
||||
@test f3-f == f
|
||||
|
||||
# spatial interpolation
|
||||
N = [1.0, 2.0]
|
||||
@test isapprox(N*f, 8.0)
|
||||
|
||||
# time interpolation
|
||||
@test isapprox(f(1.0), [2.0, 3.0])
|
||||
|
||||
# spatial interpolation of vector valued variable field
|
||||
f2 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]])
|
||||
@test isapprox(f2[1], [1.0, 2.0])
|
||||
@test isapprox(f2[2], [3.0, 4.0])
|
||||
@test length(f2) == 2
|
||||
@test isapprox(N*f2, [1.0, 2.0] + [6.0, 8.0])
|
||||
|
||||
# iteration of DVTI field
|
||||
s = zeros(2)
|
||||
for j in f2
|
||||
s += j
|
||||
end
|
||||
@test isapprox(s, [4.0, 6.0])
|
||||
|
||||
@test vec(f2) == [1.0, 2.0, 3.0, 4.0]
|
||||
@test isapprox([1.0 2.0]*f, [8.0]'')
|
||||
|
||||
new_data = [2.0, 3.0, 4.0, 5.0]
|
||||
f4 = similar(f2, new_data)
|
||||
@test isa(f4, DVTI)
|
||||
@test isapprox(f4.data[1], [2.0, 3.0])
|
||||
@test isapprox(f4.data[2], [4.0, 5.0])
|
||||
end
|
||||
|
||||
@testset "discrete, constant, time-variant field" begin
|
||||
@test isa(DCTV(), DCTV)
|
||||
f = Field(0.0 => 1.0)
|
||||
@test isa(f, DCTV)
|
||||
@test last(f).time == 0.0
|
||||
@test last(f).data == 1.0
|
||||
update!(f, 0.0 => 2.0)
|
||||
@@ -20,27 +90,81 @@ Logging.configure(level=DEBUG)
|
||||
@test last(f).time == 1.0
|
||||
@test last(f).data == 3.0
|
||||
@test length(f) == 2
|
||||
|
||||
@testset "interpolation in time direction" begin
|
||||
@test isa(f(0.0), DCTI) # converts to time-invariant after time interpolation
|
||||
@test isapprox(f(-1.0), 2.0)
|
||||
@test isapprox(f(0.0), 2.0)
|
||||
@test isapprox(f(0.5), 2.5)
|
||||
@test isapprox(f(1.0), 3.0)
|
||||
@test isapprox(f(2.0), 3.0)
|
||||
end
|
||||
|
||||
# create several time steps at once
|
||||
f = DCTV(0.0 => 1.0, 1.0 => 2.0)
|
||||
@test isapprox(f(0.5), 1.5)
|
||||
|
||||
end
|
||||
|
||||
@testset "updating time invariant fields" begin
|
||||
f = Field(1.0)
|
||||
@test f.data == 1.0
|
||||
update!(f, 2.0)
|
||||
@test f.data == 2.0
|
||||
@testset "discrete, variable, time-variant field" begin
|
||||
@test isa(DVTV(), DVTV)
|
||||
f = Field(0.0 => [1.0, 2.0])
|
||||
@test isa(f, DVTV)
|
||||
@test last(f).time == 0.0
|
||||
@test last(f).data == [1.0, 2.0]
|
||||
update!(f, 0.0 => [2.0, 3.0])
|
||||
@test last(f).time == 0.0
|
||||
@test last(f).data == [2.0, 3.0]
|
||||
@test length(f) == 1
|
||||
update!(f, 1.0 => [3.0, 4.0])
|
||||
@test last(f).time == 1.0
|
||||
@test last(f).data == [3.0, 4.0]
|
||||
@test length(f) == 2
|
||||
|
||||
@testset "interpolation in time direction" begin
|
||||
@test isa(f(0.0), DVTI) # converts to time-invariant after time interpolation
|
||||
@test isapprox(f(-1.0), [2.0, 3.0])
|
||||
@test isapprox(f(0.0), [2.0, 3.0])
|
||||
@test isapprox(f(0.5), [2.5, 3.5])
|
||||
@test isapprox(f(1.0), [3.0, 4.0])
|
||||
@test isapprox(f(2.0), [3.0, 4.0])
|
||||
end
|
||||
|
||||
# create several time steps at once
|
||||
f = DVTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0])
|
||||
@test isapprox(f(0.5), [1.5, 2.5])
|
||||
end
|
||||
|
||||
@testset "field defined using function" begin
|
||||
g(xi, t) = xi[1]*t
|
||||
f = Field(g)
|
||||
v = f([1.0], 2.0)
|
||||
@test isapprox(v, 2.0)
|
||||
@testset "continuous, constant, time-invariant field" begin
|
||||
f = Field(() -> 2.0)
|
||||
@test isapprox(f([1.0], 2.0), 2.0)
|
||||
|
||||
end
|
||||
|
||||
@testset "continuous, constant, time variant field" begin
|
||||
f = Field((time::Float64) -> 2.0*time)
|
||||
@test isapprox(f([1.0], 2.0), 4.0)
|
||||
|
||||
end
|
||||
|
||||
@testset "continuous, variable, time invariant field" begin
|
||||
f = Field((xi::Vector) -> sum(xi))
|
||||
@test isapprox(f([1.0, 2.0], 2.0), 3.0)
|
||||
end
|
||||
|
||||
@testset "continuous, variable, time variant field" begin
|
||||
f = Field((xi::Vector, t::Float64) -> xi[1]*t)
|
||||
@test isapprox(f([1.0], 2.0), 2.0)
|
||||
end
|
||||
|
||||
@testset "unknown function argument for continuous field" begin
|
||||
@test_throws ErrorException Field((a, b, c) -> a*b*c)
|
||||
end
|
||||
|
||||
@testset "dictionary fields" begin
|
||||
f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0])
|
||||
f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0])
|
||||
f = Field(0.0 => f1, 1.0 => f2)
|
||||
debug("field = $f")
|
||||
@test isa(f, DVTV)
|
||||
@test isapprox(f(0.0)[1], [0.0, 0.0])
|
||||
@test isapprox(f(1.0)[2], [1.0, 1.0])
|
||||
@@ -58,5 +182,3 @@ end
|
||||
f = Field(f1)
|
||||
@test isa(f, DVTI)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
+10
-15
@@ -3,7 +3,6 @@
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
using JuliaFEM: description
|
||||
|
||||
ALL_ELEMENTS = [
|
||||
Seg2, Seg3,
|
||||
@@ -14,6 +13,16 @@ ALL_ELEMENTS = [
|
||||
Hex8, Hex20, Hex27
|
||||
]
|
||||
|
||||
info("basic data for elements implemented so far:")
|
||||
for element_type in [Poi1; ALL_ELEMENTS]
|
||||
element = Element(element_type, Int[])
|
||||
element_length = length(element)
|
||||
element_size = size(element)
|
||||
element_description = description(element)
|
||||
info("Element $element_type, description = $element_description, length = $element_length, size = $element_size")
|
||||
end
|
||||
|
||||
|
||||
ALL_ELEMENTS_NODES = [
|
||||
[1,2], [1,2,3],
|
||||
[1,2,3], [1,2,3,4,5,6], [1,2,3,4,5,6,7],
|
||||
@@ -78,17 +87,3 @@ end
|
||||
@test length(el) == length(vec)
|
||||
end
|
||||
end
|
||||
|
||||
DESC = ["2 node segment", "3 node segment", "3 node triangle",
|
||||
"6 node triangle", "7 node triangle", "4 node quadrangle",
|
||||
"8 node Serendip quadrangle", "9 node quadrangle",
|
||||
"4 node tetrahedral element", "10 node tetrahedral element",
|
||||
"6 node prismatic element (wedge)",
|
||||
"8 node hexahedral element", "20 node hexahedral element",
|
||||
"27 node hexahedral element"]
|
||||
|
||||
@testset "element description" begin
|
||||
for (T, res) in zip(ALL_ELEMENTS, DESC)
|
||||
@test description(Type(T)) == res
|
||||
end
|
||||
end
|
||||
|
||||
@@ -5,6 +5,7 @@ using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Testing
|
||||
using JuliaFEM.Abaqus: create_surface_elements
|
||||
|
||||
@testset "test that interface transfers constant field without error" begin
|
||||
meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med"
|
||||
@@ -47,3 +47,32 @@ using JuliaFEM.Testing
|
||||
info("Temperature at point X = $X is T = $T")
|
||||
@test isapprox(T, 100.0)
|
||||
end
|
||||
|
||||
@testset "problem not found from solver" begin
|
||||
s = Solver(Linear, "demo solver")
|
||||
@test_throws KeyError getindex(s, "not_found")
|
||||
end
|
||||
|
||||
@testset "automatic determination of problem dimension if not spesified" begin
|
||||
s = Solver(Linear, "demo solver")
|
||||
p = Problem(Elasticity, "demo problem", 2)
|
||||
push!(s, p)
|
||||
get_field_assembly(s)
|
||||
@test s.ndofs == 0
|
||||
add!(p.assembly.K, [4], [4], [4.0]'')
|
||||
get_field_assembly(s)
|
||||
@test s.ndofs == 4
|
||||
end
|
||||
|
||||
@testset "test for error when overdetermined system and requesting boundary assembly" begin
|
||||
s = Solver(Linear, "demo solver")
|
||||
@test_throws AssertionError get_boundary_assembly(s) # ndofs = 0
|
||||
p1 = Problem(Dirichlet, "bc1", 2, "displacement")
|
||||
p2 = Problem(Dirichlet, "bc2", 2, "displacement")
|
||||
# third dofs constrained
|
||||
add!(p1.assembly.C2, [3], [3], [1.0]'')
|
||||
add!(p2.assembly.C2, [3], [4], [1.0]'')
|
||||
s.ndofs = 4
|
||||
push!(s, p1, p2)
|
||||
@test_throws ErrorException get_boundary_assembly(s)
|
||||
end
|
||||
@@ -17,3 +17,31 @@ end
|
||||
add!(b, sparse(b2))
|
||||
@test isapprox(full(b), full(b2))
|
||||
end
|
||||
|
||||
@testset "Failure to add data to sparse vector due dimensino mismatch" begin
|
||||
b = SparseVectorCOO()
|
||||
@test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0])
|
||||
end
|
||||
|
||||
@testset "Test combining of SparseMatrixCOO" begin
|
||||
k = convert(Matrix{Float64}, reshape(collect(1:9), 3, 3))
|
||||
dofs1 = [1, 2, 3]
|
||||
dofs2 = [2, 3, 4]
|
||||
A = SparseMatrixCOO()
|
||||
add!(A, dofs1, dofs1, k)
|
||||
add!(A, dofs2, dofs2, k)
|
||||
A1 = full(A)
|
||||
optimize!(A)
|
||||
A2 = full(A)
|
||||
@test isapprox(A1, A2)
|
||||
end
|
||||
|
||||
@testset "resize of sparse matrix and sparse vector" begin
|
||||
A = sparse(rand(3, 3))
|
||||
B = resize_sparse(A, 4, 4)
|
||||
@test size(B) == (4, 4)
|
||||
a = sparse(rand(3))
|
||||
b = resize_sparsevec(a, 4)
|
||||
@test size(b) == (4, )
|
||||
end
|
||||
|
||||
|
||||
+1354
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user