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

removed Equation type from code

Browse Source
This commit is contained in:
Jukka Aho
2015-11-27 10:10:00 +02:00
parent f5afcf2057
commit ef667e8f34
23 changed files with 730 additions and 1817 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+7 -5
View File
@@ -56,7 +56,7 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n)
end
# fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
include("fields2.jl")
include("fields.jl")
#include("basis.jl") # interpolation of discrete fields
include("symbolic.jl") # a thin symbolic layer for fields
include("types.jl") # type definitions
@@ -65,24 +65,26 @@ include("types.jl") # type definitions
include("elements.jl")
include("lagrange.jl") # Lagrange elements
#include("hierarchical.jl") # P-elements
include("mortar_elements.jl") # Mortar elements
#include("mortar_elements.jl") # Mortar elements
### EQUATIONS ###
include("integrate.jl") # default integration points for elements
include("sparse.jl")
include("equations.jl")
include("problems.jl")
include("equations.jl")
### FORMULATIION ###
include("dirichlet.jl")
include("mortar.jl") # mortar projection
include("heat.jl")
include("elasticity.jl")
### ASSEMBLY + SOLVE ###
include("assembly.jl")
include("solvers.jl")
include("directsolver.jl") # parallel sparse direct solver for non-lniear problems
include("directsolver.jl") # parallel sparse direct solver for non-linear problems
### MORTAR STUFF ###
#include("mortar.jl") # mortar projection
# PRE AND POSTPROCESS
include("xdmf.jl")
src/assembly.jl
+6 -6
View File
@@ -3,19 +3,19 @@
# Functions to handle global assembly of problem
function assemble!(assembly::Assembly, problem::Problem, time::Number=0.0, empty_assembly::Bool=true)
function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empty_assembly::Bool=true)
if empty_assembly
empty!(assembly)
end
for equation in get_equations(problem)
assemble!(assembly, equation, time, problem)
for element in get_elements(problem)
assemble!(assembly, problem, element, time)
end
end
function assemble(problem::Problem, time::Number=0.0)
function assemble(problem::AllProblems, time::Float64)
assembly = Assembly()
for equation in get_equations(problem)
assemble!(assembly, equation, time, problem)
for element in get_elements(problem)
assemble!(assembly, problem, element, time)
end
return assembly
end
src/basis.jl
-263
View File
@@ -1,266 +1,3 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
type Basis <: ContinuousField
basis :: Function
dbasisdxi :: Function
end
""" Evaluate basis. """
function Base.call(basis::Basis, xi::Vector, time::Number=0.0)
basis.basis(xi) # passing time does not make much sense actually for this...
end
""" Evaluate gradient of basis. This need geometry information to calculate Jacobian. """
function Base.call(basis::Basis, geometry::Increment, xi::Vector,
::Type{Val{:grad}})
dbasis = basis.dbasisdxi(xi)
J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
grad = inv(J)*dbasis
return grad
end
### INTERPOLATION IN SPATIAL DOMAIN ###
""" Interpolate increment in spatial domain using Basis. """
function Base.call(basis::Basis, increment::Increment, xi::Vector)
basis = basis.basis(xi)
sum([basis[i]*increment[i] for i=1:length(increment)])
end
""" Return gradient of increment in spatial domain using Basis.. """
function Base.call(basis::Basis, geometry::Increment, field::Increment,
xi::Vector, ::Type{Val{:grad}})
grad = basis(geometry, xi, Val{:grad})
gradf = sum([grad[:,i]*field[i]' for i=1:length(field)])'
return gradf
end
### INTERPOLATION IN TIME DOMAIN ###
""" Interpolate discrete field in time domain. Return Increment. """
function Base.call(field::DiscreteField, time::Number,
time_extrapolation::Symbol=:linear,
time_interpolation::Symbol=:linear)
# special cases, only 1 timestep defined or time = -Inf -> return first ts
if (length(field) == 1) || (time == -Inf)
#return field[1][end]
return first(field)
end
# special case, time = +Inf -> return last ts
if time == +Inf
#return field[end][end]
return last(field)
end
# very likely we are always near some defined timestep, usually field
# defined only on t = 0.0, test neighbourhood for timesteps
for i=1:length(field)
if isapprox(field[i].time, time)
return field[i][end]
end
end
# special case: out of time domain in positive direction, very likely
# to happen in incremental constitutive models
if time > field[end].time
if time_extrapolation == :constant
# constant time extrapolation, return last field
return field[end][end]
else
# multiple fields, pick last and second last and do linear interpolation
f1 = field[end-1]
f2 = field[end]
dt = abs(f2.time - f1.time)
i1 = f1[end]
i2 = f2[end]
di = i2 - i1
increment = Increment(i2 + di./dt * (time-f2.time))
return increment
end
end
# special case: out of time domain in negative direction
if time < field[1].time
if time_extrapolation == :constant
# constant time extrapolation, return first field
return field[1][end]
else
# multiple fields, pick first and second and do linear interpolation
f1 = field[1]
f2 = field[2]
dt = abs(f2.time - f1.time)
i1 = f1[end]
i2 = f2[end]
di = i2 - i1
increment = Increment(i1 - di./dt * (f1.time - time))
return increment
end
end
# find correct bin and perform interpolation
i = length(field)
while field[i].time >= time
i -= 1
end
if time_interpolation == :linear
t1 = field[i].time
t2 = field[i+1].time
inc1 = field[i][end]
inc2 = field[i+1][end]
dt = abs(t2 - t1)
t = (time-t1)/dt
increment = Increment((1-t)*inc1 + t*inc2)
return increment
end
if time_interpolation == :constant
# nearest neightbour interpolation, i.e. pick nearest defined field
t1 = field[i].time
t2 = field[i+1].time
dt1 = abs(t1-time)
dt2 = abs(t2-time)
if dt1 < dt2
return field[i][end]
else
return field[i+1][end]
end
end
end
""" Interpolate time derivative of field in some time t. This assumes linear
interpolation in time which is then differentiated.
Parameters
----------
field
Discrete field to interpolate. Must have timesteps and increments defined
time
Time to interpolate.
derivative
set Val{:diff} to activate this function
"""
function Base.call(field::DiscreteField, time::Number, ::Type{Val{:diff}})
# FieldSet -> Field -> TimeStep -> Increment -> data
if length(field) == 1
# just one timestep, time derivative cannot be evaluated.
error("Field length = $(length(field)), cannot evaluate time derivative")
end
function eval_field(i, j)
t1 = field[i]
t2 = field[j]
J = abs(t2.time - t1.time)
t = (time-t1.time)/J
result = 1/J*(-1*t1[end] + 1*t2[end])
return Increment(result)
end
# special cases, +Inf, -Inf, ~0.0
if (time > field[end].time) || isapprox(time, field[end].time)
return eval_field(endof(field)-1, endof(field))
end
if (time < field[1].time) || isapprox(time, field[1].time)
return eval_field(1, 2)
end
# search for a correct "bin" between time steps
i = length(field)
while (field[i].time > time) && !isapprox(field[i].time, time)
i -= 1
end
if isapprox(field[i].time, time)
# This is the hard case, maybe discontinuous time
# derivative if linear approximation.
# we are on the "mid node" in time axis
field1 = eval_field(i-1,i)
field2 = eval_field(i,i+1)
return 1/2*(field1 + field2)
end
return eval_field(i, i+1)
end
### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD
#=
""" Here we add field we are wanting to interpolate with ElementBasis. """
type ElementFieldBasis <: Basis
element_basis :: ElementBasis
field :: DiscreteField
time_extrapolation :: Symbol
time_interpolation :: Symbol
end
function Basis(basis::Function, dbasisdxi::Function, field::DiscreteField,
time_extrapolation=:linear, time_interpolation=:linear)
element_basis = ElementBasis(basis, dbasisdxi)
return ElementFieldBasis(element_basis, field, time_extrapolation,
time_interpolation)
end
function Base.call(basis::ElementFieldBasis, xi::Vector, time::Number)
increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation)
return basis.element_basis(increment, xi)
end
=#
### ELEMENT GRADIENT BASIS = ELEMENT BASIS + GEOMETRY
#=
""" Gradient of ElementBasis, needs geometry information. """
type ElementGradientBasis <: Basis
element_basis :: ElementBasis
geometry :: DiscreteField
time_extrapolation :: Symbol
time_interpolation :: Symbol
end
function grad(N::ElementBasis, f::ElementFieldBasis, X::ElementFieldBasis)
f.time_extrapolation == X.time_extrapolation || error("interpolation mismatch")
f.time_interpolation == X.time_interpolation || error("interpolation mismatch")
dN = ElementGradientBasis(N, X.field, f.time_extrapolation, f.time_interpolation)
dfdX = ElementFieldGradientBasis(dN, f.field, f.time_extrapolation, f.time_interpolation)
return dfdX
end
function grad(N::ElementBasis, f::DiscreteField, X::DiscreteField)
dN = ElementGradientBasis(N, X)
dfdX = ElementFieldGradientBasis(dN, f)
end
function ElementGradientBasis(element_basis::ElementBasis, geometry::DiscreteField)
return ElementGradientBasis(element_basis, geometry, :linear, :linear)
end
=#
### ELEMENT FIELD GRADIENT BASIS = ELEMENT GRADIENT BASIS + FIELD
#=
""" Gradient of ElementFieldBasis, needs field to interpolate. """
type ElementFieldGradientBasis <: Basis
element_gradient_basis :: ElementGradientBasis
field :: DiscreteField
time_extrapolation :: Symbol
time_interpolation :: Symbol
end
function ElementFieldGradientBasis(element_gradient_basis::ElementGradientBasis,
field::DiscreteField)
return ElementFieldGradientBasis(element_gradient_basis, field, :linear, :linear)
end
=#
### INTERPOLATION IN TIME DOMAIN ###
src/directsolver.jl
+30 -29
View File
@@ -4,7 +4,7 @@
## Direct solver
type DirectSolver <: Solver
field_problems :: Vector{FieldProblem}
field_problems :: Vector{Problem}
boundary_problems :: Vector{BoundaryProblem}
parallel :: Bool
nonlinear_problem :: Bool
@@ -12,7 +12,7 @@ type DirectSolver <: Solver
tol :: Float64
end
function push!(solver::DirectSolver, problem::FieldProblem)
function push!(solver::DirectSolver, problem::Problem)
push!(solver.field_problems, problem)
end
@@ -45,24 +45,30 @@ function call(solver::DirectSolver, time::Number=0.0)
# for this increment
for field_problem in solver.field_problems
for equation in get_equations(field_problem)
element = get_element(equation)
gdofs = get_gdofs(field_problem, equation)
if !isapprox(last(element[field_name]).time, time)
last_data = copy(last(element[field_name]).data)
push!(element[field_name], time => last_data)
for element in get_elements(field_problem)
gdofs = get_gdofs(element, field_dim)
if haskey(element, field_name)
if !isapprox(last(element[field_name]).time, time)
last_data = copy(last(element[field_name]).data)
push!(element[field_name], time => last_data)
end
else
data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
element[field_name] = (time => data)
end
end
end
for boundary_problem in solver.boundary_problems
for equation in get_equations(boundary_problem)
element = get_element(equation)
gdofs = get_gdofs(boundary_problem, equation)
eqdim = size(equation)[2]
data = Vector{Float64}[zeros(field_dim) for i in 1:eqdim]
if !isapprox(last(element["reaction force"]).time, time)
push!(element["reaction force"], time => data)
for element in get_elements(boundary_problem)
gdofs = get_gdofs(element, field_dim)
data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
if haskey(element, "reaction force")
if !isapprox(last(element["reaction force"]).time, time)
push!(element["reaction force"], time => data)
end
else
element["reaction force"] = (time => data)
end
end
end
@@ -105,32 +111,27 @@ function call(solver::DirectSolver, time::Number=0.0)
# update elements in field problems
for field_problem in solver.field_problems
for equation in get_equations(field_problem)
element = get_element(equation)
gdofs = get_gdofs(field_problem, equation)
eqsize = size(equation)
for element in get_elements(field_problem)
gdofs = get_gdofs(element, field_dim)
local_sol = vec(full(sol[gdofs])) # incremental data for element
local_sol = reshape(local_sol, eqsize)
local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)]
local_sol = reshape(local_sol, field_dim, length(element))
local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
last(element[field_name]).data += local_sol # <-- added
end
end
# update elements in boundary problems
for boundary_problem in solver.boundary_problems
for equation in get_equations(boundary_problem)
element = get_element(equation)
gdofs = get_gdofs(boundary_problem, equation) + dim
eqsize = size(equation)
for element in get_elements(boundary_problem)
gdofs = get_gdofs(element, field_dim) + dim
local_sol = vec(full(sol[gdofs]))
#info("local sol = $local_sol")
local_sol = reshape(local_sol, field_dim, eqsize[2])
local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)]
local_sol = reshape(local_sol, field_dim, length(element))
local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
last(element["reaction force"]).data = local_sol # <-- replaced
end
end
info("Iteration took $(toq()) seconds")
info("Non-linear iteration took $(toq()) seconds")
if norm(sol[1:dim]) < solver.tol
return (iter, true)
src/dirichlet.jl
+12 -53
View File
@@ -1,67 +1,26 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# Dirichlet boundary conditions in weak form
abstract DirichletProblem <: AbstractProblem
abstract DirichletEquation <: BoundaryEquation
### Dirichlet problem + equations
type DirichletProblem <: BoundaryProblem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Vector{DirichletEquation}
function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[])
return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements)
end
""" Initialize new Dirichlet boundary condition.
function assemble!{E}(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element{E}, time::Number)
Parameters
----------
dimension
dimension of unknown field (scalar, vector, ...)
# get dimension and name of PARENT field
field_dim = problem.parent_field_dim
field_name = problem.parent_field_name
Examples
--------
"""
function DirichletProblem(unknown_field_name::ASCIIString, dimension::Int=1)
DirichletProblem(unknown_field_name, dimension, [])
end
""" Dirichlet boundary condition element for 2 node line segment """
type DBC2D2 <: DirichletEquation
element :: Seg2
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::DBC2D2)
return (1, 2)
end
function Base.convert(::Type{DirichletEquation}, element::Seg2)
integration_points = get_integration_points(element, Val{3})
if !haskey(element, "reaction force")
element["reaction force"] = (0.0 => Vector{Float64}[])
end
DBC2D2(element, integration_points)
end
function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing)
isa(problem, Void) && error("Dicihlet boundary condition needs problem defined")
field_dim = problem.unknown_field_dimension
field_name = problem.unknown_field_name
element = get_element(equation)
gdofs = get_gdofs(element, field_dim)
basis = get_basis(element)
detJ = det(basis)
for ip in get_integration_points(equation)
w = ip.weight * detJ(ip)
N = basis(ip, time)
for ip in get_integration_points(element)
w = ip.weight * det(element, ip, time)
N = element(ip, time)
A = w*N'*N
if haskey(element, field_name)
# add all dimensions at once
# add all dimensions at once if defined element["blaa"] = 0.0
for i=1:field_dim
g = element(field_name, ip, time)
ldofs = gdofs[i:field_dim:end]
@@ -71,8 +30,8 @@ function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number
end
for i=1:field_dim
# add per dof if defined element["blaa 1"] = 1.0, element["blaa 2"] = 0.0 etc.
if haskey(element, field_name*" $i")
# add single component
g = element(field_name*" $i", ip, time)
ldofs = gdofs[i:field_dim:end]
add!(assembly.stiffness_matrix, ldofs, ldofs, A)
src/elasticity.jl
+36 -109
View File
@@ -3,42 +3,20 @@
# Elasticity problems
abstract ElasticityProblem <: FieldProblem
abstract ElasticityEquation <: Equation
abstract ElasticityProblem <: AbstractProblem
function get_unknown_field_name(equation::ElasticityEquation)
abstract PlaneStressElasticityProblem <: ElasticityProblem
function PlaneStressElasticityProblem(dim::Int=2, elements=[])
return Problem{PlaneStressElasticityProblem}(dim, elements)
end
function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
return "displacement"
end
### Formulation ###
""" Calculate internal energy for elasticity equation.
Override this to define your own material model. By default we use
Saint Venant-Kirchhoff material model, which is simply
S(E) = λtr(E) + 2μE
"""
function get_internal_energy(equation::ElasticityEquation, ip::IntegrationPoint, time::Number, F::Matrix)
element = get_element(equation)
basis = get_basis(element)
dbasis = grad(basis)
# material parameters
young = basis("youngs modulus", ip, time)
poisson = basis("poissons ratio", ip, time)
mu = young/(2*(1+poisson))
lambda = young*poisson/((1+poisson)*(1-2*poisson))
if isa(equation, PlaneStressElasticityEquation)
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d
end
# material model
E = 1/2*(F'*F - I) # strain
S = lambda*trace(E)*I + 2*mu*E
P = F*S
return P*dbasis(ip, time)
function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
return Vector{Float64}
end
""" Elasticity equations.
@@ -71,97 +49,44 @@ https://en.wikipedia.org/wiki/Plane_stress
https://en.wikipedia.org/wiki/Hooke's_law
"""
function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint, time::Number; variation=nothing)
function get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing)
element = get_element(equation)
basis = get_basis(element)
dbasis = grad(basis)
basis = element(ip, time)
dbasis = element(ip, time, Val{:grad})
u = basis("displacement", ip, time, variation)
gradu = dbasis("displacement", ip, time, variation)
u = element("displacement", ip, time, variation)
gradu = element("displacement", ip, time, Val{:grad}, variation)
F = I + gradu # deformation gradient
#info("Deformation gradient: $F")
# residual vector - internal energy
r = get_internal_energy(equation, ip, time, F)
#info("boundary element")
# internal forces
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))
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d
E = 1/2*(F'*F - I) # strain
S = lambda*trace(E)*I + 2*mu*E
r = F*S*dbasis
# external forces - volume load
if haskey(element, "displacement load")
b = basis("displacement load", ip, time)
r -= b*basis(ip, time)
b = element("displacement load", ip, time)
r -= b*basis
end
return vec(r)
end
""" Surface load for plane stress model. """
function get_residual_vector(problem::Problem{PlaneStressElasticityProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing)
### Plane stress elasticity ###
abstract PlaneElasticityProblem <: ElasticityProblem
abstract PlaneStressElasticityEquation <: ElasticityEquation
type PlaneStressElasticityProblem <: PlaneElasticityProblem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Vector{PlaneStressElasticityEquation}
end
function PlaneStressElasticityProblem(equations=[])
return PlaneStressElasticityProblem("displacement", 2, equations)
end
### Equations ###
""" 4-node plane stress element. """
type CPS4 <: PlaneStressElasticityEquation
element :: Quad4
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::CPS4)
return (2, 4)
end
function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
integration_points = get_integration_points(element)
if !haskey(element, "displacement")
element["displacement"] = 0.0 => [zeros(2) for i=1:4]
end
CPS4(element, integration_points)
end
""" Boundary element for plane stress problem for surface loads. """
type CPS2 <: PlaneStressElasticityEquation
element :: Seg2
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::CPS2)
return (2, 2)
end
function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
integration_points = get_integration_points(element)
if !haskey(element, "displacement")
element["displacement"] = 0.0 => [zeros(2) for i=1:2]
end
CPS2(element, integration_points)
end
function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; variation=nothing)
element = get_element(equation)
basis = get_basis(element)
u = basis("displacement", ip, time, variation)
r = zeros(size(equation))
u = element("displacement", ip, time, variation)
r = zeros(problem.dim, length(element))
if haskey(element, "displacement traction force")
T = basis("displacement traction force", ip, time)
r -= T*basis(ip, time)
T = element("displacement traction force", ip, time)
r -= T*element(ip, time)
end
return vec(r)
@@ -237,3 +162,5 @@ function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number;
end
=#
src/elements.jl
+51 -120
View File
@@ -1,10 +1,23 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using FactCheck
using ForwardDiff
abstract AbstractElement
abstract Element
type Element{E<:AbstractElement}
connectivity :: Vector{Int}
# integration_points :: Vector{IntegrationPoint}
fields :: Dict{ASCIIString, Field}
end
function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
# return Element{E}(connectivity, get_integration_points(E), Dict())
return Element{E}(connectivity, Dict())
end
function get_integration_points{E}(element::Element{E})
# return element.integration_points
return get_integration_points(E)
end
"""
Test routine for element. If this passes, element interface is properly
@@ -67,9 +80,9 @@ function Base.getindex(element::Element, field_name)
return element.fields[field_name]
end
#function get_integration_points(element)
# return get_default_integration_points(element)
#end
function Base.length{E}(element::Element{E})
size(E)[2]
end
"""Add new Field to element.
@@ -89,151 +102,69 @@ function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
element.fields[name] = Field(data...)
end
#function Base.setindex!(element::Element, field_data::Tuple, field_name)
# field = Field()
# for (time, data) in field_data
# ts = TimeStep(time, Increment[Increment(data)])
# push!(field, ts)
# end
# element[field_name] = field
#end
function get_connectivity(el::Element)
return el.connectivity
end
abstract AbstractFunctionSpace
type FunctionSpace <: AbstractFunctionSpace
basis :: CVTI
fields :: FieldSet
end
type GradientFunctionSpace <: AbstractFunctionSpace
basis :: CVTI
fields :: FieldSet
end
function get_basis(element::Element)
return FunctionSpace(element.basis, element.fields)
end
function get_dbasis(element::Element)
return GradientFunctionSpace(element.basis, element.fields)
end
function grad(u::FunctionSpace)
return GradientFunctionSpace(u.basis, u.fields)
end
""" If basis is called without a field, return basis functions evaluated at that point. """
function call(u::FunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
return u.basis(xi)
end
""" If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """
function call(gradu::GradientFunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
geometry = gradu.fields["geometry"](t)
gradu.basis(geometry, xi, Val{:grad})
end
""" Evaluate field on element function space. """
function call(u::FunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
field = !isa(variation, Void) ? variation : u.fields[field_name](t)
u.basis(field, xi)
end
""" Evaluate gradient of field on element function space. """
function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
field = !isa(variation, Void) ? variation : gradu.fields[field_name](t)
geometry = gradu.fields["geometry"](t)
gradu.basis(geometry, field, xi, Val{:grad})
end
typealias VecOrIP Union{Vector, IntegrationPoint}
function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number)
return element.basis(element[field_name](time), xi)
function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, variation=nothing)
field = isa(variation, Void) ? element[field_name](time) : variation
basis = get_basis(element)
return basis(field, xi)
end
function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}})
return element.basis(element["geometry"](time), element[field_name](time), xi, Val{:grad})
function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}}, variation=nothing)
field = isa(variation, Void) ? element[field_name](time) : variation
basis = get_basis(element)
geom = element["geometry"](time)
return basis(geom, field, xi, Val{:grad})
end
function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP)
function call(element::Element, field_name::ASCIIString, xi::VecOrIP)
return element.basis(element[field_name], xi)
end
function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
function call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
return element.basis(element["geometry"], element[field_name], xi, Val{:grad})
end
function Base.call(element::Element, field_name::ASCIIString, time::Number)
function call(element::Element, field_name::ASCIIString, time::Number)
return element[field_name](time)
end
function Base.call(element::Element, xi::VecOrIP)
element.basis(xi)
function get_basis{E}(element::Element{E}, ip::IntegrationPoint)
return get_basis(E, ip.xi)
end
function Base.call(element::Element, xi::VecOrIP, ::Type{Val{:grad}})
element.basis(element["geometry"], xi, Val{:grad})
function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0)
return get_basis(element, xi)
end
function Base.call(element::Element, field_name::ASCIIString)
function get_basis{E}(element::Element{E})
basis = CVTI(
(xi::Vector) -> get_basis(E, xi),
(xi::Vector) -> get_dbasis(E, xi))
return basis
end
function call{E}(element::Element{E}, xi::VecOrIP, time::Float64, ::Type{Val{:grad}})
basis = get_basis(element)
return basis(element["geometry"], xi, Val{:grad})
end
function call(element::Element, field_name::ASCIIString)
return element[field_name]
end
# on-line functions to get api more easy to use, ip -> xi.ip
#call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
#call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
# i think these will be the most called functions.
#call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
#call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
#call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
#call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
""" Return a field from function space. """
function get_field(u::FunctionSpace, field_name, time::Number=0.0)
return u.fields[field_name](time)
end
""" Return a field from function space. """
function get_field(u::FunctionSpace, field_name, time::Number=0.0, variation=nothing)
return !isa(variation, Void) ? variation : u.fields[field_name](time)
end
""" Return a field from function space. """
function get_fieldset(u::FunctionSpace, field_name)
return u.fields[field_name]
end
""" Get a determinant of element in point ξ. """
function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
X = u.fields["geometry"](time)
dN = u.basis(xi, Val{:grad})
function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Number=0.0)
X = element("geometry", time)
dN = get_dbasis(E, ip.xi)
J = sum([dN[:,i]*X[i]' for i=1:length(X)])
m, n = size(J)
return m == n ? det(J) : norm(J)
end
function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, time::Number=0.0)
LinAlg.det(u, ip.xi, time)
end
function LinAlg.det(u::FunctionSpace)
return (args...) -> det(u, args...)
end
function LinAlg.det(element::Element)
return det(get_basis(element))
end
#Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
#Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
#Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
#Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
""" Check does field exist. """
function Base.haskey(element::Element, what)
haskey(element.fields, what)
src/equations.jl
+31 -64
View File
@@ -3,10 +3,6 @@
# Functions to handle element level things -- integration, assembly, ...
abstract Equation
abstract FieldEquation <: Equation
abstract BoundaryEquation <: Equation
type Assembly
mass_matrix :: SparseMatrixIJV
stiffness_matrix :: SparseMatrixIJV
@@ -41,58 +37,31 @@ end
function get_residual_vector
end
function has_mass_matrix(equation::Equation)
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
function has_mass_matrix(problem::Problem, element::Element)
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
return method_exists(get_mass_matrix, default_args)
end
function has_stiffness_matrix(equation::Equation)
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
function has_stiffness_matrix(problem::Problem, element::Element)
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
return method_exists(get_stiffness_matrix, default_args)
end
function has_force_vector(equation::Equation)
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
function has_force_vector(problem::Problem, element::Element)
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
return method_exists(get_force_vector, default_args)
end
function has_potential_energy(equation::Equation)
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
function has_potential_energy(problem::Problem, element::Element)
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
return method_exists(get_potential_energy, default_args)
end
function has_residual_vector(equation::Equation)
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
function has_residual_vector(problem::Problem, element::Element)
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
return method_exists(get_residual_vector, default_args)
end
function get_element(equation::Equation)
return equation.element
end
function get_integration_points(equation::Equation)
return equation.integration_points
end
function Base.size(equation::Equation, i::Int)
return size(equation)[i]
end
""" Return global degrees of freedom of element in matrix level.
Notes
-----
This is calculated from connectivity and equation dimension.
"""
function get_gdofs(equation::Equation)
element = get_element(equation)
conn = get_connectivity(element)
dim = size(equation, 1)
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
return gdofs
end
function get_gdofs(element::Element, dim::Int)
conn = get_connectivity(element)
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
@@ -100,26 +69,23 @@ function get_gdofs(element::Element, dim::Int)
end
""" Assemble element. """
function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, problem=nothing)
function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number)
element = get_element(equation)
gdofs = get_gdofs(equation)
basis = get_basis(element)
detJ = det(basis)
unknown_field_name = get_unknown_field_name(equation)
gdofs = get_gdofs(element, problem.dim)
unknown_field_name = get_unknown_field_name(problem)
# 1. if equations are defined we just integrate them, without caring how they are done
if has_mass_matrix(equation) || has_stiffness_matrix(equation) || has_force_vector(equation)
for ip in get_integration_points(equation)
s = ip.weight*detJ(ip)
if has_mass_matrix(equation)
add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(equation, ip, time))
if has_mass_matrix(problem, element) || has_stiffness_matrix(problem, element) || has_force_vector(problem, element)
for ip in get_integration_points(element)
s = ip.weight*det(element, ip, time)
if has_mass_matrix(element)
add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(problem, element, ip, time))
end
if has_stiffness_matrix(equation)
add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(equation, ip, time))
if has_stiffness_matrix(element)
add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(problem, element, ip, time))
end
if has_force_vector(equation)
add!(assembly.force_vector, gdofs, s*get_force_vector(equation, ip, time))
if has_force_vector(element)
add!(assembly.force_vector, gdofs, s*get_force_vector(problem, element, ip, time))
end
end
# external loads -- if any nodal loads is defined add to force vector
@@ -129,7 +95,7 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
end
# 2. energy form -- user has defined potential energy W -> min!
if has_potential_energy(equation)
if has_potential_energy(problem, element)
field = element[unknown_field_name](time)
""" Wrapper for potential energy for ForwardDiff. """
@@ -137,9 +103,9 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
W = 0.0
df = similar(field, data)
# integrate potential energy
for ip in get_integration_points(equation)
s = ip.weight*detJ(ip)
dw = get_potential_energy(equation, ip, time; variation=df)
for ip in get_integration_points(element)
s = ip.weight*det(element, ip, time)
dw = get_potential_energy(problem, element, ip, time; variation=df)
W += s*dw
end
# external energy -- if any nodal loads is defined, decrease from potential energy
@@ -156,7 +122,7 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
end
# 3. virtual work -- user has defined some residual r = p - f = 0
if has_residual_vector(equation)
if has_residual_vector(problem, element)
field = DVTI(last(element[unknown_field_name]).data)
@@ -165,9 +131,9 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
R = zeros(length(data))
df = similar(field, data)
# integrate residual vector
for ip in get_integration_points(equation)
s = ip.weight*detJ(ip)
dr = get_residual_vector(equation, ip, time; variation=df)
for ip in get_integration_points(element)
s = ip.weight*det(element, ip, time)
dr = get_residual_vector(problem, element, ip, time; variation=df)
R += s*dr
end
# external loads -- if any nodal loads is defined, decrease from residual
@@ -183,3 +149,4 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))
end
end
src/fields.jl
+219 -210
View File
@@ -3,301 +3,310 @@
# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
abstract Field
abstract AbstractField
abstract DiscreteField <: Field
abstract ContinuousField <: Field
abstract Discrete <: AbstractField
abstract Continuous <: AbstractField
abstract Constant <: AbstractField
abstract Variable <: AbstractField
abstract TimeVariant <: AbstractField
abstract TimeInvariant <: AbstractField
### DEFAULT DISCRETE FIELD ###
# 1. Increment
type Increment{T} <: AbstractVector{T}
data :: Vector{T}
type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
data
end
function Base.size(increment::Increment)
return size(increment.data)
### Basic data structure for discrete field
type Increment{T}
time :: Float64
data :: T
end
function Base.linearindexing(::Type{Increment})
return LinearFast()
function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
return Increment{T}(data[1], data[2])
end
function Base.getindex(increment::Increment, i::Int)
function Base.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 Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
return increment.data[i]
end
function Base.setindex!(increment::Increment, v, i::Int)
increment.data[i] = v
function Base.(:*)(d, increment::Increment)
return d*increment.data
end
function Base.dot(k::Number, increment::Increment)
return k*increment
### Basic data structure for continuous field
type Basis
basis :: Function
dbasis :: Function
end
function Base.convert(::Type{Increment}, data::Number)
return Increment([data])
function Base.call(basis::Basis, xi::Vector)
basis.basis(xi)
end
function Base.convert{T}(::Type{Increment}, data::Array{T, 2})
return Increment([data[:,i] for i=1:size(data, 2)])
function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
basis.dbasis(xi)
end
function Base.convert{T}(::Type{Increment}, data::Array{T, 3})
return Increment([data[:,:,i] for i=1:size(data, 3)])
end
### Different field combinations and other typealiases
function Base.convert{T}(::Type{Increment}, data::Array{T, 4})
return Increment([data[:,:,:,i] for i=1:size(data, 4)])
end
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}
function Base.convert{T}(::Type{Increment}, data::Array{T, 5})
return Increment([data[:,:,:,:,i] for i=1:size(data, 5)])
end
typealias ScalarIncrement{T} Increment{T}
typealias VectorIncrement{T} Increment{Vector{T}}
typealias TensorIncrement{T} Increment{Matrix{T}}
function Base.zeros(::Type{Increment}, T, dims...)
return Increment(zeros(T, dims...))
end
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}
""" Flatten increment to Vector.
Examples
--------
### Convenient functions to create fields
>>> inc = ones(Increment, 2, 4)
>>> vec(inc)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
"""
function Base.vec(increment::Increment)
return [increment...;]
end
function Base.similar{T}(increment::Increment, data::Vector{T})
return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
end
function Base.convert{T}(::Type{Vector{T}}, increment::Increment)
return Increment[increment]
end
# 2. TimeStep
type TimeStep
time :: Float64
increments :: Vector{Increment}
end
function Base.size(timestep::TimeStep)
return size(timestep.increments)
end
function Base.endof(timestep::TimeStep)
return endof(timestep.increments)
end
function Base.length(timestep::TimeStep)
return length(timestep.increments)
end
function Base.linearindexing(::Type{TimeStep})
return Base.LinearFast()
end
function Base.getindex(timestep::TimeStep, i::Int)
return timestep.increments[i]
end
#function TimeStep(data::Union{Number, Array}...)
# return TimeStep(0.0, Increment[Increment(d) for d in data])
#function Base.convert(::Type{Field}, data)
# return Field(data)
#end
function TimeStep()
return TimeStep(0.0, [])
function Field(data)
return DCTI(data)
end
function TimeStep{T}(data::T...)
return TimeStep(0.0, Increment[Increment(d) for d in data])
function Field(data::Vector)
return DVTI(data)
end
function Base.convert(::Type{TimeStep}, value::Number)
return TimeStep(0.0, Increment[Increment(value)])
function Field{T}(data::Pair{Float64, T}...)
return DCTV([Increment{T}(d[1], d[2]) for d in data])
end
function Base.push!(timestep::TimeStep, increment::Increment)
push!(timestep.increments, increment)
function Field{T}(data::Pair{Float64, Vector{T}}...)
return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
end
# FIXME: having some serious problems here to get tuple form working.
# 3. DefaultDiscreteField
immutable DefaultDiscreteField <: DiscreteField
timesteps :: Vector{TimeStep}
#=
function DefaultDiscreteField(data::Array)
if (typeof(data) == Vector{Int64}) || (typeof(data) == Vector{Float64})
new(TimeStep[TimeStep(data)])
else
new(data)
end
end
=#
function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
end
#=
type DefaultDiscreteField <: DiscreteField
timesteps :: Vector{TimeStep}
function DefaultDiscreteField(data...)
timesteps = TimeStep[]
for (i, d) in enumerate(data)
@debug("i = $i, d = $d")
if isa(d, Tuple)
# contains time vector
increments = Increment[Increment(d[2])]
push!(timesteps, TimeStep(d[1], increments))
else
increments = Increment[Increment(d)]
push!(timesteps, TimeStep(i-1.0, increments))
end
end
new(timesteps)
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 Base.size(field::DefaultDiscreteField)
return size(field.timesteps)
function CVTI(basis::Function, dbasis::Function)
return CVTI(Basis(basis, dbasis))
end
function Base.length(field::DefaultDiscreteField)
return length(field.timesteps)
function Field(basis::Function, dbasis::Function)
return CVTI(basis, dbasis)
end
function Base.start(::DefaultDiscreteField)
### Accessing and manipulating discrete fields
function Base.getindex(field::DVTV, i::Int64)
return field.data[i]
end
function Base.push!(field::DCTV, data::Pair)
push!(field.data, data)
end
function Base.push!(field::DVTV, data::Pair)
# info("field.data = \n$(field.data)")
# info("data = \n$data")
push!(field.data, data)
end
function Base.getindex(field::DVTV, i::Int64)
return field.data[i]
end
function Base.getindex(field::DVTI, i::Int64)
return field.data[i]
end
function Base.getindex(field::DCTV, i::Int64)
return field.data[i]
end
function Base.getindex(field::Field, i::Int64)
return field.data[i]
end
function Base.length(field::DVTI)
return length(field.data)
end
function Base.length(field::DCTI)
return 1
end
function Base.next(field::DefaultDiscreteField, state)
return (field[state+1], state+1)
function Base.length(field::DVTV)
return length(field.data)
end
function Base.done(field::DefaultDiscreteField, state)
return state > length(field)
function Base.length(field::DCTV)
return length(field.data)
end
function eltype(::Type{DefaultDiscreteField})
return TimeStep
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) = ($op)(field.data, k)
@eval ($op)(k, field::DCTI) = ($op)(field.data, k)
end
function Base.linearindexing(::Type{DefaultDiscreteField})
return LinearFast()
function Base.vec(field::DVTI)
return [field.data...;]
end
function Base.getindex(field::DefaultDiscreteField, i::Int)
return field.timesteps[i]
function Base.vec(field::DCTV)
info("trying to vectorize $field")
error("does not make sense")
end
function Base.endof(field::DefaultDiscreteField)
return endof(field.timesteps)
function Base.endof(field::Field)
return endof(field.data)
end
function Base.first(field::DefaultDiscreteField)
return field[1][end]
#function Base.similar{T}(field::DVTI, data::Vector{T})
# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
#end
function Base.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 Base.last(field::DefaultDiscreteField)
return field[end][end]
function Base.start(::DVTI)
return 1
end
function Base.push!(field::DefaultDiscreteField, timestep::TimeStep)
push!(field.timesteps, timestep)
function Base.next(f::DVTI, state)
return f.data[state], state+1
end
function Base.push!(field::DefaultDiscreteField, data::Union{Vector, Matrix})
push!(field[end], Increment(data))
function Base.done(f::DVTI, s)
return s > length(f.data)
end
function Base.push!(field::DefaultDiscreteField, data::Pair)
ts = TimeStep(data[1], Increment(data[2]))
push!(field, ts)
### Accessing continuous fields
function Base.call(field::CVTI, xi::Vector)
field.data(xi)
end
"""Quickly create fields.
function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}})
field.data(xi, Val{:grad})
end
Examples
--------
>>> Field([1, 2]) # creates field with one timestep and vector value [1, 2]
>>> Field(1, 2) # creates field with two timesteps, each having scalar value
>>> Field([1, 2], [3, 4]) # creates field with two timesteps, each having vector value
>>> Field( (0.0, [1, 2]), (0.5, [3, 4]) ) # like above, but give time also
"""
function Base.convert(::Type{DefaultDiscreteField}, data...)
timesteps = TimeStep[]
for (i, d) in enumerate(data)
if isa(d, Tuple)
@debug("is tuple, has time, d = $d")
# contains time vector
increments = Increment[Increment(d[2])]
push!(timesteps, TimeStep(d[1], increments))
else
# @debug("array without time, d = $d")
# @debug(typeof(d))
increments = Increment[Increment(d)]
push!(timesteps, TimeStep(i-1.0, increments))
function Base.convert(::Type{Basis}, field::CVTI)
return field.data
end
function Base.call(field::CCTV, time::Number)
return field.data(time)
end
### Interpolation
""" Interpolate time-invariant field in time direction. """
function Base.call(field::DVTI, time::Float64)
return field
end
function Base.call(field::DCTI, time::Float64)
return field
end
function Base.call(field::CVTI, time::Float64)
return field.data()
end
function Base.call(field::CCTI, time::Float64)
return field.data()
end
""" Interpolate time-variant field in time direction. """
function Base.call(field::DCTV, time::Float64)
for i=reverse(1:length(field))
if isapprox(field[i].time, time)
return DCTI(field[i].data)
end
end
field = DefaultDiscreteField(timesteps)
return field
info(field.data)
info(time)
error("interpolate DCTV: not implemented yet")
end
function Base.convert(::Type{DefaultDiscreteField}, data::Vector{TimeStep})
field = DefaultDiscreteField(data)
# @debug(field)
return field
function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
for i=reverse(1:length(field))
if isapprox(field[i].time, time)
return DVTI(field[i].data)
end
end
info(field.data)
info(time)
error("interpolate DVTV: not implemented yet")
end
### CONTINUOUS FIELDS ###
type DefaultContinuousField <: ContinuousField
field :: Function
""" Interpolate constant field in spatial dimension. """
function Base.call(basis::CVTI, field::DCTI, xi::Vector)
return field.data
end
function Base.call(field::DefaultContinuousField, xi::Vector, time::Number)
return field.field(xi, time)
""" Interpolate variable field in spatial dimension. """
function Base.call(basis::CVTI, values::DVTI, xi::Vector)
N = basis(xi)
return sum([N[i]*values[i] for i=1:length(N)])
end
function Base.convert(::Type{DefaultContinuousField}, f::Function)
return DefaultContinuousField(f)
function Base.call(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)])
invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
grad = invJ * dbasis
return grad
end
function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
grad = call(basis, geometry, xi, Val{:grad})
gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
return length(gradf) == 1 ? gradf[1] : gradf
end
function Base.call(basis::CVTI, xi::Vector, time::Number)
call(basis, xi)
end
### FIELDSET ###
typealias FieldSet Dict{ASCIIString, Field}
# 1. given numbers, arrays or tuples -> discrete field
function Base.convert(::Type{Field}, data::Union{Number, Array, Tuple}...)
return DiscreteField(data...)
end
function Base.convert(::Type{DiscreteField}, data::Union{Number, Array, Tuple}...)
return convert(DefaultDiscreteField, data...)
end
# 2. given function -> continuous field
function Base.convert(::Type{Field}, data::Function)
return ContinuousField(data)
end
function Base.convert(::Type{ContinuousField}, data::Function)
return convert(DefaultContinuousField, data)
end
function Base.length(::Field)
return 1
end
src/fields2.jl
-306
View File
@@ -3,310 +3,4 @@
# 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
### Basic data structure for discrete field
type Increment{T}
time :: Float64
data :: T
end
function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
return Increment{T}(data[1], data[2])
end
function Base.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 Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
return increment.data[i]
end
function Base.(:*)(d, increment::Increment)
return d*increment.data
end
### Basic data structure for continuous field
type Basis
basis :: Function
dbasis :: Function
end
function Base.call(basis::Basis, xi::Vector)
basis.basis(xi)
end
function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
basis.dbasis(xi)
end
### Different field combinations and other typealiases
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}
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}
### Convenient functions to create fields
#function Base.convert(::Type{Field}, data)
# return Field(data)
#end
function Field(data)
return DCTI(data)
end
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
function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
return DCTV([Increment{Vector{T}}(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 Base.getindex(field::DVTV, i::Int64)
return field.data[i]
end
function Base.push!(field::DCTV, data::Pair)
push!(field.data, data)
end
function Base.push!(field::DVTV, data::Pair)
# info("field.data = \n$(field.data)")
# info("data = \n$data")
push!(field.data, data)
end
function Base.getindex(field::DVTV, i::Int64)
return field.data[i]
end
function Base.getindex(field::DVTI, i::Int64)
return field.data[i]
end
function Base.getindex(field::DCTV, i::Int64)
return field.data[i]
end
function Base.getindex(field::Field, i::Int64)
return field.data[i]
end
function Base.length(field::DVTI)
return length(field.data)
end
function Base.length(field::DCTI)
return 1
end
function Base.length(field::DVTV)
return length(field.data)
end
function Base.length(field::DCTV)
return length(field.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) = ($op)(field.data, k)
@eval ($op)(k, field::DCTI) = ($op)(field.data, k)
end
function Base.vec(field::DVTI)
return [field.data...;]
end
function Base.vec(field::DCTV)
info("trying to vectorize $field")
error("does not make sense")
end
function Base.endof(field::Field)
return endof(field.data)
end
#function Base.similar{T}(field::DVTI, data::Vector{T})
# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
#end
function Base.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 Base.start(::DVTI)
return 1
end
function Base.next(f::DVTI, state)
return f.data[state], state+1
end
function Base.done(f::DVTI, s)
return s > length(f.data)
end
### Accessing continuous fields
function Base.call(field::CVTI, xi::Vector)
field.data(xi)
end
function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}})
field.data(xi, Val{:grad})
end
function Base.convert(::Type{Basis}, field::CVTI)
return field.data
end
function Base.call(field::CCTV, time::Number)
return field.data(time)
end
### Interpolation
""" Interpolate time-invariant field in time direction. """
function Base.call(field::DVTI, time::Float64)
return field
end
function Base.call(field::DCTI, time::Float64)
return field
end
function Base.call(field::CVTI, time::Float64)
return field.data()
end
function Base.call(field::CCTI, time::Float64)
return field.data()
end
""" Interpolate time-variant field in time direction. """
function Base.call(field::DCTV, time::Float64)
for i=reverse(1:length(field))
if isapprox(field[i].time, time)
return DCTI(field[i].data)
end
end
info(field.data)
info(time)
error("interpolate DCTV: not implemented yet")
end
function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
for i=reverse(1:length(field))
if isapprox(field[i].time, time)
return DVTI(field[i].data)
end
end
info(field.data)
info(time)
error("interpolate DVTV: not implemented yet")
end
""" Interpolate constant field in spatial dimension. """
function Base.call(basis::CVTI, field::DCTI, xi::Vector)
return field.data
end
""" Interpolate variable field in spatial dimension. """
function Base.call(basis::CVTI, values::DVTI, xi::Vector)
N = basis(xi)
return sum([N[i]*values[i] for i=1:length(N)])
end
function Base.call(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)])
invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
grad = invJ * dbasis
return grad
end
function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
grad = call(basis, geometry, xi, Val{:grad})
gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
return length(gradf) == 1 ? gradf[1] : gradf
end
function Base.call(basis::CVTI, xi::Vector, time::Number)
call(basis, xi)
end
### FIELDSET ###
typealias FieldSet Dict{ASCIIString, Field}
src/heat.jl
+20 -68
View File
@@ -3,14 +3,20 @@
# Heat problems
abstract HeatProblem <: Problem
abstract HeatEquation <: Equation
abstract HeatProblem <: AbstractProblem
function get_unknown_field_name(equation::HeatEquation)
function HeatProblem(dim::Int=1, elements=[])
return Problem{HeatProblem}(dim, elements)
end
function get_unknown_field_name{P<:HeatProblem}(::Type{P})
return "temperature"
end
### Formulation ###
function get_unknown_field_type{P<:HeatProblem}(::Type{P})
# scalar field
return Float64
end
""" Heat equations.
@@ -36,83 +42,29 @@ References
https://en.wikipedia.org/wiki/Heat_equation
"""
function assemble!(assembly::Assembly, equation::HeatEquation, time::Number=0.0, problem=nothing)
function assemble!{E<:CG}(assembly::Assembly, problem::Problem{HeatProblem}, element::Element{E}, time::Number)
element = get_element(equation)
gdofs = get_gdofs(equation)
basis = get_basis(element)
dbasis = grad(basis)
detJ = det(basis)
for ip in get_integration_points(equation)
w = ip.weight*detJ(ip)
N = basis(ip, time)
gdofs = get_gdofs(element, problem.dim)
for ip in get_integration_points(element)
w = ip.weight*det(element, ip, time)
N = element(ip, time)
if haskey(element, "density")
rho = basis("density", ip, time)
rho = element("density", ip, time)
add!(assembly.mass_matrix, gdofs, gdofs, w*rho*N'*N)
end
if haskey(element, "temperature thermal conductivity")
dN = dbasis(ip, time)
k = basis("temperature thermal conductivity", ip, time)
dN = element(ip, time, Val{:grad})
k = element("temperature thermal conductivity", ip, time)
add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN)
end
if haskey(element, "temperature load")
f = basis("temperature load", ip, time)
f = element("temperature load", ip, time)
add!(assembly.force_vector, gdofs, w*N'*f)
end
if haskey(element, "temperature flux")
info("assemble boundary flux")
g = basis("temperature flux", ip, time)
g = element("temperature flux", ip, time)
add!(assembly.force_vector, gdofs, w*N'*g)
end
end
end
### Equations ###
""" Diffusive heat transfer for 4-node bilinear element. """
type DC2D4 <: HeatEquation
element :: Quad4
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::DC2D4)
return (1, 4)
end
""" Diffusive heat transfer for 2-node linear segment. """
type DC2D2 <: HeatEquation
element :: Seg2
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::DC2D2)
return (1, 2)
end
# Conversions element -> equation
function Base.convert(::Type{HeatEquation}, element::Quad4)
integration_points = get_integration_points(element)
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
DC2D4(element, integration_points)
end
function Base.convert(::Type{HeatEquation}, element::Seg2)
integration_points = get_integration_points(element)
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
DC2D2(element, integration_points)
end
### Problems ###
type PlaneHeatProblem <: HeatProblem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Vector{HeatEquation}
end
""" Default constructor for problem takes no arguments. """
function PlaneHeatProblem(equations=[])
return PlaneHeatProblem("temperature", 1, equations)
end
src/integrate.jl
+6 -6
View File
@@ -3,7 +3,7 @@
# Let's drop here all integration schemes and some defaults for different element types
function get_integration_points(Quad4::Element)
function get_integration_points(::Type{Quad4})
[
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
@@ -14,13 +14,13 @@ end
typealias LineElement Union{Seg2, Seg3}
function get_integration_points(element::LineElement, ::Type{Val{1}})
function get_integration_points(::Type{Seg2}, ::Type{Val{1}})
[
IntegrationPoint([0.0], 2.0)
]
end
function get_integration_points(element::LineElement, ::Type{Val{2}})
function get_integration_points(::Type{Seg2}, ::Type{Val{2}})
[
IntegrationPoint([-sqrt(1/3)], 1)
IntegrationPoint([+sqrt(1/3)], 1)
@@ -54,12 +54,12 @@ function get_integration_points(element::LineElement, ::Type{Val{5}})
]
end
function get_integration_points(element::Seg2)
return get_integration_points(element, Val{1})
function get_integration_points(::Type{Seg2})
return get_integration_points(Seg2, Val{2})
end
function get_integration_points(element::Seg3)
return get_integration_points(element, Val{2})
return get_integration_points(element, Val{3})
end
### 3D elements
src/lagrange.jl
+14 -17
View File
@@ -3,7 +3,7 @@
# Lagrange (Continous Galerkin) finite elements
abstract CG <: Element
abstract CG <: AbstractElement
"""
Given polynomial P and coordinates of reference element, calculate
@@ -29,29 +29,26 @@ Examples
>>> @create_lagrange_element(Seg2, "2 node linear segment", X, P)
"""
macro create_lagrange_element(element_name, element_description, X, P)
# Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
eltype = esc(element_name)
quote
global get_element_description
basis, dbasisdxi = calculate_lagrange_basis($P, $X)
type $eltype <: CG
connectivity :: Array{Int, 1}
basis :: CVTI
fields :: FieldSet
global get_basis, get_dbasis
basis, dbasis = calculate_lagrange_basis($P, $X)
abstract $eltype <: CG
function get_basis(::Type{$eltype}, xi::Vector{Float64})
return basis(xi)
end
function $eltype(connectivity, args...)
$eltype(connectivity, CVTI(basis, dbasisdxi), FieldSet())
function get_dbasis(::Type{$eltype}, xi::Vector{Float64})
return dbasis(xi)
end
function $eltype(args...)
return Element{$eltype}(args...)
end
function Base.size(::Type{$eltype})
return Base.size($X)
end
get_element_description(el::Type{$eltype}) = $element_description
Base.size(element::Type{$eltype}) = Base.size($X)
Base.size(element::$eltype) = Base.size($X)
end
end
# 0d Lagrange element
#@create_element(Point1, CG, "1 node point element")
# 1d Lagrange elements
@create_lagrange_element(Seg2, "2 node linear line element",
src/mortar_elements.jl
+1 -1
View File
@@ -6,7 +6,7 @@
# slave element = non-mortar element where integration happens
# master element = mortar element projected to non-mortar side
abstract MortarElement <: Element
abstract MortarElement <: AbstractElement
type MSeg2 <: MortarElement
connectivity :: Vector{Int}
src/problems.jl
+22 -63
View File
@@ -1,77 +1,36 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
abstract Problem
abstract BoundaryProblem <: Problem
abstract FieldProblem <: Problem
abstract AbstractProblem
""" Return all equations beloging to this problem. """
function get_equations(problem::Problem)
problem.equations
type Problem{T<:AbstractProblem}
dim :: Int
elements :: Vector{Element}
end
type BoundaryProblem{T<:AbstractProblem}
parent_field_name :: ASCIIString
parent_field_dim :: Int
dim :: Int
elements :: Vector{Element}
end
typealias AllProblems Union{Problem, BoundaryProblem}
function get_elements(problem::AllProblems)
return problem.elements
end
""" Return the dimension of the unknown field of this problem. """
function get_unknown_field_dimension(problem::Problem)
problem.unknown_field_dimension
return problem.dim
end
""" Return the name of the unknown field of this problem. """
function get_unknown_field_name(problem::Problem)
problem.unknown_field_name
function get_unknown_field_name{P<:AbstractProblem}(problem::Problem{P})
return get_unknown_field_name(P)
end
"""
Add new equation to problem.
Parameters
----------
problem
element
Notes
-----
Equation is automatically created during process based on problem
element -> equation mapping and element type.
"""
function Base.push!(problem::Problem, element::Element, args...)
# element_type = typeof(element)
# equation_type = problem.element_mapping[element_type]
# push!(problem.equations, equation_type(element, args...))
push!(problem.equations, element)
function Base.push!(problem::AllProblems, element::Element)
push!(problem.elements, element)
end
"""
Return the size of the problem, i.e.
maximum number of connectivity × unknown field dimension.
"""
function Base.size(problem::Problem)
mc = 0
for equation in get_equations(problem)
element = get_element(equation)
mc = max(mc, get_connectivity(element)...)
end
dim = get_unknown_field_dimension(problem)
return (dim, dim*mc)
end
""" Assign new equation mapping to problem for some element.
Examples
--------
>>> p = PlaneHeatProblem()
>>> p[Seg2] = DC2D2
"""
function Base.setindex!(problem::Problem, equation, element)
problem.element_mapping[element] = equation
end
""" Return global degrees of freedom of element in matrix level. """
function get_gdofs(problem::Problem, equation::Equation)
dim = get_unknown_field_dimension(problem)
element = get_element(equation)
conn = get_connectivity(element)
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
return gdofs
end
src/solvers.jl
+43 -99
View File
@@ -5,39 +5,6 @@
abstract Solver
"""
Solve field equations for single element with some dofs fixed. This can be used
to test nonlinear element formulations.
"""
function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false, callback=nothing)
unknown_field_name = get_unknown_field_name(equation)
element = get_element(equation)
x0 = element[unknown_field_name](0.0)
x = zeros(prod(size(equation)))
dx = fill!(similar(x), 0.0)
ass = Assembly()
for i=1:max_iterations
empty!(ass)
assemble!(ass, equation)
A = full(ass.stiffness_matrix)[free_dofs, free_dofs]
b = full(ass.force_vector)[free_dofs]
if dump_matrices
dump(full(A))
dump(full(b)')
end
dx[free_dofs] = A \ b
x += dx
eqsize = size(equation)
data = eqsize[1] != 1 ? reshape(x, eqsize) : x
push!(element[unknown_field_name], time => data)
norm(dx) < tolerance && return
if !isa(callback, Void)
callback(x)
end
end
error("Did not converge in $max_iterations iterations")
end
"""
Solve field equations for a single problem with some dofs fixed. This can be used
to test nonlinear element formulations. Dirichlet boundary is assumed to be homogeneous
@@ -72,13 +39,12 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
if !(isa(callback, Void))
callback(x)
end
for equation in get_equations(problem)
element = get_element(equation)
gdofs = get_gdofs(equation)
for element in get_elements(problem)
gdofs = get_gdofs(element, problem.dim)
data = full(x[gdofs])
eqsize = size(equation)
if eqsize[1] != 1
data = reshape(data, eqsize)
if length(data) != length(element)
data = reshape(data, problem.dim, length(element))
data = [data[:,i] for i=1:size(data,2)]
end
push!(element[field_name], time => data)
end
@@ -87,24 +53,11 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
error("Did not converge in $max_iterations iterations")
end
function Base.push!(solver::Solver, problem::Problem)
push!(solver.problems, problem)
end
""" Get all problems assigned to solver. """
function get_problems(solver::Solver)
return solver.problems
end
## SimpleSolver -- tiny direct demo solver
""" Simple solver for educational purposes. """
type SimpleSolver <: Solver
problems :: Vector{Problem}
end
""" Default initializer. """
function SimpleSolver()
SimpleSolver([])
""" Simple linear solver for educational purposes. """
type LinearSolver <: Solver
field_problem :: Problem
boundary_problem :: BoundaryProblem
end
"""
@@ -113,28 +66,30 @@ Call solver to solve a set of problems.
This is a simple direct solver for demonstration purposes. It handles the
common situation, i.e., some main field problem and it's Dirichlet boundary.
Au + C'λ = f
Ku + C'λ = f
Cu = g
"""
function call(solver::SimpleSolver, time::Number=0.0)
problem1, problem2 = get_problems(solver)
function call(solver::LinearSolver, time::Float64)
assembly1 = Assembly()
assemble!(assembly1, problem1, time)
assembly2 = Assembly()
assemble!(assembly2, problem2, time)
field_name = get_unknown_field_name(solver.field_problem)
field_dim = get_unknown_field_dimension(solver.field_problem)
info("solving $field_name problem, $field_dim dofs / nodes")
field_assembly = assemble(solver.field_problem, time)
boundary_assembly = assemble(solver.boundary_problem, time)
info("Creating sparse matrices")
K = sparse(field_assembly.stiffness_matrix)
dim = size(K, 1)
f = sparse(field_assembly.force_vector, dim, 1)
C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
g = sparse(boundary_assembly.force_vector, dim, 1)
# info("Creating sparse matrices")
A1 = sparse(assembly1.stiffness_matrix)
dims = size(A1)
b1 = sparse(assembly1.force_vector, dims[1], 1)
A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
b2 = sparse(assembly2.force_vector, dims[1], 1)
# create a saddle point problem
A = [A1 A2; A2' zeros(A2)]
b = [b1; b2]
A = [K C'; C' zeros(C)]
b = [f; g]
# solve problem
nz = unique(rowvals(A)) # take only non-zero rows
@@ -142,37 +97,26 @@ function call(solver::SimpleSolver, time::Number=0.0)
x[nz] = lufact(A[nz,nz]) \ full(b[nz])
# get "problem-wise" solution vectors
x1 = x[1:length(b1)]
x2 = x[length(b1)+1:end]
u = x[1:dim]
la = x[dim+1:end]
# update field for elements in problem 1
for equation in get_equations(problem1)
element = get_element(equation)
field_name = get_unknown_field_name(problem1)
gdofs = get_gdofs(problem1, equation)
local_sol = vec(full(x1[gdofs]))
eqsize = size(equation)
if eqsize[1] != 1
local_sol = reshape(local_sol, eqsize)
for element in get_elements(solver.field_problem)
gdofs = get_gdofs(element, field_dim)
local_sol = vec(full(u[gdofs]))
# if solving vector field, modify local solution vector
# to array of vectors
if field_dim != 1
local_sol = reshape(local_sol, field_dim, length(element))
local_sol = [local_sol[:,i] for i=1:size(local_sol,2)]
end
if haskey(element, field_name)
push!(element[field_name], time => local_sol)
else
element[field_name] = (time => local_sol)
end
#info("problem1: pushing to $field_name")
push!(element[field_name], time => local_sol)
end
# update field for elements in problem 2 (Dirichlet boundary)
for equation in get_equations(problem2)
element = get_element(equation)
field_name = "reaction force" #get_unknown_field_name(problem2)
gdofs = get_gdofs(problem2, equation)
local_sol = vec(full(x1[gdofs]))
eqsize = size(equation)
if eqsize[1] != 1
local_sol = reshape(local_sol, eqsize)
end
#info("problem2: pushing to $field_name")
#push!(element[field_name], time => local_sol)
end
return norm(x1)
return norm(u)
end
test/test_directsolver.jl
+132
View File
@@ -0,0 +1,132 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
module DirectSolverTests
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: Seg2, Quad4
using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem
using JuliaFEM: DirectSolver
function test_solver_multiple_dirichlet_bc()
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
e1 = Quad4([1, 2, 4, 3])
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
e1["youngs modulus"] = 900.0
e1["poissons ratio"] = 0.25
b1 = Seg2([3, 4])
b1["geometry"] = Vector[N[3], N[4]]
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
problem = PlaneStressElasticityProblem()
push!(problem, e1)
push!(problem, b1)
# manually solve problem 1
# free_dofs = [3, 5, 6, 8]
# free_dofs = [3, 6, 7, 8]
#solve!(problem, free_dofs, 0.0; max_iterations=10)
#disp = e1("displacement", [1.0, 1.0], 0.0)
#info("displacement at tip: $disp")
#@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
# boundary elements for dirichlet dx=0
dx = Seg2([1, 3])
dx["geometry"] = Vector[N[1], N[3]]
dx["displacement 1"] = 0.0
# boundary elements for dirichlet dy=0
dy = Seg2([1, 2])
dy["geometry"] = Vector[N[1], N[2]]
dy["displacement 2"] = 0.0
problem2 = DirichletProblem("displacement", 2)
push!(problem2, dx)
problem3 = DirichletProblem("displacement", 2)
push!(problem3, dy)
solver = DirectSolver()
push!(solver, problem)
push!(solver, problem2)
push!(solver, problem3)
# launch solver
norm = solver(0.0)
disp = e1("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
end
function test_solver_multiple_bodies_multiple_dirichlet_bc()
N = Vector[
[0.0, 0.0], [1.0, 0.0],
[0.0, 1.0], [1.0, 1.0],
[0.0, 2.0], [1.0, 2.0]]
e1 = Quad4([1, 2, 4, 3])
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
e2 = Quad4([3, 4, 6, 5])
e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
for el in [e1, e2]
el["youngs modulus"] = 900.0
el["poissons ratio"] = 0.25
end
b1 = Seg2([5, 6])
b1["geometry"] = Vector[N[5], N[6]]
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
body1 = PlaneStressElasticityProblem()
push!(body1, e1)
body2 = PlaneStressElasticityProblem()
push!(body2, e2)
push!(body2, b1)
# boundary elements for dirichlet dx=0
dx1 = Seg2([1, 3])
dx1["geometry"] = Vector[N[1], N[3]]
dx2 = Seg2([3, 5])
dx2["geometry"] = Vector[N[3], N[5]]
for dx in [dx1, dx2]
dx["displacement 1"] = 0.0
end
boundary1 = DirichletProblem("displacement", 2)
push!(boundary1, dx1)
push!(boundary1, dx2)
# boundary elements for dirichlet dy=0
dy1 = Seg2([1, 2])
dy1["geometry"] = Vector[N[1], N[2]]
dy1["displacement 2"] = 0.0
boundary2 = DirichletProblem("displacement", 2)
push!(boundary2, dy1)
solver = DirectSolver()
push!(solver, body1)
push!(solver, body2)
push!(solver, boundary1)
push!(solver, boundary2)
# launch solver
norm = solver(0.0)
disp = e2("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# code aster verification, two_elements.comm
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
end
# test_solver_multiple_bodies_multiple_dirichlet_bc()
end
test/test_dirichlet.jl
+4 -7
View File
@@ -5,7 +5,7 @@ module TestDirichletBoundaryCondition
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble!
using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble
function test_dirichlet_problem_1_dim()
element = Seg2([1, 2])
@@ -13,8 +13,7 @@ function test_dirichlet_problem_1_dim()
element["temperature"] = 0.0
problem = DirichletProblem("temperature", 1)
push!(problem, element)
assembly = Assembly()
assemble!(assembly, problem)
assembly = assemble(problem, 0.0)
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
@test isapprox(A, 1/6*[2 1; 1 2])
@@ -27,8 +26,7 @@ function test_dirichlet_problem_2_dim()
element["displacement"] = 0.0
problem = DirichletProblem("displacement", 2)
push!(problem, element)
assembly = Assembly()
assemble!(assembly, problem)
assembly = assemble(problem, 0.0)
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
A_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
@@ -42,8 +40,7 @@ function test_dirichlet_problem_2_dim_single_dof_fixed()
element["displacement 2"] = 0.0
problem = DirichletProblem("displacement", 2)
push!(problem, element)
assembly = Assembly()
assemble!(assembly, problem)
assembly = assemble(problem, 0.0)
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
info(b)
test/test_elasticity.jl
+8 -7
View File
@@ -4,9 +4,7 @@
module ElasticityTests
using JuliaFEM.Test
using JuliaFEM: Seg2, Quad4, Field, FieldSet, CPS4,
get_basis, solve!,
PlaneStressElasticityProblem
using JuliaFEM: Seg2, Quad4, PlaneStressElasticityProblem, solve!
function test_elasticity_volume_load()
element = Quad4([1, 2, 3, 4])
@@ -14,11 +12,13 @@ function test_elasticity_volume_load()
element["youngs modulus"] = 500.0
element["poissons ratio"] = 0.3
element["displacement load"] = Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]
free_dofs = [3, 4, 5, 6]
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])
problem = PlaneStressElasticityProblem()
push!(problem, element)
free_dofs = [3, 4, 5, 6]
solve!(problem, free_dofs, 0.0; max_iterations=10)
disp = get_basis(element)("displacement", [1.0, 1.0], 0.0)
disp = element("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# verified using Code Aster.
@test isapprox(disp[2], -8.77303119819776)
@@ -31,9 +31,12 @@ function test_elasticity_surface_load()
element1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
element1["youngs modulus"] = 900.0
element1["poissons ratio"] = 0.25
element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])
element2 = Seg2([3, 4])
element2["geometry"] = Vector[N[3], N[4]]
element2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
#free_dofs = [3, 5, 6, 8]
free_dofs = [3, 6, 7, 8]
@@ -41,8 +44,6 @@ function test_elasticity_surface_load()
push!(problem, element1)
push!(problem, element2)
solve!(problem, free_dofs, 0.0; max_iterations=10)
#disp = get_basis(element1)("displacement", [1.0, 1.0], 1.0)[2]
info(last(element1["displacement"]))
disp = element1("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# verified using Code Aster.
test/test_heat.jl
+10 -12
View File
@@ -7,8 +7,7 @@ module HeatTests # always wrap tests to module ending with "Tests"
using JuliaFEM.Test # always use JuliaFEM.Test, not Base.Test
using JuliaFEM: HeatEquation
using JuliaFEM: Seg2, Quad4, DC2D4, DC2D2, Assembly, assemble!
using JuliaFEM: Seg2, Quad4, HeatProblem, assemble
function test_one_element() # always start test function with name test_
@@ -26,11 +25,13 @@ function test_one_element() # always start test function with name test_
# linear ramp from 0 to 6 in time 0 to 1
boundary_element["temperature flux"] = (0.0 => 0.0, 1.0 => 6.0)
problem = HeatProblem()
push!(problem, element)
push!(problem, boundary_element)
# Set constant source f=12 with k=6. Accurate solution is
# T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
equation = convert(HeatEquation, element)
assembly = Assembly()
assemble!(assembly, equation)
assembly = assemble(problem, 0.0)
fdofs = [1, 2]
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
@@ -38,18 +39,15 @@ function test_one_element() # always start test function with name test_
# Set constant flux g=6 on boundary. Accurate solution is
# u(x,y) = x which equals T=1 on boundary.
boundary_equation = convert(HeatEquation, boundary_element)
empty!(assembly)
time = 1.0
assemble!(assembly, equation, time)
assemble!(assembly, boundary_equation, time)
# at time t=1.0 all loads should be on.
assembly = assemble(problem, 1.0)
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
T = A[fdofs, fdofs] \ b[fdofs]
info("T = $T")
@test isapprox(T, [2.0, 2.0]) # always use @test to test things.
@test isapprox(T, [2.0, 2.0])
end
end
test/test_potential_energy.jl
+31 -69
View File
@@ -3,73 +3,43 @@
module ElementTests
using JuliaFEM
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: Equation, Quad4, IntegrationPoint, Assembly, assemble!,
get_element, get_basis, grad, get_unknown_field_name,
PlaneHeatProblem, Seg2, Problem, solve!,
get_default_integration_points, Equation
using JuliaFEM: AbstractProblem, Problem
using JuliaFEM: Element, Seg2, Quad4
using JuliaFEM: IntegrationPoint, solve!
abstract MyEquation <: Equation
abstract HeatProblem <: AbstractProblem
function JuliaFEM.get_unknown_field_name(equation::MyEquation)
function HeatProblem(dim::Int=1, elements=[])
return Problem{HeatProblem}(dim, elements)
end
function JuliaFEM.get_unknown_field_name{P<:HeatProblem}(::Type{P})
return "temperature"
end
""" Diffusive heat transfer for 4-node bilinear element, with a nonlinear source term. """
type DC2D4NL <: MyEquation
element :: Quad4
integration_points :: Vector{IntegrationPoint}
function JuliaFEM.get_unknown_field_type{P<:HeatProblem}(::Type{P})
return Float64
end
function Base.size(equation::DC2D4NL)
return (1, 4)
end
""" Nonlinear flux term. """
type DC2D2NL <: MyEquation
element :: Seg2
integration_points :: Vector{IntegrationPoint}
end
function Base.size(equation::DC2D2NL)
return (1, 2)
end
function Base.convert(::Type{MyEquation}, element::Quad4)
integration_points = get_default_integration_points(element)
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
DC2D4NL(element, integration_points)
end
function Base.convert(::Type{MyEquation}, element::Seg2)
integration_points = JuliaFEM.line5()
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
DC2D2NL(element, integration_points)
end
""" Calculate a potential Π = Wint - Wext of system. """
function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing)
element = get_element(equation)
basis = get_basis(element)
k = basis("temperature thermal conductivity", ip, time)
f = basis("temperature load", ip, time)
T = basis("temperature", ip, time, variation)
c = basis("temperature nonlinearity coefficient", ip, time)
gradT = grad(basis)("temperature", ip, time, variation)
function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Quad4}, ip::IntegrationPoint, time::Number; variation=nothing)
k = element("temperature thermal conductivity", ip, time)
f = element("temperature load", ip, time)
T = element("temperature", ip, time, variation)
c = element("temperature nonlinearity coefficient", ip, time)
gradT = element("temperature", ip, time, Val{:grad}, variation)
Wint = (k + c*T) * 1/2*vecdot(gradT, gradT)
Wext = f*T
return Wint - Wext
end
function JuliaFEM.get_potential_energy(equation::DC2D2NL, ip, time; variation=nothing)
element = get_element(equation)
basis = get_basis(element)
T = basis("temperature", ip, time, variation)[1]
T_ext = basis("temperature external", ip, time)[1]
coeff = basis("temperature coefficient", ip, time)[1]
function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing)
T = element("temperature", ip, time, variation)[1]
T_ext = element("temperature external", ip, time)[1]
coeff = element("temperature coefficient", ip, time)[1]
q0 = coeff*(T_ext^4 - T^4)
Wint = 0.0
Wext = q0*T
@@ -86,28 +56,19 @@ function test_potential_energy_method()
element["temperature load"] = [0.0, 0.0, 0.0, 0.0]
element["temperature nodal load"] = [3.0, 3.0, 0.0, 0.0]
element["temperature nonlinearity coefficient"] = 6.0
equation = convert(MyEquation, element)
element["temperature"] = (0.0 => zeros(Float64, 4))
problem = HeatProblem()
push!(problem, element)
# create model -- end
solve!(equation, [1, 2], 0.0)
basis = get_basis(element)
temp = basis("temperature", [0.0, -1.0], 0.0)
solve!(problem, [1, 2], 0.0)
temp = element("temperature", [0.0, -1.0], 0.0)
err = temp - 2/3
info("error: $err")
@test isapprox(err, 0.0)
end
type TestProblem <: Problem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Vector{MyEquation}
end
function TestProblem(equations=[])
TestProblem("temperature", 1, equations)
end
function test_potential_energy_method_2()
# create model -- start
@@ -117,20 +78,21 @@ function test_potential_energy_method_2()
element1["temperature thermal conductivity"] = 6.0
element1["temperature load"] = [0.0, 0.0, 0.0, 0.0]
element1["temperature nonlinearity coefficient"] = [0.0, 0.0, 0.0, 0.0]
element1["temperature"] = (0.0 => zeros(Float64, 4))
element2 = Seg2([1, 2])
element2["geometry"] = Vector[N[1], N[2]]
element2["temperature coefficient"] = 3.0e-8 # ~ 5.7e-8 * 0.5
element2["temperature external"] = 100.0
element2["temperature"] = (0.0 => zeros(Float64, 2))
# create model -- end
problem = TestProblem()
problem = HeatProblem()
push!(problem, element1)
push!(problem, element2)
solve!(problem, [1, 2], 0.0)
basis = get_basis(element1)
temp = basis("temperature", [0.0, -1.0], 0.0)
temp = element1("temperature", [0.0, -1.0], 0.0)
err = temp - 0.5
info("error: $err")
@test isapprox(err, 0.0, atol=1.0e-6)
test/test_solver.jl
+17 -267
View File
@@ -6,14 +6,11 @@ module SolverTests
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis, MortarElement, MortarProblem, PlaneStressElasticityProblem, solve!, DirectSolver
using JuliaFEM: Seg2, Quad4
using JuliaFEM: DirichletProblem, HeatProblem
using JuliaFEM: LinearSolver
""" Define Problem 1:
- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}
- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁
"""
function get_heatproblem()
function test_linearsolver()
el1 = Quad4([1, 2, 3, 4])
el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
el1["temperature thermal conductivity"] = 6.0
@@ -26,278 +23,31 @@ function get_heatproblem()
(1.0 => 600.0)
)
problem1 = PlaneHeatProblem()
push!(problem1, el1)
push!(problem1, el2)
return problem1
end
field_problem = HeatProblem()
push!(field_problem, el1)
push!(field_problem, el2)
""" Define Problem 2:
- Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂
"""
function get_boundaryproblem()
el3 = Seg2([3, 4])
el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
el3["temperature"] = 0.0
problem2 = DirichletProblem("temperature", 1)
push!(problem2, el3)
return problem2
end
function test_simplesolver()
info("construct heat problem")
problem1 = get_heatproblem()
info("construct boundary problem")
problem2 = get_boundaryproblem()
boundary_problem = DirichletProblem("temperature", 1)
push!(boundary_problem, el3)
# Create a solver for a set of problems
info("create SimpleSolver with problems.")
solver = SimpleSolver()
push!(solver, problem1)
push!(solver, problem2)
info("solve!")
solver = LinearSolver(field_problem, boundary_problem)
# Solve problem at time t=1.0 and update fields
call(solver, 1.0)
solver(1.0)
# Postprocess.
# Interpolate temperature field along boundary of Γ₁ at time t=1.0
xi = [0.0, -1.0]
el2 = get_element(problem1.equations[2])
basis = get_basis(el2)
X = basis("geometry", xi, 1.0)
T = basis("temperature", xi, 1.0)
X = el2("geometry", xi, 1.0)
T = el2("temperature", xi, 1.0)
info("Temperature at point X = $X is T = $T")
@test isapprox(T, 100.0)
end
#test_simplesolver()
function atest_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[5], 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
function test_solver_multiple_dirichlet_bc()
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
e1 = Quad4([1, 2, 4, 3])
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
e1["youngs modulus"] = 900.0
e1["poissons ratio"] = 0.25
b1 = Seg2([3, 4])
b1["geometry"] = Vector[N[3], N[4]]
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
problem = PlaneStressElasticityProblem()
push!(problem, e1)
push!(problem, b1)
# manually solve problem 1
# free_dofs = [3, 5, 6, 8]
# free_dofs = [3, 6, 7, 8]
#solve!(problem, free_dofs, 0.0; max_iterations=10)
#disp = e1("displacement", [1.0, 1.0], 0.0)
#info("displacement at tip: $disp")
#@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
# boundary elements for dirichlet dx=0
dx = Seg2([1, 3])
dx["geometry"] = Vector[N[1], N[3]]
dx["displacement 1"] = 0.0
# boundary elements for dirichlet dy=0
dy = Seg2([1, 2])
dy["geometry"] = Vector[N[1], N[2]]
dy["displacement 2"] = 0.0
problem2 = DirichletProblem("displacement", 2)
push!(problem2, dx)
problem3 = DirichletProblem("displacement", 2)
push!(problem3, dy)
solver = DirectSolver()
push!(solver, problem)
push!(solver, problem2)
push!(solver, problem3)
# launch solver
norm = solver(0.0)
# info(e1("displacement"))
# info(last(e1["displacement"]))
disp = e1("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
end
function test_solver_multiple_bodies_multiple_dirichlet_bc()
N = Vector[
[0.0, 0.0], [1.0, 0.0],
[0.0, 1.0], [1.0, 1.0],
[0.0, 2.0], [1.0, 2.0]]
e1 = Quad4([1, 2, 4, 3])
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
e2 = Quad4([3, 4, 6, 5])
e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
for el in [e1, e2]
el["youngs modulus"] = 900.0
el["poissons ratio"] = 0.25
end
b1 = Seg2([5, 6])
b1["geometry"] = Vector[N[5], N[6]]
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
body1 = PlaneStressElasticityProblem()
push!(body1, e1)
body2 = PlaneStressElasticityProblem()
push!(body2, e2)
push!(body2, b1)
# boundary elements for dirichlet dx=0
dx1 = Seg2([1, 3])
dx1["geometry"] = Vector[N[1], N[3]]
dx2 = Seg2([3, 5])
dx2["geometry"] = Vector[N[3], N[5]]
for dx in [dx1, dx2]
dx["displacement 1"] = 0.0
end
boundary1 = DirichletProblem("displacement", 2)
push!(boundary1, dx1)
push!(boundary1, dx2)
# boundary elements for dirichlet dy=0
dy1 = Seg2([1, 2])
dy1["geometry"] = Vector[N[1], N[2]]
dy1["displacement 2"] = 0.0
boundary2 = DirichletProblem("displacement", 2)
push!(boundary2, dy1)
solver = DirectSolver()
push!(solver, body1)
push!(solver, body2)
push!(solver, boundary1)
push!(solver, boundary2)
# launch solver
norm = solver(0.0)
disp = e2("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# code aster verification, two_elements.comm
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
end
# test_solver_multiple_bodies_multiple_dirichlet_bc()
end
test/test_virtual_work.jl
+30 -36
View File
@@ -5,57 +5,48 @@ module TestAutoDiffWeakForm
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: Quad4, Equation, IntegrationPoint, assemble!, Assembly,
solve!, get_field, get_element, get_basis,
grad, get_default_integration_points
""" Plane stress formulation for 4-node bilinear element. """
type CPS4 <: Equation
element :: Quad4
integration_points :: Vector{IntegrationPoint}
using JuliaFEM: Problem, AbstractProblem, CG, Element, IntegrationPoint, Quad4, solve!
abstract PlaneStressElasticityProblem <: AbstractProblem
function PlaneStressElasticityProblem(dim::Int=2, elements=[])
return Problem{PlaneStressElasticityProblem}(dim, elements)
end
function JuliaFEM.get_unknown_field_name(equation::CPS4)
function JuliaFEM.get_unknown_field_name{P<:PlaneStressElasticityProblem}(::Type{P})
return "displacement"
end
function CPS4(element::Quad4)
integration_points = get_default_integration_points(element)
if !haskey(element, "displacement")
element["displacement"] = 0.0 => Vector{Float64}[[0.0,0.0], [0.0,0.0], [0.0,0.0], [0.0,0.0]]
end
CPS4(element, integration_points)
function JuliaFEM.get_unknown_field_type{P<:PlaneStressElasticityProblem}(::Type{P})
return Vector{Float64}
end
function Base.size(eq::CPS4)
return (2, 4)
end
function JuliaFEM.get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing)
function JuliaFEM.get_residual_vector(equation::CPS4, ip, time; variation=nothing)
element = get_element(equation)
basis = get_basis(element)
dbasis = grad(basis)
basis = element(ip, time)
dbasis = element(ip, time, Val{:grad})
# material parameters
E = basis("youngs modulus", ip, time)
nu = basis("poissons ratio", ip, time)
E = element("youngs modulus", ip, time)
nu = element("poissons ratio", ip, time)
mu = E/(2*(1+nu))
la = E*nu/((1+nu)*(1-2*nu))
la = 2*la*mu/(la + 2*mu) # <- correction for 2d
# elasticity formulation
u = basis("displacement", ip, time, variation)
gradu = dbasis("displacement", ip, time, variation)
u = element("displacement", ip, time, variation)
gradu = element("displacement", ip, time, Val{:grad}, variation)
F = I + gradu
b = basis("displacement volume load", ip, time)
E = 1/2*(F'*F - I)
S = la*trace(E)*I + 2*mu*E
P = F*S
r = F*S*dbasis
b = element("displacement volume load", ip, time)
r -= b*basis
# residual vector
r_int = P*dbasis(ip,time)
r_ext = b*basis(ip,time)
r = r_int - r_ext
return vec(r)
end
@@ -66,15 +57,18 @@ function test_residual_form()
element["youngs modulus"] = 500.0
element["poissons ratio"] = 0.3
element["displacement volume load"] = Vector[[0.0,-10.0], [0.0,-10.0], [0.0,-10.0], [0.0,-10.0]]
equation = CPS4(element)
element["displacement"] = (0.0 => Vector{Float64}[zeros(2) for i=1:length(element)])
problem = PlaneStressElasticityProblem()
push!(problem, element)
# create model -- end
free_dofs = [3, 4, 5, 6]
solve!(equation, free_dofs, 0.0) # launch a newton solver for single element
disp = get_basis(element)("displacement", [1.0, 1.0], 0.0)[2]
println("displacement at tip: $disp")
solve!(problem, free_dofs, 0.0) # launch a newton solver for single element
disp = element("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# verified using Code Aster.
@test isapprox(disp, -8.77303119819776E+00)
@test isapprox(disp[2], -8.77303119819776E+00)
end
end