fields, solvers

2026-09-06 19:26:23 +00:00 · 2015-11-20 08:48:15 +02:00
parent 106f6f80b2
commit 18ae2ee5b7
5 changed files with 213 additions and 6 deletions
@@ -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

@@ -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}
@@ -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
@@ -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()
@@ -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