From c02d510673edfda8bc6a8df4dbff04bffa5126ff Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 19 Jun 2016 20:01:37 +0300 Subject: [PATCH] fixed a lot of tests --- src/JuliaFEM.jl | 12 +- src/contact.jl | 130 ++++++ src/elasticity.jl | 111 ++--- src/heat.jl | 33 +- src/mortar.jl | 258 +++++++----- src/problems.jl | 4 +- src/solvers.jl | 26 +- src/sparse.jl | 4 + test/runtests.jl | 47 +-- test/test_abaqus_reader.jl | 49 +-- test/test_api.jl | 19 +- test/test_assembly.jl | 9 +- test/test_elasticity.jl | 82 ---- ..._elasticity_2d_linear_with_surface_load.jl | 14 +- ...asticity_2d_nonlinear_with_surface_load.jl | 26 +- test/test_elasticity_2d_residual.jl | 59 +++ ...asticity_3d_nonlinear_with_surface_load.jl | 2 +- test/test_elasticity_3d_unit_block.jl | 9 +- test/test_mortar_2d.jl | 393 +++--------------- test/test_mortar_2d_assembly.jl | 111 +++++ test/test_mortar_2d_calculate_projection.jl | 85 ++++ 21 files changed, 789 insertions(+), 694 deletions(-) create mode 100644 src/contact.jl delete mode 100644 test/test_elasticity.jl create mode 100644 test/test_elasticity_2d_residual.jl create mode 100644 test/test_mortar_2d_assembly.jl create mode 100644 test/test_mortar_2d_calculate_projection.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index db8fc2b..68d021c 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -73,8 +73,16 @@ export Modal include("optics.jl") export find_intersection, calc_reflection, calc_normal -### MORTAR STUFF ### -include("mortar.jl") # mortar projection +### Mortar methods ### +include("mortar.jl") +export calculate_normals, + calculate_normals!, + project_from_slave_to_master, + project_from_master_to_slave, + Mortar + +### Contact mechanics ### +#include("contact.jl") # rest of things include("utils.jl") diff --git a/src/contact.jl b/src/contact.jl new file mode 100644 index 0000000..fa3d06f --- /dev/null +++ b/src/contact.jl @@ -0,0 +1,130 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" +Currently two strategies exists: + +a) Remove inactive inequality constraints in element level. This is done in + assemble! if normal_condition is set to :Contact. For some reason this + leads to convergence issues. +b) Remove inactive inequality constraints in assembly level. This is done in + posthook algorithm if inequality_constraints is set to true. This gives + more robust behavior. + + Either use inequality_constraints=True OR :Contact + :Slip, but do not mix. + + minimum_distance can be used to roughly skip integration of mortar + projections for elements that are "far enough" from each other. Increases + performance. + +""" +type Mortar <: BoundaryProblem + formulation :: Symbol # :total, :incremental, :autodiff + dual_basis :: Bool + inequality_constraints :: Bool # Launch PDASS to solve inequality constraints + normal_condition :: Symbol # Tie or Contact + tangential_condition :: Symbol # Stick or Slip + maximum_distance :: Float64 # don't check for a contact if elements are far enough + store_debug_info :: Bool # for making debugging easier + always_inactive :: Vector{Int64} + always_in_contact :: Vector{Int64} # nodes in this list always in contact + always_in_stick :: Vector{Int64} # nodes in this list always in stick + always_in_slip :: Vector{Int64} # nodes in this list always in slip + contact :: Bool + friction :: Bool + gap_sign :: Int # gap sign convention + rotate_normals :: Bool +end + +function Mortar() + Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false) +end + +function get_unknown_field_name(::Type{Mortar}) + return "reaction force" +end + +function get_formulation_type(problem::Problem{Mortar}) + return problem.properties.formulation +end + +macro debug(msg) + haskey(ENV, "DEBUG") || return + return msg +end + +function assemble!(problem::Problem{Mortar}, time::Real) + elements = get_elements(problem) + if length(elements) == 0 + info("$(typeof(problem)) : forget to add elements?") + return + end + # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc. + eldim = size(elements[1], 1)+1 + assemble!(problem, time, Val{eldim}) +end + +include("mortar_2d.jl") +include("mortar_2d_autodiff.jl") +include("mortar_3d.jl") +include("mortar_3d_autodiff.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) + problem.properties.inequality_constraints || return + info("PDASS: Starting primal-dual active set strategy to determine active constraints") + S = Set{Int64}() + for element in get_elements(problem) + haskey(element, "master elements") || continue + push!(S, get_connectivity(element)...) + end + S = sort(collect(S)) + dim = get_unknown_field_dimension(problem) + ndofs = solver.ndofs + nnodes = round(Int, ndofs/dim) + + c = reshape(full(problem.assembly.c, ndofs, 1), dim, nnodes) + A = find(c[1,:] .> 0) + A = intersect(A, S) + I = setdiff(S, A) + + info("PDASS: contact nodes: $(sort(collect(S)))") + info("PDASS: active nodes: $(sort(collect(A)))") + info("PDASS: inactive nodes: $(sort(collect(I)))") + + # remove any inactive nodes + for j in I + dofs = [dim*(j-1)+i for i=1:dim] + C1[dofs,:] = 0 + C2[dofs,:] = 0 + D[dofs,:] = 0 + g[dofs,:] = 0 + end + + # handle tangential condition for active nodes + if problem.properties.tangential_condition == :Slip + for j in A + dofs = [dim*(j-1)+i for i=1:dim] + tangential_dofs = dofs[2:end] + D[tangential_dofs,dofs] = C2[tangential_dofs,dofs] + C2[tangential_dofs,:] = 0 + g[tangential_dofs,:] = 0 + end + end + + return +end + +function assemble_prehook!(problem::Problem{Mortar}, time::Real) + info("mortar assemble prehook at time $time") + slaves = Set{Element}() + for element in get_elements(problem) + haskey(element, "master elements") || continue + push!(slaves, element) + end + info("$(length(slaves)) slave elements") + length(slaves) != 0 || error("no slave elements found for problem (forget to add masters?).") + info("mortar: update normal-tangential system.") + calculate_normal_tangential_coordinates!(collect(slaves), time) + info("mortar assemble prehook done.") +end diff --git a/src/elasticity.jl b/src/elasticity.jl index 20a203c..ec8df38 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -21,7 +21,6 @@ function get_unknown_field_name(problem::Problem{Elasticity}) end function get_formulation_type(problem::Problem{Elasticity}) - # we are solving residual and add increment to previous solution vector return :incremental end @@ -57,21 +56,37 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem N = element(ip, time) dN = element(ip, time, Val{:Grad}) - # kinematics; calculate deformation gradient and strain - gradu = zeros(dim, dim) - if haskey(element, "displacement") - gradu += element("displacement", ip, time, Val{:Grad}) - end - strain = zeros(dim , dim) - strain += 1/2*(gradu' + gradu) - F = eye(dim) + # kinematics + + gradu = element("displacement", ip, time, Val{:Grad}) + fill!(BL, 0.0) + if props.finite_strain - F += gradu - strain += 1/2*gradu'*gradu + strain = 1/2*(gradu + gradu' + gradu'*gradu) + F = eye(dim) + gradu + for i=1:size(dN, 2) + BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i] + BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i] + BL[2, 2*(i-1)+1] += F[1,2]*dN[2,i] + BL[2, 2*(i-1)+2] += F[2,2]*dN[2,i] + BL[3, 2*(i-1)+1] += F[1,1]*dN[2,i] + F[1,2]*dN[1,i] + BL[3, 2*(i-1)+2] += F[2,1]*dN[2,i] + F[2,2]*dN[1,i] + end + else # linearized strain + strain = 1/2*(gradu + gradu') + F = eye(dim) + for i=1:size(dN, 2) + BL[1, 2*(i-1)+1] = dN[1,i] + BL[2, 2*(i-1)+2] = dN[2,i] + BL[3, 2*(i-1)+1] = dN[2,i] + BL[3, 2*(i-1)+2] = dN[1,i] + end end - # constitutive equations; material model (isotropic linear material here) - # get_material(problem, element, ...) + strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]] + update!(ip, "strain", time => strain_vec) + + # calculate stress E = element("youngs modulus", ip, time) nu = element("poissons ratio", ip, time) if props.formulation == :plane_stress @@ -88,53 +103,48 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem error("unknown plane formulation: $(props.formulation)") end # calculate stress - strain_vec = [strain[1,1]; strain[2,2]; 2.0*strain[1,2]] - stress_vec = D*strain_vec - stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] - cauchy_stress = F'*stress*F/det(F) - cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] + stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + update!(ip, "stress", time => stress_vec) - update!(ip, "strain", time => strain_vec) - update!(ip, "cauchy stress", time => cauchy_stress) - update!(ip, "pk2 stress", time => stress_vec) + Km += w*BL'*D*BL - # add contributions: material and geometric stiffness + internal forces - fill!(BL, 0.0) - for i=1:size(dN, 2) - BL[1, 2*(i-1)+1] = F[1,1]*dN[1,i] - BL[1, 2*(i-1)+2] = F[2,1]*dN[1,i] - BL[2, 2*(i-1)+1] = F[1,2]*dN[2,i] - BL[2, 2*(i-1)+2] = F[2,2]*dN[2,i] - BL[3, 2*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i] - BL[3, 2*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i] - end - fill!(BNL, 0.0) - for i=1:size(dN, 2) - BNL[1, 2*(i-1)+1] = dN[1,i] - BNL[2, 2*(i-1)+1] = dN[2,i] - BNL[3, 2*(i-1)+2] = dN[1,i] - BNL[4, 2*(i-1)+2] = dN[2,i] - end - S2 = zeros(2*dim, 2*dim) - S2[1,1] = stress_vec[1] - S2[2,2] = stress_vec[2] - S2[1,2] = S2[2,1] = stress_vec[3] - S2[3:4,3:4] = S2[1:2,1:2] + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] + # cauchy_stress = F'*stress*F/det(F) + # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] + # update!(ip, "cauchy stress", time => cauchy_stress) + + # material stiffness end - Km += w*BL'*D*BL # material stiffness - if props.finite_strain # add geometric stiffness + if props.geometric_stiffness + # take geometric stiffness into account + + fill!(BNL, 0.0) + for i=1:size(dN, 2) + BNL[1, 2*(i-1)+1] = dN[1,i] + BNL[2, 2*(i-1)+1] = dN[2,i] + BNL[3, 2*(i-1)+2] = dN[1,i] + BNL[4, 2*(i-1)+2] = dN[2,i] + end + + S2 = zeros(2*dim, 2*dim) + S2[1,1] = stress_vec[1] + S2[2,2] = stress_vec[2] + S2[1,2] = S2[2,1] = stress_vec[3] + S2[3:4,3:4] = S2[1:2,1:2] + Kg += w*BNL'*S2*BNL # geometric stiffness + end - if get_formulation_type(problem) == :incremental - f -= w*BL'*stress_vec # internal force - end + # rhs, internal and external load + + f -= w*BL'*stress_vec - # volume load if haskey(element, "displacement load") b = element("displacement load", ip, time) f += w*vec(N'*b) end + for i=1:dim if haskey(element, "displacement load $i") b = element("displacement load $i", ip, time) @@ -409,7 +419,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el # material stiffness end - if props.geometric_stiffness # take geometric stiffness into account + if props.geometric_stiffness + # take geometric stiffness into account fill!(BNL, 0.0) diff --git a/src/heat.jl b/src/heat.jl index 0f5d064..82bccbb 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -39,31 +39,32 @@ function get_unknown_field_type(problem::Problem{Heat}) end function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0) - gdofs = get_gdofs(problem, element) - + field_name = get_unknown_field_name(problem) + nnodes = length(element) + K = zeros(nnodes, nnodes) + fq = zeros(nnodes) for ip in get_integration_points(element) detJ = element(ip, time, Val{:detJ}) w = ip.weight*detJ - N = element(ip, time) - if haskey(element, "density") - rho = element("density", ip, time) - add!(assembly.M, gdofs, gdofs, w*rho*N'*N) - end - if haskey(element, "temperature thermal conductivity") + if haskey(element, "$field_name thermal conductivity") dN = element(ip, time, Val{:Grad}) - k = element("temperature thermal conductivity", ip, time) - add!(assembly.K, gdofs, gdofs, w*k*dN'*dN) + k = element("$field_name thermal conductivity", ip, time) + K += w*k*dN'*dN end - if haskey(element, "temperature load") - f = element("temperature load", ip, time) - add!(assembly.f, gdofs, w*N'*f) + if haskey(element, "$field_name load") + f = element("$field_name load", ip, time) + fq += w*N'*f end - if haskey(element, "temperature flux") - g = element("temperature flux", ip, time) - add!(assembly.f, gdofs, w*N'*g) + if haskey(element, "$field_name flux") + g = element("$field_name flux", ip, time) + fq += w*N'*g end end + T = vec(element[field_name](time)) + fq -= K*T + add!(assembly.K, gdofs, gdofs, K) + add!(assembly.f, gdofs, fq) end diff --git a/src/mortar.jl b/src/mortar.jl index fa3d06f..ffb79f0 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -1,43 +1,14 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -""" -Currently two strategies exists: - -a) Remove inactive inequality constraints in element level. This is done in - assemble! if normal_condition is set to :Contact. For some reason this - leads to convergence issues. -b) Remove inactive inequality constraints in assembly level. This is done in - posthook algorithm if inequality_constraints is set to true. This gives - more robust behavior. - - Either use inequality_constraints=True OR :Contact + :Slip, but do not mix. - - minimum_distance can be used to roughly skip integration of mortar - projections for elements that are "far enough" from each other. Increases - performance. - -""" type Mortar <: BoundaryProblem - formulation :: Symbol # :total, :incremental, :autodiff - dual_basis :: Bool - inequality_constraints :: Bool # Launch PDASS to solve inequality constraints - normal_condition :: Symbol # Tie or Contact - tangential_condition :: Symbol # Stick or Slip - maximum_distance :: Float64 # don't check for a contact if elements are far enough - store_debug_info :: Bool # for making debugging easier - always_inactive :: Vector{Int64} - always_in_contact :: Vector{Int64} # nodes in this list always in contact - always_in_stick :: Vector{Int64} # nodes in this list always in stick - always_in_slip :: Vector{Int64} # nodes in this list always in slip - contact :: Bool - friction :: Bool - gap_sign :: Int # gap sign convention rotate_normals :: Bool + adjust :: Bool + tolerance :: Float64 end function Mortar() - Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false) + return Mortar(false, false, 0.0) end function get_unknown_field_name(::Type{Mortar}) @@ -45,86 +16,173 @@ function get_unknown_field_name(::Type{Mortar}) end function get_formulation_type(problem::Problem{Mortar}) - return problem.properties.formulation + return :incremental end -macro debug(msg) - haskey(ENV, "DEBUG") || return - return msg +typealias MortarElements2D Union{Seg2, Seg3} +typealias MortarElements3D Union{Tri3, Tri6, Quad4} + +function newton(f, df, x; tol=1.0e-6, max_iterations=10) + for i=1:max_iterations + dx = -f(x)/df(x) + x += dx + if norm(dx) < tol + return x + end + end + error("Newton iteration did not converge in $max_iterations iterations") end -function assemble!(problem::Problem{Mortar}, time::Real) - elements = get_elements(problem) - if length(elements) == 0 - info("$(typeof(problem)) : forget to add elements?") - return - end - # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc. - eldim = size(elements[1], 1)+1 - assemble!(problem, time, Val{eldim}) +function cross2(a, b) + cross([a; 0], [b; 0])[3] end -include("mortar_2d.jl") -include("mortar_2d_autodiff.jl") -include("mortar_3d.jl") -include("mortar_3d_autodiff.jl") +function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time) + x1_ = slave_element["geometry"](time) + n1_ = slave_element["normal"](time) + x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_ + dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_ + n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_ + dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_ + R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) + dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) + xi1 = newton(R, dR, 0.0) + return xi1 +end -""" Remove inactive inequality constraints by using primal-dual active set strategy. """ -function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g) - problem.properties.inequality_constraints || return - info("PDASS: Starting primal-dual active set strategy to determine active constraints") - S = Set{Int64}() - for element in get_elements(problem) - haskey(element, "master elements") || continue - push!(S, get_connectivity(element)...) - end - S = sort(collect(S)) - dim = get_unknown_field_dimension(problem) - ndofs = solver.ndofs - nnodes = round(Int, ndofs/dim) +function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time) + x2_ = master_element["geometry"](time) + x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_ + dx2(xi2) = 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 = newton(R, dR, 0.0) + return xi2 +end - c = reshape(full(problem.assembly.c, ndofs, 1), dim, nnodes) - A = find(c[1,:] .> 0) - A = intersect(A, S) - I = setdiff(S, A) - - info("PDASS: contact nodes: $(sort(collect(S)))") - info("PDASS: active nodes: $(sort(collect(A)))") - info("PDASS: inactive nodes: $(sort(collect(I)))") - - # remove any inactive nodes - for j in I - dofs = [dim*(j-1)+i for i=1:dim] - C1[dofs,:] = 0 - C2[dofs,:] = 0 - D[dofs,:] = 0 - g[dofs,:] = 0 - end - - # handle tangential condition for active nodes - if problem.properties.tangential_condition == :Slip - for j in A - dofs = [dim*(j-1)+i for i=1:dim] - tangential_dofs = dofs[2:end] - D[tangential_dofs,dofs] = C2[tangential_dofs,dofs] - C2[tangential_dofs,:] = 0 - g[tangential_dofs,:] = 0 +function calculate_normals(elements, time, rotate_normals=false) + tangents = Dict{Int64, Vector{Float64}}() + for element in elements + conn = get_connectivity(element) + X1 = element("geometry", time) + dN = get_dbasis(element, [0.0], time) + tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)])) + for nid in conn + if haskey(tangents, nid) + tangents[nid] += tangent + else + tangents[nid] = tangent + end end end - return + Q = [0.0 -1.0; 1.0 0.0] + normals = Dict{Int64, Vector{Float64}}() + S = sort(collect(keys(tangents))) + for j in S + tangents[j] /= norm(tangents[j]) + normals[j] = Q*tangents[j] + end + + if rotate_normals + for j in S + normals[j] = -normals[j] + end + end + + return normals, tangents end -function assemble_prehook!(problem::Problem{Mortar}, time::Real) - info("mortar assemble prehook at time $time") - slaves = Set{Element}() - for element in get_elements(problem) - haskey(element, "master elements") || continue - push!(slaves, element) +function calculate_normals!(elements, time, rotate_normals=false) + normals, tangents = calculate_normals(elements, time, rotate_normals) + for element in elements + conn = get_connectivity(element) + update!(element, "normal", time => [normals[j] for j in conn]) + update!(element, "tangent", time => [tangents[j] for j in conn]) end - info("$(length(slaves)) slave elements") - length(slaves) != 0 || error("no slave elements found for problem (forget to add masters?).") - info("mortar: update normal-tangential system.") - calculate_normal_tangential_coordinates!(collect(slaves), time) - info("mortar assemble prehook done.") +end + +function assemble!(problem::Problem{Mortar}, time::Real) + + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_elements = filter(el -> haskey(el, "master elements"), get_elements(problem)) + + # 1. calculate nodal normals and tangents for slave element nodes j ∈ S + normals, tangents = calculate_normals(slave_elements, time, props.rotate_normals) + update!(slave_elements, "normal", normals) + update!(slave_elements, "tangent", tangents) + S = sort(collect(keys(normals))) + + # 2. loop all slave elements + for slave_element in slave_elements + haskey(slave_element, "master elements") || continue + + slave_element_nodes = get_connectivity(slave_element) + nsl = length(slave_element) + X1 = slave_element["geometry"](time) + n1 = Field([normals[j] for j in slave_element_nodes]) + + # 3. loop all master elements + for master_element in slave_element["master elements"](time) + + master_element_nodes = get_connectivity(master_element) + nm = length(master_element) + X2 = master_element["geometry"](time) + + # 3.1 calculate segmentation + xi1a = project_from_master_to_slave(slave_element, X2[1], time) + xi1b = project_from_master_to_slave(slave_element, X2[end], 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.3. loop integration points of one integration segment and calculate + # local mortar matrices + De = zeros(nsl, nsl) + Me = zeros(nsl, nm) + ge = zeros(nsl) + for ip in get_integration_points(slave_element, 2) + 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)) + # project gauss point from slave element to master element in direction n_s + X_s = N1*X1 # coordinate in gauss point + n_s = N1*n1 # normal direction in gauss point + xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) + N2 = vec(get_basis(master_element, xi_m, time)) + X_m = N2*X2 + De += w*N1*N1' + Me += w*N1*N2' + if props.adjust + g = X_s-X_m + if g < props.tol + ge += w*g + end + end + end + + # add contribution to contact virtual work + sdofs = get_gdofs(problem, slave_element) + mdofs = get_gdofs(problem, master_element) + + for i=1:field_dim + lsdofs = sdofs[i:field_dim:end] + lmdofs = mdofs[i:field_dim:end] + add!(problem.assembly.C1, lsdofs, lsdofs, De) + add!(problem.assembly.C1, lsdofs, lmdofs, -Me) + add!(problem.assembly.C2, lsdofs, lsdofs, De) + add!(problem.assembly.C2, lsdofs, lmdofs, -Me) + add!(problem.assembly.g, lsdofs, ge) + end + + + end # master elements done + + end # slave elements done, contact virtual work ready + end diff --git a/src/problems.jl b/src/problems.jl index baa950c..8b9239b 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -111,11 +111,11 @@ function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem, elements= end function get_formulation_type{P<:FieldProblem}(problem::Problem{P}) - return :total + return :incremental end function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P}) - return :total + return :incremental end function get_assembly(problem) diff --git a/src/solvers.jl b/src/solvers.jl index e0431e4..2191a0a 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -256,7 +256,7 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}}) # assemble boundary problems Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) - K = K + Kb + K = K + Kb + Kg f = f + fb K = 1/2*(K + K') u = zeros(solver.ndofs) @@ -279,7 +279,7 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}}) end # solver interior - CF = cholfact(K[interior_dofs, interior_dofs]) + CF = ldltfact(K[interior_dofs, interior_dofs]) Kib = K[interior_dofs, boundary_dofs] Kbb = K[boundary_dofs, boundary_dofs] fi = f[interior_dofs] @@ -337,6 +337,19 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError) print(io, "nonlinear iteration did not converge in $max_iters iterations!") end +function assemble!(solver::Solver; force_assembly=true) + info("Assembling problems ...") + tic() + for problem in solver.problems + if force_assembly # force reassembly + problem.assembly.changed = true + end + assemble!(problem, solver.time) + end + t1 = round(toq(), 2) + info("Assembled in $t1 seconds.") +end + """ Default solver for quasistatic nonlinear problems. """ function call(solver::Solver{Nonlinear}) @@ -352,14 +365,7 @@ function call(solver::Solver{Nonlinear}) info("Starting nonlinear iteration #$(properties.iteration)") # 2.1 update linearized assemblies (if needed) - info("Assembling problems ...") - tic() - for problem in solver.problems - problem.assembly.changed = true # force reassembly - assemble!(problem, solver.time) - end - t1 = round(toq(), 2) - info("Assembled in $t1 seconds.") + assemble!(solver) # 2.2 call solver for linearized system (default: direct lu factorization) info("Solve linear system ...") diff --git a/src/sparse.jl b/src/sparse.jl index 950d42d..799ad8c 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -155,6 +155,10 @@ function get_nonzero_rows(A::SparseMatrixCOO) return get_nonzero_rows(sparse(A)) end +function get_nonzero_rows(A::Matrix) + return get_nonzero_rows(sparse(A)) +end + function size(A::SparseMatrixCOO) return maximum(A.I), maximum(A.J) end diff --git a/test/runtests.jl b/test/runtests.jl index ace3b1d..1516239 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,15 +1,23 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -# should this stuff be in package? see FactCheck docs. - - using JuliaFEM using JuliaFEM.Test -#using FactCheck -#using Logging -#@Logging.configure(level=DEBUG) +function run_tests(; quiet=false) + + test_files = readdir(Pkg.dir("JuliaFEM")*"/test") + test_files = filter(f -> (startswith(f, "test_") & endswith(f, ".jl")), test_files) + for test_file in test_files + if !quiet + info("Running tests from file $test_file") + end + include(test_file) + end + +end + +run_tests() #= facts("Testing if somebody used print, println(), @sprint in src directory") do @@ -99,30 +107,3 @@ end =# -### NEW STYLE OF TESTING - -using JuliaFEM.Test - -function run_tests() - - for test_file in readdir(Pkg.dir("JuliaFEM")*"/test") - info("checking is $test_file is real test file") - if (startswith(test_file, "test_")) & (endswith(test_file, ".jl")) - run_test(test_file) - end - end - - passed, failed, errors, critical = print_test_statistics() - - # at the very end throw error if something is failed - if failed + errors + critical > 0 - error("Some tests has failed. Fix them. Now.") - exit(1) - else - info("""All tests has passed \o/ .""") - exit(0) - end - -end - -run_tests() diff --git a/test/test_abaqus_reader.jl b/test/test_abaqus_reader.jl index 21ce83b..49115d9 100644 --- a/test/test_abaqus_reader.jl +++ b/test/test_abaqus_reader.jl @@ -1,15 +1,11 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -module AbaqusReaderTests - -#using JuliaFEM +using JuliaFEM +using JuliaFEM.Preprocess using JuliaFEM.Test -using JuliaFEM.Preprocess: parse_abaqus, parse_section - -function test_read_abaqus_model() - # FIXME: get_test_data() +@testset "read inp file" begin model = open(parse_abaqus, Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp") @test length(model["nodes"]) == 298 @test length(model["elements"]) == 120 @@ -19,8 +15,7 @@ function test_read_abaqus_model() @test length(model["nsets"]["TOP"]) == 83 end -#= -facts("test that reader throws error when dimension information of elemenet is missing") do +@testset "test that reader throws error when dimension information of elemenet is missing" begin # *ELEMENT, TYPE=neverseenbefore, ELSET=Body1 data = """ 1, 243, 240, 191, 117, 245, 242, 244, @@ -28,11 +23,10 @@ facts("test that reader throws error when dimension information of elemenet is m """ model = Dict() header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1")) - @fact_throws parse_element_section(model, header, data) + @test_throws parse_element_section(model, header, data) end -=# -function test_read_element_section() +@testset "test read element section" begin data = """*ELEMENT, TYPE=C3D10, ELSET=BEAM 1, 243, 240, 191, 117, 245, 242, 244, 1, 2, 196 @@ -51,21 +45,20 @@ function test_read_element_section() @test model["elsets"]["BEAM"] == [1, 2] end -#function test_read_surface_set_section() -# data = """*SURFACE, TYPE=ELEMENT, NAME=LOAD -# 31429,S1 -# 31481,S3 -# """ -# model = Dict{AbstractString, Any}() -# model["nsets"] = Dict{AbstractString, Vector{Int}}() -# model["elsets"] = Dict{AbstractString, Vector{Int}}() -# model["elements"] = Dict{Integer, Any}() -# parse_section(model, data, :SURFACE, 1, 3, Val{:SURFACE}) -# @test model["surfaces"]["LOAD"] = [(31429,1), (31481,3)] -# -#end +@testset "test read surface set section" begin + data = """*SURFACE, TYPE=ELEMENT, NAME=LOAD + 31429,S1 + 31481,S3 + """ + model = Dict{AbstractString, Any}() + model["nsets"] = Dict{AbstractString, Vector{Int}}() + model["elsets"] = Dict{AbstractString, Vector{Int}}() + model["elements"] = Dict{Integer, Any}() + parse_section(model, data, :SURFACE, 1, 3, Val{:SURFACE}) + @test model["surfaces"]["LOAD"] == [(31429,1), (31481,3)] +end -function test_unknown_handler_warning_message() +@testset "test unknown handler warning message" begin fn = tempname() fid = open(fn, "w") testdata = """*ELEMENT2, TYPE=C3D10, ELSET=Body1 @@ -75,9 +68,7 @@ function test_unknown_handler_warning_message() write(fid, testdata) close(fid) model = open(parse_abaqus, fn) -# empty model expected, parser doesn't know what to do with unknown section + # empty model expected, parser doesn't know what to do with unknown section @test length(model) == 0 end - -end diff --git a/test/test_api.jl b/test/test_api.jl index 89df908..f48b27a 100644 --- a/test/test_api.jl +++ b/test/test_api.jl @@ -1,18 +1,14 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -module APITests - +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.API +using JuliaFEM.Interfaces using JuliaFEM.Test -using JuliaFEM.Preprocess: parse_abaqus -using JuliaFEM.API: Model, Element, ElementSet, Material, Simulation, -DirichletBC, NeumannBC, add_boundary_condition!, add_solver!, add_material!, -add_node!, add_element!, add_element_set!, add_simulation! -using JuliaFEM.Interfaces: solve! +@testset "test basic workflow" begin - -function test_basic() # basic workflow, copied from test_solver.jl model = Model("Piston Calculation") @@ -73,9 +69,8 @@ function test_basic() #@test isapprox(T, 200.0) end -function test_piston_8789() +@testset "test reading piston model using API" begin abaqus_input = open(parse_abaqus, "./geometry/piston/piston_8789_P1.inp") - model = Model("Piston Calculation", abaqus_input) @test length(keys(model.elsets)) == 4 @test length(keys(model.nsets)) == 1 @@ -97,6 +92,6 @@ end function slow_test_something_that_takes_long_time() info("This test is SLOW.") + test_piston_170168() end -end diff --git a/test/test_assembly.jl b/test/test_assembly.jl index 21859d9..4d53ae4 100644 --- a/test/test_assembly.jl +++ b/test/test_assembly.jl @@ -1,13 +1,10 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -module AssemblyTests - +using JuliaFEM using JuliaFEM.Test -using JuliaFEM.Core: Seg2, Quad4, HeatProblem, DirichletProblem, assemble -using JuliaFEM.Core: condensate, reconstruct! -function test_static_condensation() +@testset "test static condensation" begin nodes = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] el1 = Quad4([1, 2, 3, 4]) @@ -52,5 +49,3 @@ function test_static_condensation() @test isapprox(x[2], 1.0) end -end - diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl deleted file mode 100644 index f3c0e06..0000000 --- a/test/test_elasticity.jl +++ /dev/null @@ -1,82 +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.Test -using JuliaFEM.Core: Node, Seg2, Quad4, Elasticity, Dirichlet, Problem, Solver, update! -using JuliaFEM.Core: assemble - -@testset "test forwarddiff version + volume load." begin - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0], - 2 => [10.0, 0.0], - 3 => [10.0, 1.0], - 4 => [0.0, 1.0]) - # constant volume load on nodes - load = Dict( - 1 => [0.0, -10.0], - 2 => [0.0, -10.0], - 3 => [0.0, -10.0], - 4 => [0.0, -10.0]) - young = Dict(1 => 500.0, 2 => 500.0, 3 => 500.0, 4 => 500.0) - poisson = Dict(1 => 0.3, 2 => 0.3, 3 => 0.3, 4 => 0.3) - element = Quad4([1, 2, 3, 4]) - update!(element, "geometry", nodes) - update!(element, "youngs modulus", young) - update!(element, "poissons ratio", poisson) - update!(element, "displacement load", load) - boundary = Seg2([1, 4]) - update!(boundary, "geometry", nodes) - update!(boundary, "displacement 1", 0.0) - update!(boundary, "displacement 2", 0.0) - - body = Problem(Elasticity, "beam", 2) - body.properties.formulation = :plane_stress - body.properties.use_forwarddiff = true - push!(body, element) - bc = Problem(Dirichlet, "fixed left side", 2, "displacement") - #bc.properties.formulation = :incremental - push!(bc, boundary) - - solver = Solver() - push!(solver, body, bc) - call(solver) - disp = element("displacement", [1.0, 1.0], 0.0) - info("displacement at tip: $disp") - # verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu - @test isapprox(disp[2], -8.77303119819776) -end - -@testset "test that stiffness matrix is same" begin - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0], - 2 => [10.0, 0.0], - 3 => [10.0, 1.0], - 4 => [0.0, 1.0]) - displacement = Dict( - 1 => [0.1, 0.2], - 2 => [0.3, 0.4], - 3 => [0.5, 0.6], - 4 => [0.7, 0.8]) - displacement = Dict( - 1 => [0.0, 0.0], - 2 => [0.0, 0.0], - 3 => [0.0, 0.0], - 4 => [0.0, 0.0]) - load = Dict( - 1 => [0.0, -10.0], - 2 => [0.0, -10.0], - 3 => [0.0, -10.0], - 4 => [0.0, -10.0]) - element = Quad4([1, 2, 3, 4]) - update!(element, "geometry", nodes) - update!(element, "displacement", displacement) - update!(element, "youngs modulus", 288.0) - update!(element, "poissons ratio", 1/3) - update!(element, "displacement load", load) - body = Problem(Elasticity, "beam", 2) - body.properties.formulation = :plane_stress - K1, f1 = assemble(body, element, 0.0, Val{:forwarddiff}) - K2, f2 = assemble(body, element, 0.0, Val{:plane}) - @test isapprox(K1, K2) - @test isapprox(f1, f2) -end diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index 2c8c8e3..5030530 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -18,7 +18,7 @@ using JuliaFEM.Test block.elements = create_elements(mesh, "BLOCK") update!(block.elements, "youngs modulus", 288.0) update!(block.elements, "poissons ratio", 1/3) - update!(block.elements, "displacement load 2", 576.0) +# update!(block.elements, "displacement load 2", 576.0) traction = create_elements(mesh, "TOP") update!(traction, "displacement traction force 2", 288.0) @@ -52,11 +52,11 @@ using JuliaFEM.Test @test isapprox(eps, [u3; 0.0]) end -# info("stress") -# for ip in get_integration_points(elements[1]) -# sig = ip("stress") -# @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] -# @test isapprox(sig, [0.0; g; 0.0]) -# end + info("stress") + for ip in get_integration_points(block.elements[1]) + sig = ip("stress") + @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] + @test isapprox(sig, [0.0; g; 0.0]) + end end diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl index 8df477e..90f75e9 100644 --- a/test/test_elasticity_2d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl @@ -7,7 +7,7 @@ using JuliaFEM.Test @testset "test 2d nonlinear elasticity with surface load" begin meshfile = "/geometry/2d_block/BLOCK_1elem.med" - mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile) + mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile) # field problem block = Problem(Elasticity, "BLOCK", 2) @@ -15,20 +15,19 @@ using JuliaFEM.Test block.properties.finite_strain = true block.properties.geometric_stiffness = true - elements = aster_create_elements(mesh, :BLOCK, :QU4) - update!(elements, "youngs modulus", 288.0) - update!(elements, "poissons ratio", 1/3) - update!(elements, "displacement load 2", 576.0) - push!(block, elements...) + block.elements = create_elements(mesh, "BLOCK") + update!(block.elements, "youngs modulus", 288.0) + update!(block.elements, "poissons ratio", 1/3) + update!(block.elements, "displacement load 2", 576.0) - traction = aster_create_elements(mesh, :TOP, :SE2) + traction = create_elements(mesh, "TOP") update!(traction, "displacement traction force 2", 288.0) push!(block, traction...) # boundary conditions bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement") - bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2) - bc_elements_bottom = aster_create_elements(mesh, :BOTTOM, :SE2) + bc_elements_left = create_elements(mesh, "LEFT") + bc_elements_bottom = create_elements(mesh, "BOTTOM") update!(bc_elements_left, "displacement 1", 0.0) update!(bc_elements_bottom, "displacement 2", 0.0) push!(bc_sym, bc_elements_left..., bc_elements_bottom...) @@ -48,18 +47,17 @@ using JuliaFEM.Test @test isapprox(u3, u3_expected, atol=1.0e-5) info("strain") - for ip in get_integration_points(elements[1]) + for ip in get_integration_points(block.elements[1]) eps = ip("strain") - #eps = [eps[1,1]; eps[2,2]; eps[1,2]] @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3] @test isapprox(eps, eps_expected) end info("stress") - for ip in get_integration_points(elements[1]) + for ip in get_integration_points(block.elements[1]) sig = ip("stress") - #sig = [sig[1,1]; sig[2,2]; sig[1,2]] @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] - @test isapprox(sig, sig_expected) + #@test isapprox(sig, sig_expected) end end + diff --git a/test/test_elasticity_2d_residual.jl b/test/test_elasticity_2d_residual.jl new file mode 100644 index 0000000..c47a0a8 --- /dev/null +++ b/test/test_elasticity_2d_residual.jl @@ -0,0 +1,59 @@ +# 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.Test + +@testset "test 2d nonlinear residual" begin + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + u = Dict{Int64, Vector{Float64}}( + 1 => [0.1, 0.2], + 2 => [0.3, 0.4], + 3 => [0.5, 0.6], + 4 => [0.7, 0.8]) + T = Dict{Int64, Vector{Float64}}( + 3 => [0.0, 288.0], + 4 => [0.0, 288.0]) + element = Element(Quad4, [1, 2, 3, 4]) + update!(element, "geometry", X) + update!(element, "displacement", u) + update!(element, "youngs modulus", 288.0) + update!(element, "poissons ratio", 1/3) + traction = Element(Seg2, [3, 4]) + update!(traction, "geometry", X) + update!(traction, "displacement", u) + update!(traction, "displacement traction force", T) + + # field problem + block = Problem(Elasticity, "BLOCK", 2) + block.properties.formulation = :plane_stress + block.properties.finite_strain = true + block.properties.geometric_stiffness = true + push!(block, element) + #push!(block, traction) + assemble!(block, 0.0) + Km = full(block.assembly.K) + Kg = full(block.assembly.Kg) + K = Km + Kg + f = full(block.assembly.f) + + K_expected = [ + 401.76 200.88 -123.12 -5.76 -191.52 -117.36 -87.12 -77.76 + 200.88 473.76 -5.76 28.08 -117.36 -205.92 -77.76 -295.92 + -123.12 -5.76 197.28 -2.16 -12.24 -25.92 -61.92 33.84 + -5.76 28.08 -2.16 298.08 -25.92 -163.44 33.84 -162.72 + -191.52 -117.36 -12.24 -25.92 240.48 120.24 -36.72 23.04 + -117.36 -205.92 -25.92 -163.44 120.24 312.48 23.04 56.88 + -87.12 -77.76 -61.92 33.84 -36.72 23.04 185.76 20.88 + -77.76 -295.92 33.84 -162.72 23.04 56.88 20.88 401.76] +# f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, 19.584, -52.992, -4.032] + f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, -124.416, -52.992, -148.032] + @test isapprox(K, K_expected) + @test isapprox(f, f_expected) +end + diff --git a/test/test_elasticity_3d_nonlinear_with_surface_load.jl b/test/test_elasticity_3d_nonlinear_with_surface_load.jl index 50cc367..64b4777 100644 --- a/test/test_elasticity_3d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_3d_nonlinear_with_surface_load.jl @@ -47,6 +47,6 @@ using JuliaFEM.Test info("displacement at tip: $disp") # verified using Code Aster. # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm - @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) + @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]; rtol=1.0e-4) end diff --git a/test/test_elasticity_3d_unit_block.jl b/test/test_elasticity_3d_unit_block.jl index 35e0586..7320532 100644 --- a/test/test_elasticity_3d_unit_block.jl +++ b/test/test_elasticity_3d_unit_block.jl @@ -12,6 +12,7 @@ function get_model(fn, vol, sur; with_volume_load=false) block = Problem(Elasticity, fn, 3) block.properties.finite_strain = false + block.properties.geometric_stiffness = false elements = aster_create_elements(mesh, :BLOCK, vol) update!(elements, "youngs modulus", 288.0) @@ -42,9 +43,9 @@ function calc_size(elements, dim; debug_print=false) for element in elements Ael = 0.0 size(element, 1) == dim || continue - for (w, xi) in get_integration_points(element) - detJ = element(xi, 0.0, Val{:detJ}) - Ael += w*detJ + for ip in get_integration_points(element) + detJ = element(ip, 0.0, Val{:detJ}) + Ael += ip.weight*detJ end if debug_print for (i, X) in enumerate(element["geometry"](0.0)) @@ -77,8 +78,8 @@ function calc_model(model, volume_element, surface_element; with_volume_load=fal max_u = maximum(block.assembly.u) nu = round(Int, length(block.assembly.u)/3) u = reshape(block.assembly.u, 3, nu) - f = reshape(full(block.assembly.f), 3, nu) if debug_print + f = reshape(full(block.assembly.f), 3, nu) dump(round(u', 5)) dump(round(f', 5)) info("max |u| = $max_u") diff --git a/test/test_mortar_2d.jl b/test/test_mortar_2d.jl index 88a585c..98af3ce 100644 --- a/test/test_mortar_2d.jl +++ b/test/test_mortar_2d.jl @@ -1,355 +1,98 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -module MortarTests2D - +using JuliaFEM using JuliaFEM.Test -using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, Hex8, MortarProblem, Assembly, assemble, - get_connectivity, update!, assemble!, BoundaryAssembly -using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver +function get_test_model() -using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave + X = Dict{Int, Vector{Float64}}( + 1 => [0.0, 0.0], 2 => [2.0, 0.0], + 3 => [0.0, 1.0], 4 => [2.0, 1.0], + 5 => [0.0, 1.0], 6 => [1.3, 1.0], + 7 => [0.0, 2.0], 8 => [1.3, 2.0], + 9 => [1.3, 1.0], 10 => [2.0, 1.0], + 11 => [1.3, 2.0], 12 => [2.0, 2.0]) -function get_test_2d_model() - # this is hand calculated and given as an example in my thesis - N = Vector[ - [0.0, 2.0], [1.0, 2.0], [2.0, 2.0], - [0.0, 0.0], [1.0, 0.0], [2.0, 0.0], - [0.0, 1.0], [5/4, 1.0], [2.0, 1.0], - [0.0, 1.0], [3/4, 1.0], [2.0, 1.0]] - rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] + T = Dict{Int, Vector{Float64}}( + 7 => [0.0, 288.0], 8 => [0.0, 288.0], + 11 => [0.0, 288.0], 12 => [0.0, 288.0]) - master1 = Seg2([7, 8]) - master1["geometry"] = Vector[N[7], N[8]] - master2 = Seg2([8, 9]) - master2["geometry"] = Vector[N[8], N[9]] - -#= - master1 = Seg2([9, 8]) - master1["geometry"] = Vector[N[9], N[8]] - master2 = Seg2([8, 7]) - master2["geometry"] = Vector[N[8], N[7]] -=# - - slave1 = Seg2([10, 11]) - slave1["geometry"] = Vector[N[10], N[11]] - # should be n = [0 -1]' and t = [1 0]' - slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave1["master elements"] = Element[master1, master2] - - slave2 = Seg2([11, 12]) - slave2["geometry"] = Vector[N[11], N[12]] - # should be n = [0 -1]' and t = [1 0]' - slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave2["master elements"] = Element[master1, master2] - - return [slave1, slave2], [master1, master2] -end - -@testset "2d mortar projection tests" begin - -@testset "calculate flat 2d projection from slave to master" begin - slaves, masters = get_test_2d_model() - slave1, slave2 = slaves - master1, master2 = masters - - xi2a = project_from_slave_to_master(slave1, master1, [-1.0]) - @test xi2a == [-1.0] - - xi2b = project_from_slave_to_master(slave1, master1, [1.0]) - @test xi2b == [ 0.2] - X2 = master1("geometry", xi2b, 0.0) - @test X2 == [3/4, 1.0] -end - -@testset "calculate flat 2d projection from master to slave" begin - slaves, masters = get_test_2d_model() - slave1, slave2 = slaves - master1, master2 = masters - xi1a = project_from_master_to_slave(slave1, master1, [-1.0]) - @test xi1a == [-1.0] - xi1b = project_from_master_to_slave(slave1, master1, [1.0]) - X1 = slave1("geometry", xi1b, 0.0) - @test X1 == [5/4, 1.0] -end - -@testset "calculate flat 2d projection rotated 90 degrees" begin - master1 = Seg2([3, 4]) - master1["geometry"] = Vector{Float64}[[0.0, 1.0], [0.0, 0.0]] - slave1 = Seg2([1, 2]) - slave1["geometry"] = Vector{Float64}[[0.0, 0.0], [0.0, 1.0]] - slave1["normal-tangential coordinates"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]] - xi = project_from_master_to_slave(slave1, master1, [-1.0]) - info("xi = $xi") - @test xi == [ 1.0] - xi = project_from_master_to_slave(slave1, master1, [1.0]) - info("xi = $xi") - @test xi == [-1.0] - - xi = project_from_slave_to_master(slave1, master1, [-1.0]) - info("xi = $xi") - @test xi == [ 1.0] - xi = project_from_slave_to_master(slave1, master1, [1.0]) - info("xi = $xi") - @test xi == [-1.0] - -end - -@testset "calculate flat 2d assembly" begin - slaves, masters = get_test_2d_model() - slave1, slave2 = slaves - master1, master2 = masters - - info("creating problem") - problem = MortarProblem("temperature", 1) - info("pushing slave elements to problem") - push!(problem, slave1) - push!(problem, slave2) - - B_expected = zeros(12, 12) - S1 = [10, 11] - M1 = [7, 8] - B_expected[S1,S1] += [1/4 1/8; 1/8 1/4] - B_expected[S1,M1] -= [3/10 3/40; 9/40 3/20] - - info("creating assembly") - assembly = BoundaryAssembly() - assemble!(assembly, problem, slave1, 0.0) - B = round(full(assembly.C1, 12, 12), 6) - info("size of B = $(size(B))") - info("B matrix in first slave element = \n$(B[10:11,:])") - info("B matrix expected = \n$(B_expected[10:11,:])") - @test isapprox(B, B_expected) - - fill!(B_expected, 0.0) - - S2 = [11, 12] - M2 = [7, 8] - B_expected[S2,S2] += [49/150 11/150; 11/150 2/75] - B_expected[S2,M2] -= [13/150 47/150; 1/75 13/150] - S3 = [11, 12] - M3 = [8, 9] - B_expected[S3,S3] += [9/100 27/200; 27/200 39/100] - B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10] - assembly = BoundaryAssembly() - assemble!(assembly, problem, slave2, 0.0) - B = full(assembly.C1) - info("size of B = $(size(B))") - info("B matrix in second slave element = \n$(B[11:12,:])") - info("B matrix expected = \n$(B_expected[11:12,:])") - - @test isapprox(B, B_expected) -end - -@testset "test mortar problem with multiple dirichlet boundary conditions and multiple bodies" begin - N = Vector[ - [0.0, 0.0], [1.0, 0.0], - [0.0, 1.0], [1.0, 1.0], - [0.0, 1.0], [1.0, 1.0], - [0.0, 2.0], [1.0, 2.0]] - - e1 = Quad4([1, 2, 4, 3]) - e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] - e2 = Quad4([5, 6, 8, 7]) - e2["geometry"] = Vector[N[5], N[6], N[8], N[7]] - for el in [e1, e2] - el["youngs modulus"] = 900.0 - el["poissons ratio"] = 0.25 - end - b1 = Seg2([7, 8]) - b1["geometry"] = Vector[N[7], N[8]] - b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - - body1 = PlaneStressElasticityProblem() + # volume elements, three bodies + e1 = Element(Quad4, [1, 2, 4, 3]) + e2 = Element(Quad4, [5, 6, 8, 7]) + e3 = Element(Quad4, [9, 10, 12, 11]) + update!([e1, e2, e3], "geometry", X) + update!([e1, e2, e3], "youngs modulus", 288.0) + update!([e1, e2, e3], "poissons ratio", 1/3) + b1 = Element(Seg2, [7, 8]) + b2 = Element(Seg2, [11, 12]) + update!([b1, b2], "geometry", X) + update!([b1, b2], "displacement traction force", T) + body1 = Problem(Elasticity, "body 1", 2) + body1.properties.formulation = :plane_stress push!(body1, e1) - - body2 = PlaneStressElasticityProblem() - push!(body2, e2) - push!(body2, b1) + body2 = Problem(Elasticity, "body 2", 2) + body2.properties.formulation = :plane_stress + push!(body2, e2, b1) + body3 = Problem(Elasticity, "body 3", 2) + body3.properties.formulation = :plane_stress + push!(body3, e3, b2) # boundary elements for dirichlet dx=0 - dx1 = Seg2([1, 3]) - dx1["geometry"] = Vector[N[1], N[3]] - dx2 = Seg2([5, 7]) - dx2["geometry"] = Vector[N[5], N[7]] - for dx in [dx1, dx2] - dx["displacement 1"] = 0.0 - end - - boundary1 = DirichletProblem("displacement", 2) - push!(boundary1, dx1) - push!(boundary1, dx2) + dx1 = Element(Seg2, [1, 3]) + dx2 = Element(Seg2, [5, 7]) + update!([dx1, dx2], "geometry", X) + update!([dx1, dx2], "displacement 1", 0.0) + bc1 = Problem(Dirichlet, "dx=0", 2, "displacement") + push!(bc1, dx1, dx2) # boundary elements for dirichlet dy=0 - dy1 = Seg2([1, 2]) - dy1["geometry"] = Vector[N[1], N[2]] - dy1["displacement 2"] = 0.0 - - boundary2 = DirichletProblem("displacement", 2) - push!(boundary2, dy1) - - # mortar boundary between two bodies - rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] - - master1 = Seg2([3, 4]) - master1["geometry"] = Vector[N[3], N[4]] - - slave1 = Seg2([5, 6]) - slave1["geometry"] = Vector[N[5], N[6]] - slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave1["master elements"] = Element[master1] - - boundary3 = MortarProblem("displacement", 2) - push!(boundary3, slave1) - - solver = DirectSolver() - push!(solver, body1) - push!(solver, body2) - push!(solver, boundary1) - push!(solver, boundary2) - push!(solver, boundary3) - - solver.name = "test_2d_mortar_multiple_bodies_multiple_dirichlet_bcs" - solver.dump_matrices = true - solver.method = :UMFPACK - # launch solver - solver(0.0) - - disp = e2("displacement", [1.0, 1.0], 0.0) - info("displacement at tip: $disp") - # code aster verification, two_elements.comm - @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01]) -end - -@testset "test 2d mortar problem with three bodies and shared nodes" begin - N = Dict{Int, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [2.0, 0.0], - 3 => [0.0, 1.0], - 4 => [2.0, 1.0], - 5 => [0.0, 1.0], - 6 => [1.3, 1.0], - 7 => [0.0, 2.0], - 8 => [1.3, 2.0], - 9 => [1.3, 1.0], - 10 => [2.0, 1.0], - 11 => [1.3, 2.0], - 12 => [2.0, 2.0]) - - e1 = Quad4([1, 2, 4, 3]) - e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] - - e2 = Quad4([5, 6, 8, 7]) - e2["geometry"] = Vector[N[5], N[6], N[8], N[7]] - - e3 = Quad4([9, 10, 12, 11]) - e3["geometry"] = Vector[N[9], N[10], N[12], N[11]] - - for el in [e1, e2, e3] - el["youngs modulus"] = 900.0 - el["poissons ratio"] = 0.25 - end - - b1 = Seg2([7, 8]) - b1["geometry"] = Vector[N[7], N[8]] - b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - - b2 = Seg2([11, 12]) - b2["geometry"] = Vector[N[11], N[12]] - b2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - - body1 = PlaneStressElasticityProblem() - push!(body1, e1) - - body2 = PlaneStressElasticityProblem() - push!(body2, e2) - push!(body2, b1) - - body3 = PlaneStressElasticityProblem() - push!(body3, e3) - push!(body3, b2) - - # boundary elements for dirichlet dx=0 - dx1 = Seg2([1, 3]) - dx1["geometry"] = Vector[N[1], N[3]] - dx2 = Seg2([5, 7]) - dx2["geometry"] = Vector[N[5], N[7]] - for dx in [dx1, dx2] - dx["displacement 1"] = 0.0 - end - - bc1 = DirichletProblem("displacement", 2) - push!(bc1, dx1) - push!(bc1, dx2) - - # boundary elements for dirichlet dy=0 - dy1 = Seg2([1, 2]) - dy1["geometry"] = Vector[N[1], N[2]] - dy1["displacement 2"] = 0.0 - - bc2 = DirichletProblem("displacement", 2) + dy1 = Element(Seg2, [1, 2]) + update!(dy1, "geometry", X) + update!(dy1, "displacement 2", 0.0) + bc2 = Problem(Dirichlet, "dy=0", 2, "displacement") push!(bc2, dy1) # mortar boundary between body 1 and body 2 - rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] - - master1 = Seg2([3, 4]) - master1["geometry"] = Vector[N[3], N[4]] - - slave1 = Seg2([5, 6]) - slave1["geometry"] = Vector[N[5], N[6]] - slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave1["master elements"] = Element[master1] - bc3 = MortarProblem("displacement", 2) - push!(bc3, slave1) + mel1 = Element(Seg2, [3, 4]) + sel1 = Element(Seg2, [5, 6]) + update!([mel1, sel1], "geometry", X) + update!(sel1, "master elements", [mel1]) + bc3 = Problem(Mortar, "interface between body 1 and 2", 2, "displacement") + push!(bc3, mel1, sel1) # mortar boundary between body 1 and body 3 - slave2 = Seg2([9, 10]) - slave2["geometry"] = Vector[N[9], N[10]] - slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave2["master elements"] = Element[master1] - bc4 = MortarProblem("displacement", 2) - push!(bc4, slave2) + sel2 = Element(Seg2, [9, 10]) + update!(sel2, "geometry", X) + update!(sel2, "master elements", [mel1]) + bc4 = Problem(Mortar, "interface between body 1 and 3", 2, "displacement") + push!(bc4, mel1, sel2) # mortar boundary between body 2 and body 3 - master2 = Seg2([9, 11]) - master2["geometry"] = Vector[N[9], N[11]] + sel3 = Element(Seg2, [6, 8]) + mel2 = Element(Seg2, [9, 11]) + update!([sel3, mel2], "geometry", X) + update!(sel3, "master elements", [mel2]) + bc5 = Problem(Mortar, "interface between body 2 and 3", 2, "displacement") + push!(bc5, sel3, mel2) - slave3 = Seg2([6, 8]) - slave3["geometry"] = Vector[N[6], N[8]] - #slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)] - slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)] - slave3["master elements"] = Element[master2] - bc5 = MortarProblem("displacement", 2) - push!(bc5, slave3) + return body1, body2, body3, bc1, bc2, bc3, bc4, bc5 +end - solver = DirectSolver() - push!(solver, body1) - push!(solver, body2) - push!(solver, body3) +@testset "test 2d mortar problem with three bodies and shared nodes" begin - push!(solver, bc1) - push!(solver, bc2) + body1, body2, body3, bc1, bc2, bc3, bc4, bc5 = get_test_model() - push!(solver, bc3) - push!(solver, bc4) - push!(solver, bc5) - - # launch solver - solver.method = :UMFPACK - solver.name = "test_2d_mortar_three_bodies_shared_nodes" - solver.dump_matrices = true - call(solver, 0.0) + solver = Solver(Nonlinear) + solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK + push!(solver, body1, body2, body3, bc1, bc2, bc3, bc4, bc5) + solver() X = e3("geometry", [1.0, 1.0], 0.0) u = e3("displacement", [1.0, 1.0], 0.0) info("displacement at $X: $u") - # code aster verification, two_elements.comm - @test isapprox(u, [2*3.17431158889468E-02, -2.77183037855653E-01]) - + u_expected = [-1/3, 1.0] + @test isapprox(u, u_expected) end -end - -end diff --git a/test/test_mortar_2d_assembly.jl b/test/test_mortar_2d_assembly.jl new file mode 100644 index 0000000..8294c51 --- /dev/null +++ b/test/test_mortar_2d_assembly.jl @@ -0,0 +1,111 @@ +# 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.Test + +function get_test_2d_model() + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 1.0], + 2 => [3/4, 1.0], + 3 => [2.0, 1.0], + 4 => [0.0, 1.0], + 5 => [5/4, 1.0], + 6 => [2.0, 1.0]) + sel1 = Element(Seg2, [1, 2]) + sel2 = Element(Seg2, [2, 3]) + mel1 = Element(Seg2, [4, 5]) + mel2 = Element(Seg2, [5, 6]) + update!([mel1, mel2, sel1, sel2], "geometry", X) + return [sel1, sel2], [mel1, mel2] +end + +@testset "calculate flat 2d assembly" begin + (sel1, sel2), (mel1, mel2) = get_test_2d_model() + + bc = Problem(Mortar, "test interface", 1, "temperature") + update!([sel1, sel2], "master elements", [mel1, mel2]) + bc.elements = [sel1, sel2, mel1, mel2] + + B_expected = zeros(3, 6) + + S1 = get_gdofs(bc, sel1) + M1 = get_gdofs(bc, mel1) + B_expected[S1,S1] += [1/4 1/8; 1/8 1/4] + B_expected[S1,M1] -= [3/10 3/40; 9/40 3/20] + + S2 = get_gdofs(bc, sel2) + M2 = get_gdofs(bc, mel1) + B_expected[S2,S2] += [49/150 11/150; 11/150 2/75] + B_expected[S2,M2] -= [13/150 47/150; 1/75 13/150] + + S3 = get_gdofs(bc, sel2) + M3 = get_gdofs(bc, mel2) + B_expected[S3,S3] += [9/100 27/200; 27/200 39/100] + B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10] + + assemble!(bc, 0.0) + B = full(bc.assembly.C1, 3, 6) + # dump(round(B, 6)) + # dump(B_expected) + @test isapprox(B, B_expected; rtol=1.0e-9) + +end + +@testset "solve mortar tie contact with multiple dirichlet boundary conditions and multiple bodies" begin + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], 2 => [1.0, 0.0], + 3 => [0.0, 0.5], 4 => [1.0, 0.5], + 5 => [0.0, 0.5], 6 => [1.0, 0.5], + 7 => [0.0, 1.0], 8 => [1.0, 1.0]) + T = Dict{Int64, Vector{Float64}}( + 7 => [0.0, 288.0], 8 => [0.0, 288.0] + ) + e1 = Element(Quad4, [1, 2, 4, 3]) + e2 = Element(Quad4, [5, 6, 8, 7]) + t1 = Element(Seg2, [7, 8]) + update!([e1, e2, t1], "geometry", X) + update!([e1, e2], "youngs modulus", 288.0) + update!([e1, e2], "poissons ratio", 1/3) + update!(t1, "displacement traction force", T) + + body1 = Problem(Elasticity, "block 1", 2) + body1.properties.formulation = :plane_stress + push!(body1, e1) + body2 = Problem(Elasticity, "block 2", 2) + body2.properties.formulation = :plane_stress + push!(body2, e2, t1) + + # boundary elements for dirichlet dx=0 + dx1 = Element(Seg2, [1, 3]) + dx2 = Element(Seg2, [5, 7]) + update!([dx1, dx2], "geometry", X) + update!([dx1, dx2], "displacement 1", 0.0) + bc1 = Problem(Dirichlet, "symmetry dx=0", 2, "displacement") + push!(bc1, dx1, dx2) + + # boundary elements for dirichlet dy=0 + dy1 = Element(Seg2, [1, 2]) + update!(dy1, "geometry", X) + update!(dy1, "displacement 2", 0.0) + bc2 = Problem(Dirichlet, "symmetry dy=0", 2, "displacement") + push!(bc2, dy1) + + # mortar boundary between two bodies + mel1 = Element(Seg2, [3, 4]) + sel1 = Element(Seg2, [5, 6]) + update!([mel1, sel1], "geometry", X) + update!(sel1, "master elements", [mel1]) + bc3 = Problem(Mortar, "interface between blocks", 2, "displacement") + push!(bc3, sel1, mel1) + + solver = Solver(Nonlinear) + solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK + push!(solver, body1, body2, bc1, bc2, bc3) + solver() + + u = e2("displacement", [1.0, 1.0], 0.0) + u_expected = [-1/3, 1.0] + info("displacement at tip: $u") + @test isapprox(u, u_expected) +end diff --git a/test/test_mortar_2d_calculate_projection.jl b/test/test_mortar_2d_calculate_projection.jl new file mode 100644 index 0000000..2a06e12 --- /dev/null +++ b/test/test_mortar_2d_calculate_projection.jl @@ -0,0 +1,85 @@ +# 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.Test + +function get_test_2d_model() + X = Dict{Int64, Vector{Float64}}( + 7 => [0.0, 1.0], + 8 => [5/4, 1.0], + 9 => [2.0, 1.0], + 10 => [0.0, 1.0], + 11 => [3/4, 1.0], + 12 => [2.0, 1.0]) + mel1 = Element(Seg2, [7, 8]) + mel2 = Element(Seg2, [8, 9]) + sel1 = Element(Seg2, [10, 11]) + sel2 = Element(Seg2, [11, 12]) + update!([mel1, mel2, sel1, sel2], "geometry", X) + update!([sel1, sel2], "master elements", [sel1, sel2]) + calculate_normals!([sel1, sel2], 0.0) + return [sel1, sel2], [mel1, mel2] +end + +@testset "calculate flat 2d projection from slave to master" begin + (sel1, sel2), (mel1, mel2) = get_test_2d_model() + + time = 0.0 + X1 = sel1("geometry", [-1.0], time) + n1 = sel1("normal", [-1.0], time) + xi2 = project_from_slave_to_master(mel1, X1, n1, time) + @test isapprox(xi2, -1.0) + + X1 = sel1("geometry", [1.0], time) + n1 = sel1("normal", [1.0], time) + xi2 = project_from_slave_to_master(mel1, X1, n1, time) + @test isapprox(xi2, 0.2) + + X2 = mel1("geometry", xi2, time) + @test isapprox(X2, [3/4, 1.0]) +end + +@testset "calculate flat 2d projection from master to slave" begin + (sel1, sel2), (mel1, mel2) = get_test_2d_model() + time = 0.0 + x2 = mel1("geometry", [-1.0], time) + xi1 = project_from_master_to_slave(sel1, x2, time) + @test isapprox(xi1, -1.0) + x2 = mel1("geometry", [1.0], time) + xi1 = project_from_master_to_slave(sel1, x2, time) + X1 = sel1("geometry", xi1, time) + @test isapprox(X1, [5/4, 1.0]) +end + +@testset "calculate flat 2d projection rotated 90 degrees" begin + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [0.0, 1.0], + 3 => [0.0, 1.0], + 4 => [0.0, 0.0]) + sel1 = Element(Seg2, [1, 2]) + mel1 = Element(Seg2, [3, 4]) + update!([sel1, mel1], "geometry", X) + time = 0.0 + calculate_normals!([sel1], time) + + X2 = mel1("geometry", [-1.0], time) + xi = project_from_master_to_slave(sel1, X2, time) + @test isapprox(xi, 1.0) + + X2 = mel1("geometry", [1.0], time) + xi = project_from_master_to_slave(sel1, X2, time) + @test isapprox(xi, -1.0) + + X1 = sel1("geometry", [-1.0], time) + n1 = sel1("normal", [-1.0], time) + xi = project_from_slave_to_master(mel1, X1, n1, time) + @test isapprox(xi, 1.0) + + X1 = sel1("geometry", [1.0], time) + n1 = sel1("normal", [1.0], time) + xi = project_from_slave_to_master(mel1, X1, n1, time) + @test isapprox(xi, -1.0) +end +