diff --git a/geometry/2d_curved_interface/MESH_UPPER_BLOCK_SMALLER.med b/geometry/2d_curved_interface/MESH_UPPER_BLOCK_SMALLER.med new file mode 100644 index 0000000..f40ad60 Binary files /dev/null and b/geometry/2d_curved_interface/MESH_UPPER_BLOCK_SMALLER.med differ diff --git a/geometry/2d_curved_interface/PUNCH_AND_BLOCK.med b/geometry/2d_curved_interface/PUNCH_AND_BLOCK.med new file mode 100644 index 0000000..101ec9a Binary files /dev/null and b/geometry/2d_curved_interface/PUNCH_AND_BLOCK.med differ diff --git a/geometry/2d_curved_interface/UPPER_SMALLER.med b/geometry/2d_curved_interface/UPPER_SMALLER.med new file mode 100644 index 0000000..b0ec884 Binary files /dev/null and b/geometry/2d_curved_interface/UPPER_SMALLER.med differ diff --git a/geometry/2d_curved_interface/model.hdf b/geometry/2d_curved_interface/model.hdf index 322b7ed..fcbccc6 100644 Binary files a/geometry/2d_curved_interface/model.hdf and b/geometry/2d_curved_interface/model.hdf differ diff --git a/geometry/2d_finite_sliding_test/TWOBLOCKS.med b/geometry/2d_finite_sliding_test/TWOBLOCKS.med new file mode 100644 index 0000000..ee85364 Binary files /dev/null and b/geometry/2d_finite_sliding_test/TWOBLOCKS.med differ diff --git a/geometry/2d_finite_sliding_test/model.hdf b/geometry/2d_finite_sliding_test/model.hdf new file mode 100644 index 0000000..a25ba06 Binary files /dev/null and b/geometry/2d_finite_sliding_test/model.hdf differ diff --git a/src/elements.jl b/src/elements.jl index 608cfc8..4d59834 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -132,6 +132,35 @@ function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0) return get_basis(element, xi) end +""" Given a list of elementa and nodes, find a subset of elements +containing nodes. +""" +function find_elements(elements, nodes) + s = Set{Element}() + for element in elements + conn = get_connectivity(element) + for j in nodes + if j in conn + push!(s, element) + break + end + end + end + return collect(s) +end + +function get_gdofs(element::Element) + return get_gdofs(element, 1) +end + +function get_dbasis{E}(element::Element{E}, ip::IntegrationPoint) + return get_dbasis(E, ip.xi) +end + +function get_basis{E, T<:Real}(element::Element{E}, xi::T) + return get_basis(E, xi) +end + """ Return dual basis transformation matrix Ae. """ function get_dualbasis(element::Element, time::Real) if length(element.A) == 0 diff --git a/src/fields.jl b/src/fields.jl index 66e8b52..38857e0 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -216,6 +216,15 @@ function Base.(:*){T<:Real}(c::T, field::DVTI) return DVTI(c*field.data) end +""" Multiply DVTI field with another vector T. Vector length +must match to the field length and this can be used mainly +for interpolation purposes, i.e., u = ∑ Nᵢuᵢ +""" +function Base.(:*)(T::Vector, f::DVTI) + @assert length(T) == length(f) + return sum([T[i]*f[i] for i=1:length(f)]) +end + function Base.vec(field::DVTI) return [field.data...;] end diff --git a/src/lagrange.jl b/src/lagrange.jl index 98b278a..4a33f8c 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -140,3 +140,11 @@ end (xi) -> [ 1.0, xi[1], xi[2], xi[3], xi[1]^2, xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]]) + +# some helpers to make accessing 1d basis functions more easily +function get_basis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T) + get_basis(E, [xi]) +end +function get_dbasis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T) + get_dbasis(E, [xi]) +end diff --git a/src/mortar.jl b/src/mortar.jl index c2aee87..e1c8246 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -19,7 +19,7 @@ b) Remove inactive inequality constraints in assembly level. This is done in """ type Mortar <: BoundaryProblem - formulation :: Symbol # :total or :incremental + formulation :: Symbol # :total, :incremental, :autodiff dual_basis :: Bool inequality_constraints :: Bool # Launch PDASS to solve inequality constraints normal_condition :: Symbol # Tie or Contact @@ -31,10 +31,11 @@ type Mortar <: BoundaryProblem always_in_slip :: Vector{Int64} # nodes in this list always in slip contact :: Bool friction :: Bool + gap_sign :: Int # gap sign convention end function Mortar() - Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false) + Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false, -1) end function get_unknown_field_name(::Type{Mortar}) @@ -51,719 +52,8 @@ macro debug(msg) end include("mortar_2d.jl") - - -### Mortar projection calculation for 3d cases - -""" -Construct auxiliary plane for surface. - -Parameters ----------- -x::Array{Float64, 2} - Node coordinates -ximp::Array{Float64, 1} - Element mid-point in dimensionless mother element coordinates ξ -normals::Array{Float64, 2} - Normal directions in nodes - -Returns -------- -x0, Q - x0::Array{Float64, 1} - origo of auxiliary plane - Q::Array{Float64, 2} - orthogonal basis, first vector is normal direction - and two rest vectors create orthonormal right-handed basis. - -Examples --------- -Calculate auxiliary plane given nodal coordinates, midpoint of mother element, -node normals and suitable function space: - -julia> xquad = [ -... -2.5 2.5 2.0 -2.0 -... -2.0 -2.0 2.3 2.0 -... 1.0 0.7 0.0 1.0] -julia> m_midpoint = [0.0, 0.0] -julia> normals = [ -... 0.05989060 0.0590504 0.225612 0.2445800 -... -0.00748633 0.1670810 0.182034 -0.0305725 -... 0.99817700 0.9841730 0.957059 0.9691470] -julia> basis(xi) = [ -... (1-xi[1])(1-xi[2])/4 -... (1+xi[1])(1-xi[2])/4 -... (1+xi[1])(1+xi[2])/4 -... (1-xi[1])(1+xi[2])/4]' -julia> x0, Q = create_auxiliary_plane(xquad, mmidpoint, normals, basis) -julia> x0 -3-element Array{Float64,1}: - 0.0 - 0.075 - 0.675 -julia> Q -3x3 Array{Float64,2}: - 0.148586 0.988899 0.0 - 0.0784519 -0.0117877 0.996848 - 0.985783 -0.148118 -0.0793325 - -Notes ------ -- Midpoint in mother element typically (0, 0) for quadrangles and (1/3, 1/3) - for triangles. -- Uses Gram-Schmidt process to find orthogonal basis -- [1](http://www.math.umn.edu/~olver/aims_/qr.pdf) -- [2](http://www.ecs.umass.edu/ece/ece313/Online_help/gram.pdf) -- [3](http://www.terathon.com/code/tangent.html) - -""" -# function create_auxiliary_plane(x, ximp, normals, basis) -function create_auxiliary_plane{E}(element::Element{E}, time::Real) -# proj(u, v) = dot(v, u) / dot(u, u) * u -# xi = [1.0/3.0, 1.0/3.0] - - xi = get_reference_element_midpoint(E) - x0 = element("geometry", xi, time) - ntbasis = element("normal-tangential coordinates", xi, time) - return x0, ntbasis -#= - n = element("normal-tangential coordinates", xi, time)[:, 1] - n /= norm(n) - # gram-schmidt - u1 = n - j = indmax(abs(u1)) - v2 = zeros(3) - v2[mod(j,3)+1] = 1.0 - u2 = v2 - proj(u1, v2) - u3 = cross(u1, u2) - t1 = u2/norm(u2) - t2 = u3/norm(u3) - new_basis = [n t1 t2] - return x0, new_basis -=# -end - -""" -Project point q onto a plane given by a point p and normal n. - -Parameters ----------- -q::Array{Float64, 2} - point to project (row vector) -x0::Array{Float64, 2} - origo of plane -n::Array{Float64, 2} - normal vector of plane - -Returns -------- -y::Array{Float64, 2} - projected point - -Examples --------- -julia> p = [-0.5 -1.0 4.0]' -julia> x0 = [0.0 0.075 0.675]' -julia> n = [0.1485860 0.0784519 0.9857830]' -julia> project_node_to_auxiliary_plane(p, x0, n) -3-element Array{Float64,1}: - 0.963455 - -1.2447 - 0.925247 - -Notes ------ -[1](http://stackoverflow.com/questions/8942950/how-do-i-find-the-orthogonal-projection-of-a-point-onto-a-plane) - -""" -function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix) - n = Q[:,1] - ph = p - dot(p-x0, n)*n - qproj = Q'*(ph-x0) - if !isapprox(qproj[1], 0.0; atol=1.0e-12) - info("project_point_to_auxiliary_plane(): point not projected correctly.") - info("p: $p") - info("x0: $x0") - info("Q: \n$Q") - info("qproj: $qproj") - error("Failed to project point to auxiliary plane.") - end - return qproj[2:3] -end - -""" -Find edge intersections of two planar arbitrary shape polygons. - -Parameters ----------- -S::Array{Float64,2} -M::Array{Float64,2} - -Matrices with size (2, n) where n is number of vertices of each polygon. - -Returns -------- -P::Array{Float64,2} - Intersection points of polygons -n::Array{Float64,2} - Neighbour info matrix with size (ns, mn). This keeps information which - edges of polygons are intersecting. See further explanation in example - below. - -Examples --------- -Find intersection points of two triangles: - -julia> S = [0 0; 3 0; 0 3]' -julia> M = [-1 1; 2 -1/2; 1 3/2]' -julia> P, n = get_edge_intersections(S, M) -julia> P -2x4 Array{Float64,2}: - 1.0 1.75 0.0 0.0 - 0.0 0.0 0.5 1.25 -julia> n -3x3 Array{Int64,2}: - 1 1 0 - 0 0 0 - 1 0 1) - -So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25). -"Neighbour matrix" can be interpreted as following: - - 1 1 0 <--> First edge of S intersects edges 1 and 2 of M - 0 0 0 <--> Second edge of S doesn't intersect at all - 1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M - -""" -function get_edge_intersections(S::Matrix, M::Matrix) - ns = size(S, 2) - nm = size(M, 2) - P = zeros(2, 0) - n = zeros(Int64, ns, nm) - k = 0 - for i=1:ns - for j=1:nm - b = M[:,j]-S[:,i] - A = [S[:,mod(i,ns)+1]-S[:,i] -M[:,mod(j,nm)+1]+M[:,j]] - if rank(A) == 2 - r = A\b - if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found - k += 1 - f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i]) - f = f'' - P = hcat(P, f) - n[i, j] = 1 - end - end - end - end - return P, n -end - - -""" -Find any points laying inside or border of triangle. - -Parameters ----------- -Y::Array{Float64, 2} - Triangle coordinates in 2×3 matrix -X::Array{Float64, 2} - List of points to test in 2×n matrix - -Returns -------- -P::Array{Float64, 2} - List of points in triangle in 2×m matrix, where m is number of points inside triangle - -Examples --------- -julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points -julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests -julia> points_in_triangle(S, pts) -2x2 Array{Float64,2}: - 1.0 0.5 - 1.5 1.5 -""" -function get_points_inside_triangle(Y::Matrix, X::Matrix) - @assert size(Y, 2) == 3 # "Point in TRIANGLE..." - P = zeros(2, 0) - v0 = Y[:,2] - Y[:,1] - v1 = Y[:,3] - Y[:,1] # find interior points of X in Y - d00 = (v0'*v0)[1] - d01 = (v0'*v1)[1] - d11 = (v1'*v1)[1] # using baricentric coordinates - id = 1/(d00*d11 - d01*d01) - for i=1:size(X, 2) - v2 = X[:,i] - Y[:,1] - d02 = (v0'*v2)[1] - d12 = (v1'*v2)[1] - u = (d11*d02-d01*d12)*id - v = (d00*d12-d01*d02)*id - if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary - P = hcat(P, X[:,i]'') - end - end - return P -end - -""" -Determine is point P inside or on boudary of polygon X. - -http://paulbourke.net/geometry/polygonmesh/#insidepoly -""" -function is_point_inside_convex_polygon(P, X) - x, y = P - for i=1:length(X) - x0, y0 = X[i] - x1, y1 = X[mod(i, length(X))+1] - if (y-y0)*(x1-x0) - (x-x0)*(y1-y0) < 0 - return false - end - end - return true -end - -function get_points_inside_convex_polygon(pts, X) - # TODO: Make more readable - X2 = [X[:,i] for i=1:size(X,2)] - c = filter(P->is_point_inside_convex_polygon(P, X2), [pts[:,i] for i=1:size(pts, 2)]) - return length(c) == 0 ? zeros(2, 0) : hcat(c...) -end - -""" Return unique objects with some given tolerance. This is used in next function - because traditional unique() command returns row vectors as non-unique if they - differs only a "little". -""" -function uniquetol(P, dim::Int; args...) - @assert dim == 2 - items = Vector{Float64}[P[:,i] for i=1:size(P,dim)] - new_items = Vector{Float64}[] - for item in items - has_found = false - for new_item in new_items - if isapprox(item, new_item; args...) - has_found = true - break - end - end - if !has_found - push!(new_items, item) - end - end - return reshape([new_items...;], length(new_items[]), length(new_items)) -end - - -""" -Make polygon clipping of shapes S and M. - -Parameters ----------- -S::Array{Float64, 2} -M::Array{Float64, 2} - Shapes to clip. Needs to be triangles at the moment. - -Returns -------- -Array{Float64, 2}, Array{Float64, 2} -- Polygon vertices in 2×n matrix, sorted in counter-clockwise order. -- 3×3 "neighbouring" matrix, see example. - -Examples --------- -julia> S = [0 0; 3 0; 0 3]' -julia> M = [-1 1; 2 -1/2; 2 2]' -julia> P, n = clip_polygon(S, M) -julia> P -2x6 Array{Float64,2}: - 0.0 1.0 2.0 2.0 1.25 0.0 - 0.5 0.0 0.0 1.0 1.75 1.33333, -julia> n -3x3 Array{Int64,2}: - 1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]') - 1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S - 0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S - -""" -function clip_polygon(S::Matrix, M::Matrix) - P1, neighbours = get_edge_intersections(M, S) - #P2 = get_points_inside_triangle(M, S) - #P3 = get_points_inside_triangle(S, M) - P2 = get_points_inside_convex_polygon(M, S) - P3 = get_points_inside_convex_polygon(S, M) -# info("polygon clipping: P1 = $P1") -# info("polygon clipping: P2 = $P2") -# info("polygon clipping: P3 = $P3") -# info("hcat P = $P") - P = hcat(P1, P2, P3) - if length(P) == 0 - return nothing, nothing - end - P = uniquetol(P, 2) - meanval = mean(P, 2) - tmp = P .- meanval - angles = atan2(tmp[2,:], tmp[1,:]) - angles = reshape(angles, length(angles)) - order = sortperm(angles) - return P[:, order], neighbours -end - - -""" -Calculate polygon geometric center point - -Parameters ----------- -P::Array{Float64, 2} - Polygon vertices in 2×n matrix - -Returns -------- -Array{Float63, 2} - Center point - -Examples --------- -julia> P -2x6 Array{Float64,2}: - 0.0 1.0 2.0 2.0 1.25 0.0 - 0.5 0.0 0.0 1.0 1.75 1.33333, -julia> C = get_polygon_cp(P) -2x1 Array{Float64,2}: - 1.039740 - 0.804701 - -""" -function calculate_polygon_centerpoint(P::Matrix) - n = size(P, 2) - A = 0.0 - for i=1:n - A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i]) - end - Cx = 0.0 - Cy = 0.0 - for i=1:n - inext = mod(i, n)+1 - Cx += 1/(6*A)*(P[1,i] + P[1,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) - Cy += 1/(6*A)*(P[2,i] + P[2,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) - end - return Float64[Cx, Cy] -end - -""" -Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂) - -Parameters ----------- -p::Array{Float64,1} - point in auxiliary plane, in (n,t1,t2) coordinate system -x0::Array{Float64,1} - origo of auxiliary plane cs -Q::Array{Float64,2} - basis of auxiliary plane cs -x::Array{Float64,2} - surface node coords -basis::Array{Float64,2} - surface basis functions -dbasis::Array{Float64,2} - partial derivatives of surface basis functions - -Returns -------- -Array{Float64,2} - solution vector (d, ξ₁, ξ₂) where d is distance to surface - -Examples --------- -Define surface with node points, basis + dbasis - -julia> xquad = [ -... -2.5 -2.0 1.0 -... 2.5 -2.0 0.7 -... 2.0 2.3 0.0 -... -2.0 2.0 1.0]' -julia> basis(xi) = [ -... (1-xi[1])(1-xi[2])/4 -... (1+xi[1])(1-xi[2])/4 -... (1+xi[1])(1+xi[2])/4 -... (1-xi[1])(1+xi[2])/4] -julia> dbasis(xi) = [ -... -(1-xi[2])/4 -(1-xi[1])/4 -... (1-xi[2])/4 -(1+xi[1])/4 -... (1+xi[2])/4 (1+xi[1])/4 -... -(1+xi[2])/4 (1-xi[1])/4] - -We aim to find point p, which we first project to auxiliary plane defined as following -julia> p = [-2.5 -2.0 1.0]' -julia> x0 = [0.0 0.075 0.675]' -julia> Q = [ -... 0.1485860 0.9888990 0.0000000 -... 0.0784519 -0.0117877 0.9968480 -... 0.9857830 -0.1481180 -0.0793325] - -Our projected point is therefore -julia> n = Q[:,1] # first component is normal direction -julia> ph = project_node_to_auxiliary_plane(p, x0, n) -julia> ph = Q'(ph-x0) -julia> ph -3x1 Array{Float64,2}: - 1.33264e-7 - -2.49593 - -2.09424 - -Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we -project it back to surface defined by xquad*basis - -julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis) -julia> theta -3x1 Array{Float64,2}: - -0.213874 - -0.999999 - -1.0 - -We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad -[-2.5 -2.0 1.0]' correctly. - -julia> xquad*basis(theta[2:3]) -3-element Array{Float64,1}: - -2.5 - -2.0 - 1.0 - -""" -function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) - basis(xi) = get_basis(E, xi) - dbasis(xi) = get_dbasis(E, xi) - x = element("geometry", time) - ph = Q*[0; p] + x0 - theta = Float64[0.0, 0.0, 0.0] - n = Q[:,1] - for i=1:max_iterations - b = ph + theta[1]*n - basis(theta[2:3])*x - J = [n -dbasis(theta[2:3])*x] - dtheta = J \ -b - theta += dtheta - if norm(dtheta) < iter_tol - return theta - end - end - begin - info("projecting point from auxiliary plane back to surface didn't go very well.") - info("element type: $E") - info("element connectivity: $(get_connectivity(element))") - info("auxiliary plane: x0 = $x0, Q = $Q") - info("point coordinates on plane: $p") - info("element geometry: $x") - info("ph: $ph") - info("normal direction: $n") - info("parameter vector before giving up: $theta") - end - error("project_point_to_surface: did not converge in $max_iterations iterations!") -end - -typealias MortarElements3D Union{Tri3, Quad4} - -function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar}, - slave_element::Element{E}, time::Real) - haskey(slave_element, "master elements") || return - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - slave_dofs = get_gdofs(slave_element, field_dim) - - props = problem.properties - if props.formulation == :Standard && props.normal_condition == :Contact - error("for contact choose Dual formulation.""") - end - - # create auxiliary plane and project slave nodes to it - # x0 = origo, Q = local basis - x0, Q = create_auxiliary_plane(slave_element, time) - - # 1. project slave nodes to auxiliary plane - Sl = Vector{Float64}[] - for p in slave_element("geometry", time) - push!(Sl, project_point_to_auxiliary_plane(p, x0, Q)) - end - S = hcat(Sl...) - - for master_element in slave_element["master elements"] - - # if distance between elements is "far enough" cannot expect contact - if (props.normal_condition == :Contact) || props.inequality_constraints - slave_midpoint = slave_element("geometry", [0.0, 0.0], time) - master_midpoint = master_element("geometry", [0.0, 0.0], time) - if norm(slave_midpoint - master_midpoint) > props.minimum_distance - continue - end - end - - master_dofs = get_gdofs(master_element, field_dim) - - # 2. project master nodes to auxiliary plane - M = Vector{Float64}[] - for p in master_element("geometry", time) - push!(M, project_point_to_auxiliary_plane(p, x0, Q)) - end - M = hcat(M...) - - # 3. create polygon clipping on auxiliary plane - P = nothing - neighbours = nothing - try - P, neighbours = clip_polygon(S, M) - catch - info("polygon clipping failed") - info("S = ") - dump(S) - info("M = ") - dump(M) - info("original Sl = ") - info(Sl) - error("cannot continue") - end - isa(P, Void) && continue # no clipping - - # shared edge but no shared volume. skipping - size(P, 2) < 3 && continue - - C = calculate_polygon_centerpoint(P) - npts = size(P, 2) # number of vertices in polygon - - # loop vertices and create temporary integrate cells - # TODO: basically when npts == 3 or npts == 4 we could integrate without splitting to cells. - nnodes = size(slave_element, 2) - C1S3 = zeros(3*nnodes, 3*nnodes) - C1M3 = zeros(3*nnodes, 3*nnodes) - - for pnt=1:npts # integration of mortar matrices begin - cell = Field(Vector{Float64}[C, P[:,pnt], P[:,mod(pnt,npts)+1]]) - - # calculate slave side projection matrix D - # construct dual basis - Ae = zeros(nnodes, nnodes) - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - if problem.properties.formulation == :Dual # Construct dual basis - for ip in get_integration_points(Tri3, Val{5}) - N = get_basis(Tri3, ip.xi) - xi = vec(N*cell) - theta = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time) - xi_slave = theta[2:3] - N1 = slave_element(xi_slave, time) - # jacobian determinant on integration cell - dNC = get_dbasis(Tri3, ip.xi) - JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)]) - wC = ip.weight*det(JC) - De += wC*diagm(vec(N1)) - Me += wC*N1'*N1 - end - Ae = De*inv(Me) - end - for i=1:field_dim - C1S3[i:field_dim:end,i:field_dim:end] += De - end - - # Calculate master side projection matrix M - for ip in get_integration_points(Tri3, Val{5}) - # gauss point in auxiliary plane - #N = get_basis(E, ip.xi) - N = get_basis(Tri3, ip.xi) - xi = vec(N*cell) # xi defined in auxilary plane - - # find projection of gauss point to master and slave elements - theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time) - theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time) - xi_slave = theta1[2:3] - xi_master = theta2[2:3] - - # evaluate shape functions values in gauss point and add contribution to matrices - N1 = slave_element(xi_slave, time) - N2 = master_element(xi_master, time) - - # jacobian determinant on integration cell - dNC = get_dbasis(Tri3, ip.xi) - JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)]) - wC = ip.weight*det(JC) - - # extend matrices according to the problem dimension (3) - @assert length(slave_dofs) == length(master_dofs) - Me = wC*Ae*N1'*N2 - for k=1:field_dim - C1M3[k:field_dim:end,k:field_dim:end] += Me - end - end - end # integration of mortar matrices done. - - # constraints in normal-tangential direction and initial weighted gap - X1 = vec(slave_element("geometry", time)) - X2 = vec(master_element("geometry", time)) - Q_ = slave_element("normal-tangential coordinates", time) - Z = zeros(3, 3) - if nnodes == 3 - Q3 = [Q Z Z; Z Q Z; Z Z Q] - elseif nnodes == 4 - Q3 = [Q Z Z Z; Z Q Z Z; Z Z Q Z; Z Z Z Q] - end - D3 = zeros(3*nnodes, 3*nnodes) - C2S3 = Q3'*C1S3 - C2M3 = Q3'*C1M3 - G = -(C2S3*X1 - C2M3*X2) - - # complementarity condition - if haskey(slave_element, "displacement") - u1 = vec(slave_element("displacement", time)) - else - u1 = zeros(3*nnodes) - end - if haskey(master_element, "displacement") - u2 = vec(master_element("displacement", time)) - else - u2 = zeros(3*nnodes) - end - x1 = X1 + u1 - x2 = X2 + u2 - if haskey(slave_element, "reaction force") - la = vec(slave_element("reaction force", time)) - else - la = zeros(3*nnodes) - end - g = -(C2S3*x1 - C2M3*x2) - c = Q3'*la - g - inactive_nodes = find(c[1:field_dim:end] .<= 0) - active_nodes = find(c[1:field_dim:end] .> 0) - - # normal constraint: remove inactive nodes if normal condition is set to contact - if problem.properties.normal_condition == :Contact - for j in inactive_nodes - dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] - G[dofs] = 0 - C1S3[dofs,:] = 0 - C1M3[dofs,:] = 0 - C2S3[dofs,:] = 0 - C2M3[dofs,:] = 0 - end - end - - # tangential constraint: stick or slip - if problem.properties.tangential_condition == :Slip - D3 = copy(C2S3) - D3[1:field_dim:end, :] = 0 - C2S3[2:field_dim:end, :] = 0 - C2M3[2:field_dim:end, :] = 0 - C2S3[3:field_dim:end, :] = 0 - C2M3[3:field_dim:end, :] = 0 - end - - # add contributions - add!(assembly.C1, slave_dofs, slave_dofs, C1S3) - add!(assembly.C1, slave_dofs, master_dofs, -C1M3) - add!(assembly.C2, slave_dofs, slave_dofs, C2S3) - add!(assembly.C2, slave_dofs, master_dofs, -C2M3) - add!(assembly.D, slave_dofs, slave_dofs, D3) - add!(assembly.c, slave_dofs, c) - add!(assembly.g, slave_dofs, G) - end -end - +include("mortar_2d_autodiff.jl") +include("mortar_3d.jl") """ Remove inactive inequality constraints by using primal-dual active set strategy. """ function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g) diff --git a/src/mortar_2d.jl b/src/mortar_2d.jl index 64845eb..74da3dc 100644 --- a/src/mortar_2d.jl +++ b/src/mortar_2d.jl @@ -217,7 +217,6 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M} error("find projection from master to slave: did not converge") end - # Mortar assembly 2d # quadratic not tested yet @@ -546,8 +545,8 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor # Calculate normal-tangential constraints and weighted gap C2S2 = Q2'*C1S2 C2M2 = Q2'*C1M2 - G += -(C2S2*X1 - C2M2*X2) - g += -(C2S2*x1 - C2M2*x2) + G += props.gap_sign*(C2S2*X1 - C2M2*X2) + g += props.gap_sign*(C2S2*x1 - C2M2*x2) # Add contributions add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2) diff --git a/src/mortar_2d_autodiff.jl b/src/mortar_2d_autodiff.jl new file mode 100644 index 0000000..b4d2cbb --- /dev/null +++ b/src/mortar_2d_autodiff.jl @@ -0,0 +1,674 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" Find segment from slave element corresponding to master element nodes. +x1_, n1_ +slave element geometry and normal direction + +x2_ master element nodes to project onto slave +""" +function project_from_master_to_slave{E<:MortarElements2D}( + slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector; + tol=1.0e-10, max_iterations=20) + + function x1(xi1) + N = get_basis(E, xi1) + return vec(N)*x1_ + end + + function dx1(xi1) + dN = get_dbasis(E, xi1) + return vec(dN)*x1_ + end + + function n1(xi1) + N = get_basis(E, xi1) + return vec(N)*n1_ + end + + function dn1(xi1) + dN = get_dbasis(E, xi1) + return vec(dN)*n1_ + end + + cross2(a, b) = cross([a; 0], [b; 0])[3] + R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) + dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) + xi1 = 0.0 + dxi1 = 0.0 + for i=1:max_iterations + dxi1 = -R(xi1)/dR(xi1) + xi1 += dxi1 + if norm(dxi1) < tol + return xi1 + end + end + + info("x1 = $(ForwardDiff.get_value(x1_.data))") + info("n1 = $(ForwardDiff.get_value(n1_.data))") + info("x2 = $(ForwardDiff.get_value(x2))") + info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))") + info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))") + info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))") + error("find projection from master to slave: did not converge") + +end + +function project_from_slave_to_master{E<:MortarElements2D}( + master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI; + tol=1.0e-10, max_iterations=20) + + function x2(xi2) + N = get_basis(E, xi2) + return vec(N)*x2_ + end + + function dx2(xi2) + dN = get_dbasis(E, xi2) + return vec(dN)*x2_ + end + + cross2(a, b) = cross([a; 0], [b; 0])[3] + R(xi2) = cross2(x2(xi2)-x1, n1) + dR(xi2) = cross2(dx2(xi2), n1) + xi2 = 0.0 + dxi2 = 0.0 + for i=1:max_iterations + dxi2 = -R(xi2) / dR(xi2) + xi2 += dxi2 + if norm(dxi2) < tol + return xi2 + end + end + + error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2") + +end + + +function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real, + ::Type{Val{:forwarddiff_old}}) + haskey(slave_element, "master elements") || return + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_dofs = get_gdofs(slave_element, field_dim) + nnodes = size(slave_element, 2) + X1 = slave_element("geometry", time) + #u1 = slave_element("displacement", time) + #x1 = X1 + u1 + slave_element_nodes = get_connectivity(slave_element) + adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) + adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements + Q = [0.0 -1.0; 1.0 0.0] + + X = spzeros(10000, 1) + for element in get_elements(problem) + conn = get_connectivity(element) + geom = element("geometry", time) + for (c, g) in zip(conn, geom) + dofs = [field_dim*(c-1)+1, field_dim*(c-1)+2] + X[dofs] = g + end + end + + # here x does not mean deformed configuration + x = [problem.assembly.u; problem.assembly.la] + if length(x) == 0 + info("mortar_2d_autodiff: length(x) == 0") + # resize solution vectors according to initial configuration of this problem + X = vec(full(sparse(findnz(X)...))) + x = zeros(length(X)*2) + else + # resize initial configuration to match real dimension + I, J, V = findnz(X) + X = vec(full(sparse(I, J, V, length(problem.assembly.u), 1))) + end + ndofs = round(Int, length(x)/2) + # at the end we should have +# info("mortar_2d_autodiff: size of x = $(size(x))") +# info("mortar_2d_autodiff: size of X = $(size(X))") +# info("mortar_2d_autodiff: ndofs = $ndofs") + + """ Calculate normal vector for slave element nodes in current configuration. """ + function calculate_normals(u::Matrix) + normals = zeros(u) + # 1. update nodal normals + for element in adjacent_elements + conn = get_connectivity(element) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:, i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + normals[:, conn] += ip.weight*Q*t'*N + end + end + slave_normals = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) + return slave_normals + end + + function calculate_mortar_projection(u::Matrix, n1::DVTI) + B = SparseMatrixCOO{Real}([], [], []) + + u1 = Field([u[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + + for master_element in slave_element["master elements"] + X2 = master_element("geometry", time) + master_element_nodes = get_connectivity(master_element) + u2 = Field([u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + #info("master element coordinate 1 = $(ForwardDiff.get_value(x2[1]))") + #info("master element coordinate 2 = $(ForwardDiff.get_value(x2[2]))") + + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1 = clamp([xi1a; xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution + + # integrate slave side + D = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + dN = get_dbasis(slave_element, ip) + # jacobian of slave element in deformed state + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N = get_basis(slave_element, xi_s) + D += w*diagm(vec(N)) + Me += w*N'*N + end + Ae = D*inv(Me) + + # integrate master side + M = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + dN = get_dbasis(slave_element, ip) + # jacobian of slave element in deformed state + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_g = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_g) + x_g = vec(N1)*x1 + n_g = vec(N1)*n1 + #info("slave gauss point coordinates $(ForwardDiff.get_value(x_g))") + #info("slave gauss point normal direction $(ForwardDiff.get_value(n_g))") + xi_m = project_from_slave_to_master(master_element, x_g, n_g, x2) + N2 = get_basis(master_element, xi_m) + M += w*kron(Ae*N1', N2) + end + + slave_dofs = get_gdofs(slave_element, field_dim) + master_dofs = get_gdofs(master_element, field_dim) + for i=1:field_dim + add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D) + add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M) + end + + end + + return B + end + + function calculate_contact_rhs(x::Vector) + ndofs = round(Int, length(x)/2) + u = x[1:ndofs] + la = x[ndofs+1:end] +# info("calculate_contact_rhs: size of u = $(size(u))") +# info("calculate_contact_rhs: size of X = $(size(X))") +# info("calculate_contact_rhs: size of la = $(size(la))") +# info("calculate_contact_rhs: ndofs = $ndofs") + u2 = reshape(u, field_dim, round(Int, length(u)/field_dim)) + + normals = calculate_normals(u2) + B = calculate_mortar_projection(u2, normals) + B = sparse(B, ndofs, ndofs) + fc = B' * la # contact force residual for r = fint + fc - fext = 0 + + N = SparseMatrixCOO{Real}([], [], []) + T = SparseMatrixCOO{Real}([], [], []) + for (i, j) in enumerate(slave_element_nodes) + dofs = [2*(j-1)+1, 2*(j-1)+2] + add!(N, dofs, [dofs[1]], reshape(normals[i], 2, 1)) + add!(T, dofs, [dofs[2]], reshape(Q'*normals[i], 2, 1)) + end + N = sparse(N, ndofs, ndofs) + T = sparse(T, ndofs, ndofs) + gn = -N*B*(X+u) + gt = -T*B*(X+u) +# gn = -N*B*u + lan = N*la + lat = T*la + cn = 1.0e3 + C = lan - max(0, lan - cn*gn) + lat + # C = lan - gn + lat - gt <-- ihan viturallensa + # C = gn+gt <-- not working + # C = B*(X+u) + # C = B*u <- pitää kiinni, "tie". + # C = N*B*u + T*la <- palikat menee väärään suuntaan + # C = -N*B*(X+u) + T*la <- toimii suht hyvin mut kääntyy väärään suuntaan (t-suunnassa) + # C = -N*B*(X+u) - T*la <- sama + # C = -N*B*(X+u) - T*B*la <- sama + # C = N*la + T*la - max(0, N*la + N*B*(X+u)) + cond = lan[1:field_dim:end] - cn*gn[1:field_dim:end] + all_nodes = slave_element_nodes + inactive_nodes = find(cond .<= 0) + active_nodes = find(cond .> 0) + inactive_nodes = setdiff(all_nodes, inactive_nodes) + active_nodes = setdiff(all_nodes, active_nodes) + info("S = $all_nodes, I = $inactive_nodes, A = $active_nodes") + info("lambda = $(ForwardDiff.get_value(lan[slave_dofs]))") + info("gn = $(ForwardDiff.get_value(gn[slave_dofs]))") + #for j in active_nodes + # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] + # C[dofs] = 0 + #end + +# C = -(N+T)*B*(X+u) +# C = -B*(X+u) + + return [fc; C] + + end + + A, allresults = ForwardDiff.jacobian(calculate_contact_rhs, x, ForwardDiff.AllResults) + b = -ForwardDiff.value(allresults) + A = sparse(A) + b = sparse(b) + K = A[1:ndofs,1:ndofs] + C1 = A[1:ndofs,ndofs+1:end]' + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + C2[2:field_dim:end] = 0 + g[2:field_dim:end] = 0 + +# C2[2:field_dim:end] = 0 +# g[2:field_dim:end] = 0 +# inactives = find(g[1:field_dim:end] .<= 0) +# actives = find(g[1:field_dim:end] .> 0) +# inactives = setdiff(slave_element_nodes, inactives) +# actives = setdiff(slave_element_nodes, actives) +# info("all nodes = $slave_element_nodes, inactives = $inactives, actives = $actives") +# info("g = $g") + +# for j in inactives + # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] + # K[dofs,:] = 0 + # C1[dofs,:] = 0 + # C2[dofs,:] = 0 + # D[dofs,:] = 0 + # g[dofs,:] = 0 + #end + #for j in actives + # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] + # C2[dofs[1],:] = 0 + #end + add!(assembly.K, K) + add!(assembly.C1, C1) + add!(assembly.C2, C2) + add!(assembly.D, D) + add!(assembly.f, f) + add!(assembly.g, g) +end + + +function assemble!{E<:MortarElements2D}(assembly::Assembly, + problem::Problem{Mortar}, slave_element::Element{E}, + time::Real, ::Type{Val{:forwarddiff_old2}}) + haskey(slave_element, "master elements") || return + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + + function calculate_interface(x::Vector) + + ndofs = round(Int, length(x)/2) + nnodes = round(Int, ndofs/field_dim) + u = reshape(x[1:ndofs], field_dim, nnodes) + la = reshape(x[ndofs+1:end], field_dim, nnodes) + fc = zeros(u) + C = zeros(la) + + slave_element_nodes = get_connectivity(slave_element) + X1 = slave_element("geometry", time) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + + # 1. update nodal normals for this element. average nodes from adjacent elements + adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) + adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements + Q = [0.0 -1.0; 1.0 0.0] + normals = zeros(u) + for element in adjacent_elements + conn = get_connectivity(element) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:, i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + normals[:, conn] += ip.weight*Q*t'*N + end + end + # --> slave side normals in deformed state + n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) + + for master_element in slave_element["master elements"] + + master_element_nodes = get_connectivity(master_element) + X2 = master_element("geometry", time) + u2 = Field(Vector[u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1 = clamp([xi1a; xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution in this master element + + nnodes = size(slave_element, 2) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_s) + De += w*diagm(vec(N1)) + Me += w*N1'*N1 + end + Ae = De*inv(Me) + + slave_dofs = get_gdofs(slave_element, field_dim) + master_dofs = get_gdofs(master_element, field_dim) + + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + + # project gauss point from slave element to master element + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = vec(get_basis(slave_element, xi_s)) + x_s = N1*x1 # coordinate in gauss point + n_s = N1*n1 # normal direction in gauss point + t_s = Q'*n_s # tangent direction in gauss point + R_s = [n_s t_s] + xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) + N2 = vec(get_basis(master_element, xi_m)) + x_m = N2*x2 + Phi = Ae*N1 + + u_s = N1*u1 + u_m = N2*u2 + + la_s = Phi*la1 # traction force in gauss point + #lan = dot(n_s, la_s) # normal component + #lat = dot(t_s, la_s) # tangential component + la_nt = R_s*la_s + g = x_s-x_m + gn = props.gap_sign*dot(n_s, g) # normal gap + #gu = dot(n_s, u_s - u_m) # normal displacement gap + #gt = dot(t_s, x_s - x_m) # tangential gap + + fc[:,slave_element_nodes] += w*la_s*N1' + fc[:,master_element_nodes] -= w*la_s*N2' + C[1,slave_element_nodes] += w*gn*Phi' + #C[2,slave_element_nodes] += w*la_nt[2]*Phi' + #C[2,slave_element_nodes] += w*la_nt[2,:]*Phi' + #C[2,slave_element_nodes] += w*dot(t_s, u_s - u_m)*Phi' # <-- for tie + + #R += w*R_s + + end + + end # master elements done + + for (i, j) in enumerate(slave_element_nodes) + n = n1[i] + t = Q'*n + R = [n t] + la_nt = R*la[:,j] + C[2,j] += la_nt[2] + end + + return vec([fc C]) + + end + + # x doesn't mean deformed configuration here + x = [problem.assembly.u; problem.assembly.la] + ndofs = round(Int, length(x)/2) + #if ndofs == 0 + # info("INITIALIZING THINGS") + # problem.assembly.u = zeros(16) + # problem.assembly.la = zeros(16) + # x = [problem.assembly.u; problem.assembly.la] + # ndofs = round(Int, length(x)/2) + #end + + A, allresults = ForwardDiff.jacobian(calculate_interface, x, ForwardDiff.AllResults) + b = -ForwardDiff.value(allresults) + #b = -calculate_interface(x) + #info("PE = $(ForwardDiff.value(allresults))") + A = sparse(A) + b = sparse(b) + SparseMatrix.droptol!(A, 1.0e-12) + SparseMatrix.droptol!(b, 1.0e-12) + #println(A) + K = A[1:ndofs,1:ndofs] + C1 = transpose(A[1:ndofs,ndofs+1:end]) + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + add!(assembly.K, K) + add!(assembly.C1, C1) + add!(assembly.C2, C2) + add!(assembly.D, D) + add!(assembly.f, f) + add!(assembly.g, g) + + return + +end + + +function assemble!{E<:MortarElements2D}(assembly::Assembly, + problem::Problem{Mortar}, slave_element::Element{E}, + time::Real, ::Type{Val{:forwarddiff}}) + haskey(slave_element, "master elements") || return + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + + function calculate_interface(x::Vector) + + ndofs = round(Int, length(x)/2) + nnodes = round(Int, ndofs/field_dim) + u = reshape(x[1:ndofs], field_dim, nnodes) + la = reshape(x[ndofs+1:end], field_dim, nnodes) + fc = zeros(u) + C = zeros(la) + + slave_element_nodes = get_connectivity(slave_element) + X1 = slave_element("geometry", time) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + + # 1. update nodal normals for this element. average nodes from adjacent elements + adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) + adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements + Q = [0.0 -1.0; 1.0 0.0] + normals = zeros(u) + for element in adjacent_elements + conn = get_connectivity(element) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:, i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + normals[:, conn] += ip.weight*Q*t'*N + end + end + # --> slave side normals in deformed state + n1 = Field(Vector[ForwardDiff.get_value(normals[:,i]/norm(normals[:,i])) for i in slave_element_nodes]) + + + nnodes = size(slave_element, 2) + lan_tot = zeros(nnodes) # normal pressure + gap_tot = zeros(nnodes) # weighted normal gap + + for master_element in slave_element["master elements"] + + master_element_nodes = get_connectivity(master_element) + X2 = master_element("geometry", time) + u2 = Field(Vector[u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1 = clamp([xi1a; xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution in this master element + + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_s) + De += w*diagm(vec(N1)) + Me += w*N1'*N1 + end + Ae = De*inv(Me) + + slave_dofs = get_gdofs(slave_element, field_dim) + master_dofs = get_gdofs(master_element, field_dim) + + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + + # project gauss point from slave element to master element + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = vec(get_basis(slave_element, xi_s)) + x_s = N1*x1 # coordinate in gauss point + n_s = N1*n1 # normal direction in gauss point + t_s = Q'*n_s # tangent direction in gauss point + R_s = [n_s t_s] + xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) + N2 = vec(get_basis(master_element, xi_m)) + x_m = N2*x2 + Phi = Ae*N1 + + u_s = N1*u1 + u_m = N2*u2 + + la_s = Phi*la1 # traction force in gauss point + la_nt = R_s*la_s + gn = -dot(n_s, x_s - x_m) # normal gap + + fc[:,slave_element_nodes] += w*la_s*N1' + fc[:,master_element_nodes] -= w*la_s*N2' + #C[1,slave_element_nodes] += w*gn*Phi' + + lan_tot += w*la_nt[1]*Phi + gap_tot += w*gn*Phi + + end + + end # master elements done + +# ncf = lan_tot - max(0, lan_tot - gap_tot) +# info("pressure in nodes: $(ForwardDiff.get_value(lan_tot))") +# info("weighted gap in nodes: $(ForwardDiff.get_value(gap_tot))") +# info("ncf: $(ForwardDiff.get_value(ncf))") +# cond = +lan_tot + gap_tot +# cond = +lan_tot - gap_tot # singular +# cond = -lan_tot + gap_tot +# cond = -lan_tot - gap_tot + + for (i, j) in enumerate(slave_element_nodes) + n = n1[i] + t = Q'*n + R = [n t] + la_nt = R*la[:,j] + info("node $j, n=$(ForwardDiff.get_value(n)) lan = $(ForwardDiff.get_value(la_nt[1])) gap = $(ForwardDiff.get_value(gap_tot[i]))") + +# if -lan_tot[i] + gap_tot[i] < 0 + if -la_nt[1] + gap_tot[i] < 0 + #if j in [31, 32, 33, 34, 35, 36] + info("set node $j active") + C[1,j] -= gap_tot[i] +# C[1,j] += la_nt[1] - max(0, la_nt[1] - gap_tot[i]) + C[2,j] += la_nt[2] + else + info("set node $j inactive") + C[1,j] += la1[i][1] + C[2,j] += la1[i][2] + end + end + + return vec([fc C]) + + end + + # x doesn't mean deformed configuration here + x = [problem.assembly.u; problem.assembly.la] + ndofs = round(Int, length(x)/2) + A, allresults = ForwardDiff.jacobian(calculate_interface, x, ForwardDiff.AllResults) + b = -ForwardDiff.value(allresults) + #b = -calculate_interface(x) + #info("PE = $(ForwardDiff.value(allresults))") + A = sparse(A) + b = sparse(b) + SparseMatrix.droptol!(A, 1.0e-12) + SparseMatrix.droptol!(b, 1.0e-12) + #println(A) + K = A[1:ndofs,1:ndofs] + C1 = transpose(A[1:ndofs,ndofs+1:end]) + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + add!(assembly.K, K) + add!(assembly.C1, C1) + add!(assembly.C2, C2) + add!(assembly.D, D) + add!(assembly.f, f) + add!(assembly.g, g) + + return + +end + diff --git a/src/mortar_3d.jl b/src/mortar_3d.jl new file mode 100644 index 0000000..b4194e4 --- /dev/null +++ b/src/mortar_3d.jl @@ -0,0 +1,643 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# Mortar projection calculation for 3d cases + +""" Construct auxiliary plane for surface. """ +function create_auxiliary_plane{E}(element::Element{E}, time::Real) + xi = get_reference_element_midpoint(E) + x0 = element("geometry", xi, time) + ntbasis = element("normal-tangential coordinates", xi, time) + return x0, ntbasis +end + + +""" +Project point q onto a plane given by a point p and normal n. + +Parameters +---------- +q::Array{Float64, 2} + point to project (row vector) +x0::Array{Float64, 2} + origo of plane +n::Array{Float64, 2} + normal vector of plane + +Returns +------- +y::Array{Float64, 2} + projected point + +Examples +-------- +julia> p = [-0.5 -1.0 4.0]' +julia> x0 = [0.0 0.075 0.675]' +julia> n = [0.1485860 0.0784519 0.9857830]' +julia> project_node_to_auxiliary_plane(p, x0, n) +3-element Array{Float64,1}: + 0.963455 + -1.2447 + 0.925247 + +Notes +----- +[1](http://stackoverflow.com/questions/8942950/how-do-i-find-the-orthogonal-projection-of-a-point-onto-a-plane) + +""" +function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix) + n = Q[:,1] + ph = p - dot(p-x0, n)*n + qproj = Q'*(ph-x0) + if !isapprox(qproj[1], 0.0; atol=1.0e-12) + info("project_point_to_auxiliary_plane(): point not projected correctly.") + info("p: $p") + info("x0: $x0") + info("Q: \n$Q") + info("qproj: $qproj") + error("Failed to project point to auxiliary plane.") + end + return qproj[2:3] +end + +""" +Find edge intersections of two planar arbitrary shape polygons. + +Parameters +---------- +S::Array{Float64,2} +M::Array{Float64,2} + +Matrices with size (2, n) where n is number of vertices of each polygon. + +Returns +------- +P::Array{Float64,2} + Intersection points of polygons +n::Array{Float64,2} + Neighbour info matrix with size (ns, mn). This keeps information which + edges of polygons are intersecting. See further explanation in example + below. + +Examples +-------- +Find intersection points of two triangles: + +julia> S = [0 0; 3 0; 0 3]' +julia> M = [-1 1; 2 -1/2; 1 3/2]' +julia> P, n = get_edge_intersections(S, M) +julia> P +2x4 Array{Float64,2}: + 1.0 1.75 0.0 0.0 + 0.0 0.0 0.5 1.25 +julia> n +3x3 Array{Int64,2}: + 1 1 0 + 0 0 0 + 1 0 1) + +So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25). +"Neighbour matrix" can be interpreted as following: + + 1 1 0 <--> First edge of S intersects edges 1 and 2 of M + 0 0 0 <--> Second edge of S doesn't intersect at all + 1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M + +""" +function get_edge_intersections(S::Matrix, M::Matrix) + ns = size(S, 2) + nm = size(M, 2) + P = zeros(2, 0) + n = zeros(Int64, ns, nm) + k = 0 + for i=1:ns + for j=1:nm + b = M[:,j]-S[:,i] + A = [S[:,mod(i,ns)+1]-S[:,i] -M[:,mod(j,nm)+1]+M[:,j]] + if rank(A) == 2 + r = A\b + if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found + k += 1 + f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i]) + f = f'' + P = hcat(P, f) + n[i, j] = 1 + end + end + end + end + return P, n +end + + +""" +Find any points laying inside or border of triangle. + +Parameters +---------- +Y::Array{Float64, 2} + Triangle coordinates in 2×3 matrix +X::Array{Float64, 2} + List of points to test in 2×n matrix + +Returns +------- +P::Array{Float64, 2} + List of points in triangle in 2×m matrix, where m is number of points inside triangle + +Examples +-------- +julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points +julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests +julia> points_in_triangle(S, pts) +2x2 Array{Float64,2}: + 1.0 0.5 + 1.5 1.5 +""" +function get_points_inside_triangle(Y::Matrix, X::Matrix) + @assert size(Y, 2) == 3 # "Point in TRIANGLE..." + P = zeros(2, 0) + v0 = Y[:,2] - Y[:,1] + v1 = Y[:,3] - Y[:,1] # find interior points of X in Y + d00 = (v0'*v0)[1] + d01 = (v0'*v1)[1] + d11 = (v1'*v1)[1] # using baricentric coordinates + id = 1/(d00*d11 - d01*d01) + for i=1:size(X, 2) + v2 = X[:,i] - Y[:,1] + d02 = (v0'*v2)[1] + d12 = (v1'*v2)[1] + u = (d11*d02-d01*d12)*id + v = (d00*d12-d01*d02)*id + if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary + P = hcat(P, X[:,i]'') + end + end + return P +end + +""" +Determine is point P inside or on boudary of polygon X. + +http://paulbourke.net/geometry/polygonmesh/#insidepoly +""" +function is_point_inside_convex_polygon(P, X) + x, y = P + for i=1:length(X) + x0, y0 = X[i] + x1, y1 = X[mod(i, length(X))+1] + if (y-y0)*(x1-x0) - (x-x0)*(y1-y0) < 0 + return false + end + end + return true +end + +function get_points_inside_convex_polygon(pts, X) + # TODO: Make more readable + X2 = [X[:,i] for i=1:size(X,2)] + c = filter(P->is_point_inside_convex_polygon(P, X2), [pts[:,i] for i=1:size(pts, 2)]) + return length(c) == 0 ? zeros(2, 0) : hcat(c...) +end + +""" Return unique objects with some given tolerance. This is used in next function + because traditional unique() command returns row vectors as non-unique if they + differs only a "little". +""" +function uniquetol(P, dim::Int; args...) + @assert dim == 2 + items = Vector{Float64}[P[:,i] for i=1:size(P,dim)] + new_items = Vector{Float64}[] + for item in items + has_found = false + for new_item in new_items + if isapprox(item, new_item; args...) + has_found = true + break + end + end + if !has_found + push!(new_items, item) + end + end + return reshape([new_items...;], length(new_items[]), length(new_items)) +end + + +""" +Make polygon clipping of shapes S and M. + +Parameters +---------- +S::Array{Float64, 2} +M::Array{Float64, 2} + Shapes to clip. Needs to be triangles at the moment. + +Returns +------- +Array{Float64, 2}, Array{Float64, 2} +- Polygon vertices in 2×n matrix, sorted in counter-clockwise order. +- 3×3 "neighbouring" matrix, see example. + +Examples +-------- +julia> S = [0 0; 3 0; 0 3]' +julia> M = [-1 1; 2 -1/2; 2 2]' +julia> P, n = clip_polygon(S, M) +julia> P +2x6 Array{Float64,2}: + 0.0 1.0 2.0 2.0 1.25 0.0 + 0.5 0.0 0.0 1.0 1.75 1.33333, +julia> n +3x3 Array{Int64,2}: + 1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]') + 1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S + 0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S + +""" +function clip_polygon(S::Matrix, M::Matrix) + P1, neighbours = get_edge_intersections(M, S) + #P2 = get_points_inside_triangle(M, S) + #P3 = get_points_inside_triangle(S, M) + P2 = get_points_inside_convex_polygon(M, S) + P3 = get_points_inside_convex_polygon(S, M) +# info("polygon clipping: P1 = $P1") +# info("polygon clipping: P2 = $P2") +# info("polygon clipping: P3 = $P3") +# info("hcat P = $P") + P = hcat(P1, P2, P3) + if length(P) == 0 + return nothing, nothing + end + P = uniquetol(P, 2) + meanval = mean(P, 2) + tmp = P .- meanval + angles = atan2(tmp[2,:], tmp[1,:]) + angles = reshape(angles, length(angles)) + order = sortperm(angles) + return P[:, order], neighbours +end + + +""" +Calculate polygon geometric center point + +Parameters +---------- +P::Array{Float64, 2} + Polygon vertices in 2×n matrix + +Returns +------- +Array{Float63, 2} + Center point + +Examples +-------- +julia> P +2x6 Array{Float64,2}: + 0.0 1.0 2.0 2.0 1.25 0.0 + 0.5 0.0 0.0 1.0 1.75 1.33333, +julia> C = get_polygon_cp(P) +2x1 Array{Float64,2}: + 1.039740 + 0.804701 + +""" +function calculate_polygon_centerpoint(P::Matrix) + n = size(P, 2) + A = 0.0 + for i=1:n + A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i]) + end + Cx = 0.0 + Cy = 0.0 + for i=1:n + inext = mod(i, n)+1 + Cx += 1/(6*A)*(P[1,i] + P[1,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) + Cy += 1/(6*A)*(P[2,i] + P[2,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) + end + return Float64[Cx, Cy] +end + +""" +Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂) + +Parameters +---------- +p::Array{Float64,1} + point in auxiliary plane, in (n,t1,t2) coordinate system +x0::Array{Float64,1} + origo of auxiliary plane cs +Q::Array{Float64,2} + basis of auxiliary plane cs +x::Array{Float64,2} + surface node coords +basis::Array{Float64,2} + surface basis functions +dbasis::Array{Float64,2} + partial derivatives of surface basis functions + +Returns +------- +Array{Float64,2} + solution vector (d, ξ₁, ξ₂) where d is distance to surface + +Examples +-------- +Define surface with node points, basis + dbasis + +julia> xquad = [ +... -2.5 -2.0 1.0 +... 2.5 -2.0 0.7 +... 2.0 2.3 0.0 +... -2.0 2.0 1.0]' +julia> basis(xi) = [ +... (1-xi[1])(1-xi[2])/4 +... (1+xi[1])(1-xi[2])/4 +... (1+xi[1])(1+xi[2])/4 +... (1-xi[1])(1+xi[2])/4] +julia> dbasis(xi) = [ +... -(1-xi[2])/4 -(1-xi[1])/4 +... (1-xi[2])/4 -(1+xi[1])/4 +... (1+xi[2])/4 (1+xi[1])/4 +... -(1+xi[2])/4 (1-xi[1])/4] + +We aim to find point p, which we first project to auxiliary plane defined as following +julia> p = [-2.5 -2.0 1.0]' +julia> x0 = [0.0 0.075 0.675]' +julia> Q = [ +... 0.1485860 0.9888990 0.0000000 +... 0.0784519 -0.0117877 0.9968480 +... 0.9857830 -0.1481180 -0.0793325] + +Our projected point is therefore +julia> n = Q[:,1] # first component is normal direction +julia> ph = project_node_to_auxiliary_plane(p, x0, n) +julia> ph = Q'(ph-x0) +julia> ph +3x1 Array{Float64,2}: + 1.33264e-7 + -2.49593 + -2.09424 + +Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we +project it back to surface defined by xquad*basis + +julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis) +julia> theta +3x1 Array{Float64,2}: + -0.213874 + -0.999999 + -1.0 + +We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad +[-2.5 -2.0 1.0]' correctly. + +julia> xquad*basis(theta[2:3]) +3-element Array{Float64,1}: + -2.5 + -2.0 + 1.0 + +""" +function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) + basis(xi) = get_basis(E, xi) + dbasis(xi) = get_dbasis(E, xi) + x = element("geometry", time) + ph = Q*[0; p] + x0 + theta = Float64[0.0, 0.0, 0.0] + n = Q[:,1] + for i=1:max_iterations + b = ph + theta[1]*n - basis(theta[2:3])*x + J = [n -dbasis(theta[2:3])*x] + dtheta = J \ -b + theta += dtheta + if norm(dtheta) < iter_tol + return theta + end + end + begin + info("projecting point from auxiliary plane back to surface didn't go very well.") + info("element type: $E") + info("element connectivity: $(get_connectivity(element))") + info("auxiliary plane: x0 = $x0, Q = $Q") + info("point coordinates on plane: $p") + info("element geometry: $x") + info("ph: $ph") + info("normal direction: $n") + info("parameter vector before giving up: $theta") + end + error("project_point_to_surface: did not converge in $max_iterations iterations!") +end + +typealias MortarElements3D Union{Tri3, Quad4} + +function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar}, + slave_element::Element{E}, time::Real, ::Type{Val{:total}}) + assemble!(assembly, problem, slave_element, time, Val{problem.properties.formulation}) +end + +function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar}, + slave_element::Element{E}, time::Real, ::Type{Val{:total}}) + haskey(slave_element, "master elements") || return + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_dofs = get_gdofs(slave_element, field_dim) + + props = problem.properties + if props.formulation == :Standard && props.normal_condition == :Contact + error("for contact choose Dual formulation.""") + end + + # create auxiliary plane and project slave nodes to it + # x0 = origo, Q = local basis + x0, Q = create_auxiliary_plane(slave_element, time) + + # 1. project slave nodes to auxiliary plane + Sl = Vector{Float64}[] + for p in slave_element("geometry", time) + push!(Sl, project_point_to_auxiliary_plane(p, x0, Q)) + end + S = hcat(Sl...) + + for master_element in slave_element["master elements"] + + # if distance between elements is "far enough" cannot expect contact + if (props.normal_condition == :Contact) || props.inequality_constraints + slave_midpoint = slave_element("geometry", [0.0, 0.0], time) + master_midpoint = master_element("geometry", [0.0, 0.0], time) + if norm(slave_midpoint - master_midpoint) > props.minimum_distance + continue + end + end + + master_dofs = get_gdofs(master_element, field_dim) + + # 2. project master nodes to auxiliary plane + M = Vector{Float64}[] + for p in master_element("geometry", time) + push!(M, project_point_to_auxiliary_plane(p, x0, Q)) + end + M = hcat(M...) + + # 3. create polygon clipping on auxiliary plane + P = nothing + neighbours = nothing + try + P, neighbours = clip_polygon(S, M) + catch + info("polygon clipping failed") + info("S = ") + dump(S) + info("M = ") + dump(M) + info("original Sl = ") + info(Sl) + error("cannot continue") + end + isa(P, Void) && continue # no clipping + + # shared edge but no shared volume. skipping + size(P, 2) < 3 && continue + + C = calculate_polygon_centerpoint(P) + npts = size(P, 2) # number of vertices in polygon + + # loop vertices and create temporary integrate cells + # TODO: basically when npts == 3 or npts == 4 we could integrate without splitting to cells. + nnodes = size(slave_element, 2) + C1S3 = zeros(3*nnodes, 3*nnodes) + C1M3 = zeros(3*nnodes, 3*nnodes) + + for pnt=1:npts # integration of mortar matrices begin + cell = Field(Vector{Float64}[C, P[:,pnt], P[:,mod(pnt,npts)+1]]) + + # calculate slave side projection matrix D + # construct dual basis + Ae = zeros(nnodes, nnodes) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + if problem.properties.formulation == :Dual # Construct dual basis + for ip in get_integration_points(Tri3, Val{5}) + N = get_basis(Tri3, ip.xi) + xi = vec(N*cell) + theta = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time) + xi_slave = theta[2:3] + N1 = slave_element(xi_slave, time) + # jacobian determinant on integration cell + dNC = get_dbasis(Tri3, ip.xi) + JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)]) + wC = ip.weight*det(JC) + De += wC*diagm(vec(N1)) + Me += wC*N1'*N1 + end + Ae = De*inv(Me) + end + for i=1:field_dim + C1S3[i:field_dim:end,i:field_dim:end] += De + end + + # Calculate master side projection matrix M + for ip in get_integration_points(Tri3, Val{5}) + # gauss point in auxiliary plane + #N = get_basis(E, ip.xi) + N = get_basis(Tri3, ip.xi) + xi = vec(N*cell) # xi defined in auxilary plane + + # find projection of gauss point to master and slave elements + theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time) + theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time) + xi_slave = theta1[2:3] + xi_master = theta2[2:3] + + # evaluate shape functions values in gauss point and add contribution to matrices + N1 = slave_element(xi_slave, time) + N2 = master_element(xi_master, time) + + # jacobian determinant on integration cell + dNC = get_dbasis(Tri3, ip.xi) + JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)]) + wC = ip.weight*det(JC) + + # extend matrices according to the problem dimension (3) + @assert length(slave_dofs) == length(master_dofs) + Me = wC*Ae*N1'*N2 + for k=1:field_dim + C1M3[k:field_dim:end,k:field_dim:end] += Me + end + end + end # integration of mortar matrices done. + + # constraints in normal-tangential direction and initial weighted gap + X1 = vec(slave_element("geometry", time)) + X2 = vec(master_element("geometry", time)) + Q_ = slave_element("normal-tangential coordinates", time) + Z = zeros(3, 3) + if nnodes == 3 + Q3 = [Q Z Z; Z Q Z; Z Z Q] + elseif nnodes == 4 + Q3 = [Q Z Z Z; Z Q Z Z; Z Z Q Z; Z Z Z Q] + end + D3 = zeros(3*nnodes, 3*nnodes) + C2S3 = Q3'*C1S3 + C2M3 = Q3'*C1M3 + G = -(C2S3*X1 - C2M3*X2) + + # complementarity condition + if haskey(slave_element, "displacement") + u1 = vec(slave_element("displacement", time)) + else + u1 = zeros(3*nnodes) + end + if haskey(master_element, "displacement") + u2 = vec(master_element("displacement", time)) + else + u2 = zeros(3*nnodes) + end + x1 = X1 + u1 + x2 = X2 + u2 + if haskey(slave_element, "reaction force") + la = vec(slave_element("reaction force", time)) + else + la = zeros(3*nnodes) + end + g = -(C2S3*x1 - C2M3*x2) + c = Q3'*la - g + inactive_nodes = find(c[1:field_dim:end] .<= 0) + active_nodes = find(c[1:field_dim:end] .> 0) + + # normal constraint: remove inactive nodes if normal condition is set to contact + if problem.properties.normal_condition == :Contact + for j in inactive_nodes + dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + G[dofs] = 0 + C1S3[dofs,:] = 0 + C1M3[dofs,:] = 0 + C2S3[dofs,:] = 0 + C2M3[dofs,:] = 0 + end + end + + # tangential constraint: stick or slip + if problem.properties.tangential_condition == :Slip + D3 = copy(C2S3) + D3[1:field_dim:end, :] = 0 + C2S3[2:field_dim:end, :] = 0 + C2M3[2:field_dim:end, :] = 0 + C2S3[3:field_dim:end, :] = 0 + C2M3[3:field_dim:end, :] = 0 + end + + # add contributions + add!(assembly.C1, slave_dofs, slave_dofs, C1S3) + add!(assembly.C1, slave_dofs, master_dofs, -C1M3) + add!(assembly.C2, slave_dofs, slave_dofs, C2S3) + add!(assembly.C2, slave_dofs, master_dofs, -C2M3) + add!(assembly.D, slave_dofs, slave_dofs, D3) + add!(assembly.c, slave_dofs, c) + add!(assembly.g, slave_dofs, G) + end +end + + diff --git a/src/mortar_forwarddiff_pe.jl b/src/mortar_forwarddiff_pe.jl new file mode 100644 index 0000000..f0117a8 --- /dev/null +++ b/src/mortar_forwarddiff_pe.jl @@ -0,0 +1,446 @@ +using JuliaFEM.Core: MortarElements2D, DVTI, Assembly + +import JuliaFEM.Core: project_from_master_to_slave, project_from_slave_to_master, assemble!, +get_unknown_field_dimension, get_parent_field_name, get_gdofs, find_elements, get_nodes, Field, +get_integration_points, get_basis, get_dbasis, add! + +""" Find segment from slave element corresponding to master element nodes. +x1_, n1_ +slave element geometry and normal direction + +x2_ master element nodes to project onto slave +""" +function project_from_master_to_slave{E<:MortarElements2D}( + slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector) + + function x1(xi1) + N = get_basis(E, xi1) + return vec(N)*x1_ + end + + function dx1(xi1) + dN = get_dbasis(E, xi1) + return vec(dN)*x1_ + end + + function n1(xi1) + N = get_basis(E, xi1) + return vec(N)*n1_ + end + + function dn1(xi1) + dN = get_dbasis(E, xi1) + return vec(dN)*n1_ + end + + cross2(a, b) = cross([a; 0], [b; 0])[3] + R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) + dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) + xi1 = 0.0 + for i=1:5 + dxi1 = -R(xi1)/dR(xi1) + xi1 += dxi1 + if norm(dxi1) < 1.0e-10 + return xi1 + end + end + + error("find projection from master to slave: did not converge") + +end + +function project_from_slave_to_master{E<:MortarElements2D}( + master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI) + + function x2(xi2) + N = get_basis(E, xi2) + return vec(N)*x2_ + end + + function dx2(xi2) + dN = get_dbasis(E, xi2) + return vec(dN)*x2_ + end + + cross2(a, b) = cross([a; 0], [b; 0])[3] + R(xi2) = cross2(x2(xi2)-x1, n1) + dR(xi2) = cross2(dx2(xi2), n1) + xi2 = 0.0 + dxi2 = 0.0 + for i=1:5 + dxi2 = -R(xi2) / dR(xi2) + xi2 += dxi2 + if norm(dxi2) < 1.0e-10 + return xi2 + end + end + + error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2") + +end + + +function assemble!{E<:MortarElements2D}(assembly::Assembly, + problem::Problem{Mortar}, slave_element::Element{E}, + time::Real, ::Type{Val{:forwarddiff}}) + haskey(slave_element, "master elements") || return + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + + function calculate_interface(u::Matrix, la::Matrix) + + X1 = slave_element("geometry", time) + slave_element_nodes = get_connectivity(slave_element) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + + adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) + adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements + Q = [0.0 -1.0; 1.0 0.0] + # 1. update nodal normals for this element + normals = zeros(u) + for element in adjacent_elements + conn = get_connectivity(element) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:, i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + normals[:, conn] += ip.weight*Q*t'*N + end + end + # --> slave side normals in deformed state + n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) + + fc = SparseMatrixCOO{Real}([], [], []) # interface virtual work + C = SparseMatrixCOO{Real}([], [], []) # constraints + B = SparseMatrixCOO{Real}([], [], []) + + #info("u1.data = ", ForwardDiff.get_value(u1.data)) + info("normal calculations done. looping master elements.") + for master_element in slave_element["master elements"] + X2 = master_element("geometry", time) + master_element_nodes = get_connectivity(master_element) + u2 = Field(Vector[u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + info("master element ready.") + + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + info("calculating segmentation.") + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1 = clamp([xi1a; xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution + info("xi1 = $xi1") + + info("create bi-orthogonal basis") + nnodes = size(slave_element, 2) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_s) + De += w*diagm(vec(N1)) + Me += w*N1'*N1 + end + Ae = De*inv(Me) + info("bi-orthogonal basis done. integrating fc.") + + slave_dofs = get_gdofs(slave_element, field_dim) + master_dofs = get_gdofs(master_element, field_dim) + + info("integrate fc") + D = zeros(nnodes, nnodes) + M = zeros(nnodes, nnodes) + gn = zeros(nnodes) + lan = zeros(nnodes) + lat = zeros(nnodes) + + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_s) + # project gauss point to master element to evaluate shape function there + x_s = vec(N1)*x1 # coordinate in gauss point + n_s = vec(N1)*n1 # normal direction in gauss point + xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) + N2 = get_basis(master_element, xi_m) + x_m = vec(N2)*x2 + Phi = vec(Ae*N1') + la_s = Phi*la1 # traction force in gauss point + u_s = vec(N1)*u1 + u_m = vec(N2)*u2 + #info("la_s = $(ForwardDiff.get_value(la_s))") + SM = [slave_dofs; master_dofs] + #N1N2 = [N1 -N2] + #info("all_dofs = $(SM)") + #info("shape functions = $(ForwardDiff.get_value(N1N2))") + #for i=1:field_dim + # #info("add to slave dofs $(slave_dofs[i:field_dim:end])") + # #info("add to master dofs $(master_dofs[i:field_dim:end])") + # add!(fc, slave_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*N1) + # add!(fc, master_dofs[i:field_dim:end], [1, 1], +w*la_s[i]*N2) + #add!(fc, slave_dofs[i:field_dim:end], [1, 1], w*la_s'*u_s[i]) + #add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s'*u_m[i]) + # add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*u_m) + #end + #add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1, 1, 1, 1, 1], w*la_s*[u_s' -u_m']) + D += w*kron(Ae*N1', N1) + M += w*kron(Ae*N1', N2) + gn += -w*dot(n_s, x_s-x_m)*Phi + lan += w*dot(n_s, la_s)*Phi + t_s = Q'*n_s + lat += w*dot(t_s, la_s)*Phi + end + + #D2 = zeros(2*nnodes, 2*nnodes) + #M2 = zeros(2*nnodes, 2*nnodes) + #for i=1:field_dim + # D2[i:field_dim:end, i:field_dim:end] += D + # M2[i:field_dim:end, i:field_dim:end] += M + #end + #info("size of D2 = $(size(D2))") + #fco = [D2 -M2]*vec(la1) + #fco = [D -M]*la[:,slave_element_nodes] + #info("fco = $(ForwardDiff.get_value(fco))") + #add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1], fco) + + for i=1:field_dim + add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D) + add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M) + end + info("gn = $gn") + #Cj = lan - max(0, lan - gn) + lat + add!(C, slave_dofs[1:field_dim:end], [1, 1], gn') + + end # master elements done + + ndofs = prod(size(la)) + N = SparseMatrixCOO{Real}([], [], []) + T = SparseMatrixCOO{Real}([], [], []) + for (i, j) in enumerate(slave_element_nodes) + dofs = [2*(j-1)+1, 2*(j-1)+2] + add!(N, [dofs[1]], dofs, reshape(n1[i], 1, 2)) + add!(T, [dofs[2]], dofs, reshape(Q'*n1[i], 1, 2)) + end + N = sparse(N, ndofs, ndofs) + T = sparse(T, ndofs, ndofs) + B = sparse(B, ndofs, ndofs) + fc = B'*vec(la) + #println(sparse(fc)) + #fc = sparse(fc, ndofs, 1) + #println(fc) + #dump(full(fc)) + #C = sparse(C, ndofs, 1) + C = N*B*vec(u) + T*vec(la) + return fc, C + + end + + + function calculate_interface_PE(x::Vector) + + ndofs = round(Int, length(x)/2) + nnodes = round(Int, ndofs/field_dim) + u = reshape(x[1:ndofs], field_dim, nnodes) + la = reshape(x[ndofs+1:end], field_dim, nnodes) + #fixed_la = ForwardDiff.get_value(la) + #fixed_u = ForwardDiff.get_value(u) + #u = ForwardDiff.get_value(u) + + X1 = slave_element("geometry", time) + slave_element_nodes = get_connectivity(slave_element) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + #fixed_la1 = Field(Vector[fixed_la[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + + # 1. update nodal normals for this element + adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) + adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements + Q = [0.0 -1.0; 1.0 0.0] + normals = zeros(u) + for element in adjacent_elements + conn = get_connectivity(element) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:, i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + normals[:, conn] += ip.weight*Q*t'*N + end + end + # --> slave side normals in deformed state + n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) + + Wco = 0.0 + Wla = 0.0 + + for master_element in slave_element["master elements"] + X2 = master_element("geometry", time) + master_element_nodes = get_connectivity(master_element) + u2 = Field(Vector[u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1 = clamp([xi1a; xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution + + nnodes = size(slave_element, 2) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = get_basis(slave_element, xi_s) + De += w*diagm(vec(N1)) + Me += w*N1'*N1 + end + Ae = De*inv(Me) + + for ip in get_integration_points(slave_element, Val{5}) + # jacobian of slave element in deformed state + dN = get_dbasis(slave_element, ip) + j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) + w = ip.weight*norm(j)*l + xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) + N1 = vec(get_basis(slave_element, xi_s)) + + # project gauss point to master element to evaluate shape function there + x_s = N1*x1 # coordinate in gauss point + n_s = N1*n1 # normal direction in gauss point + t_s = Q'*n_s + xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) + N2 = vec(get_basis(master_element, xi_m)) + x_m = N2*x2 + Phi = Ae*N1 + + gn = -dot(n_s, x_s - x_m) + gt = dot(t_s, x_s - x_m) + lan = dot(n_s, Phi*la1) + lat = dot(t_s, Phi*la1) + u_s = N1*u1 + u_m = N2*u2 + gu = dot(n_s, u_s - u_m) + + Wco += w*dot(Phi*la1, N1*u1 - N2*u2) + #Wla += w*(lan*gn + lat*gt) + #gn = min(0, gn) + Wla += 1/2*w*1e6*gn*gn + #info("gn = $(ForwardDiff.get_value(gn))") + #Wla += w*dot(dot(n_s, Phi*la1), dot(n_s, N1*u1 - N2*u2)) + end + + end + return Wco, Wla + end + + + function calculate_contact_rhs(x::Vector) + ndofs = round(Int, length(x)/2) + nnodes = round(Int, ndofs/field_dim) + u = reshape(x[1:ndofs], field_dim, nnodes) + la = reshape(x[ndofs+1:end], field_dim, nnodes) + fc, C = calculate_interface(u, la) + info("interface vector calculated.") + return vec(full([fc; C])) + end + + # x doesn't mean deformed configuration here + x = [problem.assembly.u; problem.assembly.la] + ndofs = round(Int, length(x)/2) + if ndofs == 0 + info("INITIALIZING THINGS") + problem.assembly.u = zeros(16) + problem.assembly.la = zeros(16) + x = [problem.assembly.u; problem.assembly.la] + ndofs = round(Int, length(x)/2) + end + + function add_fco!() + get_PI(x::Vector) = calculate_interface_PE(x)[1] + A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults) + b = -ForwardDiff.gradient(allresults) + info("PE = $(ForwardDiff.value(allresults))") + A = sparse(A) + b = sparse(b) + SparseMatrix.droptol!(A, 1.0e-12) + SparseMatrix.droptol!(b, 1.0e-12) + K = A[1:ndofs,1:ndofs] + C1 = transpose(A[1:ndofs,ndofs+1:end]) + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + add!(assembly.K, K) + add!(assembly.C1, C1) + add!(assembly.C2, C2) + add!(assembly.D, D) + add!(assembly.f, f) + add!(assembly.g, g) + end + #add_fco!() + + function add_wla!() + get_PI(x::Vector) = calculate_interface_PE(x)[2] + A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults) + b = -ForwardDiff.gradient(allresults) + info("PE = $(ForwardDiff.value(allresults))") + A = sparse(A) + b = sparse(b) + SparseMatrix.droptol!(A, 1.0e-12) + SparseMatrix.droptol!(b, 1.0e-12) + #info("A") + #println(full(A)) + K = A[1:ndofs,1:ndofs] + C1 = transpose(A[1:ndofs,ndofs+1:end]) + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + + add!(assembly.K, K) + add!(assembly.C1, C1) + add!(assembly.C2, C2) + add!(assembly.D, D) + add!(assembly.f, f) + add!(assembly.g, g) + end + add_wla!() + + return + +end + +mesh, body1, body2, bc_top, bc_bottom, contact = divided_block_problem() +bc_top.properties.formulation = :incremental +bc_bottom.properties.formulation = :incremental +contact.properties.formulation = :forwarddiff +contact.assembly.u = zeros(16) +contact.assembly.la = zeros(16) +assemble!(contact.assembly, contact, contact.elements[1], 0.0, Val{:forwarddiff}) + diff --git a/src/problems.jl b/src/problems.jl index 2cd5371..3097436 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -166,16 +166,22 @@ function update_assembly!(problem, u, la) # copy current solutions to previous ones and add/replace new solution assembly.u_prev = copy(assembly.u) assembly.la_prev = copy(assembly.la) - if get_formulation_type(problem) == :incremental - info("$(problem.name): incremental formulation, adding increment to solution vector") - #info("solution vector:") - #dump(round(u, 3)') - assembly.u += u - else + if get_formulation_type(problem) == :total info("$(problem.name): total formulation, replacing solution vector with new values") assembly.u = u + assembly.la = la + elseif get_formulation_type(problem) == :incremental + info("$(problem.name): incremental formulation, adding increment to solution vector") + assembly.u += u + assembly.la = la + elseif get_formulation_type(problem) == :forwarddiff + info("$(problem.name): forwarddiff formulation, adding increment to solution vector") + assembly.u += u + assembly.la += la + else + info("$(problem.name): unknown formulation type, don't know what to do with results") + error("serious failure with problem formulation: $(get_formulation_type(problem))") end - assembly.la = la # calculate change of norm assembly.u_norm_change = norm(assembly.u - assembly.u_prev) diff --git a/src/solvers.jl b/src/solvers.jl index 22aa4af..31a172a 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -141,9 +141,11 @@ type Solver name :: ASCIIString # some descriptive name for problem time :: Real # current time iteration :: Int # iteration counter + norms :: Vector{Tuple} # solution norms for convergence studies ndofs :: Int # total dimension of global stiffness matrix, i.e., dim*nnodes problems :: Vector{Problem} is_linear_system :: Bool # setting this to true makes assumption of one step convergence + nonlinear_system_min_iterations :: Int64 nonlinear_system_max_iterations :: Int64 nonlinear_system_convergence_tolerance :: Float64 nonlinear_system_error_if_no_convergence :: Bool @@ -155,9 +157,11 @@ function Solver(name::ASCIIString="default solver", time::Real=0.0) name, time, 0, # iteration # + [], # solution norms in (norm(u), norm(la)) tuples 0, # ndofs [], # array of problems false, # is_linear_system + 1, # min nonlinear iterations 10, # max nonlinear iterations 5.0e-5, # nonlinear iteration convergence tolerance true, # throw error if no convergence @@ -277,12 +281,14 @@ crosspoints. function get_boundary_assembly(solver::Solver) ndofs = solver.ndofs @assert ndofs != 0 + Kc = spzeros(ndofs, ndofs) C1 = spzeros(ndofs, ndofs) C2 = spzeros(ndofs, ndofs) D = spzeros(ndofs, ndofs) g = spzeros(ndofs, 1) for problem in get_boundary_problems(solver) assembly = problem.assembly + Kc_ = sparse(assembly.K, ndofs, ndofs) C1_ = sparse(assembly.C1, ndofs, ndofs) C2_ = sparse(assembly.C2, ndofs, ndofs) D_ = sparse(assembly.D, ndofs, ndofs) @@ -297,12 +303,13 @@ function get_boundary_assembly(solver::Solver) handle_overconstraint_error!(problem, overconstrained_nodes, overconstrained_dofs, C1, C1_, C2, C2_, D, D_, g, g_) end + Kc += Kc_ C1 += C1_ C2 += C2_ D += D_ g += g_ end - return C1, C2, D, g + return Kc, C1, C2, D, g end @@ -316,10 +323,10 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}}) K, f = get_field_assembly(solver) # assemble boundary problems - C1, C2, D, g = get_boundary_assembly(solver) + Kc, C1, C2, D, g = get_boundary_assembly(solver) # construct global system Ax=b and solve using lu factorization - A = [K C1'; C2 D] + A = [K+Kc C1'; C2 D] b = [f; g] nz = get_nonzero_rows(A) @@ -379,6 +386,7 @@ end """ Main solver loop. """ function call(solver::Solver) + # 1. initialize each problem so that we can start nonlinear iterations for problem in solver.problems initialize!(problem, solver.time) @@ -394,6 +402,7 @@ function call(solver::Solver) # 2.2 call solver for linearized system (default: direct lu factorization) u, la = solve_linear_system(solver, Val{solver.linear_system_solver}) + push!(solver.norms, (norm(u), norm(la))) # 2.3 update solution back to elements for problem in solver.problems @@ -404,7 +413,11 @@ function call(solver::Solver) # 2.4 check convergence if has_converged(solver) info("Converged in $(solver.iteration) iterations.") - return true + if solver.iteration < solver.nonlinear_system_min_iterations + info("Converged but continuing") + else + return true + end end end diff --git a/src/sparse.jl b/src/sparse.jl index ac9527f..37250d3 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -4,18 +4,22 @@ # Sparse utils to make assembly of local and global matrices easier. # Unoptimized but should do all necessary stuff for at start. -type SparseMatrixCOO +type SparseMatrixCOO{T<:Real} I :: Vector{Int} J :: Vector{Int} - V :: Vector{Float64} + V :: Vector{T} end typealias SparseMatrixIJV SparseMatrixCOO function SparseMatrixCOO() - SparseMatrixCOO([], [], []) + SparseMatrixCOO{Float64}([], [], []) end +#function SparseMatrixCOO{T}() +# SparseMatrixCOO{T}([], [], []) +#end + function Base.convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC) return SparseMatrixCOO(findnz(A)...) end @@ -80,7 +84,7 @@ function Base.(:+)(A::SparseMatrixIJV, B::SparseMatrixIJV) return C end -function Base.full(A::SparseMatrixIJV, args...) +function Base.full(A::SparseMatrixCOO, args...) return full(sparse(A.I, A.J, A.V, args...)) end @@ -94,7 +98,7 @@ Example >>> S = [3, 4] >>> M = [6, 7, 8] >>> data = Float64[5 6 7; 8 9 10] ->>> A = SparseMatrixIJV() +>>> A = SparseMatrixCOO() >>> add!(A, S, M, data) >>> full(A) 4x8 Array{Float64,2}: @@ -104,7 +108,7 @@ Example 0.0 0.0 0.0 0.0 0.0 8.0 9.0 10.0 """ -function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix{Float64}) +function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix) n, m = size(data) for j=1:m for i=1:n @@ -112,11 +116,16 @@ function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data:: push!(A.J, dofs2[j]) end end -# append!(A.I, repeat(dofs1, outer=[m])) -# append!(A.J, repeat(dofs2, inner=[n])) append!(A.V, vec(data)) end +""" Add sparse matrix of CSC to COO. """ +function add!(A::SparseMatrixCOO, B::SparseMatrixCSC) + I, J, V = findnz(B) + C = SparseMatrixCOO(I, J, V) + append!(A, C) +end + """ Add new data to COO Sparse vector. """ function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1) if length(dofs) != length(data) diff --git a/src/xdmf.jl b/src/xdmf.jl index b6db5f2..578b9c9 100644 --- a/src/xdmf.jl +++ b/src/xdmf.jl @@ -73,11 +73,11 @@ function xdmf_new_grid(temporal_collection; time=0) return grid end -function xdmf_new_mesh!(grid, nodes, elements) +function xdmf_new_mesh!(grid, nodes, elements; datatype="XYZ") # 1. write nodes geometry = new_child(grid, "Geometry") - set_attribute(geometry, "Type", "XYZ") + set_attribute(geometry, "Type", datatype) dataitem = new_child(geometry, "DataItem") set_attribute(dataitem, "DataType", "Float") ndim = sum([length(node) for node in nodes])