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

Merge branch 'master' of git://github.com/JuliaFEM/JuliaFEM.jl

Browse Source
This commit is contained in:
Olli Väinölä
2015-12-17 15:36:46 +02:00
parent 90caf8deb2 7d1ece1368
commit ab34acedf4
16 changed files with 993 additions and 371 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+6
View File
@@ -6,6 +6,8 @@ This is JuliaFEM -- Finite Element Package
"""
module JuliaFEM
# include("common.jl")
""" JuliaFEM Core module. """
module Core
include("core.jl")
@@ -19,6 +21,10 @@ 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")
end
src/common.jl
+52
View File
@@ -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
# common routines
"""
A very simple debugging macro. It executes commands if environment variable DEBUG is set.
Usage
-----
Instead of starting session `julia file.jl`, do `DEBUG=1 julia file.jl`.
Or set `export DEBUG=1` for your `.bashrc`.
Running inside code
-------------------
julia> @debug info("moimoi heihei")
will get executed iff environment variable DEBUG is set.
Examples
--------
julia> @debug info("moimoi")
(empty)
julia> ENV["DEBUG"] = 1
julia> @debug info("moimoi")
INFO: moimoi
julia> @debug begin
... info("matrix is")
... dump([1 2; 3 4])
... end
INFO: matrix is
Array(Int64(2,2)) 2x2 Array{Int64,2}:
1 2
3 4
"""
macro debug(msg)
haskey(ENV, "DEBUG") || return
return msg
end
function set_debug_on!()
ENV["DEBUG"] = 1;
end
function set_debug_off!()
pop!(ENV, "DEBUG");
end
export @debug, set_debug_on!, set_debug_off!
src/core.jl
+2 -23
View File
@@ -1,31 +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
import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof, append!
"""
A very simple debugging macro. It prints debug message if environment variable
JULIAFEM_DEBUG is found.
Usage: instead of starting session `julia file.jl` do `DEBUG=1 julia file.jl`.
Or set `export DEBUG=1` for your `.bashrc`.
"""
macro debug(msg)
haskey(ENV, "DEBUG") || return
# return :( println("DEBUG: ", $msg) )
return msg
end
function set_debug_on!()
ENV["DEBUG"] = 1;
end
function set_debug_off!()
pop!(ENV, "DEBUG");
end
export @debug, set_debug_on!, set_debug_off!
using ForwardDiff
autodiffcache = ForwardDiffCache()
# export derivative, jacobian, hessian
src/directsolver.jl
+1
View File
@@ -143,6 +143,7 @@ end
function call(solver::DirectSolver, time::Number=0.0)
info("# of field problems: $(length(solver.field_problems))")
info("# of boundary problems: $(length(solver.boundary_problems))")
(length(solver.field_problems) != 0) || error("no field problems defined.")
timing = Dict{ASCIIString, Float64}()
tic(timing, "solver")
src/dirichlet.jl
+5 -2
View File
@@ -3,8 +3,11 @@
abstract DirichletProblem <: AbstractProblem
function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[])
return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements)
function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
return BoundaryProblem{DirichletProblem}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements)
end
function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
return BoundaryProblem{DirichletProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
end
function assemble!(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Number)
src/elements.jl
+10
View File
@@ -264,3 +264,13 @@ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::
element["normal-tangential coordinates"] = ntcoords
end
""" Pick values from nodes and set to element according to connectivity. """
function update(element::Element, field_name::ASCIIString, data::Union{Vector, Dict})
element[field_name] = [data[i] for i in get_connectivity(element)]
end
function update(elements::Vector{Element}, field_name::ASCIIString, data::Union{Vector, Dict})
# info("update $field_name for $(length(elements)) elements.")
for element in elements
update(element, field_name, data)
end
end
src/integrate.jl
+13
View File
@@ -78,6 +78,19 @@ function get_integration_points(::TriangularElement, ::Type{Val{2}})
]
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 ..?
src/linear_elasticity.jl
+13 -4
View File
@@ -5,8 +5,8 @@
abstract LinearElasticityProblem <: ElasticityProblem
function LinearElasticityProblem(dim::Int=3, elements=[])
return Problem{LinearElasticityProblem}(dim, elements)
function LinearElasticityProblem(name="linear elasticity", dim::Int=3, elements=[])
return Problem{LinearElasticityProblem}(name, dim, elements)
end
""" Elasticity equations, general 3D case. """
@@ -57,13 +57,22 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
add!(assembly.force_vector, gdofs, vec(L))
end
for dim in 1:problem.dim
if haskey(element, "displacement traction force $dim")
T = element("displacement traction force $dim", ip, time)
ldofs = gdofs[dim:problem.dim:end]
JT = transpose(J)
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
add!(assembly.force_vector, ldofs, vec(L))
end
end
end
end
abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem
function PlaneStressLinearElasticityProblem(dim::Int=2, elements=[])
return Problem{PlaneStressLinearElasticityProblem}(dim, elements)
function PlaneStressLinearElasticityProblem(name="plane stress linear elasticity", dim::Int=2, elements=[])
return Problem{PlaneStressLinearElasticityProblem}(name, dim, elements)
end
""" Elasticity equations, plane stress. """
src/mortar.jl
+151 -35
View File
@@ -3,6 +3,13 @@
# Mortar projection calculation for 2d
macro debug(msg)
haskey(ENV, "DEBUG") || return
return msg
end
""" Find projection from slave nodes to master element, i.e. find xi2 from
master element corresponding to the xi1.
"""
@@ -617,7 +624,18 @@ function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix
return theta
end
end
error("project_point_to_auxiliary_plane: did not converge in $max_iterations iterations!")
begin
info("projecting point from auxiliary plane back to surface didn't go very well.")
info("element type: $E")
info("element connectivity: $(get_connectivity(element))")
info("auxiliary plane: x0 = $x0, Q = $Q")
info("point coordinates on plane: $p")
info("element geometry: $x")
info("ph: $ph")
info("normal direction: $n")
info("parameter vector before giving up: $theta")
end
error("project_point_to_surface: did not converge in $max_iterations iterations!")
end
@@ -632,8 +650,12 @@ node_csys
"""
abstract MortarProblem <: AbstractProblem
function MortarProblem(parent_field_name, parent_field_dim, dim=1, elements=[])
return BoundaryProblem{MortarProblem}(parent_field_name, parent_field_dim, dim, elements)
function MortarProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[])
return BoundaryProblem{MortarProblem}("mortar problem", parent_field_name, parent_field_dim, dim, elements)
end
function MortarProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[])
return BoundaryProblem{MortarProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
end
# Mortar assembly
@@ -688,6 +710,31 @@ end
typealias MortarElements3D Union{Tri3, Quad4}
""" Find master elements from list of potential master elements. """
function find_master_elements(slave_element::Element, time::Real)
x0, Q = create_auxiliary_plane(slave_element, time)
Sl = Vector{Float64}[]
for p in slave_element("geometry", time)
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
end
S = hcat(Sl...)
master_elements = Element[]
for master_element in slave_element["master elements"]
M = Vector{Float64}[]
for p in master_element("geometry", time)
push!(M, project_point_to_auxiliary_plane(p, x0, Q))
end
M = hcat(M...)
P, neighbours = clip_polygon(S, M)
isa(P, Void) && continue # no clipping
size(P, 2) < 3 && continue # shared edge, no contribution
push!(master_elements, master_element)
end
return master_elements
end
function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element{E}, time::Real)
field_dim = problem.parent_field_dim
field_name = problem.parent_field_name
@@ -702,11 +749,14 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
for p in slave_element("geometry", time)
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
end
@debug info("auxiliary plane coords and basis: origo = $x0")
@debug info("basis:")
@debug dump(round(Q, 3))
@debug begin
info("auxiliary plane coords and basis: origo = $x0")
info("basis:")
dump(round(Q, 3))
end
#S = reshape([S...;], 2, size(slave_element)[2])
S = hcat(Sl...)
slave_geom = Field(Vector{Float64}[S[:,j] for j=1:size(S,2)])
for master_element in slave_element["master elements"]
master_dofs = get_gdofs(master_element, field_dim)
@@ -717,11 +767,16 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
end
#M = reshape([M...;], 2, size(master_element)[2])
M = hcat(M...)
master_geom = Field(Vector{Float64}[M[:,j] for j=1:size(M,2)])
P = nothing
neighbours = nothing
@debug info("applying polygon clip algorithm, S & M = ")
@debug dump(round(S, 3))
@debug dump(round(M, 3))
@debug begin
info("applying polygon clip algorithm, S & M = ")
dump(round(S, 3))
dump(round(M, 3))
end
try
P, neighbours = clip_polygon(S, M)
catch
@@ -735,53 +790,114 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
error("cannot continue")
end
isa(P, Void) && continue # no clipping
@debug info("polygon on auxilyary plane: ")
@debug dump(round(P, 3))
C = calculate_polygon_centerpoint(P)
@debug info("center point = $C")
@debug begin
info("polygon coords on auxilyary plane: ")
dump(round(P, 3))
end
if size(P, 2) < 3
# shared edge but no shared volume. skipping
continue
info("this is not polygon at all.")
info("clipping S")
dump(S)
info("clipping M")
dump(M)
error("size(P, 2) < 3")
end
C = calculate_polygon_centerpoint(P)
npts = size(P, 2) # number of vertices in polygon
@debug info("number of vectices in polygon: $npts")
# S = zeros(3, 3)
# M = zeros(3, 3)
@debug begin
info("clip polygon info")
theta = project_point_from_plane_to_surface(C, x0, Q, slave_element, time)
CC = slave_element("geometry", theta[2:3], time)
info("center point on slave: $CC")
info("number of vectices in polygon: $npts")
on_slave = zeros(3, 0)
on_master = zeros(3, 0)
for i=1:size(P, 2)
theta = project_point_from_plane_to_surface(P[:,i], x0, Q, slave_element, time)
on_slave = [on_slave slave_element("geometry", theta[2:3], time)]
theta = project_point_from_plane_to_surface(P[:,i], x0, Q, master_element, time)
on_master = [on_master master_element("geometry", theta[2:3], time)]
end
info("polygon coords projected to slave element")
dump(round(on_slave, 3))
info("polygon coords projected to master element")
dump(round(on_master, 3))
end
for i=1:npts # loop vertices and create temporary integrate cells
xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]]
yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]]
X = hcat(xvec, yvec)'
@debug info("cell $i, coords = ")
@debug dump(round(X, 3))
geom = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)])
@debug begin
on_slave = zeros(3, 0)
on_master = zeros(3, 0)
for j=1:size(X, 2)
theta = project_point_from_plane_to_surface(X[:,j], x0, Q, slave_element, time)
on_slave = [on_slave slave_element("geometry", theta[2:3], time)]
theta = project_point_from_plane_to_surface(X[:,j], x0, Q, master_element, time)
on_master = [on_master master_element("geometry", theta[2:3], time)]
end
info("cell $i coords projected to slave element")
dump(round(on_slave, 3))
info("cell $i coords projected to master element")
dump(round(on_master, 3))
end
# integration cell geometry, i.e., Tri3
cell = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)])
# info("geom = $geom")
for ip in get_integration_points(Tri3, Val{5})
# calculate determiant of jacobian
#dN = get_dbasis(E, ip.xi)
dN = get_dbasis(Tri3, ip.xi)
J = sum([kron(dN[:,j], geom[j]') for j=1:length(geom)])
w = ip.weight*det(J)
# gauss point in auxiliary plane
#N = get_basis(E, ip.xi)
N = get_basis(Tri3, ip.xi)
x = vec(N*geom)
xi = vec(N*cell) # xi defined in auxilary plane
#xi = ip.xi
# info("x = $x")
# find projection of gauss point to master and slave elements
theta1 = project_point_from_plane_to_surface(x, x0, Q, slave_element, time)
theta2 = project_point_from_plane_to_surface(x, x0, Q, master_element, time)
theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time)
xi_slave = theta1[2:3]
xi_master = theta2[2:3]
@debug begin
X_slave = slave_element("geometry", xi_slave, time)
X_master = master_element("geometry", xi_master, time)
info("integration point on slave: $xi_slave => $X_slave")
info("integration point on master: $xi_master => $X_master")
end
# evaluate shape functions values in gauss point and add contribution to matrices
N1 = slave_element(theta1[2:3], time)
N2 = master_element(theta2[2:3], time)
Sm = w*N1'*N1
N1 = slave_element(xi_slave, time)
N2 = master_element(xi_master, time)
# calculate determiant of jacobian
dNC = get_dbasis(Tri3, ip.xi)
dNS = get_dbasis(Quad4, xi_slave)
dNM = get_dbasis(Quad4, xi_master)
JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)])
JN = sum([kron(dNS[:,j], slave_geom[j]') for j=1:length(slave_geom)])
JM = sum([kron(dNM[:,j], master_geom[j]') for j=1:length(master_geom)])
wS = det(JN)
wM = det(JM)
wC = det(JC)
@debug info("weight S = $wS, weight M = $wM, weight C = $wC")
Sm = ip.weight*N1'*N1*wC
# FIXME: master side transpose -- why?
Mm = w*(N1'*N2)'
Mm = ip.weight*(N1'*N2)'*wC
for k=1:field_dim
sd = slave_dofs[k:field_dim:end]
md = master_dofs[k:field_dim:end]
add!(assembly.stiffness_matrix, sd, sd, Sm)
add!(assembly.stiffness_matrix, sd, md, -Mm)
# info("sd = $sd")
# info("md = $md")
end
end
# info("breaking on first")
# break
end
# info("S = \n$S")
# info("M = \n$M")
end
end
src/preprocess_aster_reader.jl
+140
View File
@@ -1,6 +1,8 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using HDF5
function aster_parse_nodes(section::ASCIIString; strip_characters=true)
@@ -53,6 +55,144 @@ function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}})
return model
end
function aster_renumber_nodes_!(mesh, node_numbering)
old_nodes = mesh["nodes"]
new_nodes = typeof(old_nodes)()
for (node_id, node_coords) in old_nodes
new_node_id = node_numbering[node_id]
new_nodes[new_node_id] = node_coords
end
mesh["nodes"] = new_nodes
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
new_elcon = [node_numbering[node_id] for node_id in elcon]
mesh["connectivity"][elid] = (eltype, elset, new_elcon)
end
end
function aster_renumber_nodes!(mesh1, mesh2)
reserved_node_ids = Set(collect(keys(mesh1["nodes"])))
@debug info("already reserved node ids: $reserved_node_ids")
mesh2_node_numbering = Dict{Int64, Int64}()
# find new node ids assigned for mesh 2
k = 1
for node_id in sort(collect(keys(mesh2["nodes"])))
@debug info("mesh2: processing node $node_id")
# if node id is reserved in mesh 1, find new number
if node_id in reserved_node_ids
@debug info("node id conflict, $node_id already defined in mesh 1, renumbering")
while k in reserved_node_ids
k += 1
end
@debug info("mesh2: node $node_id -> $k")
mesh2_node_numbering[node_id] = k
push!(reserved_node_ids, k)
else
mesh2_node_numbering[node_id] = node_id
end
end
@debug info("new node numering:")
@debug println(mesh2_node_numbering)
aster_renumber_nodes_!(mesh2, mesh2_node_numbering)
#=
# create new nodes
mesh2_old_nodes = mesh2["nodes"]
mesh2_new_nodes = typeof(mesh2_old_nodes)()
for (node_id, node_coords) in mesh2_old_nodes
new_node_id = mesh2_node_numbering[node_id]
mesh2_new_nodes[new_node_id] = node_coords
end
mesh2["nodes"] = mesh2_new_nodes
# update connectivity
for (elid, (eltype, elset, elcon)) in mesh2["connectivity"]
new_elcon = [mesh2_node_numbering[node_id] for node_id in elcon]
mesh2["connectivity"][elid] = (eltype, elset, new_elcon)
end
=#
end
function aster_renumber_elements!(mesh1, mesh2)
reserved_element_ids = Set(collect(keys(mesh1["connectivity"])))
@debug info("already reserved element ids: $reserved_element_ids")
mesh2_element_numbering = Dict{Int64, Int64}()
# find new element ids assigned for mesh 2
k = 1
for element_id in sort(collect(keys(mesh2["connectivity"])))
@debug info("mesh2: processing element $element_id")
# if node id is reserved in mesh 1, find new number
if element_id in reserved_element_ids
@debug info("element id conflict, $element_id already defined in mesh 1, renumbering")
while k in reserved_element_ids
k += 1
end
@debug info("mesh2: element $element_id -> $k")
mesh2_element_numbering[element_id] = k
push!(reserved_element_ids, k)
else
mesh2_element_numbering[element_id] = element_id
end
end
@debug info("element numbering for mesh 2:")
@debug info(mesh2_element_numbering)
# create new elements
mesh2_old_elements = mesh2["connectivity"]
mesh2_new_elements = typeof(mesh2_old_elements)()
for (element_id, element_data) in mesh2_old_elements
new_element_id = mesh2_element_numbering[element_id]
mesh2_new_elements[new_element_id] = element_data
end
mesh2["connectivity"] = mesh2_new_elements
end
function aster_combine_meshes(mesh1, mesh2)
# check that meshes are ready to be combined
node_ids_mesh_1 = collect(keys(mesh1["nodes"]))
node_ids_mesh_2 = collect(keys(mesh2["nodes"]))
if length(intersect(node_ids_mesh_1, node_ids_mesh_2)) != 0
error("nodes with same id number in both meshes, failed.")
end
element_ids_mesh_1 = collect(keys(mesh1["connectivity"]))
element_ids_mesh_2 = collect(keys(mesh2["connectivity"]))
if length(intersect(element_ids_mesh_1, element_ids_mesh_2)) != 0
error("elements with same id number in both meshes, failed.")
end
@assert similar(mesh1) == similar(mesh2)
@assert similar(mesh1["nodes"]) == similar(mesh2["nodes"])
@assert similar(mesh1["connectivity"]) == similar(mesh2["connectivity"])
new_mesh = similar(mesh1)
new_mesh["nodes"] = similar(mesh1["nodes"])
new_mesh["connectivity"] = similar(mesh1["connectivity"])
for (node_id, node_coords) in mesh1["nodes"]
new_mesh["nodes"][node_id] = node_coords
end
for (node_id, node_coords) in mesh2["nodes"]
new_mesh["nodes"][node_id] = node_coords
end
for (element_id, element_data) in mesh1["connectivity"]
new_mesh["connectivity"][element_id] = element_data
end
for (element_id, element_data) in mesh2["connectivity"]
new_mesh["connectivity"][element_id] = element_data
end
return new_mesh
end
"""
Code Aster binary file (.med), which is exported from SALOME.
"""
src/problems.jl
+3
View File
@@ -4,11 +4,13 @@
abstract AbstractProblem
type Problem{T<:AbstractProblem}
name :: ASCIIString
dim :: Int
elements :: Vector{Element}
end
type BoundaryProblem{T<:AbstractProblem}
name :: ASCIIString
parent_field_name :: ASCIIString
parent_field_dim :: Int
dim :: Int
@@ -36,3 +38,4 @@ end
function Base.push!(problem::AllProblems, element::Element)
push!(problem.elements, element)
end
src/test.jl
+2 -3
View File
@@ -1,7 +1,7 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using Base.Test
using BaseTestNext
abstract TestResult
@@ -159,5 +159,4 @@ function print_test_statistics()
return passed, failed, errors, critical
end
export @test, run_test, print_test_statistics
export @test, @testset, @test_throws, run_test, print_test_statistics