From b39d88ecca04a2146a1413021bf6a0e6c5961c48 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 14 Dec 2015 17:48:35 +0200 Subject: [PATCH] mortar tests, med file node ordering, ... --- src/core.jl | 15 ++-- src/directsolver.jl | 5 +- src/mortar.jl | 25 ++++--- src/preprocess_aster_reader.jl | 19 ++++- test/test_mortar.jl | 128 ++++++++++++++++++++++++++++----- 5 files changed, 152 insertions(+), 40 deletions(-) diff --git a/src/core.jl b/src/core.jl index ece2689..b88fbbd 100644 --- a/src/core.jl +++ b/src/core.jl @@ -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! diff --git a/src/directsolver.jl b/src/directsolver.jl index 6f463a7..5a961bd 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -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) diff --git a/src/mortar.jl b/src/mortar.jl index 1206bad..5834cf6 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -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] diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index c833dad..003a6b5 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -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 diff --git a/test/test_mortar.jl b/test/test_mortar.jl index a856083..e7b1434 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -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