diff --git a/REQUIRE b/REQUIRE index 08eff8e..aa17f2b 100644 --- a/REQUIRE +++ b/REQUIRE @@ -13,3 +13,4 @@ FEMQuad Reexport HeatTransfer MortarContact2D +MortarContact2DAD diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 1fc5a14..c7d3aa7 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -47,10 +47,10 @@ export assemble!, postprocess! ### Mortar methods ### @reexport using MortarContact2D +@reexport using MortarContact2DAD include("problems_mortar.jl") include("problems_mortar_3d.jl") -include("problems_mortar_2d_autodiff.jl") export calculate_normals, calculate_normals!, project_from_slave_to_master, project_from_master_to_slave, Mortar, get_slave_elements, get_polygon_clip @@ -67,7 +67,7 @@ include("solvers_modal.jl") export Modal include("problems_contact.jl") include("problems_contact_3d.jl") -include("problems_contact_2d_autodiff.jl") +#include("problems_contact_3d_autodiff.jl") export Contact # Preprocess module diff --git a/src/problems_contact_2d_autodiff.jl b/src/problems_contact_2d_autodiff.jl deleted file mode 100644 index 4466304..0000000 --- a/src/problems_contact_2d_autodiff.jl +++ /dev/null @@ -1,380 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using ForwardDiff - -""" Find segment from slave element corresponding to master element nodes. - -Parameters ----------- -x1_, n1_ - slave element geometry and normal direction -x2 - master element node to project onto slave - -Returns -------- -xi - dimensionless coordinate on slave corresponding to - projected master - -""" -function project_from_master_to_slave_ad{E<:MortarElements2D}( - slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector; - tol=1.0e-10, max_iterations=20, debug=false) - - """ Multiply basis / dbasis at `xi` with field. """ - function mul(func, xi, field) - B = func(slave_element, [xi], time) - return sum(B[i]*field[i] for i=1:length(B)) - end - - x1(xi1) = mul(get_basis, xi1, x1_) - dx1(xi1) = mul(get_dbasis, xi1, x1_) - n1(xi1) = mul(get_basis, xi1, n1_) - dn1(xi1) = mul(get_dbasis, xi1, n1_) - 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 - xi1_next = 0.0 - dxi1 = 0.0 - for i=1:max_iterations - dxi1 = -R(xi1)/dR(xi1) - dxi1 = clamp.(dxi1, -0.3, 0.3) - xi1_next = clamp.(xi1 + dxi1, -1.0, 1.0) - if norm(xi1_next - xi1) < tol - return xi1_next - end - if debug - info("xi1 = $xi1") - info("R(xi1) = $(R(xi1))") - info("dR(xi1) = $(dR(xi1))") - info("dxi1 = $dxi1") - info("norm = $(norm(xi1_next - xi1))") - info("xi1_next = $xi1_next") - end - xi1 = xi1_next - end - - info("x1 = $x1_") - info("n1 = $n1_") - info("x2 = $x2") - info("xi1 = $xi1, dxi1 = $dxi1") - info("-R(xi1) = $(-R(xi1))") - info("dR(xi1) = $(dR(xi1))") - error("find projection from master to slave: did not converge") - -end - -function project_from_slave_to_master_ad{E<:MortarElements2D}( - master_element::Element{E}, x1, n1, x2_; - tol=1.0e-10, max_iterations=20) - - x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) - dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) - 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 - -""" -Frictionless 2d finite sliding contact with forwarddiff. - -true/false flags: finite_sliding, friction, use_forwarddiff -""" -function assemble!(problem::Problem{Contact}, time::Float64, - ::Type{Val{1}}, ::Type{Val{true}}, - ::Type{Val{false}}, ::Type{Val{true}}) - - props = problem.properties - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - slave_elements = get_slave_elements(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) - gap = zeros(u) - C = zeros(la) - S = Set{Int64}() - - # 1. update nodal normals for slave elements - Q = [0.0 -1.0; 1.0 0.0] - normals = zeros(u) - for element in slave_elements - conn = get_connectivity(element) - push!(S, conn...) - gdofs = get_gdofs(problem, element) - X_el = element("geometry", time) - x_el = tuple( (X_el[i] + u[:,j] for (i,j) in enumerate(conn))... ) - #= - for ip in get_integration_points(element, 3) - dN = get_dbasis(element, ip, time) - 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 - =# - dN = get_dbasis(element, [0.0], time) - t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) - n = Q*t' - n /= norm(n) - for c in conn - normals[:,c] += n - end - end - for i in 1:size(normals,2) - normals[:,i] /= norm(normals[:,i]) - end - # swap element normals in 2d if they point to inside of body - if props.rotate_normals - for i=1:size(normals,2) - normals[:,i] = -normals[:,i] - end - end - normals2 = Dict() - for j in S - normals2[j] = normals[:,j] - end - update!(slave_elements, "normal", time => normals2) - - # 2. loop all slave elements - for slave_element in slave_elements - - slave_element_nodes = get_connectivity(slave_element) - X1 = slave_element("geometry", time) - u1 = ((u[:,i] for i in slave_element_nodes)...) - x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...) - la1 = ((la[:,i] for i in slave_element_nodes)...) - n1 = ((normals[:,i] for i in slave_element_nodes)...) - nnodes = size(slave_element, 2) - - # construct dual basis - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for master_element in slave_element("master elements", time) - - master_element_nodes = get_connectivity(master_element) - X2 = master_element("geometry", time) - u2 = ((u[:,i] for i in master_element_nodes)...) - x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) - - # calculate segmentation: we care only about endpoints - xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1]) - xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2]) - 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 - - for ip in get_integration_points(slave_element, 3) - # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip, time) - j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - w = ip.weight*norm(j)*l - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = get_basis(slave_element, xi_s, time) - De += w*diagm(vec(N1)) - Me += w*N1'*N1 - end - end - - Ae = De*inv(Me) - - # 3. loop all master elements - for master_element in slave_element("master elements", time) - - master_element_nodes = get_connectivity(master_element) - X2 = master_element("geometry", time) - u2 = ((u[:,i] for i in master_element_nodes)...) - x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) - - #x1_midpoint = 1/2*(x1[1]+x1[2]) - #x2_midpoint = 1/2*(x2[1]+x2[2]) - #distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint)) - #distance > props.maximum_distance && continue - - # calculate segmentation: we care only about endpoints - xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1]) - xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2]) - 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 - - slave_dofs = get_gdofs(problem, slave_element) - master_dofs = get_gdofs(problem, master_element) - - # 4. loop integration points of segment - for ip in get_integration_points(slave_element, 3) - # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip, time) - 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 = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - x_s = interpolate(N1, x1) # coordinate in gauss point - n_s = interpolate(N1, n1) # normal direction in gauss point - t_s = Q'*n_s # tangent direction in gauss point - xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2) - N2 = vec(get_basis(master_element, xi_m, time)) - x_m = interpolate(N2, x2) - Phi = Ae*N1 - - la_s = interpolate(Phi, la1) # traction force in gauss point - 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' - gap[1,slave_element_nodes] += w*gn*Phi - #gap[1,slave_element_nodes] += w*gn*N1' - - end # done integrating segment - - end # master elements done - - end # slave elements done - - # at this point we have calculated contact force fc and gap for all slave elements. - # next task is to find out are they in contact or not and remove inactive nodes - - state = problem.properties.contact_state_in_first_iteration - if problem.properties.iteration == 1 - info("First contact iteration, initial contact state = $state") - if state == :AUTO - avg_gap = ForwardDiff.value(mean([gap[1, j] for j in S])) - std_gap = ForwardDiff.value(std([gap[1, j] for j in S])) - if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12) - state = :ACTIVE - else - state = :UNKNOWN - end - info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state") - end - end - - is_active = Dict{Int, Bool}() - condition = Dict() - - for j in S - if j in props.always_in_contact - is_active[j] = true - continue - end - lan = dot(normals[:,j], la[:,j]) - condition[j] = ForwardDiff.value(lan - gap[1, j]) - is_active[j] = condition[j] > 0 - end - - if problem.properties.iteration == 1 && state == :ACTIVE - for j in S - is_active[j] = true - end - end - - if problem.properties.iteration == 1 && state == :INACTIVE - for j in S - is_active[j] = false - end - end - - if Logging._root.level == DEBUG - debug("Summary of nodes") - for j in sort(collect(keys(is_active))) - n = map(ForwardDiff.value, normals[:,j]) - debug("$j, c=$(condition[j]), s=$(is_active[j]), n=$n") - end - end - - for j in S - - if is_active[j] - n = normals[:,j] - t = Q'*n - lan = dot(n, la[:,j]) - lat = dot(t, la[:,j]) - C[1,j] += gap[1, j] - C[2,j] += lat - else - C[:,j] = la[:,j] - end - - end - - return vec([fc C]) - - end - - # x doesn't mean deformed configuration here - x = [problem.assembly.u; problem.assembly.la] - if length(x) == 0 - error("2d autodiff contact problem: initialize problem.assembly.u & la before solution") - end - - A = ForwardDiff.jacobian(calculate_interface, x) - b = calculate_interface(x) - A = sparse(A) - b = sparse(b) - SparseArrays.droptol!(A, 1.0e-9) - SparseArrays.droptol!(b, 1.0e-9) - - ndofs = round(Int, length(x)/2) - 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] - - f += C1*problem.assembly.la - g += D*problem.assembly.la - - #= - if !haskey(problem, "contact force") - problem.fields["contact force"] = Field(time => f) - else - update!(problem.fields["contact force"], time => f) - end - - fc = problem.fields["contact force"] - - if length(fc) > 1 - # kick in generalized alpha rule for time integration - alpha = 0.5 - info("Applying Generalized alpha time integration") - K = (1-alpha)*K - C1 = (1-alpha)*C1 - f = alpha*fc[end-1].data - end - =# - - problem.assembly.K = K - problem.assembly.C1 = transpose(C1) - problem.assembly.C2 = C2 - problem.assembly.D = D - problem.assembly.f = sparse(f) - problem.assembly.g = sparse(g) - -end - diff --git a/src/problems_mortar_2d_autodiff.jl b/src/problems_mortar_2d_autodiff.jl deleted file mode 100644 index 9975e9c..0000000 --- a/src/problems_mortar_2d_autodiff.jl +++ /dev/null @@ -1,247 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using ForwardDiff - -const MortarElements2D = Union{Seg2,Seg3} - -# forwarddiff version of mesh tying in 2d - -function project_from_master_to_slave_ad{E<:MortarElements2D}( - slave_element::Element{E}, x1_, n1_, x2, time; - tol=1.0e-10, max_iterations=20) - - x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_) - dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_) - n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_) - dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_) - 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_ad{E<:MortarElements2D}( - master_element::Element{E}, x1, n1, x2_, time; - tol=1.0e-10, max_iterations=20) - - x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) - dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) - 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 - -""" 2d mesh tie using ForwardDiff. - -Construct .. + fc*la and C(d,la)=0 - -""" -function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{true}}) - - props = problem.properties - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - slave_elements = get_slave_elements(problem) - if field_name != "displacement" - error("mortar forwarddiff assembly: only displacement field with adjust=yes supported") - end - - 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) - gap = zeros(u) - C = zeros(la) - - S = Set{Int64}() - # 1. update nodal normals for slave elements - tangents = zeros(u) - for element in slave_elements - conn = get_connectivity(element) - push!(S, conn...) - X1 = element("geometry", time) - u1 = ((u[:,i] for i in conn)...) - x1 = map(+, X1, u1) - dN = get_dbasis(element, [0.0], time) - tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - for nid in conn - tangents[:,nid] += tangent[:] - end - end - - Q = [0.0 -1.0; 1.0 0.0] - normals = zeros(u) - for j in S - tangents[:,j] /= norm(tangents[:,j]) - normals[:,j] = Q*tangents[:,j] - end - - if props.rotate_normals - for j in S - normals[:,j] = -normals[:,j] - end - end - - #update!(slave_elements, "normal", time => Dict(j => normals[:,j] for j in S)) - #update!(slave_elements, "tangent", time => Dict(j => tangents[:,j] for j in S)) - - # 2. loop all slave elements - for slave_element in slave_elements - - nsl = length(slave_element) - slave_element_nodes = get_connectivity(slave_element) - X1 = slave_element("geometry", time) - u1 = ((u[:,i] for i in slave_element_nodes)...) - x1 = map(+, X1, u1) - la1 = ((la[:,i] for i in slave_element_nodes)...) - n1 = ((normals[:,i] for i in slave_element_nodes)...) - - - # 3. loop all master elements - for master_element in slave_element("master elements", time) - - nm = length(master_element) - master_element_nodes = get_connectivity(master_element) - X2 = master_element("geometry", time) - u2 = ((u[:,i] for i in master_element_nodes)...) - x2 = map(+, X2, u2) - - # 3.1 calculate segmentation - xi1a = project_from_master_to_slave_ad(slave_element, x1, n1, x2[1], time) - xi1b = project_from_master_to_slave_ad(slave_element, x1, n1, x2[2], time) -# xi1a = project_from_master_to_slave(slave_element, X2[1], time) -# xi1b = project_from_master_to_slave(slave_element, X2[2], time) - 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 - - # 3.2. bi-orthogonal basis - De = zeros(nsl, nsl) - Me = zeros(nsl, nsl) - Ae = zeros(nsl, nsl) - if props.dual_basis - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - De += w*diagm(N1) - Me += w*N1*N1' - end - Ae = De*inv(Me) - else - Ae = eye(nsl) - end - - # 3.3. loop integration points of one integration segment and calculate - # local mortar matrices - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - #dN = get_dbasis(slave_element, ip, time) - #j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - #w = ip.weight*norm(j)*l - - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - Phi = Ae*N1 - # project gauss point from slave element to master element in direction n_s - x_s = interpolate(N1, x1) # coordinate in gauss point - n_s = interpolate(N1, n1) # normal direction in gauss point - #xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) - xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2, time) - N2 = vec(get_basis(master_element, xi_m, time)) - x_m = interpolate(N2, x2) - - la_s = interpolate(Phi, la1) - gn = dot(n_s, x_s-x_m) - - u_s = interpolate(N1, u1) - u_m = interpolate(N2, u2) - X_s = interpolate(N1, X1) - X_m = interpolate(N2, X2) - - fc[:,slave_element_nodes] += w*la_s*N1' - fc[:,master_element_nodes] -= w*la_s*N2' - #gap[1,slave_element_nodes] += w*gn*Phi' - gap[:,slave_element_nodes] += w*(u_s-u_m)*Phi' - if props.adjust - G = w*(X_s-X_m)*Phi' - gap[:,slave_element_nodes] += G - end - end - - end # master elements done - - end # slave elements done, contact virtual work ready - - C = gap - - 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 = ForwardDiff.jacobian(calculate_interface, x) - b = -calculate_interface(x) - - A = sparse(A) - b = sparse(b) - SparseArrays.droptol!(A, 1.0e-12) - SparseArrays.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] - - empty!(problem.assembly) - problem.assembly.K = K - problem.assembly.C1 = C1 - problem.assembly.C2 = C2 - problem.assembly.D = D - problem.assembly.f = f - problem.assembly.g = g - -end diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index 9911e59..826d85f 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -33,9 +33,9 @@ function Modal(nev=10, which=:SM) end """ Eliminate Dirichlet boundary condition from matrices K, M. """ -function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, - M_red::SparseMatrixCSC, - problem::Problem{Dirichlet}, ndim::Int) +function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC, + M_red::SparseMatrixCSC, + problem::Problem{Dirichlet}, ndim::Int) K = sparse(problem.assembly.K, ndim, ndim) C1 = sparse(problem.assembly.C1, ndim, ndim) C2 = sparse(problem.assembly.C2, ndim, ndim) @@ -100,10 +100,10 @@ end """ Eliminate mesh tie constraints from matrices K, M. """ -function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, - M_red::SparseMatrixCSC, - problem::Union{Problem{Mortar}, Problem{Mortar2D}}, - ndim::Int) +function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC, + M_red::SparseMatrixCSC, + problem::Union{Problem{Mortar}, Problem{Mortar2D}}, + ndim::Int) C1 = sparse(problem.assembly.C1, ndim, ndim) C2 = sparse(problem.assembly.C2, ndim, ndim) diff --git a/test/test_contact_2d_finite_sliding.jl b/test/test_contact_2d_finite_sliding.jl index 9a31a18..8a07d2f 100644 --- a/test/test_contact_2d_finite_sliding.jl +++ b/test/test_contact_2d_finite_sliding.jl @@ -36,15 +36,13 @@ using JuliaFEM.Testing update!(bc_lower, "displacement 1", 0.0) update!(bc_lower, "displacement 2", 0.0) - contact = Problem(Contact, "contact between upper and lower block", 2, "displacement") + contact = Problem(Contact2DAD, "contact between upper and lower block", 2, "displacement") contact.properties.rotate_normals = true - contact.properties.finite_sliding = true - contact.properties.friction = false - contact.properties.use_forwarddiff = true contact_slave_elements = create_elements(mesh, "LOWER_TOP") contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") - update!(contact_slave_elements, "master elements", contact_master_elements) - contact.elements = [contact_master_elements; contact_slave_elements] + add_slave_elements!(contact, contact_slave_elements) + add_master_elements!(contact, contact_master_elements) + nnodes = length(mesh.nodes) contact.assembly.u = zeros(2*nnodes) contact.assembly.la = zeros(2*nnodes) diff --git a/test/test_mortar_2d_mesh_tie_forwarddiff.jl b/test/test_mortar_2d_mesh_tie_forwarddiff.jl deleted file mode 100644 index b92e9f6..0000000 --- a/test/test_mortar_2d_mesh_tie_forwarddiff.jl +++ /dev/null @@ -1,178 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Preprocess -using JuliaFEM.Postprocess -using JuliaFEM.Testing - -function get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}}; - dy=0.0, adjust=false, tolerance=0.0, rotate_normals=false, swap=false, - dual_basis=false, use_forwarddiff=true, finite_strain=false, - geometric_stiffness=false) - - meshfile = @__DIR__() * "/testdata/block_2d_curved.med" - mesh = aster_read_mesh(meshfile) - - upper = Problem(Elasticity, "upper", 2) - upper.properties.formulation = :plane_stress - upper.properties.finite_strain = finite_strain - upper.properties.geometric_stiffness = geometric_stiffness - upper.elements = create_elements(mesh, "UPPER") - update!(upper.elements, "youngs modulus", 96.0) - update!(upper.elements, "poissons ratio", 1/3) - - lower = Problem(Elasticity, "lower", 2) - lower.properties.formulation = :plane_stress - lower.properties.finite_strain = finite_strain - lower.properties.geometric_stiffness = geometric_stiffness - lower.elements = create_elements(mesh, "LOWER") - update!(lower.elements, "youngs modulus", 96.0) - update!(lower.elements, "poissons ratio", 1/3) - - bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement") - bc_upper.elements = create_elements(mesh, "UPPER_TOP") - update!(bc_upper.elements, "displacement 1", 0.0) - update!(bc_upper.elements, "displacement 2", dy) - - bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement") - bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") - update!(bc_lower.elements, "displacement 1", 0.0) - update!(bc_lower.elements, "displacement 2", 0.0) - - interface = Problem(Mortar, "interface between upper and lower block", 2, "displacement") - interface_slave_elements = create_elements(mesh, "LOWER_TOP") - interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") - if swap - interface_slave_elements, interface_master_elements = interface_master_elements, interface_slave_elements - end - update!(interface_slave_elements, "master elements", interface_master_elements) - interface.elements = [interface_master_elements; interface_slave_elements] - interface.properties.adjust = adjust - interface.properties.distval = tolerance - interface.properties.rotate_normals = rotate_normals - interface.properties.dual_basis = dual_basis - interface.properties.use_forwarddiff = use_forwarddiff - interface.assembly.u = zeros(2*length(mesh.nodes)) - interface.assembly.la = zeros(2*length(mesh.nodes)) - - solver = Solver(Linear) - push!(solver, upper, lower, bc_upper, bc_lower, interface) - - return solver - -end - -#= - -@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=0.0" begin - # TODO: analytical solution now known, verify using other fem software - solver = get_model("mesh tie with curved 2d block"; - adjust=false, tolerance=10, dy=-0.1, rotate_normals=true, - dual_basis=true, use_forwarddiff=true, finite_strain=false, - geometric_stiffness=false) - solver() - interface = solver["interface between upper and lower block"] - @test isapprox(norm(interface.assembly.u), 0.11339715157447851) -end - - -@testset "curved surface with adjust=true, dual lagrange, slave=lower surface, dy=0.0" begin - # TODO: analytical solution now known, verify using other fem software - solver = get_model("mesh tie with curved 2d block"; - adjust=true, tolerance=10, dy=0.0, rotate_normals=true, - dual_basis=true, use_forwarddiff=true) - solver() - interface = solver["interface between upper and lower block"] - @test solver.properties.iteration == 2 - # differs -- why? - @test isapprox(norm(interface.assembly.u), 0.11660422877751599) -end - -@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=-0.1" begin - # TODO: analytical solution now known, verify using other fem software - solver = get_model("mesh tie with curved 2d block"; - adjust=true, tolerance=10, dy=-0.1, rotate_normals=true, - dual_basis=false, use_forwarddiff=true) - solver() - interface = solver["interface between upper and lower block"] - @test solver.properties.iteration == 2 - @test isapprox(norm(interface.assembly.u), 0.34230262165505887) -end - -@testset "curved surface, adjust=true, dual basis, slave=lower surface, dy=-0.1" begin - # TODO: analytical solution now known, verify using other fem software - solver = get_model("mesh tie with curved 2d block"; - adjust=true, tolerance=10, dy=-0.1, rotate_normals=true, - dual_basis=true, use_forwarddiff=true) - solver() - interface = solver["interface between upper and lower block"] - @test solver.properties.iteration == 2 - @test isapprox(norm(interface.assembly.u), 0.34318800698017704) -end - -=# - - - -@testset "compare forwarddiff solution to normal" begin - X = Dict( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [0.0, 1.0], - 4 => [1.0, 1.0]) - u = Dict( - 1 => [0.0, 0.0], - 2 => [0.0, 0.0], - 3 => [0.0, 0.0], - 4 => [0.0, 0.0]) - sel1 = Element(Seg2, [1, 2]) - mel1 = Element(Seg2, [3, 4]) - update!([sel1, mel1], "geometry", X) - update!([sel1, mel1], "displacement", u) - update!(sel1, "master elements", [mel1]) - - p1 = Problem(Mortar2D, "test 1", 2, "displacement") - add_slave_elements!(p1, [sel1]) - add_master_elements!(p1, [mel1]) - assemble!(p1, 0.0) - - p2 = Problem(Mortar, "test 2", 2, "displacement") - push!(p2, sel1, mel1) - p2.properties.use_forwarddiff = true - p2.assembly.u = zeros(8) - p2.assembly.la = zeros(8) - assemble!(p2, 0.0) - - @test isapprox(p1.assembly, p2.assembly) - - #= - empty!(p1.assembly) - empty!(p2.assembly) - p1.properties.adjust = true - p2.properties.adjust = true - assemble!(p1, 0.0) - assemble!(p2, 0.0) - C11 = full(p1.assembly.C1, 4, 8) - C12 = full(p2.assembly.C1, 4, 8) - C21 = full(p1.assembly.C2, 4, 8) - C22 = full(p2.assembly.C2, 4, 8) - D1 = full(p1.assembly.D) - D2 = full(p2.assembly.D) - g1 = full(p1.assembly.g, 4, 1) - g2 = full(p2.assembly.g, 4, 1) - println("C1") - dump(C11) - dump(C12) - println("C2") - dump(C21) - dump(C22) - println("D") - dump(D1) - dump(D2) - println("g") - dump(g1) - dump(g2) - @test isapprox(p1.assembly, p2.assembly) - =# -end diff --git a/test/test_problems_contact_2d_autodiff.jl b/test/test_problems_contact_2d_autodiff.jl index 599032d..95fbc57 100644 --- a/test/test_problems_contact_2d_autodiff.jl +++ b/test/test_problems_contact_2d_autodiff.jl @@ -6,88 +6,77 @@ using JuliaFEM.Preprocess using JuliaFEM.Postprocess using JuliaFEM.Testing -datadir = first(splitext(basename(@__FILE__))) +pkg_dir = Pkg.dir("JuliaFEM") +datadir = joinpath(pkg_dir, "test", first(splitext(basename(@__FILE__)))) -function get_model() - meshfile = joinpath(datadir, "block_2d.med") - mesh = aster_read_mesh(meshfile) - #error("mesh has $(length(mesh.nodes)) nodes") +meshfile = joinpath(datadir, "block_2d.med") +mesh = aster_read_mesh(meshfile) - upper = Problem(mesh, Elasticity, "UPPER", 2) - lower = Problem(mesh, Elasticity, "LOWER", 2) +upper = Problem(mesh, Elasticity, "UPPER", 2) +lower = Problem(mesh, Elasticity, "LOWER", 2) - for body in [upper, lower] - body.properties.formulation = :plane_stress - update!(body, "youngs modulus", 288.0) - update!(body, "poissons ratio", 1/3) - end - - load = Problem(mesh, Elasticity, "UPPER_TOP", 2) - load.properties.formulation = :plane_stress - update!(load, "displacement traction force 2", 0.0 => 0.0) - update!(load, "displacement traction force 2", 1.0 => -28.8) - bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement") - update!(bc1, "displacement 2", 0.0) - bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement") - update!(bc2, "displacement 1", 0.0) - bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement") - update!(bc3, "displacement 1", 0.0) - - interface = Problem(Contact, "interface", 2, "displacement") - interface_slave_elements = create_elements(mesh, "LOWER_TOP") - interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") - update!(interface_slave_elements, "master elements", interface_master_elements) - interface.elements = [interface_master_elements; interface_slave_elements] - interface.properties.rotate_normals = true - - # in LOWER_LEFT we have node belonging also to contact interface - # let's remove it from dirichlet bc - create_node_set_from_element_set!(mesh, "LOWER_LEFT") - nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT") - coords = mesh.nodes[nid] - info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords") - dofs = [2*(nid-1)+1, 2*(nid-1)+2] - info("removing nid $nid, dofs $dofs from LOWER_LEFT") - push!(bc2.assembly.removed_dofs, dofs...) - - solver = Solver(Nonlinear) - push!(solver, upper, lower, load, bc1, bc2, bc3, interface) - return solver +for body in [upper, lower] + body.properties.formulation = :plane_stress + update!(body, "youngs modulus", 288.0) + update!(body, "poissons ratio", 1/3) end -@testset "finite sliding 2d patch test, linear Seg2 elements, standard basis" begin +load = Problem(mesh, Elasticity, "UPPER_TOP", 2) +load.properties.formulation = :plane_stress +update!(load, "displacement traction force 2", 0.0 => 0.0) +update!(load, "displacement traction force 2", 1.0 => -28.8) +bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement") +update!(bc1, "displacement 2", 0.0) +bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement") +update!(bc2, "displacement 1", 0.0) +bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement") +update!(bc3, "displacement 1", 0.0) - solver = get_model() - interface = solver["interface"] - interface.assembly.u = zeros(48) - interface.assembly.la = zeros(48) - upper = solver["UPPER"] - lower = solver["LOWER"] - for body in [upper, lower] - body.properties.geometric_stiffness = true - body.properties.finite_strain = true - end - interface.properties.finite_sliding = true - interface.properties.use_forwarddiff = true +interface = Problem(Contact2DAD, "interface", 2, "displacement") +interface_slave_elements = create_elements(mesh, "LOWER_TOP") +interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") +add_slave_elements!(interface, interface_slave_elements) +add_master_elements!(interface, interface_master_elements) +interface.properties.rotate_normals = true - for time in [0.0, 1/3, 2/3, 1.0] - interface.properties.iteration = 1 - solve!(solver, time) - end +# in LOWER_LEFT we have node belonging also to contact interface +# let's remove it from dirichlet bc +create_node_set_from_element_set!(mesh, "LOWER_LEFT") +nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT") +coords = mesh.nodes[nid] +info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords") +dofs = [2*(nid-1)+1, 2*(nid-1)+2] +info("removing nid $nid, dofs $dofs from LOWER_LEFT") +push!(bc2.assembly.removed_dofs, dofs...) - node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0) - node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0) - node_ids, la = get_nodal_vector(get_slave_elements(interface), "lambda", 1.0) - u2 = [u[2] for u in displacement] - f2 = [f[2] for f in la] - maxabsu2 = maximum(abs.(u2)) - stdabsu2 = std(abs.(u2)) - info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2") - @test isapprox(stdabsu2, 0.0; atol=1.0e-12) - maxabsf2 = maximum(abs.(f2)) - stdabsf2 = std(abs.(f2)) - info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2") - @test isapprox(stdabsf2, 0.0; atol=1.0e-12) - # for linear case pressure 28.8 - @test isapprox(mean(abs.(f2)), 27.76616800689944; rtol=1.0e-3) +solver = Solver(Nonlinear) +push!(solver, upper, lower, load, bc1, bc2, bc3, interface) + +interface.assembly.u = zeros(48) +interface.assembly.la = zeros(48) + +for body in [upper, lower] + body.properties.geometric_stiffness = true + body.properties.finite_strain = true end + +for time in [0.0, 1/3, 2/3, 1.0] + interface.properties.iteration = 0 + solve!(solver, time) +end + +node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0) +node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0) +node_ids, la = get_nodal_vector(interface.elements, "lambda", 1.0) +u2 = [u[2] for u in displacement] +f2 = [f[2] for f in la] +maxabsu2 = maximum(abs.(u2)) +stdabsu2 = std(abs.(u2)) +info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2") +@test isapprox(stdabsu2, 0.0; atol=1.0e-12) +maxabsf2 = maximum(abs.(f2)) +stdabsf2 = std(abs.(f2)) +info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2") +@test isapprox(stdabsf2, 0.0; atol=1.0e-12) +# for linear case pressure 28.8 +@test isapprox(mean(abs.(f2)), 27.76616800689944; rtol=1.0e-3)