diff --git a/src/fields.jl b/src/fields.jl index bffdd4f..e1d26e4 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -4,6 +4,7 @@ # https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb abstract Field + abstract DiscreteField <: Field abstract ContinuousField <: Field diff --git a/src/fields2.jl b/src/fields2.jl new file mode 100644 index 0000000..14bbfe5 --- /dev/null +++ b/src/fields2.jl @@ -0,0 +1,96 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb + +abstract AbstractField + +abstract Discrete <: AbstractField +abstract Continuous <: AbstractField +abstract Constant <: AbstractField +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 + +# Different field combinations +typealias DCTI Field{Discrete, Constant, TimeInvariant} +typealias DVTI Field{Discrete, Variable, TimeInvariant} +typealias DCTV Field{Discrete, Constant, TimeVariant} +typealias DVTV Field{Discrete, Variable, TimeVariant} +typealias CCTI Field{Continuous, Constant, TimeInvariant} +typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to interpolate in spatial dimension +typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time +typealias CVTV Field{Continuous, Variable, TimeVariant} + +# Basic data structure for discrete field +type Increment{T} + time :: Float64 + data :: T +end + +typealias VectorIncrement Increment{Vector} + +function Base.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 + +# Functions simplifying definition of fields. + +""" +All other data than vectors are considered as constant time invariant fields. +""" +function Field(data) + DCTI(data) +end + +""" +Vector data is considered as variable field time invariant field. +""" +function Field(data::Vector) + DVTI(data) +end + +""" +Data given in (time, value) pairs, where value is not vector, is considered as +constant time variant field. +""" +function Field{T}(data::Pair{Float64, T}...) + increments = [Increment{T}(d[1], d[2]) for d in data] + DCTV(increments) +end + +""" +Data given in (time, value) pairs, where value is a vector, is considered as +variable time variant field. +""" +function Field{T}(data::Pair{Float64, Vector{T}}...) + increments = [Increment{Vector{T}}(d[1], d[2]) for d in data] + DVTV(increments) +end + +""" Special case, constant time-variant vector, converted automatically. """ +function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...) + increments = [Increment(d[1], d[2]) for d in data] + DCTV(increments) +end + +## Other field related functions + +function Base.getindex(field::DVTV, i::Int64) + return field.data[i] +end + + +### FIELDSET ### + +typealias FieldSet Dict{ASCIIString, Field} diff --git a/src/solvers.jl b/src/solvers.jl index ac37c5e..f0cda11 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -86,8 +86,8 @@ function Base.push!(solver::Solver, problem::Problem) end """ Get all problems assigned to solver. """ -function get_problems(s::Solver) - return s.problems +function get_problems(solver::Solver) + return solver.problems end ## SimpleSolver -- tiny direct demo solver diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 528b9aa..cb99e6e 100644 --- a/test/test_elasticity.jl +++ b/test/test_elasticity.jl @@ -8,7 +8,6 @@ using JuliaFEM: Seg2, Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem - function test_elasticity_volume_load() element = Quad4([1, 2, 3, 4]) element["geometry"] = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]] @@ -19,10 +18,10 @@ function test_elasticity_volume_load() problem = PlaneStressElasticityProblem() push!(problem, element) solve!(problem, free_dofs; max_iterations=10) - disp = get_basis(element)("displacement", [1.0, 1.0])[2] + disp = get_basis(element)("displacement", [1.0, 1.0]) info("displacement at tip: $disp") # verified using Code Aster. - @test isapprox(disp, -8.77303119819776) + @test isapprox(disp[2], -8.77303119819776) end function test_elasticity_surface_load() diff --git a/test/test_solver.jl b/test/test_solver.jl index 732d415..567385f 100644 --- a/test/test_solver.jl +++ b/test/test_solver.jl @@ -6,7 +6,7 @@ module SolverTests using JuliaFEM.Test using JuliaFEM -using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis +using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis, MortarElement, MortarProblem, DirectSolver, PlaneStressElasticityProblem """ Define Problem 1: @@ -67,4 +67,115 @@ function test_simplesolver() @test isapprox(T, 100.0) end +function test_direct_solver() + + N = Dict{Int, Vector}( + 1 => [0.0, 0.0], + 2 => [2.0, 0.0], + 3 => [4.0, 0.0], + 4 => [0.0, 1.0], + 5 => [2.0, 1.0], + 6 => [4.0, 1.0], + 7 => [0.0, 1.0], + 8 => [1.0, 1.0], + 9 => [3.0, 1.0], + 10 => [4.0, 1.0], + 11 => [0.0, 2.0], + 12 => [1.0, 2.0], + 13 => [3.0, 2.0], + 13 => [4.0, 1.0]) + + # volume elements + e1 = Quad4([1, 2, 5, 4]) + e1["geometry"] = Vector[N[1], N[2], N[3], N[4]] + e2 = Quad4([2, 3, 6, 5]) + e2["geometry"] = Vector[N[2], N[3], N[6], N[5]] + e3 = Quad4([7, 8, 12, 11]) + e3["geometry"] = Vector[N[7], N[8], N[12], N[11]] + e4 = Quad4([8, 9, 13, 12]) + e4["geometry"] = Vector[N[8], N[9], N[13], N[12]] + e5 = Quad4([9, 10, 14, 13]) + e5["geometry"] = Vector[N[9], N[10], N[14], N[13]] + + # boundary elements for boundary load + b1 = Seg2([11, 12]) + b1["geometry"] = Vector[N[11], N[12]] + b1["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] + b2 = Seg2([12, 13]) + b2["geometry"] = Vector[N[12], N[13]] + b2["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] + b3 = Seg3([13, 14]) + b3["geometry"] = Vector[N[13], N[14]] + b3["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] + + # boundary elements for dirichlet dy=0 + d1 = Seg2([1, 2]) + d1["geometry"] = Vector[N[1], N[2]] + d1["displacement 2"] = 0.0 + d2 = Seg2([2, 3]) + d2["geometry"] = Vector[N[2], N[3]] + d2["displacement 2"] = 0.0 + + # boundary elements for dirichlet dx=0 + d3 = Seg2([1, 4]) + d3["geometry"] = Vector[N[1], N[4]] + d3["displacement 1"] = 0.0 + d4 = Seg2([4, 11]) + d4["geometry"] = Vector[N[4], N[11]] + d4["displacmeent 1"] = 0.0 + + # mortar elements to tie meshes -- masters + m1 = MSeg2([4, 5]) + m1["geometry"] = Vector[N[4], N[5]] + m2 = MSeg2([5, 6]) + m2["geometry"] = Vector[N[5], N[6]] + + # mortar elements to tie meshes -- slaves + rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] + phi = rotation_matrix(-pi/2) + m3 = MSeg2([7, 8]) + m3["geometry"] = Vector[N[7], N[8]] + m3["nodal ntsys"] = Matrix[phi, phi] + m3["master elements"] = MortarElement[m1, m2] + m4 = MSeg2([8, 9]) + m4["geometry"] = Vector[N[8], N[9]] + m4["nodal ntsys"] = Matrix[phi, phi] + m4["master elements"] = MortarElement[m1, m2] + m5 = MSeg2([9, 10]) + m5["geometry"] = Vector[N[9], N[10]] + m5["nodal ntsys"] = Matrix[phi, phi] + m5["master elements"] = MortarElement[m1, m2] + + problem1 = PlaneStressElasticityProblem() + push!(problem1, e1) + push!(problem1, e2) + push!(problem1, e3) + push!(problem1, e4) + push!(problem1, e5) + push!(problem1, b1) + push!(problem1, b2) + push!(problem1, b3) + + problem2 = DirichletProblem() + push!(problem2, d1) + push!(problem2, d2) + push!(problem2, d3) + push!(problem2, d4) + + problem3 = MortarProblem() + push!(problem3, m1) + push!(problem3, m2) + push!(problem3, m3) + push!(problem3, m4) + push!(problem3, m5) + + solver = DirectSolver() + push!(solver, problem1) + push!(solver, problem2) + push!(solver, problem3) + + call(solver) + +end + end