mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-30 08:02:50 +00:00
removed some obsolete code
This commit is contained in:
Jukka Aho
@@ -1,3 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
@@ -1,2 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
@@ -1,132 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
function get_mass_matrix
|
||||
end
|
||||
|
||||
function get_stiffness_matrix
|
||||
end
|
||||
|
||||
function get_force_vector
|
||||
end
|
||||
|
||||
function get_potential_energy
|
||||
end
|
||||
|
||||
function get_residual_vector
|
||||
end
|
||||
|
||||
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(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(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(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(problem::Problem, element::Element)
|
||||
default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
|
||||
return method_exists(get_residual_vector, default_args)
|
||||
end
|
||||
|
||||
|
||||
""" Assemble element. """
|
||||
function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number)
|
||||
|
||||
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(problem, element) || has_stiffness_matrix(problem, element) || has_force_vector(problem, element)
|
||||
for ip in get_integration_points(element)
|
||||
w = ip.weight*det(J)
|
||||
if has_mass_matrix(element)
|
||||
add!(assembly.mass_matrix, gdofs, gdofs, w*get_mass_matrix(problem, element, ip, time))
|
||||
end
|
||||
if has_stiffness_matrix(element)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*get_stiffness_matrix(problem, element, ip, time))
|
||||
end
|
||||
if has_force_vector(element)
|
||||
add!(assembly.force_vector, gdofs, w*get_force_vector(problem, element, ip, time))
|
||||
end
|
||||
end
|
||||
# external loads -- if any nodal loads is defined add to force vector
|
||||
if haskey(element, "$unknown_field_name nodal load")
|
||||
add!(assembly.force_vector, gdofs, vec(element["$unknown_field_name nodal load"](time)))
|
||||
end
|
||||
end
|
||||
|
||||
# 2. energy form -- user has defined potential energy W -> min!
|
||||
if has_potential_energy(problem, element) && haskey(element, unknown_field_name)
|
||||
field = element[unknown_field_name](time)
|
||||
|
||||
""" 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(element)
|
||||
dw = get_potential_energy(problem, element, ip, time; variation=df)
|
||||
W += ip.weight*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)
|
||||
W -= dot(vec(P), vec(df))
|
||||
end
|
||||
return W[1]
|
||||
end
|
||||
|
||||
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 -- user has defined some residual r = p - f = 0
|
||||
if has_residual_vector(problem, element) && haskey(element, unknown_field_name)
|
||||
|
||||
field = DVTI(last(element[unknown_field_name]).data)
|
||||
|
||||
""" Wrapper for virtual work for ForwardDiff. """
|
||||
function calc_R(data::Vector)
|
||||
R = zeros(length(data))
|
||||
df = similar(field, data)
|
||||
gauss_fields = IntegrationPoint[]
|
||||
# integrate residual vector
|
||||
for ip in get_integration_points(element)
|
||||
dr = get_residual_vector(problem, element, ip, time; variation=df)
|
||||
R += ip.weight*dr
|
||||
if ip.changed
|
||||
push!(gauss_fields, ip)
|
||||
end
|
||||
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
|
||||
#info("return = $R")
|
||||
if length(gauss_fields) != 0
|
||||
update_gauss_fields!(element, gauss_fields, time)
|
||||
end
|
||||
return R
|
||||
end
|
||||
|
||||
jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
|
||||
residual_vector = -ForwardDiff.value(allresults)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
|
||||
add!(assembly.force_vector, gdofs, residual_vector)
|
||||
end
|
||||
end
|
||||
|
||||
-116
@@ -1,116 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
|
||||
immutable SymbolicField <: AbstractField
|
||||
name :: ASCIIString
|
||||
end
|
||||
|
||||
immutable Expression
|
||||
expr :: Expr
|
||||
end
|
||||
|
||||
function Expression(expression::ASCIIString)
|
||||
return Expression(parse(expression))
|
||||
end
|
||||
|
||||
function Base.convert(::Type{Field}, name::ASCIIString)
|
||||
return SymbolicField(name)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{Symbol}, field::SymbolicField)
|
||||
return Symbol(field.name)
|
||||
end
|
||||
|
||||
function Base.(:*)(n::Number, field::SymbolicField)
|
||||
expr = Expr(:call, :*, n, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function grad(field::SymbolicField)
|
||||
expr = Expr(:call, :grad, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function diff(field::SymbolicField)
|
||||
expr = Expr(:call, :diff, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function Base.call(basis::Basis, expr::Expression, fieldset::FieldSet,
|
||||
xi::Vector, time::Number)
|
||||
return replace(expr.expr, basis, fieldset, xi, time)
|
||||
end
|
||||
|
||||
# Unbelievable code. I don't know why or how this works.
|
||||
""" Replace symbolic fields with real arrays. """
|
||||
function Base.replace(expression::Expr, basis::Basis, fieldset::FieldSet,
|
||||
xi::Vector, time::Number, data=Dict())
|
||||
|
||||
info("expression = $expression")
|
||||
if expression.head == symbol("'")
|
||||
info("transpose")
|
||||
expr = Expr(:call, :transpose, expression.args...)
|
||||
return replace(expr, basis, fieldset, xi, time, data)
|
||||
end
|
||||
operator = expression.args[1]
|
||||
|
||||
if operator == :diff
|
||||
info("inside diff operator")
|
||||
if !haskey(data, expression)
|
||||
data[expression] = fieldset[string(expression.args[2])](time, Val{:diff})
|
||||
end
|
||||
info("data = $(data[expression])")
|
||||
#return basis(data[expression], xi)
|
||||
return data[expression]
|
||||
end
|
||||
|
||||
if operator == :grad
|
||||
info("inside gradient operator")
|
||||
info("gradient args: $(expression.args)")
|
||||
field_name = expression.args[2]
|
||||
if isa(field_name, Expr)
|
||||
info("expression inside gradient")
|
||||
end
|
||||
#if startswith(string(field_name), "diff")
|
||||
if !haskey(data, field_name)
|
||||
info("grad: evaluate field $field_name")
|
||||
if isa(field_name, Symbol)
|
||||
data[field_name] = fieldset[string(field_name)](time)
|
||||
elseif isa(field_name, Expr) && (field_name.args[1] == :diff)
|
||||
info("taking time derivative of $(field_name.args[2])")
|
||||
#data[field_name] = fieldset[string(field_name.args[2])](time, Val{:diff})
|
||||
data[field_name] = replace(field_name, basis, fieldset, xi, time, data)
|
||||
end
|
||||
end
|
||||
if !haskey(data, :geometry)
|
||||
info("evaluate geometry")
|
||||
data[:geometry] = fieldset["geometry"](time)
|
||||
end
|
||||
#return Expr(:call, basis, data[:geometry], data[field_name], xi, Val{:gradient})
|
||||
#return :(basis($(data[:geometry]), $(data[field_name]), $xi, Val{:gradient}))
|
||||
info("evaluate gradient")
|
||||
return basis(data[:geometry], data[field_name], xi, Val{:grad})
|
||||
end
|
||||
|
||||
for i in 2:length(expression.args)
|
||||
arg = expression.args[i]
|
||||
if isa(arg, Expr)
|
||||
expression.args[i] = replace(arg, basis, fieldset, xi, time, data)
|
||||
end
|
||||
if isa(arg, Symbol) && haskey(fieldset, string(arg))
|
||||
if !haskey(data, arg)
|
||||
info("evaluate field $arg")
|
||||
data[arg] = fieldset[string(arg)](time)
|
||||
end
|
||||
#expression.args[i] = Expr(:call, basis, data[arg], xi)
|
||||
#return :(basis($(data[arg]), $xi))
|
||||
expression.args[i] = basis(data[arg], xi)
|
||||
end
|
||||
end
|
||||
|
||||
info("new expression: $expression")
|
||||
return expression
|
||||
end
|
||||
|
||||
@@ -1,119 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
|
||||
function get_basis()
|
||||
basis(xi) = 1/4*[
|
||||
(1-xi[1])*(1-xi[2])
|
||||
(1+xi[1])*(1-xi[2])
|
||||
(1+xi[1])*(1+xi[2])
|
||||
(1-xi[1])*(1+xi[2])]'
|
||||
dbasis(xi) = 1/4*[
|
||||
-(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2])
|
||||
-(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])]
|
||||
return Basis(basis, dbasis)
|
||||
end
|
||||
|
||||
function get_fieldset()
|
||||
X = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
|
||||
T = Field(
|
||||
(0.0, [0, 0, 0, 0]),
|
||||
(1.0, [1, 2, 3, 4]))
|
||||
u = Field(
|
||||
(0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
|
||||
(1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
|
||||
# FIXME: this is not working
|
||||
#fieldset = FieldSet("geometry" => X, "temperature" => T, "displacement" => u)
|
||||
fieldset = FieldSet()
|
||||
fieldset["geometry"] = X
|
||||
fieldset["temperature"] = T
|
||||
fieldset["displacement"] = u
|
||||
return fieldset
|
||||
end
|
||||
|
||||
function test_create_symbolic_field()
|
||||
f = Field("temperature")
|
||||
@test isa(f, Field)
|
||||
end
|
||||
|
||||
function test_evaluate_symbolic_field()
|
||||
T = Field("temperature")
|
||||
expr = Symbol(T)
|
||||
@test expr == :(temperature)
|
||||
end
|
||||
|
||||
function test_simple_math()
|
||||
T = Field("temperature")
|
||||
eq = 1/2*T
|
||||
@test eq.expr == :(0.5*temperature)
|
||||
end
|
||||
|
||||
function test_evaluate_expression()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = 1/2*T
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), 1/2*mean([1, 2, 3, 4]))
|
||||
end
|
||||
|
||||
function test_evaluate_gradient_of_scalar_field()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = grad(T)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [0.0 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_gradient_of_vector_field()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
expr = grad(u)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_strain_rate()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
#expr = grad(u)
|
||||
expr = Expression("1/2*(grad(diff(displacement)) + grad(diff(displacement))')")
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.75; 0.75 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_grad_diff()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
#expr = grad(u)
|
||||
expr = Expression("grad(diff(displacement))")
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
|
||||
end
|
||||
|
||||
function test_grad_diff_simplification()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
expr1 = diff(grad(u))
|
||||
expr2 = Expression("grad(diff(displacement))")
|
||||
info("expr1 = $expr1")
|
||||
info("expr2 = $expr2")
|
||||
@test expr1 == expr2
|
||||
end
|
||||
|
||||
function test_evaluate_time_derivative()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = diff(T)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), mean([1.0, 2.0, 3.0, 4.0]))
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user