diff --git a/src/integrate.jl b/src/integrate.jl index 47a4d43..14219f0 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -148,23 +148,56 @@ function get_integration_points(element::TriangularElement, ::Type{Val{4}}) return zip(weights, points) end +""" 7 point integration rule for triangular elements. + +References +---------- +Code Aster documentation, http://code-aster.org/doc/default/fr/man_r/r3/r3.01.01.pdf +""" function get_integration_points(element::TriangularElement, ::Type{Val{5}}) - weights = 0.5*[ - 0.22500000000000, - 0.13239415278851, - 0.13239415278851, - 0.13239415278851, - 0.12593918054483, - 0.12593918054483, - 0.12593918054483] + A = 0.470142064105115 + B = 0.101286507323456 + P1 = 0.066197076394253 + P2 = 0.062969590272413 + weights = [9/80, P1, P1, P1, P2, P2, P2] points = Vector{Float64}[ - [0.33333333333333, 0.33333333333333], - [0.47014206410511, 0.47014206410511], - [0.47014206410511, 0.05971587178977], - [0.05971587178977, 0.47014206410511], - [0.10128650732346, 0.10128650732346], - [0.10128650732346, 0.79742698535309], - [0.79742698535309, 0.10128650732346]] + [1/3, 1/3], + [A, A], + [1-2A, A], + [A, 1-2A], + [B, B], + [1-2B, B], + [B, 1-2B]] + return zip(weights, points) +end + +""" 12 point integration fule for triangular elements. +References +---------- +Code Aster documentation, http://code-aster.org/doc/default/fr/man_r/r3/r3.01.01.pdf +""" +function get_integration_points{E<:TriangularElement}(element::Element{E}, ::Type{Val{:FPG12}}) + A = 0.063089014491502 + B = 0.249286745170910 + C = 0.310352451033785 + D = 0.053145049844816 + P1 = 0.025422453185103 + P2 = 0.058393137863189 + P3 = 0.041425537809187 + weights = [P1, P1, P1, P2, P2, P2, P3, P3, P3, P3, P3, P3] + points = Vector{Float64}[ + [A, A], + [1-2A, A], + [A, 1-2A], + [B, B], + [1-2B, B], + [B, 1-2B], + [C, D], + [D, C], + [1-C-D, C], + [1-C,D, D], + [C, 1-C-D], + [D, 1-C-D]] return zip(weights, points) end diff --git a/src/preprocess.jl b/src/preprocess.jl index 58561a1..0a8ef78 100644 --- a/src/preprocess.jl +++ b/src/preprocess.jl @@ -17,13 +17,13 @@ import Base: copy using JuliaFEM type Mesh - nodes :: Dict{Int64, Vector{Float64}} - node_sets :: Dict{Symbol, Set{Int64}} - elements :: Dict{Int64, Vector{Int64}} - element_types :: Dict{Int64, Symbol} - element_codes :: Dict{Int64, Symbol} - element_sets :: Dict{Symbol, Set{Int64}} - surface_sets :: Dict{Symbol, Vector{Tuple{Int64, Symbol}}} + nodes :: Dict{Int, Vector{Float64}} + node_sets :: Dict{Symbol, Set{Int}} + elements :: Dict{Int, Vector{Int}} + element_types :: Dict{Int, Symbol} + element_codes :: Dict{Int, Symbol} + element_sets :: Dict{Symbol, Set{Int}} + surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}} surface_types :: Dict{Symbol, Symbol} end @@ -35,7 +35,7 @@ function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64}) mesh.nodes[nid] = ncoords end -function add_nodes!(mesh::Mesh, nodes::Dict{Int64, Vector{Float64}}) +function add_nodes!(mesh::Mesh, nodes::Dict{Int, Vector{Float64}}) for (nid, ncoords) in nodes add_node!(mesh, nid, ncoords) end @@ -43,17 +43,28 @@ end function add_node_to_node_set!(mesh::Mesh, set_name, nids...) if !haskey(mesh.node_sets, set_name) - mesh.node_sets[set_name] = Set{Int64}() + mesh.node_sets[set_name] = Set{Int}() end push!(mesh.node_sets[set_name], nids...) + return end -function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int64}) +""" Create a new node set from nodes in element set. """ +function create_node_set_from_element_set!(mesh::Mesh, set_name) + node_ids = Set{Int}() + for elid in mesh.element_sets[set_name] + push!(node_ids, mesh.elements[elid]...) + end + mesh.node_sets[set_name] = node_ids + return +end + +function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int}) mesh.elements[elid] = connectivity mesh.element_types[elid] = eltype end -function add_elements!(mesh::Mesh, elements::Dict{Int64, Tuple{Symbol, Vector{Int64}}}) +function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}}) for (elid, (eltype, elcon)) in elements add_element!(mesh, elid, eltype, elcon) end @@ -61,7 +72,7 @@ end function add_element_to_element_set!(mesh::Mesh, set_name, elids...) if !haskey(mesh.element_sets, set_name) - mesh.element_sets[set_name] = Set{Int64}() + mesh.element_sets[set_name] = Set{Int}() end push!(mesh.element_sets[set_name], elids...) end @@ -76,7 +87,7 @@ function copy(mesh::Mesh) return mesh2 end -function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int64}) +function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int}) mesh2 = copy(mesh) mesh2.elements = Dict() for elid in element_ids @@ -113,7 +124,7 @@ function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothi if isempty(element_sets) element_ids = collect(keys(mesh.elements)) else - element_ids = Set{Int64}() + element_ids = Set{Int}() for set_name in element_sets element_ids = union(element_ids, mesh.element_sets[set_name]) end @@ -135,7 +146,7 @@ end """ find npts nearest nodes from mesh and return id numbers as list. """ function find_nearest_nodes(mesh::Mesh, coords::Vector, npts=1) - dist = Dict{Int64, Float64}() + dist = Dict{Int, Float64}() for (nid, c) in mesh.nodes dist[nid] = norm(coords-c) end @@ -166,7 +177,7 @@ function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{ end end -function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int64) +function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int) problem = Problem(P, name, dimension) problem.elements = create_elements(mesh, name) return problem @@ -177,18 +188,3 @@ function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimen problem.elements = create_elements(mesh, name) return problem end - -""" -Swap surface element connectivity s.t. normals point outward -""" -function check_orientation! - # TODO -end - -""" -Partition model using METIS -""" -function partition_model! - # TODO -end - diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl index 8aad738..351c12a 100644 --- a/src/problems_mortar.jl +++ b/src/problems_mortar.jl @@ -42,12 +42,14 @@ type Mortar <: BoundaryProblem linear_surface_elements :: Bool split_quadratic_slave_elements :: Bool split_quadratic_master_elements :: Bool + alpha :: Float64 + drop_tolerance :: Float64 store_fields :: Vector{Symbol} end function Mortar() default_fields = [] - return Mortar(-1, false, false, false, false, Inf, true, true, true, default_fields) + return Mortar(-1, false, false, false, false, Inf, true, true, true, 0.0, 1.0e-9, default_fields) end function get_unknown_field_name(problem::Problem{Mortar}) @@ -63,13 +65,202 @@ function get_formulation_type(problem::Problem{Mortar}) end function assemble!(problem::Problem{Mortar}, time::Float64) + if length(problem.elements) == 0 + warn("No elements defined in interface $(problem.name), this will result empty assembly!") + return + end if problem.properties.dimension == -1 problem.properties.dimension = dim = size(first(problem.elements), 1) - info("assuming dimension of mesh tie surface is $dim") - info("if this is wrong set is manually using problem.properties.dimension") + info("Assuming dimension of mesh tie surface is $dim. If this is wrong set is manually using problem.properties.dimension") end dimension = Val{problem.properties.dimension} use_forwarddiff = Val{problem.properties.use_forwarddiff} assemble!(problem, time, dimension, use_forwarddiff) end +""" Given a CCW ordered set of vertices, calculate area of polygon. + +Examples +-------- + +julia> P = Vector[[1/3, 5/12, 1/2], [1/3, 1/2, 1/2], [1/2, 1/2, 1/2], [1/2, 1/3, 1/2], [5/12, 1/3, 1/2]] +5-element Array{Array{T,1},1}: + [0.333333,0.416667,0.5] + [0.333333,0.5,0.5] + [0.5,0.5,0.5] + [0.5,0.333333,0.5] + [0.416667,0.333333,0.5] + +julia> A = calculate_polygon_area(P) +0.02430555555555556 + +julia> isapprox(A, 7/288) +true + +""" +function calculate_polygon_area(P) + N_P = length(P) + A = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + return A +end + +""" Function to print useful debug information from interface to find bugs. """ +function diagnose_interface(problem::Problem{Mortar}, time::Float64) + info("Diagnosing Mortar interface...") + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_elements = get_slave_elements(problem) + + I_area = 0.0 + + if props.split_quadratic_slave_elements + info("props.split_quadratic_slave_elements = true") + if !props.linear_surface_elements + warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior") + end + slave_elements = split_quadratic_elements(slave_elements, time) + end + info("Number of slave elements in interface: $(length(slave_elements))") + + # 1. calculate nodal normals and tangents for slave element nodes j ∈ S + normals = calculate_normals(slave_elements, time, Val{2}; + rotate_normals=props.rotate_normals) + update!(slave_elements, "normal", time => normals) + + S_areas = [] + C_areas = [] + P_areas = [] + + for slave_element in slave_elements + + info(repeat("-", 80)) + info("Processing slave element $(slave_element.id), type = $(get_element_type(slave_element))") + info(repeat("-", 80)) + + S_area = 0.0 + S_area_in_contact = 0.0 + for ip in get_integration_points(slave_element) + S_area += ip.weight*slave_element(ip, time, Val{:detJ}) + end + info("Total area of slave element = $S_area") + + + if props.linear_surface_elements + info("Converting slave element to linear surface element") + slave_element = convert_to_linear_element(slave_element) + end + + slave_element_nodes = get_connectivity(slave_element) + info("Slave element connectivity = $slave_element_nodes") + nsl = length(slave_element) + X1 = slave_element("geometry", time) + n1 = Field([normals[j] for j in slave_element_nodes]) + + # project slave nodes to auxiliary plane (x0, Q) + xi = mean(get_reference_coordinates(slave_element)) + N = vec(get_basis(slave_element, xi, time)) + x0 = N*X1 + n0 = N*n1 + info("Auxiliary plane x0 = $x0, n0 = $n0") + S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] + check_orientation!(S, n0) + info("Slave element $(slave_element.id) vertices in auxiliary plane: $S") + + # 3. loop all master elements + master_elements = slave_element("master elements", time) + if props.split_quadratic_master_elements + master_elements = split_quadratic_elements(master_elements, time) + end + + for master_element in master_elements + + if props.linear_surface_elements + master_element = convert_to_linear_element(master_element) + end + + master_element_nodes = get_connectivity(master_element) + nm = length(master_element) + X2 = master_element("geometry", time) + + if norm(mean(X1) - mean(X2)) > problem.properties.distval + # elements are "far enough" + continue + end + + # 3.1 project master nodes to auxiliary plane and create polygon clipping + M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm] + check_orientation!(M, n0) + P = get_polygon_clip(S, M, n0) + if length(P) < 3 + if length(P) == 0 + continue + end + if length(P) == 1 + info("length(P) == 1, shared vertex") + end + if length(P) == 2 + info("length(P) == 2, shared edge") + end + continue + end + info("Master element $(master_element.id) vertices in auxiliary plane = $M") + check_orientation!(P, n0) + P_area_ = calculate_polygon_area(P) + info("Polygon clip found, P=$P, N_P = $(length(P)), area of polygon = $P_area_") + if isapprox(P_area_, 0.0) + error("Polygon P has zero area: $P_area_") + end + + P_area = 0.0 + + C0 = calculate_centroid(P) + info("Centroid of polygon = $C0") + + # 4. loop integration cells + all_cells = get_cells(P, C0) + info("Polygon is splitted to $(length(all_cells)) integration cells.") + for (cell_id, cell) in enumerate(all_cells) + C_area = 0.0 + virtual_element = Element(Tri3, Int[]) + update!(virtual_element, "geometry", cell) + + # 5. loop integration point of integration cell + for ip in get_integration_points(virtual_element, 3) + N = vec(get_basis(virtual_element, ip, time)) + detJ = virtual_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + # project gauss point from auxiliary plane to master and slave element + x_gauss = virtual_element("geometry", ip, time) + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) + xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) + C_area += w + end # integration points done + info("Cell $cell_id has area of $C_area") + P_area += C_area + push!(C_areas, C_area) + end # integration cells done + + if !isapprox(P_area, P_area_) + error("P_area = $P_area, should be $P_area_") + end + + S_area_in_contact += P_area + push!(P_areas, P_area) + + end # master elements done + + S_perc = S_area_in_contact / S_area * 100.0 + push!(S_areas, S_area_in_contact) + info("Area of slave element in contact: $S_area_in_contact, it's $S_perc % of total element area") + + I_area += S_area_in_contact + + end # slave elements done, contact virtual work ready + + info("Area of interface: $I_area") + info("Smallest cell area: $(minimum(C_areas))") + info("Smallest polygon area: $(minimum(P_areas))") + info("Smallest slave element area in contact: $(minimum(S_areas))") + +end diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl index 92a6b2a..d0affd7 100644 --- a/src/problems_mortar_3d.jl +++ b/src/problems_mortar_3d.jl @@ -112,7 +112,7 @@ function get_polygon_clip(xs, xm, n) # 1. test is master point inside slave, if yes, add to clip for i=1:nm if vertex_inside_polygon(xm[i], xs) - debug("1. $(xm[i]) inside S -> push") + # debug("1. $(xm[i]) inside S -> push") push!(P, xm[i]) end end @@ -121,7 +121,7 @@ function get_polygon_clip(xs, xm, n) for i=1:ns if vertex_inside_polygon(xs[i], xm) approx_in(xs[i], P) && continue - debug("2. $(xs[i]) inside M -> push") + # debug("2. $(xs[i]) inside M -> push") push!(P, xs[i]) end end @@ -144,7 +144,7 @@ function get_polygon_clip(xs, xm, n) #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))") if vertex_inside_polygon(q, xm) approx_in(q, P) && continue - debug("3. $q inside M -> push") + # debug("3. $q inside M -> push") push!(P, q) end end @@ -228,12 +228,39 @@ function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false) return normals end +""" Given polygon P and normal direction n, check that polygon vertices are +ordered in counter clock wise direction with respect to surface normal and +sort if necessary. It is assumed that polygon is convex. + +Examples +-------- +Unit triangle, normal in z-direction: + +julia> P = Vector[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]] +3-element Array{Array{T,1},1}: + [0.0,0.0,0.0] + [0.0,1.0,0.0] + [1.0,0.0,0.0] + +julia> n = [0.0, 0.0, 1.0] +3-element Array{Float64,1}: + 0.0 + 0.0 + 1.0 + +julia> check_orientation!(P, n) +3-element Array{Array{T,1},1}: + [1.0,0.0,0.0] + [0.0,0.0,0.0] + [0.0,1.0,0.0] + +""" function check_orientation!(P, n) C = mean(P) np = length(P) s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np] all(s .< 0) && return - debug("polygon not in ccw order, fixing") + # debug("polygon not in ccw order, fixing") # project points to new orthogonal basis Q and sort there t1 = (P[1]-C)/norm(P[1]-C) t2 = cross(n, t1) @@ -274,10 +301,9 @@ function split_quadratic_element(element::Element{Tri6}, time::Float64) u = element("displacement", time) update!(new_element, "displacement", time => u[elmap]) end - #n = element("normal", time) - #update!(new_element, "normal", time => n[elmap]) - if haskey(element, "master elements") - update!(new_element, "master elements", time => element("master elements", time)) + if haskey(element, "normal") + n = element("normal", time) + update!(new_element, "normal", time => n[elmap]) end push!(new_elements, new_element) end @@ -298,70 +324,62 @@ function split_quadratic_elements(elements::Vector, time::Float64) end n1 = length(elements) n2 = length(new_elements) - info("Splitted $n1 (maybe quadratic) elements to $n2 (linear) sub-elements") + if n1 != n2 + info("Splitted $n1 elements to $n2 (linear) sub-elements") + end return new_elements end -function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}) +""" Assemble linear surface element to problem. +Dual basis is constructed such that partially integrated slave segments are taken into account in a proper way. + +Notes +----- +For full integrated slave element, coefficient matrix for Tri3 is +Ae = [3.0 -1.0 -1.0; -1.0 3.0 -1.0; -1.0 -1.0 3.0] + +References +---------- + +[Popp2013] Popp, Alexander, et al. "Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach." Computer Methods in Applied Mechanics and Engineering 264 (2013): 67-80. +""" +function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) + props = problem.properties field_dim = get_unknown_field_dimension(problem) field_name = get_parent_field_name(problem) - slave_elements = get_slave_elements(problem) area = 0.0 - if props.split_quadratic_slave_elements - if !props.linear_surface_elements - warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior") - end - slave_elements = split_quadratic_elements(slave_elements, time) - end + slave_element_nodes = get_connectivity(slave_element) + nsl = length(slave_element) + X1 = slave_element("geometry", time) + n1 = slave_element("normal", time) - # 1. calculate nodal normals and tangents for slave element nodes j ∈ S - normals = calculate_normals(slave_elements, time, Val{2}; - rotate_normals=props.rotate_normals) - update!(slave_elements, "normal", time => normals) + # project slave nodes to auxiliary plane (x0, Q) + xi = mean(get_reference_coordinates(slave_element)) + first_slave_element && debug("midpoint xi = $xi") + N = vec(get_basis(slave_element, xi, time)) + x0 = N*X1 + n0 = N*n1 + S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] - # 2. loop all slave elements - first_slave_element = true + master_elements = slave_element("master elements", time) - for slave_element in slave_elements + if props.dual_basis - if props.linear_surface_elements - slave_element = convert_to_linear_element(slave_element) - end - - slave_element_nodes = get_connectivity(slave_element) - nsl = length(slave_element) - X1 = slave_element("geometry", time) - n1 = Field([normals[j] for j in slave_element_nodes]) - - # project slave nodes to auxiliary plane (x0, Q) - xi = mean(get_reference_coordinates(slave_element)) - first_slave_element && debug("midpoint xi = $xi") - N = vec(get_basis(slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 - S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] - - # 3. loop all master elements - master_elements = slave_element("master elements", time) - if props.split_quadratic_master_elements - master_elements = split_quadratic_elements(master_elements, time) - end + debug("Creating dual basis for element $(slave_element.id)") + De = zeros(nsl, nsl) + Me = zeros(nsl, nsl) + for master_element in master_elements - if props.linear_surface_elements - master_element = convert_to_linear_element(master_element) - end - master_element_nodes = get_connectivity(master_element) nm = length(master_element) X2 = master_element("geometry", time) if norm(mean(X1) - mean(X2)) > problem.properties.distval - # elements are "far enough" continue end @@ -373,170 +391,486 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ N_P = length(P) P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) - if first_slave_element - debug("Polygon clip info for first slave element:") - debug("S = $S") - debug("M = $M") - debug("P = $P") - debug("N_P = $N_P") - debug("P_area = $P_area") - end - if isapprox(P_area, 0.0) info("Polygon P has zero area: $P_area") continue end - C0 = calculate_centroid(P) - - #= - if isnan(C0[1]) - info("C0 = $C0") - info("P = $P") - info("S = $S") - info("M = $M") - info("n0 = $n0") - error("Calculation of centroid of polygon clip P failed.") - end - =# - - De = zeros(nsl, nsl) - Me = zeros(nsl, nm) - ge = zeros(field_dim*nsl) - # 4. loop integration cells + C0 = calculate_centroid(P) all_cells = get_cells(P, C0) for cell in all_cells virtual_element = Element(Tri3, Int[]) update!(virtual_element, "geometry", cell) - #x_cell = Field(cell) - - # construct bi-orthogonal basis - nnodes = length(slave_element) - if props.dual_basis - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(virtual_element, 3) - x_gauss = nothing - #try - x_gauss = virtual_element("geometry", ip, time) - xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) - detJ = virtual_element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N1 = vec(get_basis(slave_element, xi_s, time)) - De += w*diagm(vec(N1)) - Me += w*N1*N1' - #catch - # info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.") - # info("x_gauss = $x_gauss") - # info("cell = $cell") - # info("C0 = $C0") - # info("P = $P") - # info("S = $S") - # info("M = $M") - # info("n0 = $n0") - # rethrow() - #end - end - Ae = De*inv(Me) - else - Ae = eye(nnodes) - end - - # 5. loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) - N = vec(get_basis(virtual_element, ip, time)) - #dN = vec(get_dbasis(virtual_element, ip, time)) - #JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)])) - #wC = ip.weight*norm(cross(JC[:,1], JC[:,2])) detJ = virtual_element(ip, time, Val{:detJ}) w = ip.weight*detJ - - # project gauss point from auxiliary plane to master and slave element - #x_gauss = N*x_cell x_gauss = virtual_element("geometry", ip, time) - #= - if isnan(x_gauss[1]) - info("is nan") - info("x_gauss = $x_gauss") - info("cell = $cell") - info("C0 = $C0") - info("P = $P") - info("S = $S") - info("M = $M") - info("n0 = $n0") - error("nan, unable to continue") - end - =# - - xi_s = nothing - xi_m = nothing - alpha = nothing - - #try - xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) - xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) - #catch - # info("projecting vertex back to surface has failed.") - # info("x_gauss = $x_gauss") - # info("cell = $cell") - # info("C0 = $C0") - # info("P = $P") - # info("S = $S") - # info("M = $M") - # info("n0 = $n0") - # rethrow() - #end - - # add contributions - N1 = vec(get_basis(slave_element, xi_s, time)) - N2 = vec(get_basis(master_element, xi_m, time)) - Phi = Ae*N1 - De += w*Phi*N1' - Me += w*Phi*N2' - if props.adjust - u1 = slave_element("displacement", time) - u2 = master_element("displacement", time) - x_s = N1*(X1+u1) - x_m = N2*(X2+u2) - ge += w*vec((x_m-x_s)*Phi') - end - area += w - end # integration points done - + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) + N1 = slave_element(xi_s, time) + De += w*diagm(vec(N1)) + Me += w*N1'*N1 + end end # integration cells done - # 6. add contribution to contact virtual work - sdofs = get_gdofs(problem, slave_element) - mdofs = get_gdofs(problem, master_element) - - for i=1:field_dim - lsdofs = sdofs[i:field_dim:end] - lmdofs = mdofs[i:field_dim:end] - add!(problem.assembly.C1, lsdofs, lsdofs, De) - add!(problem.assembly.C1, lsdofs, lmdofs, -Me) - add!(problem.assembly.C2, lsdofs, lsdofs, De) - add!(problem.assembly.C2, lsdofs, lmdofs, -Me) + end # master elements done + + Ae = De*inv(Me) + + info("Dual basis coefficient matrix: $Ae") + + else + Ae = eye(nsl) + end + + for master_element in master_elements + + master_element_nodes = get_connectivity(master_element) + nm = length(master_element) + X2 = master_element("geometry", time) + + if norm(mean(X1) - mean(X2)) > problem.properties.distval + continue + end + + # 3.1 project master nodes to auxiliary plane and create polygon clipping + M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm] + P = get_polygon_clip(S, M, n0) + length(P) < 3 && continue # no clipping or shared edge (no volume) + check_orientation!(P, n0) + N_P = length(P) + P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + + if first_slave_element + debug("Polygon clip info for first slave element:") + debug("S = $S") + debug("M = $M") + debug("P = $P") + debug("N_P = $N_P") + debug("P_area = $P_area") + end + + if isapprox(P_area, 0.0) + info("Polygon P has zero area: $P_area") + continue + end + + C0 = calculate_centroid(P) + + De = zeros(nsl, nsl) + Me = zeros(nsl, nm) + ge = zeros(field_dim*nsl) + + # 4. loop integration cells + all_cells = get_cells(P, C0) + for cell in all_cells + virtual_element = Element(Tri3, Int[]) + update!(virtual_element, "geometry", cell) + + # 5. loop integration point of integration cell + for ip in get_integration_points(virtual_element, 3) + detJ = virtual_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + + # project gauss point from auxiliary plane to master and slave element + x_gauss = virtual_element("geometry", ip, time) + + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) + xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) + + # add contributions + N1 = vec(get_basis(slave_element, xi_s, time)) + N2 = vec(get_basis(master_element, xi_m, time)) + Phi = Ae*N1 + # Phi = [3.0-4.0*xi_s[1]-4.0*xi_s[2], 4.0*xi_s[1]-1.0, 4.0*xi_s[2]-1.0] + De += w*Phi*N1' + Me += w*Phi*N2' + if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement") + u1 = slave_element("displacement", time) + u2 = master_element("displacement", time) + x_s = N1*(X1+u1) + x_m = N2*(X2+u2) + ge += w*vec((x_m-x_s)*Phi') + end + area += w + end # integration points done + + end # integration cells done + + # 6. add contribution to contact virtual work + sdofs = get_gdofs(problem, slave_element) + mdofs = get_gdofs(problem, master_element) + + for i=1:field_dim + lsdofs = sdofs[i:field_dim:end] + lmdofs = mdofs[i:field_dim:end] + add!(problem.assembly.C1, lsdofs, lsdofs, De) + add!(problem.assembly.C1, lsdofs, lmdofs, -Me) + add!(problem.assembly.C2, lsdofs, lsdofs, De) + add!(problem.assembly.C2, lsdofs, lmdofs, -Me) + end + add!(problem.assembly.g, sdofs, ge) + + end # master elements done + + return area +end + + +""" Assemble quadratic surface element to problem. + +In polygon clipping element is divided to linear sub-elements proposed in [Puso2008]. + +References +---------- + +[Puso2008] Puso, Michael A., T. A. Laursen, and Jerome Solberg. "A segment-to-segment mortar contact method for quadratic elements and large deformations." Computer Methods in Applied Mechanics and Engineering 197.6 (2008): 555-566. + +[Popp1012] Popp, Alexander, et al. "Dual quadratic mortar finite element methods for 3D finite deformation contact." SIAM Journal on Scientific Computing 34.4 (2012): B421-B446. + +""" +function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) + + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + area = 0.0 + + Xs = slave_element("geometry", time) + + alp = props.alpha + + if alp != 0.0 + T = [ + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + alp alp 0.0 1.0-2*alp 0.0 0.0 + 0.0 alp alp 0.0 1.0-2*alp 0.0 + alp 0.0 alp 0.0 0.0 1.0-2*alp + ] + else + T = eye(6) + end + + #= + invT = [ + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0 + 0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 + -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp) + ] + =# + + if props.dual_basis + # info("Creating dual basis for element $(slave_element.id)") + nsl = length(slave_element) + De = zeros(nsl, nsl) + Me = zeros(nsl, nsl) + + # split slave element to linear sub-elements and loop + for sub_slave_element in split_quadratic_element(slave_element, time) + + slave_element_nodes = get_connectivity(sub_slave_element) + nsl = length(sub_slave_element) + X1 = sub_slave_element("geometry", time) + n1 = sub_slave_element("normal", time) + + # create auxiliary plane + xi = mean(get_reference_coordinates(sub_slave_element)) + first_slave_element && debug("midpoint xi = $xi") + N = vec(get_basis(sub_slave_element, xi, time)) + x0 = N*X1 + n0 = N*n1 + + # project slave nodes to auxiliary plane + S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] + + # 3. loop all master elements + master_elements = slave_element("master elements", time) + + for master_element in master_elements + + Xm = master_element("geometry", time) + + if norm(mean(Xs) - mean(Xm)) > problem.properties.distval + continue + end + + # split master element to linear sub-elements and loop + for sub_master_element in split_quadratic_element(master_element, time) + + master_element_nodes = get_connectivity(sub_master_element) + nm = length(sub_master_element) + X2 = sub_master_element("geometry", time) + + # 3.1 project master nodes to auxiliary plane + M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm] + + # create polygon clipping P + P = get_polygon_clip(S, M, n0) + length(P) < 3 && continue # no clipping or shared edge (no volume) + check_orientation!(P, n0) + N_P = length(P) + P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + + C0 = calculate_centroid(P) + + # 4. loop integration cells + all_cells = get_cells(P, C0) + for cell in all_cells + virtual_element = Element(Tri3, Int[]) + update!(virtual_element, "geometry", cell) + for ip in get_integration_points(virtual_element, 3) + x_gauss = virtual_element("geometry", ip, time) + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time) + detJ = virtual_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + N1 = vec(slave_element(xi_s, time)*T) + De += w*diagm(N1) + Me += w*N1*N1' + end + + end # integration cells done + + end # sub aster elements done + + end # master elements done + + end # sub slave elements done + + Ae = De*inv(Me) + # info("Dual basis construction finished.") + # info("Slave element geometry = $Xs") + # info("De = $De") + # info("Me = $Me") + # info("Dual basis coefficient matrix: $Ae") + + else + nsl = length(slave_element) + Ae = eye(nsl) + end + + # split slave element to linear sub-elements and loop + for sub_slave_element in split_quadratic_element(slave_element, time) + + slave_element_nodes = get_connectivity(sub_slave_element) + nsl = length(sub_slave_element) + X1 = sub_slave_element("geometry", time) + n1 = sub_slave_element("normal", time) + + # create auxiliary plane + xi = mean(get_reference_coordinates(sub_slave_element)) + first_slave_element && debug("midpoint xi = $xi") + N = vec(get_basis(sub_slave_element, xi, time)) + x0 = N*X1 + n0 = N*n1 + + # project slave nodes to auxiliary plane + S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] + + # 3. loop all master elements + master_elements = slave_element("master elements", time) + + for master_element in master_elements + + Xm = master_element("geometry", time) + + if norm(mean(Xs) - mean(Xm)) > problem.properties.distval + continue end - add!(problem.assembly.g, sdofs, ge) + + # split master element to linear sub-elements and loop + for sub_master_element in split_quadratic_element(master_element, time) + + master_element_nodes = get_connectivity(sub_master_element) + nm = length(sub_master_element) + X2 = sub_master_element("geometry", time) + + # 3.1 project master nodes to auxiliary plane + M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm] + + # create polygon clipping P + P = get_polygon_clip(S, M, n0) + length(P) < 3 && continue # no clipping or shared edge (no volume) + check_orientation!(P, n0) + N_P = length(P) + P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + + if first_slave_element + debug("Polygon clip info for first slave element:") + debug("S = $S") + debug("M = $M") + debug("P = $P") + debug("N_P = $N_P") + debug("P_area = $P_area") + end + + if isapprox(P_area, 0.0) + warn("Polygon P has zero area: $P_area") + continue + end + + C0 = calculate_centroid(P) + + # while our polygon clipping algorithm is working in linear sub elements + # contributions is calculated using quadratic shape functions + De = zeros(length(slave_element), length(slave_element)) + Me = zeros(length(slave_element), length(master_element)) + ge = zeros(field_dim*length(slave_element)) + + # 4. loop integration cells + all_cells = get_cells(P, C0) + for cell in all_cells + virtual_element = Element(Tri3, Int[]) + update!(virtual_element, "geometry", cell) + + # 5. loop integration point of integration cell + for ip in get_integration_points(virtual_element, 3) + + x_gauss = virtual_element("geometry", ip, time) + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time) + xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time) + + # add contributions + N1 = vec(slave_element(xi_s, time)*T) + N2 = vec(master_element(xi_m, time)) + Phi = Ae*N1 + + detJ = virtual_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + + De += w*Phi*N1' + Me += w*Phi*N2' + if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement") + u1 = slave_element("displacement", time) + u2 = master_element("displacement", time) + xs = N1*(Xs+u1) + xm = N2*(Xm+u2) + ge += w*vec((xm-xs)*Phi') + end + area += w + end # integration points done + + end # integration cells done + + # 6. add contribution to contact virtual work + sdofs = get_gdofs(problem, slave_element) + mdofs = get_gdofs(problem, master_element) + + for i=1:field_dim + lsdofs = sdofs[i:field_dim:end] + lmdofs = mdofs[i:field_dim:end] + add!(problem.assembly.C1, lsdofs, lsdofs, De) + add!(problem.assembly.C1, lsdofs, lmdofs, -Me) + add!(problem.assembly.C2, lsdofs, lsdofs, De) + add!(problem.assembly.C2, lsdofs, lmdofs, -Me) + end + add!(problem.assembly.g, sdofs, ge) + + end # sub aster elements done end # master elements done + + end # sub slave elements done + return area +end + + +function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}) + + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_elements = get_slave_elements(problem) + area = 0.0 + + #= + if props.split_quadratic_slave_elements + if !props.linear_surface_elements + warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior") + end + slave_elements = split_quadratic_elements(slave_elements, time) + end + =# + + # 1. calculate nodal normals and tangents for slave element nodes j ∈ S + normals = calculate_normals(slave_elements, time, Val{2}; + rotate_normals=props.rotate_normals) + + update!(slave_elements, "normal", time => normals) + + # 2. loop all slave elements + first_slave_element = true + + for slave_element in slave_elements + + area += assemble!(problem, slave_element, time; first_slave_element=first_slave_element) first_slave_element = false end # slave elements done, contact virtual work ready - - if problem.properties.dual_basis - tol = 1.0e-9 - debug("Dual basis is used, dropping small values for C1 & C2, tol = $tol") - C1 = sparse(problem.assembly.C1) - C2 = sparse(problem.assembly.C2) - SparseArrays.droptol!(C1, tol) - SparseArrays.droptol!(C2, tol) - problem.assembly.C1 = C1 - problem.assembly.C2 = C2 + + C1 = sparse(problem.assembly.C1) + C2 = sparse(problem.assembly.C2) + + maxdim = maximum(size(C1)) + if problem.properties.alpha != 0.0 + debug("mortar_3d: size C1 = ", size(C1), " max dim = $maxdim") + debug("alpha != 0.0, applying transformation D = Dh*T^-1") + alp = problem.properties.alpha + Te = [ + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + alp alp 0.0 1.0-2*alp 0.0 0.0 + 0.0 alp alp 0.0 1.0-2*alp 0.0 + alp 0.0 alp 0.0 0.0 1.0-2*alp + ] + invTe = [ + 1.0 0.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 + -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0 + 0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 + -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp) + ] + # construct global transformation matrices T and invT + T = SparseMatrixCOO() + invT = SparseMatrixCOO() + for element in slave_elements + dofs = get_gdofs(problem, element) + for i=1:field_dim + ldofs = dofs[i:field_dim:end] + add!(T, ldofs, ldofs, Te) + add!(invT, ldofs, ldofs, invTe) + end + end + T = sparse(T, maxdim, maxdim, (a, b) -> b) + invT = sparse(invT, maxdim, maxdim, (a, b) -> b) + # fill diagonal + d = ones(size(T, 1)) + d[get_nonzero_rows(T)] = 0.0 + T += spdiagm(d) + invT += spdiagm(d) + #invT2 = sparse(inv(full(T))) + #info("invT == invT2? ", invT == invT2) + #maxabsdiff = maximum(abs(invT - invT2)) + #info("max diff = $maxabsdiff") + C1 = C1*invT + C2 = C2*invT end + tol = problem.properties.drop_tolerance + debug("Dropping small values from C1 & C2, tolerace = $tol") + SparseArrays.droptol!(C1, tol) + SparseArrays.droptol!(C2, tol) + + problem.assembly.C1 = C1 + problem.assembly.C2 = C2 + debug("area of interface: $area") end diff --git a/src/solvers.jl b/src/solvers.jl index 918e504..5781eba 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -109,11 +109,13 @@ end function check_for_overconstrained_dofs(solver::Solver) overdetermined = false constrained_dofs = Set{Int}() + all_overconstrained_dofs = Set{Int}() boundary_problems = get_boundary_problems(solver) for problem in boundary_problems new_constraints = Set(problem.assembly.C2.I) new_constraints = setdiff(new_constraints, problem.assembly.removed_dofs) overconstrained_dofs = intersect(constrained_dofs, new_constraints) + all_overconstrained_dofs = union(all_overconstrained_dofs, overconstrained_dofs) if length(overconstrained_dofs) != 0 warn("problem is overconstrained, finding overconstrained dofs... ") overdetermined = true @@ -133,6 +135,8 @@ function check_for_overconstrained_dofs(solver::Solver) constrained_dofs = union(constrained_dofs, new_constraints) end if overdetermined + warn("List of all overconstrained dofs:") + warn(sort(collect(all_overconstrained_dofs))) error("problem is overconstrained, not continuing to solution.") end return true @@ -276,6 +280,7 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}}) u[:] = x[1:solver.ndofs] la[:] = x[solver.ndofs+1:end] + return true end diff --git a/test/test_mortar_3d_polygon_clip.jl b/test/test_mortar_3d_polygon_clip.jl index 508a06d..cb51b9c 100644 --- a/test/test_mortar_3d_polygon_clip.jl +++ b/test/test_mortar_3d_polygon_clip.jl @@ -1,10 +1,11 @@ # 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: get_polygon_clip, calculate_polygon_area using JuliaFEM.Testing -@testset "polygon clip case 1" begin +@testset "polygon clipping" begin + S = Vector[ [0.375, 0.0, 0.5], [0.6, 0.0, 0.5], @@ -15,17 +16,9 @@ using JuliaFEM.Testing [0.375, 0.25, 0.5]] n0 = [0.0, 0.0, 1.0] P = get_polygon_clip(S, M, n0) - P_expected = Vector{Float64}[ - [0.500, 0.0, 0.5], - [0.375, 0.0, 0.5], - [0.4375, 0.125, 0.5]] - @test length(P) == length(P_expected) - for (Pi, Pj) in zip(P, P_expected) - @test isapprox(Pi, Pj) - end -end - -@testset "polygon clip case 2" begin + @test length(P) == 3 + @test isapprox(calculate_polygon_area(P), 1/128) + S = Vector[ [0.25, 0.0, 0.5], [0.75, 0.0, 0.5], @@ -37,5 +30,13 @@ end n0 = [0.0, 0.0, 1.0] P = get_polygon_clip(S, M, n0) @test length(P) == 3 -end + @test isapprox(calculate_polygon_area(P), 1/48) + # visually inspected + Xs = Vector[[0.0, 0.0, 0.5], [1.0, 0.0, 0.5], [0.0, 1.0, 0.5]] + Xm = Vector[[-0.25, 0.50, 0.5], [0.50, -0.25, 0.5], [0.75,0.75, 0.5]] + P_ = Vector[[0.65,0.35,0.0], [0.5625,0.0,0.0], [0.25,0.0,0.0], + [0.0,0.25,0.0], [0.0,0.5625,0.0], [0.35,0.65,0.0]] + P = get_polygon_clip(Xs, Xm, [0.0, 0.0, 1.0]) + @test length(P) == length(P_) +end diff --git a/test/test_problems_mortar_3d.jl b/test/test_problems_mortar_3d.jl index 0e36c2d..db2325a 100644 --- a/test/test_problems_mortar_3d.jl +++ b/test/test_problems_mortar_3d.jl @@ -7,32 +7,40 @@ using JuliaFEM.Postprocess using JuliaFEM.Testing using JuliaFEM.Abaqus: create_surface_elements -@testset "test that interface transfers constant field without error" begin - meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med" - mesh = aster_read_mesh(meshfile) +### temperature patch tests, sl tet4, dl tet4, sl tet10, dl tet 10 - upper = Problem(Heat, "upper", 1) +tet4_meshfile = "test_problems_mortar_3d/tet4.inp" +tet10_meshfile = "test_problems_mortar_3d/tet10.inp" + +@testset "patch test temperature + abaqus inp + tet4" begin + mesh = abaqus_read_mesh(tet4_meshfile) + + upper = Problem(Heat, "UPPER", 1) upper.elements = create_elements(mesh, "UPPER") update!(upper, "temperature thermal conductivity", 1.0) - lower = Problem(Heat, "lower", 1) + + lower = Problem(Heat, "LOWER", 1) lower.elements = create_elements(mesh, "LOWER") update!(lower, "temperature thermal conductivity", 1.0) - bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature") - bc_upper.elements = create_elements(mesh, "UPPER_TOP") + bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") update!(bc_upper, "temperature 1", 0.0) - bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature") - bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") update!(bc_lower, "temperature 1", 1.0) interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature") - interface_slave_elements = create_elements(mesh, "LOWER_TOP") - interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") update!(interface_slave_elements, "master elements", interface_master_elements) interface.elements = [interface_master_elements; interface_slave_elements] + JuliaFEM.diagnose_interface(interface, 0.0) solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) + solver.xdmf = Xdmf("sl_lin_temp_results") + solver() node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) @@ -42,4 +50,413 @@ using JuliaFEM.Abaqus: create_surface_elements info("minT = $minT, maxT = $maxT") @test isapprox(minT, 0.5) @test isapprox(maxT, 0.5) + + #= + initialize!(solver) + assemble!(solver) + M, K, Kg, f, fg = get_field_assembly(solver) + Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) + K = K + Kg + Kb + f = f + fg + fb + K = 1/2*(K + K') + M = 1/2*(M + M') + =# + end + +@testset "patch test temperature + abaqus inp + tet4 + dual basis + adjust" begin + mesh = abaqus_read_mesh(tet4_meshfile) + + upper = Problem(Heat, "UPPER", 1) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "temperature thermal conductivity", 1.0) + + lower = Problem(Heat, "LOWER", 1) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "temperature thermal conductivity", 1.0) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "temperature 1", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "temperature 1", 1.0) + + interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_master_elements; interface_slave_elements] + interface.properties.dual_basis = true + #interface.properties.adjust = true + + JuliaFEM.diagnose_interface(interface, 0.0) + solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) + solver.xdmf = Xdmf("dl_lin_temp_results") + + solver() + + node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) + T = [t[1] for t in temperature] + minT = minimum(T) + maxT = maximum(T) + info("minT = $minT, maxT = $maxT") + @test isapprox(minT, 0.5) + @test isapprox(maxT, 0.5) + +end + +@testset "patch test temperature + abaqus inp + tet10, quadratic surface elements" begin + mesh = abaqus_read_mesh(tet10_meshfile) + + upper = Problem(Heat, "UPPER", 1) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "temperature thermal conductivity", 1.0) + + lower = Problem(Heat, "LOWER", 1) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "temperature thermal conductivity", 1.0) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "temperature 1", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "temperature 1", 1.0) + + interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_master_elements; interface_slave_elements] + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.alpha = 0.0 +# JuliaFEM.diagnose_interface(interface, 0.0) + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) + solver.xdmf = Xdmf("sl_quad_temp_results") + solver() + + node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) + T = [t[1] for t in temperature] + + minT = minimum(T) + maxT = maximum(T) + stdT = std(T) + info("minT = $minT, maxT = $maxT, stdT = $stdT") + @test maxT - minT < 1.0e-10 + @test isapprox(stdT, 0.0; atol=1.0e-10) + +end + +@testset "patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2" begin + mesh = abaqus_read_mesh(tet10_meshfile) + + upper = Problem(Heat, "UPPER", 1) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "temperature thermal conductivity", 1.0) + + lower = Problem(Heat, "LOWER", 1) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "temperature thermal conductivity", 1.0) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "temperature 1", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "temperature 1", 1.0) + + interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_master_elements; interface_slave_elements] + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.dual_basis = true + interface.properties.alpha = 0.2 + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) + solver.xdmf = Xdmf("dl_quad_temp_results") + solver() + + node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) + #node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0) + #= + for (j, (nid, T)) in enumerate(zip(node_ids, temperature)) + info("$j: $nid -> $(T[1])") + j == 10 && break + end + =# + T = [t[1] for t in temperature] + minT = minimum(T) + maxT = maximum(T) + stdT = std(T) + info("minT = $minT, maxT = $maxT, stdT = $stdT") + @test maxT - minT < 1.0e-10 + @test isapprox(stdT, 0.0; atol=1.0e-10) +end + +### displacement patch tests, sl tet4, dl tet4, sl tet10, dl tet 10 + +@testset "patch test displacement + abaqus inp + tet4 + adjust" begin + + mesh = abaqus_read_mesh(tet4_meshfile) + # modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust + JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER) + for nid in mesh.node_sets[:UPPER] + mesh.nodes[nid][3] += 0.2 + end + + upper = Problem(Elasticity, "UPPER", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 288.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "LOWER", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 288.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 3", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 3", 0.0) + + bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement") + bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")] + update!(bc_sym13, "displacement 2", 0.0) + + bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement") + bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")] + update!(bc_sym23, "displacement 1", 0.0) + + interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_slave_elements; interface_master_elements] + + append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627]) + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.adjust = true + interface.properties.dual_basis = false + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) + solver.xdmf = Xdmf("sl_lin_disp_results") + solver() + + node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) + node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) + u3 = [u[3] for u in displacement] + maxabsu3 = maximum(abs(u3)) + stdabsu3 = std(abs(u3)) + info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + @test isapprox(stdabsu3, 0.0; atol=1.0e-10) +end + +@testset "patch test displacement + abaqus inp + tet4 + adjust + dual basis" begin + + mesh = abaqus_read_mesh(tet4_meshfile) + # modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust + JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER) + for nid in mesh.node_sets[:UPPER] + mesh.nodes[nid][3] += 0.2 + end + + upper = Problem(Elasticity, "UPPER", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 288.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "LOWER", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 288.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 3", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 3", 0.0) + + bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement") + bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")] + update!(bc_sym13, "displacement 2", 0.0) + + bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement") + bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")] + update!(bc_sym23, "displacement 1", 0.0) + + interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_slave_elements; interface_master_elements] + + append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627]) + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.adjust = true + interface.properties.dual_basis = true + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) + solver.xdmf = Xdmf("dl_lin_disp_results") + solver() + + node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) + node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) + u3 = [u[3] for u in displacement] + maxabsu3 = maximum(abs(u3)) + stdabsu3 = std(abs(u3)) + info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + @test isapprox(stdabsu3, 0.0; atol=1.0e-10) +end + +@testset "patch test displacement + abaqus inp + tet10 + adjust" begin + + mesh = abaqus_read_mesh(tet10_meshfile) + # modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust + JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER) + for nid in mesh.node_sets[:UPPER] + mesh.nodes[nid][3] += 0.2 + end + + upper = Problem(Elasticity, "UPPER", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 288.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "LOWER", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 288.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 3", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 3", 0.0) + + bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement") + bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")] + update!(bc_sym13, "displacement 2", 0.0) + + bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement") + bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")] + update!(bc_sym23, "displacement 1", 0.0) + + interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_slave_elements; interface_master_elements] + + removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657] + append!(interface.assembly.removed_dofs, removed_dofs) + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.adjust = true + interface.properties.dual_basis = false + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) + solver.xdmf = Xdmf("sl_quad_disp_results") + solver() + + node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) + node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) + u3 = [u[3] for u in displacement] + maxabsu3 = maximum(abs(u3)) + stdabsu3 = std(abs(u3)) + info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + @test isapprox(stdabsu3, 0.0; atol=1.0e-6) +end + + +@testset "patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2" begin + + mesh = abaqus_read_mesh(tet10_meshfile) + # modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust + JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER) + for nid in mesh.node_sets[:UPPER] + mesh.nodes[nid][3] += 0.2 + end + + upper = Problem(Elasticity, "UPPER", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 288.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "LOWER", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 288.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 3", 0.0) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 3", 0.0) + + bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement") + bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")] + update!(bc_sym13, "displacement 2", 0.0) + + bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement") + bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")] + update!(bc_sym23, "displacement 1", 0.0) + + interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_slave_elements; interface_master_elements] + + removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657] + append!(interface.assembly.removed_dofs, removed_dofs) + + interface.properties.linear_surface_elements = false + interface.properties.split_quadratic_slave_elements = false + interface.properties.split_quadratic_master_elements = false + interface.properties.adjust = true + interface.properties.dual_basis = true + interface.properties.alpha = 0.2 + + solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) + solver.xdmf = Xdmf("dl_quad_disp_results") + solver() + + node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) + node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) + u3 = [u[3] for u in displacement] + maxabsu3 = maximum(abs(u3)) + stdabsu3 = std(abs(u3)) + info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + @test isapprox(stdabsu3, 0.0; atol=1.0e-6) +end + diff --git a/test/testdata/test_problems_mortar_3d_tet10.inp b/test/test_problems_mortar_3d/tet10.inp old mode 100755 new mode 100644 similarity index 100% rename from test/testdata/test_problems_mortar_3d_tet10.inp rename to test/test_problems_mortar_3d/tet10.inp diff --git a/test/test_problems_mortar_3d/tet4.inp b/test/test_problems_mortar_3d/tet4.inp new file mode 100644 index 0000000..fabfbb8 --- /dev/null +++ b/test/test_problems_mortar_3d/tet4.inp @@ -0,0 +1,589 @@ +**NSET COUNT = 116 +*NODE +127, 1.00000, 0.00000, 0.75000 +128, 1.00000, 0.00000, 1.00000 +129, 0.75000, 0.00000, 1.00000 +133, 0.75000, 0.00000, 0.75000 +136, 1.00000, 0.00000, 0.50000 +139, 0.75000, 0.00000, 0.50000 +142, 0.50000, 0.00000, 1.00000 +145, 0.50000, 0.00000, 0.75000 +149, 0.50000, 0.00000, 0.50000 +152, 0.25000, 0.00000, 1.00000 +155, 0.25000, 0.00000, 0.75000 +159, 0.25000, 0.00000, 0.50000 +162, 0.00000, 0.00000, 1.00000 +165, 0.00000, 0.00000, 0.75000 +169, 0.00000, 0.00000, 0.50000 +172, 0.75000, 1.00000, 1.00000 +173, 1.00000, 1.00000, 1.00000 +174, 1.00000, 1.00000, 0.75000 +178, 0.75000, 1.00000, 0.75000 +181, 1.00000, 1.00000, 0.50000 +184, 0.75000, 1.00000, 0.50000 +187, 0.50000, 1.00000, 1.00000 +190, 0.50000, 1.00000, 0.75000 +194, 0.50000, 1.00000, 0.50000 +197, 0.25000, 1.00000, 1.00000 +200, 0.25000, 1.00000, 0.75000 +204, 0.25000, 1.00000, 0.50000 +207, 0.00000, 1.00000, 1.00000 +210, 0.00000, 1.00000, 0.75000 +214, 0.00000, 1.00000, 0.50000 +217, 1.00000, 0.75000, 1.00000 +220, 1.00000, 0.75000, 0.75000 +224, 1.00000, 0.75000, 0.50000 +227, 1.00000, 0.50000, 1.00000 +230, 1.00000, 0.50000, 0.75000 +234, 1.00000, 0.50000, 0.50000 +237, 1.00000, 0.25000, 1.00000 +240, 1.00000, 0.25000, 0.75000 +244, 1.00000, 0.25000, 0.50000 +252, 0.00000, 0.75000, 1.00000 +255, 0.00000, 0.75000, 0.75000 +259, 0.00000, 0.75000, 0.50000 +262, 0.00000, 0.50000, 1.00000 +265, 0.00000, 0.50000, 0.75000 +269, 0.00000, 0.50000, 0.50000 +272, 0.00000, 0.25000, 1.00000 +275, 0.00000, 0.25000, 0.75000 +279, 0.00000, 0.25000, 0.50000 +288, 0.75000, 0.75000, 1.00000 +292, 0.50000, 0.75000, 1.00000 +296, 0.25000, 0.75000, 1.00000 +302, 0.75000, 0.50000, 1.00000 +306, 0.50000, 0.50000, 1.00000 +310, 0.25000, 0.50000, 1.00000 +316, 0.75000, 0.25000, 1.00000 +320, 0.50000, 0.25000, 1.00000 +324, 0.25000, 0.25000, 1.00000 +337, 0.75000, 0.75000, 0.50000 +341, 0.50000, 0.75000, 0.50000 +345, 0.25000, 0.75000, 0.50000 +351, 0.75000, 0.50000, 0.50000 +355, 0.50000, 0.50000, 0.50000 +359, 0.25000, 0.50000, 0.50000 +365, 0.75000, 0.25000, 0.50000 +369, 0.50000, 0.25000, 0.50000 +373, 0.25000, 0.25000, 0.50000 +438, 1.00000, 0.00000, 0.25000 +439, 1.00000, 0.00000, 0.50000 +440, 0.66667, 0.00000, 0.50000 +444, 0.66667, 0.00000, 0.25000 +447, 1.00000, 0.00000, 0.00000 +450, 0.66667, 0.00000, 0.00000 +453, 0.33333, 0.00000, 0.50000 +456, 0.33333, 0.00000, 0.25000 +460, 0.33333, 0.00000, 0.00000 +463, 0.00000, 0.00000, 0.50000 +466, 0.00000, 0.00000, 0.25000 +470, 0.00000, 0.00000, 0.00000 +473, 0.66667, 1.00000, 0.50000 +474, 1.00000, 1.00000, 0.50000 +475, 1.00000, 1.00000, 0.25000 +479, 0.66667, 1.00000, 0.25000 +482, 1.00000, 1.00000, 0.00000 +485, 0.66667, 1.00000, 0.00000 +488, 0.33333, 1.00000, 0.50000 +491, 0.33333, 1.00000, 0.25000 +495, 0.33333, 1.00000, 0.00000 +498, 0.00000, 1.00000, 0.50000 +501, 0.00000, 1.00000, 0.25000 +505, 0.00000, 1.00000, 0.00000 +508, 1.00000, 0.66667, 0.50000 +511, 1.00000, 0.66667, 0.25000 +515, 1.00000, 0.66667, 0.00000 +518, 1.00000, 0.33333, 0.50000 +521, 1.00000, 0.33333, 0.25000 +525, 1.00000, 0.33333, 0.00000 +533, 0.00000, 0.66667, 0.50000 +536, 0.00000, 0.66667, 0.25000 +540, 0.00000, 0.66667, 0.00000 +543, 0.00000, 0.33333, 0.50000 +546, 0.00000, 0.33333, 0.25000 +550, 0.00000, 0.33333, 0.00000 +559, 0.66667, 0.66667, 0.50000 +563, 0.33333, 0.66667, 0.50000 +569, 0.66667, 0.33333, 0.50000 +573, 0.33333, 0.33333, 0.50000 +584, 0.66667, 0.66667, 0.00000 +588, 0.33333, 0.66667, 0.00000 +594, 0.66667, 0.33333, 0.00000 +598, 0.33333, 0.33333, 0.00000 +608, 0.66006, 0.47128, 0.70078 +609, 0.19762, 0.64074, 0.81889 +610, 0.65771, 0.82813, 0.74829 +611, 0.59993, 0.15953, 0.76106 +612, 0.40915, 0.16311, 0.74937 +613, 0.50000, 0.50000, 0.25000 +** +**ELSET COUNT = 172 +**HWCOLOR COMP 54 0 +*ELEMENT, TYPE=C3D4, ELSET=UPPER + 570, 252, 262, 255, 609 + 571, 252, 310, 262, 609 + 572, 252, 296, 310, 609 + 573, 200, 252, 255, 609 + 574, 296, 306, 310, 609 + 575, 262, 310, 265, 609 + 576, 259, 265, 269, 359 + 577, 184, 337, 220, 610 + 578, 259, 345, 265, 359 + 579, 288, 608, 292, 610 + 580, 187, 288, 292, 610 + 581, 172, 288, 187, 610 + 582, 178, 172, 187, 610 + 583, 178, 220, 172, 610 + 584, 172, 220, 288, 610 + 585, 230, 240, 237, 316 + 586, 178, 184, 220, 610 + 587, 178, 190, 184, 610 + 588, 178, 187, 190, 610 + 589, 187, 292, 190, 610 + 590, 190, 292, 609, 610 + 591, 190, 609, 341, 610 + 592, 190, 341, 194, 610 + 593, 184, 190, 194, 610 + 594, 184, 194, 337, 610 + 595, 194, 341, 337, 610 + 596, 337, 341, 608, 610 + 597, 230, 337, 608, 610 + 598, 230, 608, 288, 610 + 599, 220, 230, 288, 610 + 600, 220, 337, 230, 610 + 601, 341, 609, 608, 610 + 602, 292, 608, 609, 610 + 603, 306, 316, 320, 611 + 604, 306, 608, 316, 611 + 605, 365, 608, 369, 611 + 606, 240, 608, 365, 611 + 607, 240, 316, 608, 611 + 608, 240, 133, 316, 611 + 609, 240, 365, 133, 611 + 610, 365, 139, 133, 611 + 611, 365, 369, 139, 611 + 612, 369, 149, 139, 611 + 613, 149, 145, 611, 612 + 614, 240, 127, 128, 316 + 615, 320, 129, 142, 611 + 616, 316, 129, 320, 611 + 617, 316, 133, 129, 611 + 618, 129, 133, 142, 611 + 619, 142, 133, 145, 611 + 620, 139, 145, 133, 611 + 621, 139, 149, 145, 611 + 622, 310, 320, 324, 612 + 623, 306, 320, 310, 612 + 624, 306, 611, 320, 612 + 625, 306, 608, 611, 612 + 626, 306, 609, 608, 612 + 627, 306, 310, 609, 612 + 628, 265, 609, 310, 612 + 629, 265, 310, 324, 612 + 630, 265, 324, 275, 612 + 631, 265, 275, 373, 612 + 632, 265, 373, 359, 612 + 633, 265, 359, 609, 612 + 634, 355, 609, 359, 612 + 635, 355, 608, 609, 612 + 636, 355, 369, 608, 612 + 637, 369, 611, 608, 612 + 638, 214, 255, 259, 345 + 639, 142, 611, 145, 612 + 640, 320, 611, 142, 612 + 641, 320, 142, 324, 612 + 642, 324, 142, 152, 612 + 643, 324, 152, 155, 612 + 644, 275, 324, 155, 612 + 645, 275, 155, 373, 612 + 646, 373, 155, 159, 612 + 647, 373, 159, 149, 612 + 648, 369, 373, 149, 612 + 649, 359, 373, 369, 612 + 650, 355, 359, 369, 612 + 651, 369, 149, 611, 612 + 652, 149, 155, 145, 612 + 653, 149, 159, 155, 612 + 654, 152, 145, 155, 612 + 655, 142, 145, 152, 612 + 656, 237, 240, 128, 316 + 657, 244, 136, 127, 365 + 658, 240, 365, 127, 133 + 659, 240, 351, 244, 365 + 660, 220, 224, 230, 337 + 661, 227, 288, 230, 302 + 662, 200, 207, 252, 296 + 663, 240, 127, 316, 133 + 664, 230, 337, 234, 351 + 665, 275, 162, 165, 155 + 666, 172, 220, 217, 288 + 667, 200, 252, 210, 255 + 668, 190, 292, 197, 296 + 669, 230, 234, 240, 351 + 670, 204, 210, 214, 345 + 671, 210, 255, 214, 345 + 672, 275, 165, 373, 155 + 673, 172, 174, 217, 220 + 674, 200, 204, 341, 345 + 675, 240, 244, 127, 365 + 676, 220, 230, 227, 288 + 677, 272, 275, 324, 155 + 678, 272, 152, 162, 155 + 679, 200, 255, 210, 345 + 680, 190, 200, 194, 341 + 681, 197, 207, 200, 296 + 682, 190, 197, 200, 296 + 683, 279, 165, 169, 159 + 684, 200, 207, 210, 252 + 685, 224, 234, 230, 337 + 686, 187, 197, 190, 292 + 687, 275, 165, 279, 373 + 688, 265, 275, 269, 373 + 689, 272, 324, 152, 155 + 690, 262, 265, 310, 324 + 691, 265, 272, 275, 324 + 692, 184, 224, 220, 337 + 693, 272, 162, 275, 155 + 694, 262, 272, 265, 324 + 695, 178, 174, 172, 220 + 696, 136, 365, 139, 133 + 697, 227, 230, 237, 302 + 698, 178, 184, 181, 220 + 699, 172, 174, 173, 217 + 700, 194, 200, 204, 341 + 701, 234, 244, 240, 351 + 702, 200, 210, 204, 345 + 703, 230, 237, 302, 316 + 704, 128, 316, 127, 129 + 705, 217, 220, 227, 288 + 706, 127, 365, 136, 133 + 707, 265, 269, 359, 373 + 708, 279, 373, 165, 159 + 709, 165, 159, 373, 155 + 710, 269, 275, 279, 373 + 711, 127, 129, 316, 133 + 712, 181, 220, 184, 224 + 713, 178, 181, 174, 220 + 714, 292, 302, 306, 608 + 715, 288, 302, 292, 608 + 716, 230, 302, 288, 608 + 717, 230, 316, 302, 608 + 718, 230, 240, 316, 608 + 719, 230, 351, 240, 608 + 720, 230, 337, 351, 608 + 721, 302, 316, 306, 608 + 722, 240, 351, 365, 608 + 723, 351, 355, 365, 608 + 724, 341, 355, 351, 608 + 725, 337, 341, 351, 608 + 726, 355, 369, 365, 608 + 727, 341, 355, 608, 609 + 728, 292, 608, 306, 609 + 729, 292, 306, 296, 609 + 730, 190, 292, 296, 609 + 731, 190, 200, 341, 609 + 732, 190, 296, 200, 609 + 733, 200, 296, 252, 609 + 734, 200, 255, 345, 609 + 735, 200, 345, 341, 609 + 736, 341, 345, 355, 609 + 737, 345, 359, 355, 609 + 738, 265, 359, 345, 609 + 739, 255, 265, 345, 609 + 740, 255, 265, 259, 345 + 741, 255, 262, 265, 609 +** +**ELSET COUNT = 92 +**HWCOLOR COMP 1 0 +*ELEMENT, TYPE=C3D4, ELSET=LOWER + 742, 456, 460, 466, 598 + 743, 444, 450, 456, 598 + 744, 444, 525, 447, 594 + 745, 536, 540, 588, 598 + 746, 453, 456, 463, 543 + 747, 533, 543, 536, 573 + 748, 444, 456, 453, 573 + 749, 466, 550, 546, 598 + 750, 536, 543, 546, 573 + 751, 491, 536, 495, 588 + 752, 438, 447, 444, 525 + 753, 444, 594, 450, 598 + 754, 444, 521, 525, 594 + 755, 438, 444, 440, 521 + 756, 540, 546, 550, 598 + 757, 453, 456, 543, 573 + 758, 488, 491, 559, 563 + 759, 505, 536, 540, 588 + 760, 440, 444, 453, 573 + 761, 491, 495, 584, 588 + 762, 444, 447, 450, 594 + 763, 485, 491, 495, 584 + 764, 440, 569, 444, 573 + 765, 495, 501, 505, 536 + 766, 473, 479, 511, 559 + 767, 438, 440, 439, 521 + 768, 491, 533, 536, 563 + 769, 444, 569, 521, 594 + 770, 473, 479, 475, 511 + 771, 491, 498, 501, 533 + 772, 495, 536, 505, 588 + 773, 438, 444, 521, 525 + 774, 473, 511, 508, 559 + 775, 440, 521, 444, 569 + 776, 511, 515, 521, 584 + 777, 491, 498, 533, 563 + 778, 515, 525, 521, 584 + 779, 460, 550, 466, 598 + 780, 439, 521, 440, 569 + 781, 456, 463, 543, 546 + 782, 479, 482, 511, 515 + 783, 533, 536, 563, 573 + 784, 479, 515, 511, 584 + 785, 450, 460, 456, 598 + 786, 521, 569, 559, 594 + 787, 559, 594, 569, 613 + 788, 559, 584, 594, 613 + 789, 521, 559, 584, 594 + 790, 536, 546, 540, 598 + 791, 511, 521, 559, 584 + 792, 479, 559, 491, 584 + 793, 491, 559, 563, 584 + 794, 511, 521, 518, 559 + 795, 491, 584, 563, 588 + 796, 563, 588, 584, 613 + 797, 559, 563, 584, 613 + 798, 559, 569, 563, 613 + 799, 563, 569, 573, 613 + 800, 536, 563, 573, 588 + 801, 563, 573, 588, 613 + 802, 573, 598, 588, 613 + 803, 491, 501, 495, 536 + 804, 536, 573, 546, 598 + 805, 536, 588, 573, 598 + 806, 588, 598, 594, 613 + 807, 584, 588, 594, 613 + 808, 444, 569, 594, 598 + 809, 521, 584, 525, 594 + 810, 444, 573, 569, 598 + 811, 569, 598, 573, 613 + 812, 569, 594, 598, 613 + 813, 444, 456, 573, 598 + 814, 439, 518, 521, 569 + 815, 456, 546, 543, 573 + 816, 479, 485, 482, 515 + 817, 508, 511, 518, 559 + 818, 456, 546, 573, 598 + 819, 488, 498, 491, 563 + 820, 479, 485, 515, 584 + 821, 491, 563, 536, 588 + 822, 473, 488, 479, 559 + 823, 479, 491, 485, 584 + 824, 491, 533, 501, 536 + 825, 479, 511, 559, 584 + 826, 479, 488, 491, 559 + 827, 475, 479, 482, 511 + 828, 518, 559, 521, 569 + 829, 473, 475, 508, 511 + 830, 473, 475, 474, 508 + 831, 460, 470, 466, 550 + 832, 456, 466, 463, 546 + 833, 456, 466, 546, 598 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER +799,S1 +798,S1 +764,S2 +814,S2 +747,S2 +828,S2 +819,S2 +822,S2 +777,S3 +746,S4 +757,S4 +760,S4 +780,S4 +783,S4 +817,S4 +758,S4 +774,S4 +830,S4 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM +806,S1 +807,S1 +831,S2 +779,S2 +785,S2 +778,S2 +753,S3 +762,S3 +744,S3 +745,S3 +809,S3 +761,S3 +820,S3 +816,S3 +756,S4 +759,S4 +772,S4 +763,S4 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13 +831,S1 +742,S1 +832,S1 +746,S1 +785,S1 +743,S1 +748,S1 +760,S1 +762,S1 +752,S1 +755,S1 +767,S1 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23 +749,S1 +756,S1 +790,S1 +750,S1 +747,S1 +759,S1 +831,S3 +832,S3 +781,S3 +765,S3 +824,S3 +771,S3 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP +678,S1 +689,S1 +642,S1 +641,S1 +615,S1 +616,S1 +622,S1 +623,S1 +603,S1 +721,S1 +571,S1 +572,S1 +574,S1 +729,S1 +714,S1 +715,S1 +580,S1 +581,S1 +704,S2 +694,S2 +661,S2 +681,S2 +686,S2 +703,S3 +662,S3 +668,S3 +656,S4 +690,S4 +697,S4 +705,S4 +666,S4 +699,S4 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER +647,S1 +648,S1 +612,S1 +611,S1 +696,S1 +649,S1 +650,S1 +726,S1 +723,S1 +737,S1 +736,S1 +724,S1 +725,S1 +595,S1 +594,S1 +708,S2 +657,S2 +701,S2 +578,S2 +685,S2 +692,S2 +707,S3 +659,S3 +664,S3 +674,S3 +683,S4 +710,S4 +576,S4 +638,S4 +670,S4 +700,S4 +712,S4 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13 +653,S1 +652,S1 +654,S1 +655,S1 +621,S1 +620,S1 +619,S1 +618,S1 +709,S2 +711,S2 +683,S3 +665,S3 +678,S3 +696,S4 +706,S4 +704,S4 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23 +683,S1 +687,S1 +665,S1 +693,S1 +710,S1 +688,S1 +691,S1 +694,S1 +576,S1 +740,S1 +741,S1 +570,S1 +638,S1 +671,S1 +667,S3 +684,S3 +** +**Property Definitions +** +*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material +*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material +** +**Material Definitions +** +**Material:Def_Material +*MATERIAL,NAME=Def_Material +*ELASTIC,TYPE=ISO +2.08000e+005,3.00000e-001 +*DENSITY +7.80000e-009, +*SPECIFIC HEAT +5.00000e-001 +*CONDUCTIVITY +4.98100e-002 +** \ No newline at end of file diff --git a/test/test_problems_mortar_3d_lowlevel.jl b/test/test_problems_mortar_3d_lowlevel.jl new file mode 100644 index 0000000..81c535a --- /dev/null +++ b/test/test_problems_mortar_3d_lowlevel.jl @@ -0,0 +1,196 @@ +# 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.Postprocess +using JuliaFEM.Testing +using JuliaFEM.Abaqus: create_surface_elements + +@testset "forget to add elements to problem" begin + X = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [0.0, 1.0, 0.0], + 4 => [-0.25, 0.50, 0.00], + 5 => [0.50, -0.25, 0.00], + 6 => [0.75, 0.75, 0.00]) + s = Element(Tri3, [1, 2, 3]) + m = Element(Tri3, [4, 5, 6]) + update!([s, m], "geometry", X) + update!(s, "master elements", [m]) + p = Problem(Mortar, "two elements", 1, "temperature") + initialize!(p) + assemble!(p) + @test true +end + +""" Calculate mortar projection matrix P = D^-1*M from mortar assembly. """ +function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int) + + C1 = sparse(problem.assembly.C1, ndim, ndim) + C2 = sparse(problem.assembly.C2, ndim, ndim) + + @assert nnz(sparse(problem.assembly.K)) == 0 + @assert nnz(sparse(problem.assembly.D)) == 0 + @assert nnz(sparse(problem.assembly.Kg)) == 0 + @assert nnz(sparse(problem.assembly.fg)) == 0 + @assert nnz(sparse(problem.assembly.f)) == 0 + @assert nnz(sparse(problem.assembly.g)) == 0 + + @assert C1 == C2 + #@assert problem.properties.dual_basis == true + @assert problem.properties.adjust == false + + S = get_nonzero_rows(C2) + M = setdiff(get_nonzero_columns(C2), S) + + # Construct matrix P = D^-1*M + D_ = C2[S,S] + M_ = -C2[S,M] + + #= + P = nothing + if !isdiag(D_) + warn("D is not diagonal, is dual basis used? This might take a long time.") + P = ldltfact(1/2*(D_ + D_')) \ M_ + else + P = D_ \ M_ + end + =# + + P = lufact(D_) \ full(M_) + + return S, M, P +end + +@testset "two linear element clipping, calculation of projection matrix P for standard and dual basis" begin + X = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [0.0, 1.0, 0.0], + 4 => [-0.25, 0.50, 0.00], + 5 => [0.50, -0.25, 0.00], + 6 => [0.75, 0.75, 0.00]) + s = Element(Tri3, [1, 2, 3]) + m = Element(Tri3, [4, 5, 6]) + update!([s, m], "geometry", X) + update!(s, "master elements", [m]) + p = Problem(Mortar, "two elements", 1, "temperature") + p.properties.dual_basis = false + p.elements = [s; m] + initialize!(p) + assemble!(p) + C1 = sparse(p.assembly.C1) + C2 = sparse(p.assembly.C2) + D = sparse(p.assembly.D) + @test length(D) == 0 + @test C1 == C2 + S, M, P = calculate_mortar_projection_matrix(p, 6) + @test S == [1, 2, 3] + @test M == [4, 5, 6] + # visually inspected to be ok result + P_expected = 1/15*[9 9 -3; -7 13 9; 13 -7 9] + @test isapprox(P, P_expected) + um = [7.5, 15.0, 22.5] + @test isapprox(P*um, [9.0, 23.0, 13.0]) + + empty!(p.assembly) + p.properties.dual_basis = true + assemble!(p) + C1 = sparse(p.assembly.C1) + C2 = sparse(p.assembly.C2) + D = sparse(p.assembly.D) + @test length(D) == 0 + @test C1 == C2 + S, M, P = calculate_mortar_projection_matrix(p, 6) + @test S == [1, 2, 3] + @test M == [4, 5, 6] + @test isapprox(P, P_expected) + um = [7.5, 15.0, 22.5] + @test isapprox(P*um, [9.0, 23.0, 13.0]) +end + + +@testset "two quadratic element clipping, calculation of projection matrix P for standard basis" begin + X = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [0.0, 1.0, 0.0], + 7 => [-0.25, 0.50, 0.00], + 8 => [0.50, -0.25, 0.00], + 9 => [0.75, 0.75, 0.00]) + # middle nodes + X[4] = 1/2*(X[1] + X[2]) + X[5] = 1/2*(X[2] + X[3]) + X[6] = 1/2*(X[3] + X[1]) + X[10] = 1/2*(X[7] + X[8]) + X[11] = 1/2*(X[8] + X[9]) + X[12] = 1/2*(X[9] + X[7]) + s = Element(Tri6, [1, 2, 3, 4, 5, 6]) + m = Element(Tri6, [7, 8, 9, 10, 11, 12]) + update!([s, m], "geometry", X) + update!(s, "master elements", [m]) + p = Problem(Mortar, "two elements", 1, "temperature") + p.properties.dual_basis = false + p.properties.alpha = 0.2 + p.elements = [s; m] + initialize!(p) + assemble!(p) + C1 = sparse(p.assembly.C1) + C2 = sparse(p.assembly.C2) + D = sparse(p.assembly.D) + @test length(D) == 0 + @test C1 == C2 + S, M, P = calculate_mortar_projection_matrix(p, 12) + @test S == [1, 2, 3, 4, 5, 6] + @test M == [7, 8, 9, 10, 11, 12] + println(full(P)) + # visually inspected to be ok result + P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396] + @test isapprox(P, P_expected) + um = 15/2*[1, 2, 3] + um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])] + us = P*um + @test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0]) +end + +@testset "two quadratic element clipping, calculation of projection matrix P for dual lagrange basis" begin + X = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [0.0, 1.0, 0.0], + 7 => [-0.25, 0.50, 0.00], + 8 => [0.50, -0.25, 0.00], + 9 => [0.75, 0.75, 0.00]) + # middle nodes + X[4] = 1/2*(X[1] + X[2]) + X[5] = 1/2*(X[2] + X[3]) + X[6] = 1/2*(X[3] + X[1]) + X[10] = 1/2*(X[7] + X[8]) + X[11] = 1/2*(X[8] + X[9]) + X[12] = 1/2*(X[9] + X[7]) + s = Element(Tri6, [1, 2, 3, 4, 5, 6]) + m = Element(Tri6, [7, 8, 9, 10, 11, 12]) + update!([s, m], "geometry", X) + update!(s, "master elements", [m]) + p = Problem(Mortar, "two elements", 1, "temperature") + p.properties.dual_basis = true + p.properties.alpha = 0.2 + p.elements = [s; m] + initialize!(p) + assemble!(p) + C1 = sparse(p.assembly.C1) + C2 = sparse(p.assembly.C2) + D = sparse(p.assembly.D) + @test length(D) == 0 + @test C1 == C2 + S, M, P = calculate_mortar_projection_matrix(p, 12) + P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396] + @test isapprox(P, P_expected) + um = 15/2*[1, 2, 3] + um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])] + us = P*um + @test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0]) +end +