mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
first verification.. failed!
This commit is contained in:
Jukka Aho
@@ -1499,7 +1499,7 @@
|
||||
" # solve problem\n",
|
||||
" nz = unique(rowvals(A))\n",
|
||||
" x = zeros(b)\n",
|
||||
" x[nz] = lufact(Atot[nz,nz]) \\ full(b[nz])\n",
|
||||
" x[nz] = lufact(A[nz,nz]) \\ full(b[nz])\n",
|
||||
"\n",
|
||||
" # get \"problem-wise\" solution vectors\n",
|
||||
" x1 = x[1:length(b1)]\n",
|
||||
|
||||
+18
-8
@@ -7,6 +7,14 @@ using Lexicon
|
||||
using Logging
|
||||
@Logging.configure(level=DEBUG)
|
||||
|
||||
"""
|
||||
Simple linspace extension to multidimensional values. Contribute to julialang?
|
||||
"""
|
||||
function Base.linspace(X1, X2, n)
|
||||
[1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
|
||||
end
|
||||
|
||||
|
||||
include("types.jl") # type definitions
|
||||
|
||||
include("interpolate.jl") # interpolation routines
|
||||
@@ -16,16 +24,18 @@ include("elements.jl")
|
||||
include("lagrange.jl") # Lagrange elements
|
||||
#include("hierarchical.jl") # P-elements
|
||||
|
||||
include("equations.jl")
|
||||
include("problems.jl")
|
||||
include("solvers.jl")
|
||||
|
||||
include("equations.jl") # formulations
|
||||
include("problems.jl") # problems
|
||||
|
||||
|
||||
|
||||
include("math.jl") # basic mathematical operations -- obsolete ..?
|
||||
include("elasticity_solver.jl")
|
||||
#include("math.jl") # basic mathematical operations -- obsolete ..?
|
||||
# pre- and postprocess
|
||||
include("xdmf.jl")
|
||||
include("abaqus_reader.jl")
|
||||
include("interfaces.jl")
|
||||
#include("interfaces.jl")
|
||||
|
||||
include("dirichlet.jl")
|
||||
include("heat.jl")
|
||||
#include("elasticity_solver.jl")
|
||||
|
||||
end # module
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
# 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 DirichletEquation <: Equation
|
||||
|
||||
get_unknown_field_name(eq::DirichletEquation) = symbol("reaction force")
|
||||
|
||||
### Dirichlet problem + equations
|
||||
|
||||
type DirichletProblem <: BoundaryProblem
|
||||
equations :: Array{DirichletEquation, 1}
|
||||
end
|
||||
function DirichletProblem()
|
||||
DirichletProblem([])
|
||||
end
|
||||
get_dimension(pr::Type{DirichletProblem}) = 1 # ..?
|
||||
get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2
|
||||
|
||||
"""
|
||||
Dirichlet boundary condition element for 2 node line segment
|
||||
"""
|
||||
type DBC2D2 <: DirichletEquation
|
||||
element :: Seg2
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
global_dofs :: Array{Int64, 1}
|
||||
fieldval :: Function
|
||||
end
|
||||
function DBC2D2(el::Seg2)
|
||||
integration_points = [
|
||||
IntegrationPoint([-sqrt(1/3)], 1.0),
|
||||
IntegrationPoint([+sqrt(1/3)], 1.0)]
|
||||
new_fieldset!(el, "reaction force")
|
||||
fieldval(X, t) = 0.0
|
||||
DBC2D2(el, integration_points, [], fieldval)
|
||||
end
|
||||
function get_lhs(eq::DBC2D2, ip, t)
|
||||
el = get_element(eq)
|
||||
h = get_basis(el)(ip.xi)
|
||||
return h*h'
|
||||
end
|
||||
function get_rhs(eq::DBC2D2, ip, t)
|
||||
el = get_element(eq)
|
||||
h = get_basis(el, ip.xi)
|
||||
f = eq.fieldval
|
||||
X = interpolate(el, "geometry", ip.xi, t)
|
||||
return h*f(X, t)
|
||||
end
|
||||
has_lhs(eq::DBC2D2) = true
|
||||
has_rhs(eq::DBC2D2) = true
|
||||
|
||||
+75
@@ -0,0 +1,75 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Heat problems
|
||||
|
||||
abstract HeatProblem <: Problem
|
||||
abstract HeatEquation <: Equation
|
||||
|
||||
get_unknown_field_name(eq::HeatEquation) = symbol("temperature")
|
||||
|
||||
|
||||
### Plane heat problem + equations ###
|
||||
|
||||
type PlaneHeatProblem <: HeatProblem
|
||||
equations :: Array{HeatEquation, 1}
|
||||
end
|
||||
|
||||
""" Default constructor for problem takes no arguments. """
|
||||
function PlaneHeatProblem()
|
||||
return PlaneHeatProblem([])
|
||||
end
|
||||
|
||||
""" Return dimension of unknown field variable, temperature is scalar field. """
|
||||
get_dimension(pr::Type{PlaneHeatProblem}) = 1
|
||||
|
||||
""" Map Lagrange element Quad4 to equation DC2D4 """
|
||||
get_equation(pr::Type{PlaneHeatProblem}, el::Type{Quad4}) = DC2D4
|
||||
|
||||
""" Map Lagrange element Seg2 to equation DC2D2 """
|
||||
get_equation(pr::Type{PlaneHeatProblem}, el::Type{Seg2}) = DC2D2
|
||||
|
||||
|
||||
""" Diffusive heat transfer for 4-node bilinear element. """
|
||||
type DC2D4 <: HeatEquation
|
||||
element :: Quad4
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
global_dofs :: Array{Int64, 1}
|
||||
end
|
||||
function DC2D4(el::Quad4)
|
||||
integration_points = [
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]
|
||||
new_fieldset!(el, "temperature")
|
||||
DC2D4(el, integration_points, [])
|
||||
end
|
||||
function get_lhs(eq::DC2D4, ip, t)
|
||||
el = get_element(eq)
|
||||
dNdX = get_dbasisdX(el, ip.xi, t)
|
||||
k = interpolate(el, "temperature thermal conductivity", ip.xi, t)
|
||||
return dNdX*k*dNdX'
|
||||
end
|
||||
JuliaFEM.has_lhs(eq::DC2D4) = true
|
||||
|
||||
""" Diffusive heat transfer for 2-node linear segment. """
|
||||
type DC2D2 <: HeatEquation
|
||||
element :: Seg2
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
global_dofs :: Array{Int64, 1}
|
||||
end
|
||||
function DC2D2(el::Seg2)
|
||||
integration_points = [IntegrationPoint([0.0], 1.0)]
|
||||
new_fieldset!(el, "temperature")
|
||||
DC2D2(el, integration_points, [])
|
||||
end
|
||||
function get_rhs(eq::DC2D2, ip, t)
|
||||
el = get_element(eq)
|
||||
h = get_basis(el, ip.xi)
|
||||
f = interpolate(el, "temperature flux", ip.xi, t)
|
||||
println("h = $h")
|
||||
println("f = $f")
|
||||
return h*f
|
||||
end
|
||||
JuliaFEM.has_rhs(eq::DC2D2) = true
|
||||
@@ -1,58 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract Heat <: Equation
|
||||
|
||||
"""
|
||||
Diffusive heat transfer for 4-node bilinear element.
|
||||
"""
|
||||
type DC2D4 <: Heat
|
||||
element :: Quad4
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
global_dofs :: Array{Int64, 1}
|
||||
end
|
||||
function DC2D4(el::Quad4)
|
||||
integration_points = [
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]
|
||||
set_field(el, :temperature, zeros(2, 4))
|
||||
DC2D4(el, integration_points, [])
|
||||
end
|
||||
function get_lhs(eq::DC2D4)
|
||||
function get_lhs_(eq, ip)
|
||||
el = get_element(eq)
|
||||
xi = ip.xi
|
||||
dNdX = get_dbasisdX(el, xi)'
|
||||
hc = interpolate(el, :"temperature heat coefficient", xi)
|
||||
return dNdX'*hc*dNdX
|
||||
end
|
||||
integrate(eq, get_lhs_)
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Diffusive heat transfer for 2-node linear segment.
|
||||
"""
|
||||
type DC2D2 <: Heat
|
||||
element :: Seg2
|
||||
integration_points :: Array{IntegrationPoint, 1}
|
||||
global_dofs :: Array{Int64, 1}
|
||||
end
|
||||
function DC2D2(el::Seg2)
|
||||
integration_points = [IntegrationPoint([0.0], 1.0)]
|
||||
set_field(el, :temperature, zeros(2, 1))
|
||||
set_field(el, :"temperature flux", zeros(2, 1))
|
||||
DC2D2(el, integration_points, [])
|
||||
end
|
||||
function get_rhs(eq::DC2D2)
|
||||
function get_rhs_(eq, ip)
|
||||
el = get_element(eq)
|
||||
xi = ip.xi
|
||||
N = get_basis(el, xi)
|
||||
f = interpolate(el, :"temperature flux", xi)
|
||||
return f*N
|
||||
end
|
||||
integrate(eq, get_rhs_)
|
||||
end
|
||||
+90
-2
@@ -2,10 +2,18 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract Problem
|
||||
abstract BoundaryProblem <: Problem
|
||||
abstract FieldProblem <: Problem
|
||||
|
||||
get_equations(pr::Problem) = pr.equations
|
||||
get_dimension(pr::Type{Problem}) = nothing
|
||||
get_equation(pr::Type{Problem}, el::Type{Element}) = nothing
|
||||
|
||||
function get_dimension(pr::Type{Problem})
|
||||
throw("Unable to determine problem dimension for problem $pr")
|
||||
end
|
||||
|
||||
function get_equation(pr::Type{Problem}, el::Type{Element})
|
||||
throw("Could not find corresponding equation for element $el in problem $pr")
|
||||
end
|
||||
|
||||
"""
|
||||
Add new element to problem
|
||||
@@ -50,3 +58,83 @@ function set_global_dofs!(pr::Problem)
|
||||
set_global_dofs!(eq, gconn)
|
||||
end
|
||||
end
|
||||
|
||||
function get_connectivity(pr::Problem)
|
||||
conn = Int[]
|
||||
for eq in get_equations(pr)
|
||||
el = get_element(eq)
|
||||
append!(conn, get_connectivity(el))
|
||||
end
|
||||
conn = unique(conn)
|
||||
return conn
|
||||
end
|
||||
|
||||
"""
|
||||
Calculate global dofs for equations, maybe using some bandwidth
|
||||
minimizing or fill reducing algorithm
|
||||
"""
|
||||
function calculate_global_dofs(pr::Problem)
|
||||
conn = get_connectivity(pr)
|
||||
dim = get_dimension(typeof(pr))
|
||||
ndofs = dim*length(conn)
|
||||
Logging.debug("total dofs: $ndofs")
|
||||
|
||||
mconn = maximum(conn)
|
||||
gdofs = reshape(collect(1:mconn), dim, mconn)
|
||||
dofmap = Dict{Int64, Array{Int64, 1}}()
|
||||
for (i, c) in enumerate(conn)
|
||||
dofmap[c] = gdofs[:, i]
|
||||
end
|
||||
return dofmap
|
||||
end
|
||||
|
||||
"""
|
||||
Assign global dofs for equations.
|
||||
"""
|
||||
function assign_global_dofs!(pr::Problem, dofmap)
|
||||
for eq in get_equations(pr)
|
||||
el = get_element(eq)
|
||||
c = get_connectivity(el)
|
||||
#gdofs = [dofmap[ci] for ci in c]
|
||||
gdofs = Int64[]
|
||||
for ci in c
|
||||
append!(gdofs, dofmap[ci])
|
||||
end
|
||||
set_global_dofs!(eq, gdofs)
|
||||
end
|
||||
end
|
||||
|
||||
function get_lhs(pr::Problem, t::Float64)
|
||||
I = Int64[]
|
||||
J = Int64[]
|
||||
V = Float64[]
|
||||
dim = get_dimension(typeof(pr))
|
||||
for eq in filter(has_lhs, get_equations(pr))
|
||||
dofs = get_global_dofs(eq)
|
||||
lhs = integrate_lhs(eq, t)
|
||||
for (li, i) in enumerate(dofs)
|
||||
for (lj, j) in enumerate(dofs)
|
||||
push!(I, i)
|
||||
push!(J, j)
|
||||
push!(V, lhs[li, lj])
|
||||
end
|
||||
end
|
||||
end
|
||||
return I, J, V
|
||||
end
|
||||
|
||||
function get_rhs(pr::Problem, t::Float64)
|
||||
I = Int64[]
|
||||
V = Float64[]
|
||||
dim = get_dimension(typeof(pr))
|
||||
for eq in filter(has_rhs, get_equations(pr))
|
||||
dofs = get_global_dofs(eq)
|
||||
rhs = integrate_rhs(eq, t)
|
||||
for (li, i) in enumerate(dofs)
|
||||
push!(I, i)
|
||||
push!(V, rhs[li])
|
||||
end
|
||||
end
|
||||
return I, V
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,95 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Solver stuff
|
||||
|
||||
abstract Solver
|
||||
|
||||
"""
|
||||
Add new problem to solver
|
||||
"""
|
||||
function add_problem!(solver::Solver, problem::Problem)
|
||||
push!(solver.problems, problem)
|
||||
end
|
||||
|
||||
"""
|
||||
Get all problems assigned to solver
|
||||
"""
|
||||
function get_problems(s::Solver)
|
||||
return s.problems
|
||||
end
|
||||
|
||||
## SimpleSolver -- tiny direct demo solver
|
||||
""" Simple solver for educational purposes. """
|
||||
type SimpleSolver <: Solver
|
||||
problems
|
||||
end
|
||||
""" Default initializer. """
|
||||
function SimpleSolver()
|
||||
SimpleSolver(Problem[])
|
||||
end
|
||||
|
||||
"""
|
||||
Call solver to solve a set of problems.
|
||||
|
||||
This is simple serial solver for demonstration purposes. It handles the most
|
||||
common situation, i.e., some main field problem and it's Dirichlet boundary.
|
||||
"""
|
||||
function call(solver::SimpleSolver, t)
|
||||
problems = get_problems(solver)
|
||||
problem1 = problems[1]
|
||||
problem2 = problems[2]
|
||||
|
||||
# calculate order of degrees of freedom in global matrix
|
||||
# and set the ordering to problems
|
||||
dofmap = calculate_global_dofs(problem1)
|
||||
assign_global_dofs!(problem1, dofmap)
|
||||
assign_global_dofs!(problem2, dofmap)
|
||||
|
||||
# assemble problem 1
|
||||
A1 = sparse(get_lhs(problem1, t)...)
|
||||
b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))
|
||||
|
||||
# assemble problem 2
|
||||
A2 = sparse(get_lhs(problem2, t)...)
|
||||
b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))
|
||||
|
||||
# make one monolithic assembly
|
||||
A = [A1 A2; A2' zeros(A2)]
|
||||
b = [b1; b2]
|
||||
dump(full(A))
|
||||
dump(full(b))
|
||||
|
||||
# solve problem
|
||||
nz = unique(rowvals(A))
|
||||
x = zeros(b)
|
||||
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]
|
||||
|
||||
# check residual
|
||||
R1 = A1*x1 - b1
|
||||
R2 = A2*x2 - b2
|
||||
println("Residual norm: $(norm(R1+R2))")
|
||||
|
||||
# update field for elements in problem 1
|
||||
for equation in get_equations(problem1)
|
||||
gdofs = get_global_dofs(equation)
|
||||
element = get_element(equation)
|
||||
field_name = get_unknown_field_name(equation) # field we are solving
|
||||
field = Field(t, full(x1[gdofs])[:])
|
||||
add_field!(element, field_name, field)
|
||||
end
|
||||
|
||||
# update field for elements in problem 2
|
||||
for equation in get_equations(problem2)
|
||||
gdofs = get_global_dofs(equation)
|
||||
element = get_element(equation)
|
||||
field_name = get_unknown_field_name(equation)
|
||||
field = Field(t, full(x2[gdofs]))
|
||||
add_field!(element, field_name, field)
|
||||
end
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user