From de3383288a0e917aa874256e30a2d778d89fb7c0 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 10 Feb 2016 22:21:30 +0200 Subject: [PATCH] little changes --- src/elements.jl | 10 ++--- src/fields.jl | 63 ++++++++++++++++++++++++------- src/mortar.jl | 90 ++++++++++++++------------------------------- src/solver_utils.jl | 47 +++++++++++++++-------- src/sparse.jl | 8 ++++ 5 files changed, 120 insertions(+), 98 deletions(-) diff --git a/src/elements.jl b/src/elements.jl index 90fca47..eb7f18b 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -427,19 +427,19 @@ function update!(element::Element, field_name::ASCIIString, data::Dict) element[field_name] = [data[i] for i in get_connectivity(element)] end -function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair}) +function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair}...) element[field_name] = data end """ Update values for several elements at once. """ # FIXME: with or without {T} ? -function update!{T}(elements::Vector{Element{T}}, field_name::ASCIIString, data) +function update!{T}(elements::Vector{Element{T}}, field_name::ASCIIString, data...) for element in elements - update!(element, field_name, data) + update!(element, field_name, data...) end end -function update!(elements::Vector{Element}, field_name::ASCIIString, data) +function update!(elements::Vector{Element}, field_name::ASCIIString, data...) for element in elements - update!(element, field_name, data) + update!(element, field_name, data...) end end diff --git a/src/fields.jl b/src/fields.jl index 1d1dc3e..66e8b52 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -102,10 +102,25 @@ function Field{T}(data::Pair{Float64, Vector{T}}...) return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) end -function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...) +function Base.convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...) return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) end +""" Create new discrete, constant, time variant field. + +Examples +-------- +julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0 +julia> f = DCTV(t0 => y0, t1 => y1) + +""" +function Base.convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...) + return DCTV([Increment(d[1],d[2]) for d in data]) +end +#function Base.convert(::Type{DCTV}, data::Pair{Real, Any}...) +# return DCTV([Increment{Vector}(d[1], d[2]) for d in data]) +#end + function Field(func::Function) if method_exists(func, Tuple{}) return CCTI(func) @@ -176,6 +191,14 @@ function Base.length(field::DCTV) return length(field.data) end +function Base.first(field::Union{DCTV, DVTV}) + return field[1] +end + +function Base.isapprox(f1::DCTI, f2::DCTI) + isapprox(f1.data, f2.data) +end + for op = (:+, :*, :/, :-) @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data) @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data) @@ -263,27 +286,39 @@ function Base.call(field::CCTI, time::Float64) return field.data() end -""" Interpolate time-variant field in time direction. """ -function Base.call(field::DCTV, time::Float64) +""" Interpolate constant time-variant field in time direction. """ +function Base.call(field::DCTV, time::Real) + time < first(field).time && return DCTI(first(field).data) + time > last(field).time && return DCTI(last(field).data) for i=reverse(1:length(field)) - if isapprox(field[i].time, time) - return DCTI(field[i].data) + isapprox(field[i].time, time) && return DCTI(field[i].data) + end + for i=reverse(2:length(field)) + t0 = field[i-1].time + t1 = field[i].time + if t0 < time < t1 + new_data = field[i-1].data + (time-t0)/(t1-t0)*field[i].data + return DCTI(new_data) end end - info(field.data) - info(time) - error("interpolate DCTV: not implemented yet") + error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time") end -function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear) +function Base.call(field::DVTV, time::Float64) + time < first(field).time && return DVTI(first(field).data) + time > last(field).time && return DVTI(last(field).data) for i=reverse(1:length(field)) - if isapprox(field[i].time, time) - return DVTI(field[i].data) + isapprox(field[i].time, time) && return DVTI(field[i].data) + end + for i=reverse(2:length(field)) + t0 = field[i-1].time + t1 = field[i].time + if t0 < time < t1 + new_data = field[i-1].data + (time-t0)/(t1-t0)*field[i].data + return DVTI(new_data) end end - info(field.data) - info(time) - error("interpolate DVTV: not implemented yet") + error("interpolate DVTV: unknown failure when interpolating $(field.data) for time $time") end """ Interpolate constant field in spatial dimension. """ diff --git a/src/mortar.jl b/src/mortar.jl index 756990b..353ea36 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -742,8 +742,6 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor # if distance between elements is "far enough" cannot expect contact if props.contact && (props.minimum_distance < Inf) - #slave_midpoint = slave_element("geometry", [0.0], time) - #master_midpoint = master_element("geometry", [0.0], time) slave_midpoint = Float64[mean(x1[1:field_dim:2]), mean(x1[2:field_dim:2])] master_midpoint = Float64[mean(x2[1:field_dim:2]), mean(x2[2:field_dim:2])] if norm(slave_midpoint - master_midpoint) > props.minimum_distance @@ -772,7 +770,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor Me += w*N'*N end Ae = De*inv(Me) - else + else # Standard Lagrange basis for ip in get_integration_points(slave_element, Val{5}) J = get_jacobian(slave_element, ip, time) w = ip.weight*norm(J)*l @@ -815,21 +813,12 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor G += -(C2S2*X1 - C2M2*X2) # change in weighted gap caused by deformation u += -(C2S2*u1 - C2M2*u2) -# g = G + u - # complementarity condition -# c += la - (G + u) # Add contributions add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2) add!(local_assembly.C1, slave_dofs, master_dofs, -C1M2) add!(local_assembly.C2, slave_dofs, slave_dofs, C2S2) add!(local_assembly.C2, slave_dofs, master_dofs, -C2M2) -# add!(local_assembly.g, slave_dofs, G) -# add!(local_assembly.D, slave_dofs, slave_dofs, D2) -# add!(c_, slave_dofs, c) -# add!(u_, slave_dofs, u) -# add!(la_, slave_dofs, la) -# add!(g_, slave_dofs, g) end # all master elements are done @@ -843,17 +832,22 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor C1 = sparse(local_assembly.C1) C2 = sparse(local_assembly.C2) - D = spzeros(size(C2)...)#copy(C2) + D = spzeros(size(C2)...) g = sparse(local_assembly.g) # complementarity condition -# g = G + u - c = la - (G + u) + lan = la[1:field_dim:end] + lat = la[2:field_dim:end] + Gn = G[1:field_dim:end] + un = u[1:field_dim:end] + cn = lan - (Gn + un) + #Cn = lan - max(0, lan - (Gn+un)) # normal condition - cn = c[1:field_dim:end] inactive_nodes = find(cn .<= 0) active_nodes = find(cn .> 0) + #inactive_nodes = find(Cn .>= 0) + #active_nodes = find(Cn .== 0) # inactive element if length(active_nodes) == 0 @@ -874,55 +868,33 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor g[gdofs] = 0 end - # frictionless contact - if !props.friction - for j in node_ids[active_nodes] - gdofs = [2*(j-1)+1, 2*(j-1)+2] - #D[gdofs[1],:] = 0 - #D[gdofs[2],:] = 0 - D[gdofs[2],gdofs] = C2[gdofs[2],gdofs] - C2[gdofs[2],:] = 0 - g[gdofs[2]] = 0 - end - local_assembly.C1 = C1 - local_assembly.C2 = C2 - local_assembly.D = D - local_assembly.g = g - append!(assembly, local_assembly) - return - end - # frictional contact, see Gitterle2010 mu = 0.3 - lan = la[1:field_dim:end] - lat = la[2:field_dim:end] -# ut = c[2:field_dim:end] ct = lat + c[2:field_dim:end] - #println("cn, ct, lat") - #println(cn) - #println(ct) - #println(lat) - @eval begin - global la = $la - global c = $c - global lat = $lat - global lan = $lan - global ct = $ct - end + C = max(mu*cn, abs(ct)).*lat - mu*max(0, cn).*ct stick_nodes = find(abs(ct) - mu*cn .< 0) slip_nodes = find(abs(ct) - mu*cn .>= 0) stick_nodes = setdiff(stick_nodes, inactive_nodes) slip_nodes = setdiff(slip_nodes, inactive_nodes) - #println("C = ") - #println(C) for (i, j) in enumerate(node_ids[active_nodes]) - gdofs = [2*(j-1)+1, 2*(j-1)+2] - D[gdofs[2],gdofs] = C2[gdofs[2],gdofs] - C2[gdofs[2],:] = 0 - g[gdofs[2]] = C[i] - end + gdofs = [2*(j-1)+1, 2*(j-1)+2] + #D[gdofs[2],gdofs] = C2[gdofs[2],gdofs] + D[gdofs[2],gdofs] = Q[i][:,2]' + C2[gdofs[2],:] = 0 + if props.friction + g[gdofs[2]] = C[i] + else + g[gdofs[2]] = 0.0 + end + end + + local_assembly.C1 = C1 + local_assembly.C2 = C2 + local_assembly.D = D + local_assembly.g = g + append!(assembly, local_assembly) if props.store_debug_info slave_element["G"] = G @@ -934,14 +906,6 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor slave_element["active nodes"] = active_nodes end - local_assembly.C1 = C1 - local_assembly.C2 = C2 - local_assembly.D = D - local_assembly.g = g - #local_assembly.c = c - - append!(assembly, local_assembly) - end typealias MortarElements3D Union{Tri3, Quad4} diff --git a/src/solver_utils.jl b/src/solver_utils.jl index 46c10f7..61ad64f 100644 --- a/src/solver_utils.jl +++ b/src/solver_utils.jl @@ -45,7 +45,7 @@ function pretty_print_C1_row(r) return s end -function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2, D_, D, g_, g) +function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2, D_, D, g_, g; show_info=false) # old, new, old, new... #= herzian contact with symmetry boundary @@ -69,7 +69,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 # INFO: algorithm 2 solved issue? true # INFO: fixed: new setting is # INFO: dof 1109: 0.0*u₁₅₃ - 0.0*u₁₅₄ - 0.0*u₁₅₅ + 0.15*u₁₅₆ + 0.0*u₁₁₀₉ - 0.15*u₁₁₁₀ = -0.0 - +=# if 555 in nodes info("overconstraint DIRTY HACK") # It is possible to selectively remove mortar constraints and the associated @@ -79,15 +79,30 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 # new configuration is mortar constraint # 1. remove mortar constrains and associated Lagrange multiplier components # in dof 1109, that is, x direction of node 555. - C1[1109,:] = 0 - C2[1109,:] = 0 - D[1109,:] = 0 - g[1109,:] = 0 -# D[1110,1110] = 0 -# g[1110] = 0 + + # works quite well +# C1_[1109,:] = 0 +# C2_[1110,:] = C2_[1109,:] +# C2[1110,:] = 0 +# D[1110,:] = 0 + + #C1_[1110,:] = C1_[1109,:] + C2_[1110,:] = C2_[1109,:] + #C1_[1109,:] = 0 + C2_[1109,:] = 0 + + #C1[1110,:] = 0 + #C2[1110,:] = 0 + #C2[:,1109] = 0 + D[1110,:] = 0 + #D[:,1109] = 0 + g[1110,:] = 0 + #D[:,1109] = 0 + #C2[:,1109] = 0 + #C1[1109,:] = 0 + #C1_[1110,:] = 0 return end -=# """ Return all other dofs which connects to overconstrained dofs. """ function get_related_dofs(dofs) @@ -297,7 +312,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 for node_id in nodes dofs = [2*(node_id-1)+1, 2*(node_id-1)+2] - print_summary(node_id, dofs) + show_info && print_summary(node_id, dofs) # try to resolve issue automatically resolved = false @@ -311,16 +326,16 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 if resolved info("fixed: new setting is") - show_rows_in_constraint_matrix(dofs, C2, D; show_status=false) - show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false) - #show_related_equations(dofs, C2, C2_, D, D_) - info() + show_info && show_rows_in_constraint_matrix(dofs, C2, D; show_status=false) + show_info && show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false) + show_info && show_related_equations(dofs, C2, C2_, D, D_) + show_info && info() continue end - info("unable to resolve overconstrained situation, not continuing") + show_info && info("unable to resolve overconstrained situation, not continuing") throw("failed to resolve overconstraint situation") - info() + show_info && info() end end diff --git a/src/sparse.jl b/src/sparse.jl index 4445258..ac9527f 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -20,6 +20,14 @@ function Base.convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC) return SparseMatrixCOO(findnz(A)...) end +function Base.convert(::Type{SparseMatrixCOO}, A::Matrix) + return SparseMatrixCOO(findnz(A)...) +end + +function Base.convert(::Type{SparseMatrixCOO}, A::Vector) + return SparseMatrixCOO(findnz(sparse(A))...) +end + """ Convert from COO format to CSC. Parameters