Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity

mortar tests, med file node ordering, ...

Browse Source
This commit is contained in:
Jukka Aho
2015-12-14 17:48:35 +02:00
parent 69bd7ce58c
commit a3f84b7292
5 changed files with 152 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/core.jl
+7 -8
View File
@@ -7,22 +7,21 @@ import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, c
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 `JULIAFEM_DEBUG=1 julia file.jl`.
Or set `export JULIAFEM_DEBUG=1` for your `.bashrc`.
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)
if !haskey(ENV, "JULIAFEM_DEBUG")
return
end
return :( println("DEBUG: ", $msg) )
haskey(ENV, "DEBUG") || return
# return :( println("DEBUG: ", $msg) )
return msg
end
function set_debug_on!()
ENV["JULIAFEM_DEBUG"] = 1;
ENV["DEBUG"] = 1;
end
function set_debug_off!()
pop!(ENV, "JULIAFEM_DEBUG");
pop!(ENV, "DEBUG");
end
export @debug, set_debug_on!, set_debug_off!
src/directsolver.jl
+2 -3
View File
@@ -143,7 +143,6 @@ 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))")
@assert solver.nonlinear_problem == true
timing = Dict{ASCIIString, Float64}()
tic(timing, "solver")
@@ -275,7 +274,7 @@ function call(solver::DirectSolver, time::Number=0.0)
toc(timing, "non-linear iteration")
if true
info("timing info for non-linear iteration:")
info("timing info for iteration:")
info("boundary assembly : ", time_elapsed(timing, "boundary assembly"))
info("field assembly : ", time_elapsed(timing, "field assembly"))
info("dump matrices to disk : ", time_elapsed(timing, "dump matrices to disk"))
@@ -284,7 +283,7 @@ function call(solver::DirectSolver, time::Number=0.0)
info("non-linear iteration : ", time_elapsed(timing, "non-linear iteration"))
end
if norm(sol) < solver.tol
if (norm(sol) < solver.tol) || !solver.nonlinear_problem
toc(timing, "solver")
info("solver finished in ", time_elapsed(timing, "solver"), " seconds.")
return (iter, true)
src/mortar.jl
+15 -10
View File
@@ -702,11 +702,12 @@ 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))
#S = reshape([S...;], 2, size(slave_element)[2])
S = hcat(Sl...)
integration_points = get_integration_points(Tri3, Val{5})
for master_element in slave_element["master elements"]
master_dofs = get_gdofs(master_element, field_dim)
# project master nodes to auxiliary plane and create polygon clipping
@@ -718,6 +719,9 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
M = hcat(M...)
P = nothing
neighbours = nothing
@debug info("applying polygon clip algorithm, S & M = ")
@debug dump(round(S, 3))
@debug dump(round(M, 3))
try
P, neighbours = clip_polygon(S, M)
catch
@@ -731,23 +735,23 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
error("cannot continue")
end
isa(P, Void) && continue # no clipping
# info("polygon on auxilyary plane: ")
# dump(round(P, 3))
@debug info("polygon on auxilyary plane: ")
@debug dump(round(P, 3))
C = calculate_polygon_centerpoint(P)
# info("center point = $C")
@debug info("center point = $C")
npts = size(P, 2) # number of vertices in polygon
# info("number of vectices in polygon: $npts")
@debug info("number of vectices in polygon: $npts")
# S = zeros(3, 3)
# M = zeros(3, 3)
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)'
# info("cell $i, coords = ")
# dump(round(X, 3))
@debug info("cell $i, coords = ")
@debug dump(round(X, 3))
geom = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)])
for ip in integration_points
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)
@@ -764,7 +768,8 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
N1 = slave_element(theta1[2:3], time)
N2 = master_element(theta2[2:3], time)
Sm = w*N1'*N1
Mm = w*N1'*N2
# FIXME: master side transpose -- why?
Mm = w*(N1'*N2)'
for k=1:field_dim
sd = slave_dofs[k:field_dim:end]
md = master_dofs[k:field_dim:end]
src/preprocess_aster_reader.jl
+18 -1
View File
@@ -100,6 +100,14 @@ function get_element_sets(med::MEDFile, mesh_name)
return es
end
# hex8 nodes rotating cw first in yz plane then x+1
global const med_elmap = Dict{Symbol, Vector{Int}}(
:HE8 => [4, 8, 7, 3, 1, 5, 6, 2],
:QU4 => [4, 3, 2, 1],
:SE2 => [2, 1]
)
function get_connectivity(med::MEDFile, elsets, mesh_name)
elsets[0] = :OTHER
increments = keys(med.data["ENS_MAA"][mesh_name])
@@ -116,7 +124,16 @@ function get_connectivity(med::MEDFile, elsets, mesh_name)
element_dim = round(Int, length(element_connectivity)/nelements)
element_connectivity = reshape(element_connectivity, nelements, element_dim)'
for i=1:nelements
d[element_ids[i]] = (Symbol(eltype), Symbol(elsets[elset_ids[i]]), element_connectivity[:, i])
eltype = Symbol(eltype)
elco = element_connectivity[:, i]
elset = Symbol(elsets[elset_ids[i]])
if haskey(med_elmap, eltype)
elco = elco[med_elmap[eltype]]
else
warn("no element mapping info found for element type $eltype")
warn("consider this as a warning: element may have french nodal ordering")
end
d[element_ids[i]] = (eltype, elset, elco)
end
end
return d
test/test_mortar.jl
+110 -18
View File
@@ -562,26 +562,17 @@ end
function test_assemble_3d_problem_quad4()
info("assemble 3d problem in quad4-quad4")
nodes = Vector{Float64}[
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[1.0, 1.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 0.1],
[1.0, 0.0, 0.1],
[1.0, 1.0, 0.1],
[0.0, 1.0, 0.1]]
#=
nodes = Vector{Float64}[
[-1.0, -1.0, 0.0],
[+1.0, -1.0, 0.0],
[+1.0, +1.0, 0.0],
[-1.0, +1.0, 0.0],
[-1.0, -1.0, 0.1],
[+1.0, -1.0, 0.1],
[+1.0, +1.0, 0.1],
[-1.0, +1.0, 0.1]]
=#
[2.0, 0.0, 0.1],
[2.0, 2.0, 0.1],
[0.0, 2.0, 0.1]]
mel = Quad4([5, 6, 7, 8])
mel["geometry"] = Vector{Float64}[nodes[5], nodes[6], nodes[7], nodes[8]]
sel = Quad4([1, 2, 3, 4])
@@ -591,14 +582,115 @@ function test_assemble_3d_problem_quad4()
prob = MortarProblem("temperature", 1)
push!(prob, sel)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
info("stiffness matrix for this problem:\n$stiffness_matrix")
M = D = 1/36*[4 2 1 2; 2 4 2 1; 1 2 4 2; 2 1 2 4]
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*144
D = [16 8 4 8; 8 16 8 4; 4 8 16 8; 8 4 8 16]
M = [25 5 1 5; 20 10 2 4; 16 8 4 8; 20 4 2 10]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:\n$B")
info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
end
#test_assemble_3d_problem_quad4()
function test_assemble_3d_problem_quad4_2()
info("assemble 3d problem in quad4-quad4")
nodes = Vector{Float64}[
[0.0, 0.0, 0.0],
[1/4, 0.0, 0.0],
[1/4, 1/4, 0.0],
[0.0, 1/4, 0.0],
[0.0, 0.0, 0.0],
[1/3, 0.0, 0.0],
[1/3, 1/3, 0.0],
[0.0, 1/3, 0.0]]
mel = Quad4([5, 6, 7, 8])
mel["geometry"] = Vector{Float64}[nodes[5], nodes[6], nodes[7], nodes[8]]
sel = Quad4([1, 2, 3, 4])
sel["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3], nodes[4]]
calculate_normal_tangential_coordinates!(sel, 0.0)
sel["master elements"] = Element[mel]
prob = MortarProblem("temperature", 1)
push!(prob, sel)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*589824
D = [
4096 2048 1024 2048
2048 4096 2048 1024
1024 2048 4096 2048
2048 1024 2048 4096
]
M = [
5184 1728 576 1728
3456 3456 1152 1152
2304 2304 2304 2304
3456 1152 1152 3456
]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
end
#test_assemble_3d_problem_quad4_2()
function test_assemble_3d_problem_quad4_3()
info("assemble 3d problem in quad4-quad4")
a = 1/4
b = 1/3
nodes = Vector{Float64}[
[2*a, a, 0],
[3*a, a, 0],
[3*a, 2*a, 0],
[2*a, 2*a, 0],
[ b, 0, 0],
[2*b, 0, 0],
[2*b, b, 0],
[ b, b, 0]]
mel = Quad4([5, 6, 7, 8])
mel["geometry"] = Vector{Float64}[nodes[5], nodes[6], nodes[7], nodes[8]]
sel = Quad4([1, 2, 3, 4])
sel["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3], nodes[4]]
calculate_normal_tangential_coordinates!(sel, 0.0)
sel["master elements"] = Element[mel]
prob = MortarProblem("temperature", 1)
push!(prob, sel)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*186624*9
D = [
7904 3040 560 1456
3040 2432 448 560
560 448 128 160
1456 560 160 416
]
M = [
504 1224 7956 3276
144 720 4680 936
18 90 990 198
63 153 1683 693
]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
end
test_assemble_3d_problem_quad4_3()
end