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

some tests pass now

Browse Source
This commit is contained in:
Jukka Aho
2016-05-22 02:34:38 +03:00
parent 283621f4b1
commit 53d640f477
9 changed files with 213 additions and 217 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+8 -8
View File
@@ -23,9 +23,9 @@ export DCTI
### ELEMENTS ###
include("elements.jl") # common element routines
export Element
export Element, update!
include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
export Quad4
export Seg2, Tri3, Quad4, Hex8, Tet4
#include("hierarchical.jl") # P-elements
#include("mortar_elements.jl") # Mortar elements
@@ -41,13 +41,18 @@ include("elasticity.jl") # elasticity equations
export Elasticity
include("dirichlet.jl")
export Dirichlet
include("heat.jl")
export Heat
export assemble
### ASSEMBLY + SOLVE ###
include("assembly.jl")
include("solver_utils.jl")
include("solvers.jl")
export Solver
### MORTAR STUFF ###
include("mortar.jl") # mortar projection
@@ -64,14 +69,9 @@ include("api.jl")
end
module Preprocess
#=
macro debug(msg)
haskey(ENV, "DEBUG") || return
return msg
end
=#
include("abaqus_reader.jl")
include("preprocess_aster_reader.jl")
export aster_create_elements, parse_aster_med_file
end
module Postprocess
src/dirichlet.jl
+5 -6
View File
@@ -43,26 +43,25 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
end
# right hand side
for ip in get_integration_points(element, Val{3})
w = ip.weight
J = get_jacobian(element, ip, time)
for (w, xi) in get_integration_points(element, Val{3})
J = element(xi, time, Val{:Jacobian})
JT = transpose(J)
if size(JT, 2) == 1 # plane problem
w *= norm(JT)
else
w *= norm(cross(JT[:,1], JT[:,2]))
end
N = element(ip, time)
N = element(xi, time)
for i=1:field_dim
ldofs = gdofs[i:field_dim:end]
if haskey(element, field_name*" $i")
g = element(field_name*" $i", ip, time)
g = element(field_name*" $i", xi, time)
if get_formulation_type(problem) == :incremental
# if having incremental formulation need to add previous
# displacement to rhs (solving increment Δu !
haskey(element, "displacement") || continue
g_prev = element(field_name, ip, time)
g_prev = element(field_name, xi, time)
g -= g_prev[i]
end
add!(assembly.g, ldofs, w*g*Ae*N')
src/elasticity.jl
+30 -17
View File
@@ -44,11 +44,15 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El
end
function assemble(problem::Problem{Elasticity}, element::Element, time=0.0)
assemble(problem, element, time, Val{problem.properties.formulation})
problem.properties
if problem.properties.formulation in [:plane_stress, :plane_strain]
return assemble(problem, element, time, Val{:plane})
end
return assemble(problem, element, time, Val{problem.properties.formulation})
end
""" Elasticity equations for 2d cases. """
function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane_stress}})
function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
props = problem.properties
dim = get_unknown_field_dimension(problem)
@@ -82,10 +86,19 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
# get_material(problem, element, ...)
E = element("youngs modulus", xi, time)
nu = element("poissons ratio", xi, time)
D = E/(1.0 - nu^2) .* [
1.0 nu 0.0
nu 1.0 0.0
0.0 0.0 (1.0-nu)/2.0]
if props.formulation == :plane_stress
D = E/(1.0 - nu^2) .* [
1.0 nu 0.0
nu 1.0 0.0
0.0 0.0 (1.0-nu)/2.0]
elseif props.formulation == :plane_strain
D = E/((1+nu)*(1-2*nu)) .* [
1-nu nu 0
nu 1-nu 0
0 0 (1-2*nu)/2]
else
error("unknown plane formulation: $(props.formulation)")
end
# calculate stress
S = D*[strain[1,1]; strain[2,2]; 2*strain[1,2]]
@@ -138,30 +151,30 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E
Kt = zeros(dim*nnodes, dim*nnodes)
f = zeros(dim*nnodes)
for ip in get_integration_points(element)
for (w, xi) in get_integration_points(element)
J = get_jacobian(element, ip, time)
N = element(ip, time)
w = ip.weight*norm(J)
J = element(xi, time, Val{:Jacobian})
detJ = norm(J)
N = element(xi, time)
if haskey(element, "displacement traction force")
T = element("displacement traction force", ip, time)
f += vec(w*T*N)
T = element("displacement traction force", xi, time)
f += w*vec(T*N)*detJ
end
for i=1:dim
# traction force for ith component
if haskey(element, "displacement traction force $i")
T = element("displacement traction force $i", ip, time)
f[i:dim:end] += vec(w*T*N)
T = element("displacement traction force $i", xi, time)
f[i:dim:end] += w*vec(T*N)*detJ
end
end
if haskey(element, "nt displacement traction force")
# traction force given in normal-tangential direction
T = element("nt displacement traction force", ip, time)
Q = element("normal-tangential coordinates", ip, time)
f += vec(w*Q'*T*N)
T = element("nt displacement traction force", xi, time)
Q = element("normal-tangential coordinates", xi, time)
f += w*vec(Q'*T*N)*detJ
end
end
src/elements.jl
+41 -39
View File
@@ -74,6 +74,10 @@ function size{E}(element::Element{E})
size(element.properties)
end
function size(element::Element, dim::Int)
size(element)[dim]
end
""" Update element field based on a dictionary of nodal data and connectivity information.
Examples
@@ -89,10 +93,16 @@ function update!(element::Element, field_name::ASCIIString, data::Dict)
element[field_name] = [data[i] for i in get_connectivity(element)]
end
function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair}...)
function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair})
element[field_name] = data
end
function update!(elements::Vector, field_name::ASCIIString, data)
for element in elements
update!(element, field_name, data)
end
end
""" Evaluate partial derivatives of basis functions using ForwardDiff. """
function get_dbasis(element::Element, xi::Vector, time)
basis(xi) = vec(get_basis(element, xi, time))
@@ -104,6 +114,36 @@ function haskey(element::Element, field_name)
haskey(element.fields, field_name)
end
function get_connectivity(element::Element)
return element.connectivity
end
function get_gdofs(element::Element)
return get_gdofs(element, 1)
end
""" Return dual basis transformation matrix Ae. """
function get_dualbasis(element::Element, time)
nnodes = length(element)
De = zeros(nnodes, nnodes)
Me = zeros(nnodes, nnodes)
for (w, xi) in get_integration_points(element, Val{3})
J = element(xi, time, Val{:Jacobian})
JT = transpose(J)
if size(JT, 2) == 1 # plane problem
# || ∂X/∂ξ ||
w *= norm(JT)
else
# || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
w *= norm(cross(JT[:,1], JT[:,2]))
end
N = element(xi, time)
De += w*diagm(vec(N))
Me += w*N'*N
end
return De, Me, De*inv(Me)
end
#=
type Element{E}
@@ -172,9 +212,6 @@ function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
element.fields[name] = Field(data...)
end
function get_connectivity(el::Element)
return el.connectivity
end
typealias VecOrIP Union{Vector, IntegrationPoint}
@@ -248,9 +285,6 @@ function find_elements(elements, nodes)
return collect(s)
end
function get_gdofs(element::Element)
return get_gdofs(element, 1)
end
function get_dbasis{E}(element::Element{E}, ip::IntegrationPoint)
return get_dbasis(E, ip.xi)
@@ -260,33 +294,6 @@ function get_basis{E, T<:Real}(element::Element{E}, xi::T)
return get_basis(E, xi)
end
""" Return dual basis transformation matrix Ae. """
function get_dualbasis(element::Element, time::Real)
if length(element.A) == 0
nnodes = size(element, 2)
De = zeros(nnodes, nnodes)
Me = zeros(nnodes, nnodes)
for ip in get_integration_points(element, Val{3})
w = ip.weight
J = get_jacobian(element, ip, time)
JT = transpose(J)
if size(JT, 2) == 1 # plane problem
# || ∂X/∂ξ ||
w *= norm(JT)
else
# || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
w *= norm(cross(JT[:,1], JT[:,2]))
end
N = element(ip, time)
De += w*diagm(vec(N))
Me += w*N'*N
end
element.D = De
element.M = Me
element.A = De*inv(Me)
end
return element.D, element.M, element.A
end
function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}})
De, Me, Ae = get_dualbasis(element, time)
@@ -486,10 +493,5 @@ function update!{T}(elements::Vector{Element{T}}, field_name::ASCIIString, data.
update!(element, field_name, data...)
end
end
function update!(elements::Vector{Element}, field_name::ASCIIString, data...)
for element in elements
update!(element, field_name, data...)
end
end
=#
src/equations.jl
-13
View File
@@ -41,19 +41,6 @@ function has_residual_vector(problem::Problem, element::Element)
return method_exists(get_residual_vector, default_args)
end
function get_gdofs(element::Element, dim::Int)
conn = get_connectivity(element)
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
return gdofs
end
function get_gdofs(element::Element, problem::Problem)
return get_gdofs(element, problem.dimension)
end
function get_gdofs(problem::Problem, element::Element)
return get_gdofs(element, problem.dimension)
end
""" Assemble element. """
function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number)
src/integrate.jl
+112 -129
View File
@@ -1,123 +1,9 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# Let's drop here all integration schemes and some defaults for different element types
# maybe parse from txt file ..?
### 1d elements
typealias LineElement Union{Type{Seg2}, Type{Seg3}}
function get_integration_points(::LineElement, ::Type{Val{1}})
[
IntegrationPoint([0.0], 2.0)
]
end
function get_integration_points(::LineElement, ::Type{Val{2}})
[
IntegrationPoint([-sqrt(1/3)], 1)
IntegrationPoint([+sqrt(1/3)], 1)
]
end
function get_integration_points(::LineElement, ::Type{Val{3}})
[
IntegrationPoint([0.0], 8/9),
IntegrationPoint([-sqrt(3/5)], 5/9),
IntegrationPoint([+sqrt(3/5)], 5/9)
]
end
function get_integration_points(::LineElement, ::Type{Val{4}})
[
IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
]
end
function get_integration_points(::LineElement, ::Type{Val{5}})
[
IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
IntegrationPoint([0.0], 128/225),
IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900)
]
end
function get_integration_points(::Type{Seg2})
return get_integration_points(Seg2, Val{2})
end
function get_integration_points(::Type{Seg3})
return get_integration_points(Seg3, Val{3})
end
### 2d triangular elements
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
typealias TriangularElement Union{Type{Tri3}, Type{Tri6}}
function get_integration_points(::TriangularElement, ::Type{Val{1}})
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
[
IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5)
]
end
function get_integration_points(::TriangularElement, ::Type{Val{2}})
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
[
IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0),
IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0),
IntegrationPoint([1.0/6.0, 1.0/6.0], 1.0/6.0)
]
end
function get_integration_points(::TriangularElement, ::Type{Val{3}})
[
IntegrationPoint([1/3, 1/3], 0.5*-0.5625),
IntegrationPoint([0.2, 0.2], 0.5*0.5208333333333333),
IntegrationPoint([0.2, 0.6], 0.5*0.5208333333333333),
IntegrationPoint([0.6, 0.2], 0.5*0.5208333333333333),
]
end
function get_integration_points(::TriangularElement, ::Type{Val{4}})
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
# FIXME: something wrong here with weights ..?
[
IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801),
IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801),
IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801),
IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532),
IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532),
IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532)
]
end
function get_integration_points(::TriangularElement, ::Type{Val{5}})
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
# FIXME: something wrong here with weights ..?
[
IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000),
IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851),
IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851),
IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851),
IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483),
IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483),
IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483)
]
end
function get_integration_points(::Type{Tri3})
return get_integration_points(Tri3, Val{1})
end
# Let's drop here all integration schemes and some defaults for different element types maybe parse from txt file ..?
### Gauss quadrature rules for one dimension
function get_integration_points(::Type{Val{1}})
return [2.0], [0.0]
@@ -161,8 +47,30 @@ function get_integration_points(::Type{Val{5}})
return weights, points
end
function get_integration_points(element::Quad4)
return get_integration_points(element, Val{2})
### "cartesian" elements, integration rules comes from tensor product
### 1d elements
typealias LineElement Union{Seg2, Seg3}
function get_integration_points(element::LineElement, ::Type{Val{1}})
w, xi = get_integration_points(Val{1})
[ (w[i], [xi[i]]) for i=1:1 ]
end
function get_integration_points(element::LineElement, ::Type{Val{2}})
w, xi = get_integration_points(Val{2})
[ (w[i], [xi[i]]) for i=1:2 ]
end
function get_integration_points(element::LineElement, ::Type{Val{3}})
w, xi = get_integration_points(Val{3})
[ (w[i], [xi[i]]) for i=1:3 ]
end
function get_integration_points{E<:LineElement}(element::Element{E}, ::Type{Val{3}})
w, xi = get_integration_points(Val{3})
[ (w[i], [xi[i]]) for i=1:3 ]
end
function get_integration_points(element::Quad4, ::Type{Val{2}})
@@ -175,19 +83,94 @@ function get_integration_points(element::Quad4, ::Type{Val{3}})
[ (w[i]*w[j], [xi[i], xi[j]]) for i=1:3, j=1:3 ]
end
function get_integration_points(element::Hex8, ::Type{Val{2}})
w, xi = get_integration_points(Val{2})
[ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:2, j=1:2, k=1:2 ]
end
### default number of integration points for each element
function get_integration_points(element::Seg2)
get_integration_points(element, Val{2})
end
function get_integration_points(element::Seg3)
get_integration_points(element, Val{3})
end
function get_integration_points(element::Quad4)
get_integration_points(element, Val{2})
end
function get_integration_points(element::Hex8)
get_integration_points(element, Val{2})
end
### triangular and tetrahedral elements
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
typealias TriangularElement Union{Type{Tri3}, Type{Tri6}}
function get_integration_points(::TriangularElement, ::Type{Val{1}})
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
weights = [0.5]
points = Vector{Float64}[1.0/3.0*[1.0, 1.0]]
return weights, points
end
function get_integration_points(::TriangularElement, ::Type{Val{2}})
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
weights = 1.0/6.0*[1.0, 1.0, 1.0]
points = Vector{Float64}[
[2.0/3.0, 1.0/6.0],
[1.0/6.0, 2.0/3.0],
[1.0/6.0, 1.0/6.0]]
return weights, points
end
function get_integration_points(::TriangularElement, ::Type{Val{3}})
[
IntegrationPoint([1/3, 1/3], 0.5*-0.5625),
IntegrationPoint([0.2, 0.2], 0.5*0.5208333333333333),
IntegrationPoint([0.2, 0.6], 0.5*0.5208333333333333),
IntegrationPoint([0.6, 0.2], 0.5*0.5208333333333333),
]
end
function get_integration_points(::TriangularElement, ::Type{Val{4}})
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
# FIXME: something wrong here with weights ..?
[
IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801),
IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801),
IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801),
IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532),
IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532),
IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532)
]
end
function get_integration_points(::TriangularElement, ::Type{Val{5}})
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
# FIXME: something wrong here with weights ..?
[
IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000),
IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851),
IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851),
IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851),
IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483),
IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483),
IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483)
]
end
function get_integration_points(::Type{Tri3})
return get_integration_points(Tri3, Val{1})
end
### 3d elements
function get_integration_points(::Type{Hex8}, ::Type{Val{2}})
p = 1.0/sqrt(3.0)*[-1.0, 1.0]
w = [1.0, 1.0]
return vec([IntegrationPoint([p[i], p[j], p[k]], w[i]*w[j]) for i=1:2, j=1:2, k=1:2])
end
function get_integration_points(::Type{Hex8})
return get_integration_points(Hex8, Val{2})
end
function get_integration_points(::Type{Tet4})
# http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
[
src/preprocess_aster_reader.jl
+1 -2
View File
@@ -4,7 +4,6 @@
using HDF5
using JuliaFEM
# TODO: this should be elsewhere
function aster_create_elements(mesh, element_set, element_type=nothing; reverse_connectivity=false)
elements = Element[]
mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8, :TE4 => Tet4)
@@ -25,7 +24,7 @@ function aster_create_elements(mesh, element_set, element_type=nothing; reverse_
if reverse_connectivity
elcon = reverse(elcon)
end
element = mapping[eltype](elcon)
element = Element(mapping[eltype], elcon)
push!(elements, element)
end
update!(elements, "geometry", mesh["nodes"])
src/problems.jl
+14
View File
@@ -267,6 +267,20 @@ function push!(problem::Problem, element)
push!(problem.elements, element)
end
function get_gdofs(element::Element, dim::Int)
conn = get_connectivity(element)
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
return gdofs
end
function get_gdofs(element::Element, problem::Problem)
return get_gdofs(element, problem.dimension)
end
function get_gdofs(problem::Problem, element::Element)
return get_gdofs(element, problem.dimension)
end
""" Find dofs corresponding to nodes. """
function find_dofs_by_nodes(problem::Problem, nodes)
dim = get_unknown_field_dimension(problem)
test/test_elasticity_2d_linear_with_surface_load.jl
+2 -3
View File
@@ -1,11 +1,10 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Test
using JuliaFEM.Preprocess: aster_create_elements, parse_aster_med_file
using JuliaFEM.Core: Problem, Elasticity, Dirichlet, Solver, update!
@testset "test 2d linear elasticity with surface load" begin
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile)