From d24b1e178a92aecef52d00ae0190d8950011fb30 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 15 Dec 2016 09:18:15 +0200 Subject: [PATCH] a new strategy to deal with quadratic elements, change order of surface elements to linear --- src/problems_mortar.jl | 28 ++++++++++++++- src/problems_mortar_3d.jl | 74 +++++++++++++++++++++++++++++---------- 2 files changed, 82 insertions(+), 20 deletions(-) diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl index da2bb43..fe12447 100644 --- a/src/problems_mortar.jl +++ b/src/problems_mortar.jl @@ -1,6 +1,31 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +""" +Parameters +---------- +dimension + dimension of surface, 1 for 2d problems (plane strain, plane stress, + axisymmetric) and 2 for 3d problems. It not given, try to determine problem + dimension from first element +rotate_normals + if all surface elements are in cw order instead of ccw, this can be used to + swap normal directions so that normals point to outward of body +adjust + for elasticity problems only; closes any gaps between surfaces if found +dual_basis + use bi-orthogonal basis when interpolating Lagrange multiplier space +use_forwarddiff + use forwarddiff to linearize contact constraints directly from weighted + gap function +distval + charasteristic measure, contact pairs with distance over this value are + skipped from contact segmentation algorithm +linear_surface_elements + convert quadratic surface elements to linear elements on the fly +store_fields + not used +""" type Mortar <: BoundaryProblem dimension :: Int rotate_normals :: Bool @@ -8,12 +33,13 @@ type Mortar <: BoundaryProblem dual_basis :: Bool use_forwarddiff :: Bool distval :: Float64 + linear_surface_elements :: Bool store_fields :: Vector{Symbol} end function Mortar() default_fields = [] - return Mortar(-1, false, false, false, false, Inf, default_fields) + return Mortar(-1, false, false, false, false, Inf, true, default_fields) end function get_unknown_field_name(problem::Problem{Mortar}) diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl index e43f32d..f04f9bd 100644 --- a/src/problems_mortar_3d.jl +++ b/src/problems_mortar_3d.jl @@ -105,7 +105,7 @@ function contains{T}(P::Vector{T}, q::T; check_is_close=true, rtol=1.0e-4) return false end -function get_polygon_clip(xs, xm, n; debug=false) +function get_polygon_clip(xs, xm, n) # objective: search does line xm1 - xm2 clip xs nm = length(xm) ns = length(xs) @@ -114,7 +114,7 @@ function get_polygon_clip(xs, xm, n; debug=false) # 1. test is master point inside slave, if yes, add to clip for i=1:nm if vertex_inside_polygon(xm[i], xs) - debug && info("1. $(xm[i]) inside S -> push") + debug("1. $(xm[i]) inside S -> push") push!(P, xm[i]) end end @@ -123,7 +123,7 @@ function get_polygon_clip(xs, xm, n; debug=false) for i=1:ns if vertex_inside_polygon(xs[i], xm) contains(P, xs[i]) && continue - debug && info("2. $(xs[i]) inside M -> push") + debug("2. $(xs[i]) inside M -> push") push!(P, xs[i]) end end @@ -146,7 +146,7 @@ function get_polygon_clip(xs, xm, n; debug=false) #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))") if vertex_inside_polygon(q, xm) contains(P, q) && continue - debug && info("3. $q inside M -> push") + debug("3. $q inside M -> push") push!(P, q) end end @@ -155,18 +155,21 @@ function get_polygon_clip(xs, xm, n; debug=false) return P end +""" Project some vertex p to surface of element E using Newton's iterations. """ function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector, - element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-6) + element::Element{E}, x::DVTI, time::Real; + max_iterations::Int=10, iter_tol::Float64=1.0e-6) basis(xi) = get_basis(element, xi, time) dbasis(xi) = get_dbasis(element, xi, time) + nnodes = length(element) f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p L(theta) = inv3([dbasis(theta[1:2])*x -n0]) -# L2(theta) = inv(ForwardDiff.get_value([dbasis(theta[2:3])*x -n0])) - # FIXME: for some reason forwarddiff gives NaN's here. theta = zeros(3) dtheta = zeros(3) for i=1:max_iterations - dtheta = L(theta) * f(theta) + invA = L(theta) + b = f(theta) + dtheta = invA * b theta -= dtheta if norm(dtheta) < iter_tol return theta[1:2], theta[3] @@ -227,12 +230,12 @@ function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false) return normals end -function check_orientation!(P, n; debug=false) +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 && info("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) @@ -246,7 +249,20 @@ function check_orientation!(P, n; debug=false) end) end -function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}; debug=true) +function convert_to_linear_element{E}(element::Element{E}) + debug("No linear convert rule for element $E") + return element +end + +function convert_to_linear_element(element::Element{Tri6}) + debug("converting Tri6 to Tri3") + new_element = Element(Tri3, element.connectivity[1:3]) + new_element.id = element.id + new_element.fields = element.fields + return new_element +end + +function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}) props = problem.properties field_dim = get_unknown_field_dimension(problem) @@ -257,27 +273,37 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ # 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", normals) + update!(slave_elements, "normal", time => normals) # 2. loop all slave elements + first_slave_element = true + for slave_element in slave_elements + 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 = get_reference_element_midpoint(slave_element) - xi = [1/3, 1/3] + 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(p, x0, n0) for p in X1] + S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] # 3. loop all master elements for master_element in slave_element("master elements", time) + 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) @@ -288,12 +314,20 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ end # 3.1 project master nodes to auxiliary plane and create polygon clipping - M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2] + 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 @@ -315,7 +349,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ # 4. loop integration cells all_cells = get_cells(P, C0) for cell in all_cells - virtual_element = Element(Tri3) + virtual_element = Element(Tri3, Int[]) update!(virtual_element, "geometry", cell) #x_cell = Field(cell) @@ -428,11 +462,13 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ end # master elements done + first_slave_element = false + end # slave elements done, contact virtual work ready if problem.properties.dual_basis tol = 1.0e-9 - debug && info("Dual basis is used, dropping small values for C1 & C2, tol = $tol") + 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) @@ -441,7 +477,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ problem.assembly.C2 = C2 end - debug && info("area of interface: $area") + debug("area of interface: $area") end