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

modifications

Browse Source
This commit is contained in:
Jukka Aho
2016-05-25 02:44:29 +03:00
parent 043ba38d31
commit 070e42c7fc
5 changed files with 106 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+4 -3
View File
@@ -15,7 +15,7 @@ autodiffcache = ForwardDiffCache()
include("common.jl")
include("fields.jl")
export DCTI
export DCTI, Field
#include("basis.jl") # interpolation of discrete fields
#include("symbolic.jl") # a thin symbolic layer for fields
#include("types.jl") # type definitions
@@ -36,7 +36,8 @@ export get_integration_points
include("sparse.jl")
include("problems.jl") # common problem routines
export Problem
export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
get_unknown_field_dimension
include("elasticity.jl") # elasticity equations
export Elasticity
@@ -79,7 +80,7 @@ module Postprocess
include("xdmf.jl")
export xdmf_new_temporal_collection, xdmf_new_grid,
xdmf_new_mesh!, xdmf_new_nodal_field!,
xdmf_save_model, xdmf_new_model
xdmf_save_model, xdmf_new_model, xdmf_dump
end
""" JuliaFEM testing routines. """
src/dirichlet.jl
+19 -15
View File
@@ -22,16 +22,28 @@ function get_formulation_type(problem::Problem{Dirichlet})
return problem.properties.formulation
end
function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time::Real)
@assert problem.properties.dual_basis
function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time)
# get dimension and name of PARENT field
nnodes = length(element)
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
gdofs = get_gdofs(element, field_dim)
De, Me, Ae = get_dualbasis(element, time)
# if problem.properties.formulation == :dual_basis
De, Me, Ae = get_dualbasis(element, time)
# else
# Ae = eye(nnodes)
# De = zeros(nnodes, nnodes)
# for (w, xi) in get_integration_points(element, Val{3})
# N = element(xi, time)
# detJ = element(xi, time, Val{:detJ})
# De += w*N'*N*detJ
# end
# end
# De = Ae = eye(nnodes)
# left hand side
for i=1:field_dim
@@ -44,27 +56,19 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
# right hand side
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
detJ = element(xi, time, Val{:detJ})
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", xi, time)
if get_formulation_type(problem) == :incremental
# if having incremental formulation need to add previous
# displacement to rhs (solving increment Δu !
if true
haskey(element, "displacement") || continue
g_prev = element(field_name, xi, time)
g -= g_prev[i]
end
add!(assembly.g, ldofs, w*g*Ae*N')
add!(assembly.g, ldofs, w*g*Ae*N'*detJ)
end
end
src/elasticity.jl
+29 -12
View File
@@ -27,6 +27,7 @@ end
function get_formulation_type(problem::Problem{Elasticity})
# we are solving residual and add increment to previous solution vector
return :incremental
#return :total
end
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real)
@@ -129,13 +130,22 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
if props.finite_strain # add geometric stiffness
Kt += w*BNL'*S2*BNL*detJ # geometric stiffness
end
f -= w*BL'*S*detJ # internal force
if get_formulation_type(problem) == :incremental
f -= w*BL'*S*detJ # internal force
end
# volume load
if haskey(element, "displacement load")
b = element("displacement load", xi, time)
f += w*vec(N'*b)*detJ
end
for i=1:dim
if haskey(element, "displacement load $i")
b = element("displacement load $i", xi, time)
f[i:dim:end] += w*vec(b*N)*detJ
end
end
end
@@ -213,16 +223,13 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
E = element("youngs modulus", xi, time)
nu = element("poissons ratio", xi, time)
a = 1 - nu
b = 1 - 2*nu
c = 1 + nu
D = E/(b*c) .* [
a nu nu 0 0 0
nu a nu 0 0 0
nu nu a 0 0 0
0 0 0 b 0 0
0 0 0 0 b 0
0 0 0 0 0 b]
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
1.0-nu nu nu 0.0 0.0 0.0
nu 1.0-nu nu 0.0 0.0 0.0
nu nu 1.0-nu 0.0 0.0 0.0
0.0 0.0 0.0 0.5-nu 0.0 0.0
0.0 0.0 0.0 0.0 0.5-nu 0.0
0.0 0.0 0.0 0.0 0.0 0.5-nu]
# # PK2 stress tensor in voigt notation
S = D*[strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[2,3]; 2*strain[1,3]; 2*strain[1,2]]
@@ -249,6 +256,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
end
fill!(BNL, 0.0)
for i=1:size(dN, 2)
BNL[1, 3*(i-1)+1] = dN[1,i]
@@ -274,13 +282,22 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
if props.finite_strain
Kt += w*BNL'*S3*BNL*detJ
end
f -= w*BL'*S*detJ
if get_formulation_type(problem) == :incremental
f -= w*BL'*S*detJ
end
# volume load
if haskey(element, "displacement load")
T = element("displacement load", ip, time)
f += w*vec(T*N)*detJ
end
for i=1:dim
if haskey(element, "displacement load $i")
b = element("displacement load $i", xi, time)
f[i:dim:end] += w*vec(b*N)*detJ
end
end
end
return Kt, f
src/elements.jl
+7 -7
View File
@@ -6,13 +6,14 @@ abstract AbstractElement
typealias Node Vector{Float64}
type Element{E<:AbstractElement}
id :: Int
connectivity :: Vector{Int}
fields :: Dict{ASCIIString, Field}
properties :: E
end
function Element{E<:AbstractElement}(::Type{E}, connectivity=[], fields=Dict(), properties...)
Element{E}(connectivity, fields, E(properties...))
function Element{E<:AbstractElement}(::Type{E}, connectivity=[], id=-1, fields=Dict(), properties...)
Element{E}(id, connectivity, fields, E(properties...))
end
function getindex(element::Element, field_name::ASCIIString)
@@ -59,12 +60,8 @@ function call(element::Element, xi::Vector, time, ::Type{Val{:detJ}})
end
end
function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
element(xi, time, Val{:Jacobian})
end
function call(element::Element, xi::Vector, time, ::Type{Val{:Grad}})
J = get_jacobian(element, xi, time)
J = element(xi, time, Val{:Jacobian})
return inv(J)*get_dbasis(element, xi, time)
end
@@ -72,6 +69,9 @@ function call(element::Element, field_name, xi::Vector, time, ::Type{Val{:Grad}}
element(xi, time, Val{:Grad})*element[field_name](time)
end
#function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
# element(xi, time, Val{:Jacobian})
#end
#function get_basis(element::Element, xi::Vector, time=0.0)
# get_basis(element.properties, xi, time)
#end
src/xdmf.jl
+47
View File
@@ -2,6 +2,7 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using LightXML
using JuliaFEM
# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
# > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
@@ -124,3 +125,49 @@ function xdmf_save_model(xdoc, filename)
save_file(xdoc, filename)
end
function xdmf_dump(all_elements, eltype, elsym, time=0.0, filename="/tmp/xdmf_result.xmf")
info("$(length(all_elements)) elements.")
xdoc, xmodel = xdmf_new_model()
coll = xdmf_new_temporal_collection(xmodel)
grid = xdmf_new_grid(coll; time=time)
Xg = Dict{Int64, Vector{Float64}}()
ug = Dict{Int64, Vector{Float64}}()
nids = Dict{Int64, Int64}()
for element in all_elements
conn = get_connectivity(element)
for (i, c) in enumerate(conn)
nids[c] = c
end
X = element("geometry", time)
for (i, c) in enumerate(conn)
Xg[c] = X[i]
end
haskey(element, "displacement") || continue
u = element("displacement", time)
for (i, c) in enumerate(conn)
ug[c] = u[i]
end
end
perm = sort(collect(keys(Xg)))
nodes = Vector{Float64}[Xg[i] for i in perm]
disp = Vector{Float64}[ug[i] for i in perm]
nids = Int[nids[i] for i in perm]
inids = Dict{Int64, Int64}()
for (i, nid) in enumerate(nids)
inids[nid] = i
end
elements = []
for element in all_elements
isa(element, eltype) || continue
conn = get_connectivity(element)
nconn = [inids[i] for i in conn]
push!(elements, (elsym, nconn))
end
xdmf_new_mesh!(grid, nodes, elements)
xdmf_new_nodal_field!(grid, "displacement", disp)
xdmf_save_model(xdoc, filename)
info("model dumped to $filename")
end