little changes

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

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. """

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}

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
 

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