mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-28 20:46:58 +00:00
rewrite assembly, see #69. a lot of tests probably fail but the most important ones pass
This commit is contained in:
Jukka Aho
@@ -63,6 +63,7 @@ include("lagrange.jl") # Lagrange elements
|
||||
|
||||
### EQUATIONS ###
|
||||
include("integrate.jl") # default integration points for elements
|
||||
include("sparse.jl")
|
||||
include("equations.jl")
|
||||
include("problems.jl")
|
||||
|
||||
|
||||
+4
-56
@@ -3,61 +3,9 @@
|
||||
|
||||
# Functions to handle global assembly of problem
|
||||
|
||||
""" Global assembly. """
|
||||
type GlobalAssembly <: Assembly
|
||||
ndofs :: Int
|
||||
mass_matrix :: SparseMatrixCSC
|
||||
stiffness_matrix :: SparseMatrixCSC
|
||||
force_vector :: SparseMatrixCSC
|
||||
end
|
||||
|
||||
""" Initialize global assembly of size ndofs. """
|
||||
function initialize_global_assembly(ndofs::Int=1)
|
||||
mass_matrix = spzeros(ndofs, ndofs)
|
||||
stiffness_matrix = spzeros(ndofs, ndofs)
|
||||
force_vector = spzeros(ndofs, 1)
|
||||
return GlobalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector)
|
||||
end
|
||||
|
||||
""" Initialize global assembly, get dimension from problem. """
|
||||
function initialize_global_assembly(problem::Problem)
|
||||
dim, ndofs = size(problem)
|
||||
return initialize_global_assembly(ndofs)
|
||||
end
|
||||
|
||||
""" Initialize or empty workspace for global assembly. """
|
||||
function initialize_global_assembly!(assembly::GlobalAssembly, problem::Problem)
|
||||
ndofs = prod(size(problem))
|
||||
if ndofs != assembly.ndofs
|
||||
# if problem size changes, automatically initialize new work space
|
||||
assembly.ndofs = ndofs
|
||||
assembly.mass_matrix = spzeros(ndofs, ndofs)
|
||||
assembly.stiffness_matrix = spzeros(ndofs, ndofs)
|
||||
assembly.force_vector = spzeros(ndofs, 1)
|
||||
return
|
||||
end
|
||||
# otherwise, empty workspace ready for next iteration
|
||||
fill!(assembly.mass_matrix, 0.0)
|
||||
fill!(assembly.stiffness_matrix, 0.0)
|
||||
fill!(assembly.force_vector, 0.0)
|
||||
return
|
||||
end
|
||||
|
||||
""" Calculate global assembly for a problem. """
|
||||
function calculate_global_assembly!(assembly::GlobalAssembly, problem::Problem, time::Number=Inf)
|
||||
|
||||
unknown_field_name = get_unknown_field_name(problem)
|
||||
initialize_global_assembly!(assembly, problem) # zero all
|
||||
dim, ndofs = size(problem)
|
||||
info("assembling problem for $unknown_field_name")
|
||||
info("dimension of unknown field: $dim, problem dofs: $ndofs")
|
||||
local_assembly = initialize_local_assembly()
|
||||
for (i, equation) in enumerate(get_equations(problem))
|
||||
calculate_local_assembly!(local_assembly, equation, unknown_field_name, time)
|
||||
conn = get_connectivity(get_element(equation))
|
||||
gdofs = get_gdofs(problem, equation)
|
||||
assembly.mass_matrix[gdofs, gdofs] += local_assembly.mass_matrix
|
||||
assembly.stiffness_matrix[gdofs, gdofs] += local_assembly.stiffness_matrix
|
||||
assembly.force_vector[gdofs] += local_assembly.force_vector
|
||||
function assemble!(assembly::Assembly, problem::Problem, time::Number=0.0)
|
||||
empty!(assembly)
|
||||
for equation in get_equations(problem)
|
||||
assemble!(assembly, equation, time, problem)
|
||||
end
|
||||
end
|
||||
|
||||
+8
-4
@@ -5,6 +5,10 @@
|
||||
|
||||
abstract DirichletEquation <: Equation
|
||||
|
||||
function get_unknown_field_name(equation::DirichletEquation)
|
||||
return "reaction force"
|
||||
end
|
||||
|
||||
### Dirichlet problem + equations
|
||||
|
||||
type DirichletProblem <: BoundaryProblem
|
||||
@@ -56,19 +60,19 @@ function DBC2D2(element::Seg2)
|
||||
end
|
||||
Base.size(equation::DBC2D2) = (1, 2)
|
||||
|
||||
function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing)
|
||||
initialize_local_assembly!(assembly, equation)
|
||||
function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing)
|
||||
gdofs = get_gdofs(equation)
|
||||
element = get_element(equation)
|
||||
basis = get_basis(element)
|
||||
detJ = det(basis)
|
||||
for ip in get_integration_points(equation)
|
||||
w = ip.weight * detJ(ip)
|
||||
N = basis(ip, time)
|
||||
assembly.stiffness_matrix += w * N'*N
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*N'*N)
|
||||
if !isa(problem, Void)
|
||||
X = basis("geometry", ip, time)
|
||||
u = problem.field_value(X)
|
||||
assembly.force_vector += w * N'*u
|
||||
add!(assembly.force_vector, gdofs, w*N'*u)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -6,6 +6,10 @@
|
||||
abstract ElasticityProblem <: Problem
|
||||
abstract ElasticityEquation <: Equation
|
||||
|
||||
function get_unknown_field_name(equation::ElasticityEquation)
|
||||
return "displacement"
|
||||
end
|
||||
|
||||
### Formulation ###
|
||||
|
||||
""" Calculate internal energy for elasticity equation.
|
||||
|
||||
+115
-141
@@ -5,203 +5,177 @@
|
||||
|
||||
abstract Equation
|
||||
|
||||
abstract Assembly
|
||||
|
||||
""" Local element assembly. """
|
||||
type LocalAssembly <: Assembly
|
||||
ndofs :: Int
|
||||
mass_matrix :: Matrix
|
||||
stiffness_matrix :: Matrix
|
||||
force_vector :: Matrix
|
||||
potential_energy
|
||||
residual_vector :: Vector
|
||||
type Assembly
|
||||
mass_matrix :: SparseMatrixIJV
|
||||
stiffness_matrix :: SparseMatrixIJV
|
||||
force_vector :: SparseMatrixIJV
|
||||
lhs :: SparseMatrixIJV
|
||||
rhs :: SparseMatrixIJV
|
||||
end
|
||||
|
||||
""" Initialize workspace for local matrices for dimension ndofs. """
|
||||
function initialize_local_assembly(ndofs::Int=1)
|
||||
mass_matrix = zeros(ndofs, ndofs)
|
||||
stiffness_matrix = zeros(ndofs, ndofs)
|
||||
force_vector = zeros(ndofs, 1)
|
||||
potential_energy = 0.0
|
||||
residual_vector = zeros(ndofs)
|
||||
return LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector,
|
||||
potential_energy, residual_vector)
|
||||
function Assembly()
|
||||
return Assembly(
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV())
|
||||
end
|
||||
|
||||
""" Initialize workspace for local matrices, get dimension from equation. """
|
||||
function initialize_local_assembly(equation::Equation)
|
||||
ndofs = prod(size(equation))
|
||||
return initialize_local_assembly(ndofs)
|
||||
function Base.empty!(assembly::Assembly)
|
||||
empty!(assembly.mass_matrix)
|
||||
empty!(assembly.stiffness_matrix)
|
||||
empty!(assembly.force_vector)
|
||||
empty!(assembly.lhs)
|
||||
empty!(assembly.rhs)
|
||||
end
|
||||
|
||||
""" Initialize or zero workspace. """
|
||||
function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation)
|
||||
ndofs = prod(size(equation))
|
||||
if ndofs != assembly.ndofs
|
||||
# if problem size changes, automatically initialize new work space
|
||||
assembly.ndofs = ndofs
|
||||
assembly.mass_matrix = zeros(ndofs, ndofs)
|
||||
assembly.stiffness_matrix = zeros(ndofs, ndofs)
|
||||
assembly.force_vector = zeros(ndofs, 1)
|
||||
assembly.potential_energy = 0.0
|
||||
assembly.residual_vector = zeros(ndofs)
|
||||
return
|
||||
end
|
||||
# otherwise, empty workspace ready for next iteration
|
||||
fill!(assembly.mass_matrix, 0.0)
|
||||
fill!(assembly.stiffness_matrix, 0.0)
|
||||
fill!(assembly.force_vector, 0.0)
|
||||
assembly.potential_energy = 0.0
|
||||
fill!(assembly.residual_vector, 0.0)
|
||||
return
|
||||
function get_mass_matrix
|
||||
end
|
||||
|
||||
has_mass_matrix(equation::Equation) = false
|
||||
function get_mass_matrix(equation::Equation, ip, time=0.0, problem=nothing)
|
||||
get_mass_matrix(equation, ip, time)
|
||||
end
|
||||
function get_mass_matrix(equation::Equation, ip, time=0.0)
|
||||
get_mass_matrix(equation, ip)
|
||||
end
|
||||
function get_mass_matrix(equation::Equation, ip)
|
||||
nothing
|
||||
function get_stiffness_matrix
|
||||
end
|
||||
|
||||
has_stiffness_matrix(equation::Equation) = false
|
||||
function get_stiffness_matrix(equation::Equation, ip, time=0.0, problem=nothing)
|
||||
get_stiffness_matrix(equation, ip, time)
|
||||
end
|
||||
function get_stiffness_matrix(equation::Equation, ip, time=0.0)
|
||||
get_stiffness_matrix(equation, ip)
|
||||
end
|
||||
function get_stiffness_matrix(equation::Equation, ip)
|
||||
nothing
|
||||
function get_force_vector
|
||||
end
|
||||
|
||||
has_force_vector(equation::Equation) = false
|
||||
function get_force_vector(equation::Equation, ip, time=0.0, problem=nothing)
|
||||
get_force_vector(equation, ip, time)
|
||||
end
|
||||
function get_force_vector(equation::Equation, ip, time=0.0)
|
||||
get_force_vector(equation, ip)
|
||||
end
|
||||
function get_force_vector(equation::Equation, ip)
|
||||
nothing
|
||||
function get_potential_energy
|
||||
end
|
||||
|
||||
has_residual_vector(equation::Equation) = false
|
||||
function get_residual_vector(equation::Equation, ip, time=0.0, problem=nothing)
|
||||
get_residual_vector(equation, ip, time)
|
||||
end
|
||||
function get_residual_vector(equation::Equation, ip, time=0.0)
|
||||
get_residual_vector(equation, ip)
|
||||
end
|
||||
function get_residual_vector(equation::Equation, ip)
|
||||
nothing
|
||||
function get_residual_vector
|
||||
end
|
||||
|
||||
has_potential_energy(equation::Equation) = false
|
||||
function get_potential_energy(equation::Equation, ip, time=0.0, problem=nothing)
|
||||
get_potential_energy(equation, ip, time)
|
||||
end
|
||||
function get_potential_energy(equation::Equation, ip, time=0.0)
|
||||
get_potential_energy(equation, ip)
|
||||
end
|
||||
function get_potential_energy(equation::Equation, ip)
|
||||
nothing
|
||||
function has_mass_matrix(equation::Equation)
|
||||
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
|
||||
return method_exists(get_mass_matrix, default_args)
|
||||
end
|
||||
|
||||
get_element(equation::Equation) = equation.element
|
||||
get_integration_points(equation::Equation) = equation.integration_points
|
||||
function has_stiffness_matrix(equation::Equation)
|
||||
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
|
||||
return method_exists(get_stiffness_matrix, default_args)
|
||||
end
|
||||
|
||||
function has_force_vector(equation::Equation)
|
||||
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
|
||||
return method_exists(get_force_vector, default_args)
|
||||
end
|
||||
|
||||
""" Return a local assembly for element. """
|
||||
function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
|
||||
unknown_field_name::ASCIIString, time::Number=0.0,
|
||||
problem=nothing)
|
||||
function has_potential_energy(equation::Equation)
|
||||
default_args = Tuple{typeof(equation), IntegrationPoint, Float64}
|
||||
return method_exists(get_potential_energy, default_args)
|
||||
end
|
||||
|
||||
initialize_local_assembly!(assembly, equation) # zero all
|
||||
function has_residual_vector(equation::Equation)
|
||||
default_args = Tuple{typeof(equation), 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
|
||||
|
||||
""" Assemble element. """
|
||||
function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, problem=nothing)
|
||||
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(equation)
|
||||
basis = get_basis(element)
|
||||
detJ = det(basis)
|
||||
unknown_field_name = get_unknown_field_name(equation)
|
||||
|
||||
# 1. if equations are defined we just integrate them
|
||||
# 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)
|
||||
assembly.mass_matrix += s*get_mass_matrix(equation, ip, time, problem)
|
||||
add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(equation, ip, time))
|
||||
end
|
||||
if has_stiffness_matrix(equation)
|
||||
assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time, problem)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(equation, ip, time))
|
||||
end
|
||||
if has_force_vector(equation)
|
||||
assembly.force_vector += s*get_force_vector(equation, ip, time, problem)
|
||||
add!(assembly.force_vector, gdofs, s*get_force_vector(equation, ip, time))
|
||||
end
|
||||
end
|
||||
# external loads -- if any nodal loads is defined add to force vector
|
||||
if haskey(element, "$unknown_field_name nodal load")
|
||||
assembly.force_vector += vec(element["$unknown_field_name nodal load"](time))
|
||||
add!(assembly.force_vector, gdofs, vec(element["$unknown_field_name nodal load"](time)))
|
||||
end
|
||||
end
|
||||
|
||||
# 2. variational / energy form - user has defined some potential energy / variational form
|
||||
# 2. energy form -- user has defined potential energy W -> min!
|
||||
if has_potential_energy(equation)
|
||||
element = get_element(equation)
|
||||
field = element[unknown_field_name](time)
|
||||
function potential_energy(data::Vector)
|
||||
# calculate potential energy for some setting. this is needed by forwarddiff
|
||||
assembly.potential_energy = 0.0
|
||||
|
||||
""" Wrapper for potential energy for ForwardDiff. """
|
||||
function calc_W(data::Vector)
|
||||
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)
|
||||
assembly.potential_energy += ip.weight * dw * detJ(ip)
|
||||
W += s*dw
|
||||
end
|
||||
# external energy -- if any nodal loads is defined, decrease from potential energy
|
||||
if haskey(element, "$unknown_field_name nodal load")
|
||||
P = element["$unknown_field_name nodal load"](time)
|
||||
assembly.potential_energy -= dot(vec(P), vec(df))
|
||||
W -= dot(vec(P), vec(df))
|
||||
end
|
||||
if isa(assembly.potential_energy, Array)
|
||||
return assembly.potential_energy[1]
|
||||
end
|
||||
return assembly.potential_energy
|
||||
return isa(W, Array) ? W[1] : W
|
||||
end
|
||||
hessian, allresults = ForwardDiff.hessian(potential_energy, vec(field),
|
||||
AllResults, cache=autodiffcache)
|
||||
assembly.stiffness_matrix += hessian
|
||||
assembly.force_vector -= ForwardDiff.gradient(allresults) # <--- minus explained in tutorial
|
||||
assembly.potential_energy = ForwardDiff.value(allresults)
|
||||
#info("potential energy of system: $(assembly.potential_energy)")
|
||||
|
||||
hessian, allresults = ForwardDiff.hessian(calc_W, vec(field), AllResults, cache=autodiffcache)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, hessian)
|
||||
add!(assembly.force_vector, gdofs, -ForwardDiff.gradient(allresults))
|
||||
end
|
||||
|
||||
# 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v
|
||||
# 3. virtual work -- user has defined some residual r = p - f = 0
|
||||
if has_residual_vector(equation)
|
||||
element = get_element(equation)
|
||||
field = element[unknown_field_name](time)
|
||||
function residual_vector(data::Vector)
|
||||
fill!(assembly.residual_vector, 0.0)
|
||||
#@debug("field: $field, length = $(size(field))")
|
||||
#@debug("data: $data, size = $(size(data))")
|
||||
#df = similar(field, data)
|
||||
df = Increment(reshape(data, size(equation)...))
|
||||
# integrate W
|
||||
for ip in get_integration_points(equation)
|
||||
dr = get_residual_vector(equation, ip, time; variation=df)
|
||||
assembly.residual_vector += ip.weight*dr*detJ(ip)
|
||||
end
|
||||
# external loads -- if any nodal loads is defined, remove from residual
|
||||
if haskey(element, "$unknown_field_name nodal load")
|
||||
assembly.residual_vector -= vec(element["$unknown_field_name nodal load"](time))
|
||||
end
|
||||
return assembly.residual_vector
|
||||
end
|
||||
jacobian, allresults = ForwardDiff.jacobian(residual_vector, vec(field),
|
||||
AllResults, cache=autodiffcache)
|
||||
assembly.stiffness_matrix += jacobian
|
||||
assembly.force_vector -= ForwardDiff.value(allresults) # <-- minus explained in tutorial
|
||||
end
|
||||
|
||||
""" Wrapper for virtual work for ForwardDiff. """
|
||||
function calc_R(data::Vector)
|
||||
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)
|
||||
R += s*dr
|
||||
end
|
||||
# external loads -- if any nodal loads is defined, decrease from residual
|
||||
if haskey(element, "$unknown_field_name nodal load")
|
||||
R -= vec(element["$unknown_field_name nodal load"](time))
|
||||
end
|
||||
return R
|
||||
end
|
||||
|
||||
jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
|
||||
add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))
|
||||
end
|
||||
end
|
||||
|
||||
+30
-27
@@ -6,6 +6,10 @@
|
||||
abstract HeatProblem <: Problem
|
||||
abstract HeatEquation <: Equation
|
||||
|
||||
function get_unknown_field_name(equation::HeatEquation)
|
||||
return "temperature"
|
||||
end
|
||||
|
||||
### Formulation ###
|
||||
|
||||
""" Heat equations.
|
||||
@@ -32,55 +36,36 @@ References
|
||||
https://en.wikipedia.org/wiki/Heat_equation
|
||||
|
||||
"""
|
||||
function calculate_local_assembly!(assembly::LocalAssembly, equation::HeatEquation,
|
||||
unknown_field_name::ASCIIString, time::Number=Inf,
|
||||
problem=nothing)
|
||||
|
||||
initialize_local_assembly!(assembly, equation)
|
||||
function assemble!(assembly::Assembly, equation::HeatEquation, time::Number=0.0, problem=nothing)
|
||||
|
||||
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)
|
||||
w = ip.weight*detJ(ip)
|
||||
N = basis(ip, time)
|
||||
if haskey(element, "density")
|
||||
rho = basis("density", ip, time)
|
||||
assembly.mass_matrix += w * rho*N'*N
|
||||
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)
|
||||
assembly.stiffness_matrix += w * k*dN'*dN
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN)
|
||||
end
|
||||
if haskey(element, "temperature load")
|
||||
f = basis("temperature load", ip, time)
|
||||
assembly.force_vector += w * N'*f
|
||||
add!(assembly.force_vector, gdofs, w*N'*f)
|
||||
end
|
||||
if haskey(element, "temperature flux")
|
||||
g = basis("temperature flux", ip, time)
|
||||
assembly.force_vector += w * N'*g
|
||||
add!(assembly.force_vector, gdofs, w*N'*g)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
### Problems ###
|
||||
|
||||
type PlaneHeatProblem <: HeatProblem
|
||||
unknown_field_name :: ASCIIString
|
||||
unknown_field_dimension :: Int
|
||||
equations :: Array{HeatEquation, 1}
|
||||
element_mapping :: Dict{DataType, DataType}
|
||||
end
|
||||
|
||||
""" Default constructor for problem takes no arguments. """
|
||||
function PlaneHeatProblem()
|
||||
element_mapping = Dict(
|
||||
Quad4 => DC2D4,
|
||||
Seg2 => DC2D2)
|
||||
return PlaneHeatProblem("temperature", 1, [], element_mapping)
|
||||
end
|
||||
|
||||
### Equations ###
|
||||
|
||||
@@ -101,7 +86,7 @@ Base.size(equation::DC2D4) = (1, 4)
|
||||
""" Diffusive heat transfer for 2-node linear segment. """
|
||||
type DC2D2 <: HeatEquation
|
||||
element :: Seg2
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
integration_points :: Vector{IntegrationPoint}
|
||||
end
|
||||
function DC2D2(element::Seg2)
|
||||
integration_points = get_default_integration_points(element)
|
||||
@@ -112,3 +97,21 @@ function DC2D2(element::Seg2)
|
||||
end
|
||||
Base.size(equation::DC2D2) = (1, 2)
|
||||
|
||||
### Problems ###
|
||||
|
||||
type PlaneHeatProblem <: HeatProblem
|
||||
unknown_field_name :: ASCIIString
|
||||
unknown_field_dimension :: Int
|
||||
equations :: Vector{Equation}
|
||||
#element_mapping :: Dict{Element, Equation}
|
||||
# FIXME: Why is not working ^
|
||||
element_mapping :: Dict{Any, Any}
|
||||
end
|
||||
|
||||
""" Default constructor for problem takes no arguments. """
|
||||
function PlaneHeatProblem()
|
||||
element_mapping = Dict(
|
||||
Quad4 => DC2D4,
|
||||
Seg2 => DC2D2)
|
||||
return PlaneHeatProblem("temperature", 1, [], element_mapping)
|
||||
end
|
||||
|
||||
+23
-18
@@ -9,18 +9,19 @@ 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, unknown_field_name::ASCIIString,
|
||||
free_dofs::Array{Int, 1}, time::Number=0.0;
|
||||
function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number=0.0;
|
||||
max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
|
||||
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)
|
||||
la = initialize_local_assembly()
|
||||
ass = Assembly()
|
||||
for i=1:max_iterations
|
||||
calculate_local_assembly!(la, equation, unknown_field_name)
|
||||
A = la.stiffness_matrix[free_dofs, free_dofs]
|
||||
b = la.force_vector[free_dofs]
|
||||
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)')
|
||||
@@ -39,30 +40,34 @@ to test nonlinear element formulations. Dirichlet boundary is assumed to be homo
|
||||
and degrees of freedom are eliminated. So if boundary condition is known in nodal
|
||||
points and everything is zero this should be quite good.
|
||||
"""
|
||||
function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=1.0;
|
||||
max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
|
||||
function solve!(problem::Problem, free_dofs::Vector{Int}, time::Number=1.0; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
|
||||
info("start solver")
|
||||
ga = initialize_global_assembly(problem)
|
||||
x = zeros(ga.ndofs)
|
||||
dx = fill!(similar(x), 0.0)
|
||||
assembly = Assembly()
|
||||
# x = zeros(ga.ndofs)
|
||||
# dx = fill!(similar(x), 0.0)
|
||||
# FIXME: better.
|
||||
x = nothing
|
||||
dx = nothing
|
||||
field_name = get_unknown_field_name(problem)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
for i=1:max_iterations
|
||||
info("calculate global assembly")
|
||||
calculate_global_assembly!(ga, problem)
|
||||
info("done")
|
||||
A = ga.stiffness_matrix[free_dofs, free_dofs]
|
||||
b = ga.force_vector[free_dofs]
|
||||
assemble!(assembly, problem, time)
|
||||
A = sparse(assembly.stiffness_matrix)
|
||||
b = sparse(assembly.force_vector)
|
||||
if dump_matrices
|
||||
dump(full(A))
|
||||
dump(full(b)')
|
||||
end
|
||||
dx[free_dofs] = lufact(A) \ full(b)
|
||||
if isa(dx, Void)
|
||||
x = zeros(length(b))
|
||||
dx = zeros(length(b))
|
||||
end
|
||||
dx[free_dofs] = lufact(A[free_dofs,free_dofs]) \ full(b)[free_dofs]
|
||||
info("Difference in solution norm: $(norm(dx))")
|
||||
x += dx
|
||||
for equation in get_equations(problem)
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(problem, equation)
|
||||
gdofs = get_gdofs(equation)
|
||||
data = reshape(full(x[gdofs]), size(equation))
|
||||
push!(element[field_name], data)
|
||||
end
|
||||
|
||||
@@ -14,6 +14,10 @@ function SparseMatrixIJV()
|
||||
SparseMatrixIJV([], [], [])
|
||||
end
|
||||
|
||||
function Base.sparse(A::SparseMatrixIJV)
|
||||
return sparse(A.I, A.J, A.V)
|
||||
end
|
||||
|
||||
function Base.push!(A::SparseMatrixIJV, I::Int, J::Int, V::Float64)
|
||||
push!(A.I, I)
|
||||
push!(A.J, J)
|
||||
@@ -63,3 +67,10 @@ function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data::
|
||||
append!(A.V, vec(data))
|
||||
end
|
||||
|
||||
""" Sparse vector version. """
|
||||
function add!(A::SparseMatrixIJV, dofs::Vector{Int}, data::Array{Float64})
|
||||
append!(A.I, dofs)
|
||||
append!(A.J, ones(Int, length(dofs)))
|
||||
append!(A.V, vec(data))
|
||||
end
|
||||
|
||||
|
||||
Reference in New Issue
Block a user