diff --git a/docs/make.jl b/docs/make.jl index 8a256e9..92fcbc2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,10 +3,17 @@ using Documenter, Mortar2D +if haskey(ENV, "TRAVIS") + println("inside TRAVIS, installing PyPlot + matplotlib") + Pkg.add("PyPlot") + run(`pip install matplotlib`) +end + makedocs(modules=[Mortar2D], format = :html, sitename = "Mortar2D", pages = [ "Introduction" => "index.md", - "API" => "api.md", + "Theory" => "theory.md", + "API" => "api.md" ]) diff --git a/docs/plots.jl b/docs/plots.jl new file mode 100644 index 0000000..df13af4 --- /dev/null +++ b/docs/plots.jl @@ -0,0 +1,99 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE + +using PyPlot +using Mortar2D: calculate_normals, project_from_master_to_slave, project_from_slave_to_master, calculate_segments, calculate_mortar_matrices, calculate_mortar_assembly + +coords = Dict(1 => [8.0, 10.0], + 2 => [7.0, 7.0], + 3 => [4.0, 3.0], + 4 => [0.0, 0.0], + 5 => [-3.0, 0.0], + 6 => [12.0, 10.0], + 7 => [10.0, 4.0], + 8 => [7.0, 2.0], + 9 => [4.0, -2.0], + 10 => [0.0, -3.0], + 11 => [-4.0, -3.0]) + +elements = Dict( + 1 => [1, 2], + 2 => [2, 3], + 3 => [3, 4], + 4 => [4, 5], + 5 => [6, 7], + 6 => [7, 8], + 7 => [8, 9], + 8 => [9, 10], + 9 => [10, 11]) + +slave_element_ids = [1, 2, 3, 4] +slave_elements = Dict(i => elements[i] for i in slave_element_ids) +master_element_ids = [5, 6, 7, 8, 9] +element_types = Dict(i => :Seg2 for i=1:length(elements)) +normals = calculate_normals(slave_elements, element_types, coords) +segmentation = calculate_segments(slave_element_ids, master_element_ids, + elements, element_types, coords, normals) + +function plot1(; plot_element_normals=false, plot_nodal_normals=false, plot_segmentation=false) + + figure(figsize=(4, 4)) + + for i in slave_element_ids + x1,y1 = coords[elements[i][1]] + x2,y2 = coords[elements[i][2]] + xmid = 1/2*(x1+x2) + ymid = 1/2*(y1+y2) + n = [y1-y2, x2-x1] + n /= norm(n) + plot([x1,x2], [y1,y2], "-bo") + if plot_element_normals + arrow(x1, y1, n[1], n[2], head_width=0.1, head_length=0.2, fc="b", ec="b") + arrow(x2, y2, n[1], n[2], head_width=0.1, head_length=0.2, fc="b", ec="b") + end + end + + for i in master_element_ids + x1,y1 = coords[elements[i][1]] + x2,y2 = coords[elements[i][2]] + plot([x1,x2], [y1,y2], "-ro") + end + + if plot_nodal_normals + for i in keys(normals) + x = coords[i] + n = normals[i] + arrow(x[1], x[2], n[1], n[2], head_width=0.1, head_length=0.2, fc="k", ec="k") + end + end + + if plot_segmentation + for (sid, seg) in segmentation + scon = elements[sid] + xs1 = coords[scon[1]] + xs2 = coords[scon[2]] + ns1 = normals[scon[1]] + ns2 = normals[scon[2]] + for (mid, xi) in seg + mcon = elements[mid] + xm1 = coords[mcon[1]] + xm2 = coords[mcon[2]] + for xi_s in xi + #xi_s = (1-xi)/2*p + (1+xi)/2*p + N1 = [(1-xi_s)/2 (1+xi_s)/2] + n_s = N1[1]*ns1 + N1[2]*ns2 + x_g = N1[1]*xs1 + N1[2]*xs2 + xi_m = project_from_slave_to_master(Val{:Seg2}, x_g, n_s, xm1, xm2) + N2 = [(1-xi_m)/2 (1+xi_m)/2] + x_m = N2[1]*xm1 + N2[2]*xm2 + plot([x_g[1], x_m[1]], [x_g[2], x_m[2]], "--k") + end + end + end + end + + axis("equal") + annotate("\$\\Gamma_{\\mathrm{c}}^{\\left(1\\right)}\$", xy=(2.5, 5.0)) + annotate("\$\\Gamma_{\\mathrm{c}}^{\\left(2\\right)}\$", xy=(8.0, 0.0)) + axis("off") +end diff --git a/docs/src/api.md b/docs/src/api.md index 438bf1f..2de4d8b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,5 +1,10 @@ # API documentation +## Index + +```@index +``` + ```@meta DocTestSetup = quote using Mortar2D @@ -16,8 +21,3 @@ Mortar2D.calculate_mortar_matrices Mortar2D.calculate_mortar_assembly ``` -## Index - -```@index -``` - diff --git a/docs/src/figs/poisson_problem.png b/docs/src/figs/poisson_problem.png new file mode 100644 index 0000000..3c22101 Binary files /dev/null and b/docs/src/figs/poisson_problem.png differ diff --git a/docs/src/index.md b/docs/src/index.md index a92619e..4beabcd 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,7 +2,7 @@ ![Typical 2d segmentation](figs/contact_segmentation.png) ```@contents -Pages = ["index.md", "api.md"] +Pages = ["index.md", "theory.md", "api.md"] ``` Mortar2D.jl is a julia package to calculate discrete projections between diff --git a/docs/src/theory.md b/docs/src/theory.md new file mode 100644 index 0000000..5df1d8e --- /dev/null +++ b/docs/src/theory.md @@ -0,0 +1,167 @@ +# Theory + +Let us consider the simple following simple Laplace problem in the +domain ``\Omega=\left\{ \left(x,y\right)\in\mathbb{R}^{2}|0\leq x,y\leq2\right\}``. + +![](figs/poisson_problem.png) + +The strong form of the problem is +```math +\begin{align} +-\Delta u^{\left(i\right)} & =0\qquad\text{in }\Omega^{\left(i\right)},i=1,2\\ +u^{\left(1\right)} & =0\qquad\text{on }\Gamma_{\mathrm{D}}^{\left(1\right)}\\ +u^{\left(1\right)}-u^{\left(2\right)} & =0\qquad\text{on }\Gamma_{\mathrm{C}}\\ +\frac{\partial u^{\left(2\right)}}{\partial n} & =g\qquad\text{on }\Gamma_{\mathrm{N}}^{\left(2\right)} +\end{align} +``` + +Corresponding weak form of the problem is to find ``u^{\left(i\right)}\in\mathcal{U}$ +and $\lambda\in\mathcal{M}`` such that +```math +\begin{align} +\int_{\Omega}\nabla u\cdot\nabla v\,\mathrm{d}x+\int_{\Gamma_{\mathrm{C}}}\lambda\left(v^{\left(1\right)}-v^{\left(2\right)}\right)\,\mathrm{d}s & =\int_{\Omega}fv\,\mathrm{d}x+\int_{\Gamma_{\mathrm{N}}}gv\,\mathrm{d}s & & \forall v^{\left(i\right)}\in\mathcal{V}^{\left(i\right)}\\ +\int_{\Gamma_{\mathrm{C}}}\mu\left(u^{\left(1\right)}-u^{\left(2\right)}\right)\,\mathrm{d}s & =0 & & \forall\mu\in\mathcal{M} +\end{align} +``` + +In more general form is to find ``u^{\left(i\right)}\in\mathcal{U}`` +and ``\lambda\in\mathcal{M}`` such that +```math +\begin{align} +a\left(u^{\left(i\right)},v^{\left(i\right)}\right)+b\left(\lambda,v^{\left(i\right)}\right) & =0\qquad\forall v^{\left(i\right)}\in\mathcal{V}^{\left(i\right)}\\ +b\left(\mu,u^{\left(i\right)}\right) & =0\qquad\forall\mu\in\mathcal{M} +\end{align}, +``` +where +```math +\begin{align} +b\left(\lambda,v^{\left(i\right)}\right) & =\int_{\Gamma_{\mathrm{C}}}\lambda\left(v^{\left(1\right)}-v^{\left(2\right)}\right)\,\mathrm{d}s\\ +b\left(\mu,u^{\left(i\right)}\right) & =\int_{\Gamma_{\mathrm{C}}}\mu\left(u^{\left(1\right)}-u^{\left(2\right)}\right)\,\mathrm{d}s +\end{align} +``` + +After substituting interpolation polynomials to weak form we get so called mortar matrices ``\boldsymbol{D}`` and ``\boldsymbol{M}``: +```math +\begin{align} +\boldsymbol{D}\left[j,k\right] & =\int_{\Gamma_{\mathrm{c}}^{\left(1\right)}}N_{j}N_{k}^{\left(1\right)}\,\mathrm{d}s,\\ +\boldsymbol{M}\left[j,l\right] & =\int_{\Gamma_{\mathrm{c}}^{\left(1\right)}}N_{j}\left(N_{l}^{\left(2\right)}\circ\chi\right)\,\mathrm{d}s, +\end{align} +``` +where ``\chi`` is mapping between contacting surfaces. Let us define some contact pair: + +```@example 0 +coords = Dict(1 => [8.0, 10.0], + 2 => [7.0, 7.0], + 3 => [4.0, 3.0], + 4 => [0.0, 0.0], + 5 => [-3.0, 0.0], + 6 => [12.0, 10.0], + 7 => [10.0, 4.0], + 8 => [7.0, 2.0], + 9 => [4.0, -2.0], + 10 => [0.0, -3.0], + 11 => [-4.0, -3.0]) + +elements = Dict( + 1 => [1, 2], + 2 => [2, 3], + 3 => [3, 4], + 4 => [4, 5], + 5 => [6, 7], + 6 => [7, 8], + 7 => [8, 9], + 8 => [9, 10], + 9 => [10, 11]) + +slave_element_ids = [1, 2, 3, 4] + +slave_elements = Dict(i => elements[i] for i in slave_element_ids) + +master_element_ids = [5, 6, 7, 8, 9] + +element_types = Dict(i => :Seg2 for i=1:length(elements)); +``` + +```@example 0 +using PyPlot # hide +include("plots.jl") # hide +plot1(plot_element_normals=true) # hide +savefig("fig1.svg"); nothing # hide +``` +![](fig1.svg) + +For first order elements, normal direction is not unique. For that reason some preprocessing needs to be done to calculate unique nodal normals. + +Unique nodal normals can be calculated several different ways, more or less sophisticated. An easy solution is just to take average of the normals of adjacing elements connecting to node ``k``, i.e. +```math +\begin{equation} +\boldsymbol{n}_{k}=\frac{\sum_{e=1}^{n_{k}^{\mathrm{adj}}}\boldsymbol{n}_{k}^{\left(e\right)}}{\left\Vert \sum_{e=1}^{n_{k}^{\mathrm{adj}}}\boldsymbol{n}_{k}^{\left(e\right)}\right\Vert }, +\end{equation} +``` +where ``\boldsymbol{n}_{k}^{\left(e\right)}`` means the normal calculated +in element ``e`` in node ``k``, and adj means adjacing elements. + +This is implemented in function `calculate_normals`: + +```@example 0 +plot1(;plot_nodal_normals=true) # hide +savefig("fig2.svg"); nothing # hide +normals = calculate_normals(slave_elements, element_types, coords) +``` +![](fig2.svg) + +This package follows the idea of continuous normal field, proposed by Yang et al., where all the quantities are projected using only slave side normals. +If we wish to find the projection of a slave node ``\boldsymbol{x}_{\mathrm{s}}``, +having normal vector ``\boldsymbol{n}_{\mathrm{s}}`` onto a master +element with nodes ``\boldsymbol{x}_{\mathrm{m1}}`` and ``\boldsymbol{x}_{\mathrm{m2}}``, +we are solving ``\xi^{\left(2\right)}`` from the equation +```math +\begin{equation} +\left[N_{1}\left(\xi^{\left(2\right)}\right)\boldsymbol{x}_{\mathrm{m1}}+N_{2}\left(\xi^{\left(2\right)}\right)\boldsymbol{x}_{\mathrm{m2}}-\boldsymbol{x}_{\mathrm{s}}\right]\times\boldsymbol{n}_{\mathrm{s}}=\boldsymbol{0}. +\end{equation} +``` + +The equation to find the projection of a master node ``\boldsymbol{x}_{\mathrm{m}}`` +onto a slave element with nodes ``\boldsymbol{x}_{\mathrm{s1}}`` and +``\boldsymbol{x}_{\mathrm{s2}}`` and normals ``\boldsymbol{n}_{\mathrm{s1}}`` +and ``\boldsymbol{n}_{\mathrm{s1}}`` is +```math +\begin{equation} +\left[N_{1}\left(\xi^{\left(1\right)}\right)\boldsymbol{x}_{\mathrm{s1}}+N_{2}\left(\xi^{\left(1\right)}\right)\boldsymbol{x}_{\mathrm{s2}}-\boldsymbol{x}_{\mathrm{m}}\right]\times\left[N_{1}\left(\xi^{\left(1\right)}\right)\boldsymbol{n}_{s1}+N_{2}\left(\xi^{\left(1\right)}\right)\boldsymbol{n}_{\mathrm{s2}}\right]=\boldsymbol{0}, +\end{equation} +``` +where ``\xi^{\left(1\right)}`` is the unknown parameter. First equation +is linear and second is quadratic (in general). Second equation is +also linear if ``\boldsymbol{n}_{\mathrm{s1}}=\boldsymbol{n}_{\mathrm{s2}}``. + +These equations are solved in function `project_from_master_to_slave` and `project_from_slave_to_master`. They are used in function `calculate_segments`, which is used to calculate segmentation of interface. + +```@example 0 +plot1(;plot_segmentation=true) # hide +savefig("fig3.svg"); nothing # hide +segmentation = calculate_segments(slave_element_ids, master_element_ids, + elements, element_types, coords, normals) + +``` +![](fig3.svg) + +After segmentation is calculated, it's possible to integrate over +non-conforming surface to calculate mortar matrices ``\boldsymbol{D}`` +and ``\boldsymbol{M}`` or ``\boldsymbol{P}=\boldsymbol{D}^{-1}\boldsymbol{M}``. +Calculation projection matrix ``\boldsymbol{P}`` is implemented as function `calculate_mortar_assembly`: + +```@example 0 +P = calculate_mortar_assembly(elements, element_types, coords, + slave_element_ids, master_element_ids) +``` + +This last command combines everything above to single command to calculate +projection matrix needed for finite element codes. + +# References + +- Wikipedia contributors. "Mortar methods." Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia. +- Maday, Yvon, Cathy Mavriplis, and Anthony Patera. "Nonconforming mortar element methods: Application to spectral discretizations." (1988). +- Yang, Bin, Tod A. Laursen, and Xiaonong Meng. "Two dimensional mortar contact methods for large deformation frictional sliding." International journal for numerical methods in engineering 62.9 (2005): 1183-1225. +- Yang, Bin, and Tod A. Laursen. "A contact searching algorithm including bounding volume trees applied to finite sliding mortar formulations." Computational Mechanics 41.2 (2008): 189-205. +- Wohlmuth, Barbara I. "A mortar finite element method using dual spaces for the Lagrange multiplier." SIAM journal on numerical analysis 38.3 (2000): 989-1012. diff --git a/src/Mortar2D.jl b/src/Mortar2D.jl index 4265221..98c5ffa 100644 --- a/src/Mortar2D.jl +++ b/src/Mortar2D.jl @@ -6,5 +6,5 @@ include("calculate_normals.jl") include("calculate_projections.jl") include("calculate_segments.jl") include("calculate_mortar_matrices.jl") -include("calculate_assembly.jl") +include("calculate_mortar_assembly.jl") end diff --git a/src/calculate_assembly.jl b/src/calculate_mortar_assembly.jl similarity index 61% rename from src/calculate_assembly.jl rename to src/calculate_mortar_assembly.jl index f8ed3d6..029b1dd 100644 --- a/src/calculate_assembly.jl +++ b/src/calculate_mortar_assembly.jl @@ -1,9 +1,27 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE -""" Given data, calculate projection matrix `P`. """ -function calculate_mortar_assembly(elements, element_types, coords, - slave_element_ids, master_element_ids) +""" + calculate_mortar_assembly(elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + slave_element_ids::Vector{Int}, + master_element_ids::Vector{Int}) + +Given data, calculate projection matrix `P`. This is the main function of +package. + +Matrix ``P`` is defined as a projection between slave and master surface, +i.e. +```math + D u_s = M u_m \\Rightarrow u_s = D^{-1} M u_m = P u_m. +``` +""" +function calculate_mortar_assembly(elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + slave_element_ids::Vector{Int}, + master_element_ids::Vector{Int}) S = Int[] M = Int[] for sid in slave_element_ids @@ -16,7 +34,6 @@ function calculate_mortar_assembly(elements, element_types, coords, M = sort(unique(M)) ns = length(S) nm = length(M) - println("$ns slave nodes, $nm master nodes") normals = calculate_normals(elements, element_types, coords) segmentation = calculate_segments(slave_element_ids, master_element_ids, elements, element_types, coords, normals) @@ -28,7 +45,7 @@ function calculate_mortar_assembly(elements, element_types, coords, M_V = Float64[] for sid in slave_element_ids mids = [mid for (mid, xi) in segmentation[sid]] - De, Me = calculate_mortar_matrices(sid, mids, elements, element_types, + De, Me = calculate_mortar_matrices(sid, elements, element_types, coords, normals, segmentation) for (i,I) in enumerate(elements[sid]) for (j,J) in enumerate(elements[sid]) diff --git a/src/calculate_mortar_matrices.jl b/src/calculate_mortar_matrices.jl index 25c03e9..fb92413 100644 --- a/src/calculate_mortar_matrices.jl +++ b/src/calculate_mortar_matrices.jl @@ -7,21 +7,61 @@ const seg_integration_points = Dict( (5.0/9.0, +sqrt(3.0/5.0)))) """ -Given segmentation, calculate mortar matrices De and Me. + calculate_mortar_matrices(slave_element_id::Int, + elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + normals::Dict{Int, Vector{Float64}}, + segmentation:MortarSegmentation) + +Calculate mortar matrices De and Me for slave element. + +# Example + +```jldoctest +elements = Dict( + 1 => [1, 2], + 2 => [3, 4]) +element_types = Dict( + 1 => :Seg2, + 2 => :Seg2) +coords = Dict( + 1 => [1.0, 2.0], + 2 => [3.0, 2.0], + 3 => [2.0, 2.0], + 4 => [0.0, 2.0]) +normals = Dict( + 1 => [0.0, -1.0], + 2 => [0.0, -1.0]) +segmentation = Dict(1 => [(2, [-1.0, 0.0])]) +De, Me = calculate_mortar_matrices(1, elements, element_types, + coords, normals, segmentation) + +# output + +([0.583333 0.166667; 0.166667 0.0833333], Dict(2=>[0.541667 0.208333; 0.208333 0.0416667])) + +``` + """ -function calculate_mortar_matrices(sid, mids, elements, element_types, coords, normals, segmentation) - @assert element_types[sid] == :Seg2 +function calculate_mortar_matrices(slave_element_id::Int, + elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + normals::Dict{Int, Vector{Float64}}, + segmentation::MortarSegmentation) + @assert element_types[slave_element_id] == :Seg2 # Initialization + calculate jacobian De = zeros(2, 2) - Me = Dict() - scon = elements[sid] + Me = Dict{Int, Matrix{Float64}}() + scon = elements[slave_element_id] xs1 = coords[scon[1]] xs2 = coords[scon[2]] ns1 = normals[scon[1]] ns2 = normals[scon[2]] J = norm(xs2-xs1)/2.0 # 1. calculate De - for (mid, (xi1, xi2)) in segmentation[sid] + for (mid, (xi1, xi2)) in segmentation[slave_element_id] s = abs(xi2-xi1)/2.0 for (w, ip) in seg_integration_points[3] xi_s = (1-ip)/2*xi1 + (1+ip)/2*xi2 @@ -30,7 +70,7 @@ function calculate_mortar_matrices(sid, mids, elements, element_types, coords, n end end # 2. calculate Me - for (mid, (xi1, xi2)) in segmentation[sid] + for (mid, (xi1, xi2)) in segmentation[slave_element_id] @assert element_types[mid] == :Seg2 Me[mid] = zeros(2, 2) mcon = elements[mid] diff --git a/src/calculate_normals.jl b/src/calculate_normals.jl index 336ef9c..c84aaa9 100644 --- a/src/calculate_normals.jl +++ b/src/calculate_normals.jl @@ -10,7 +10,11 @@ Given elements, element types and node locations, calculate nodal normals by first calculating normal directions for each element and then averaging them in nodes. As a result we get unique normal direction defined to each node. +# Notes +Only linear elements supported. + # Example + ```jldoctest X = Dict(1 => [7.0, 7.0], 2 => [4.0, 3.0], 3 => [0.0, 0.0]) elements = Dict(1 => [1, 2], 2 => [2, 3]) diff --git a/src/calculate_projections.jl b/src/calculate_projections.jl index 6545d63..71ffe28 100644 --- a/src/calculate_projections.jl +++ b/src/calculate_projections.jl @@ -2,27 +2,30 @@ # License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE """ -Find the projection of a slave node `xs`, having normal vector `ns`, onto master -elements with nodes (xm1, xm2). - -## Mathematics - -## References + project_from_slave_to_master(Val{:Seg2}, xs, ns, xm1, xm2) -- Yang2005 +Find the projection of a slave node `xs`, having normal vector `ns`, onto master +elements with nodes (`xm1`, `xm2`). """ -function project_from_slave_to_master{T<:Number}(::Type{Val{:Seg2}}, xs::Vector{T}, ns::Vector{T}, xm1::Vector{T}, xm2::Vector{T}) +function project_from_slave_to_master{T<:Number}(::Type{Val{:Seg2}}, xs::Vector{T}, + ns::Vector{T}, xm1::Vector{T}, + xm2::Vector{T}) nom = ns[1]*(xm1[2] + xm2[2] - 2*xs[2]) - ns[2]*(xm1[1] + xm2[1] - 2*xs[1]) denom = ns[1]*(xm1[2] - xm2[2]) + ns[2]*(xm2[1] - xm1[1]) return nom/denom end """ + project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2) + Find the projection of a master node `xm`, to the slave surface with nodes -(xs1, xs2), in direction of slave surface normal defined by (ns1, ns2). +(`xs1`, `xs2`), in direction of slave surface normal defined by (`ns1`, `ns2`). + """ -function project_from_master_to_slave{T<:Number}(::Type{Val{:Seg2}}, xm::Vector{T}, xs1::Vector{T}, xs2::Vector{T}, ns1::Vector{T}, ns2::Vector{T}) +function project_from_master_to_slave{T<:Number}(::Type{Val{:Seg2}}, xm::Vector{T}, + xs1::Vector{T}, xs2::Vector{T}, + ns1::Vector{T}, ns2::Vector{T}) # Special case when normal is constant. Then we have # unique solution and linear equation to solve. diff --git a/src/calculate_segments.jl b/src/calculate_segments.jl index 755b8a2..c13a13e 100644 --- a/src/calculate_segments.jl +++ b/src/calculate_segments.jl @@ -1,18 +1,65 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE +const MortarSegmentation = Dict{Int, Vector{Tuple{Int, Vector{Float64}}}} + """ + calculate_segments(slave_element_ids::Vector{Int}, + master_element_ids::Vector{Int}, + elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + normals::Dict{Int, Vector{Float64}}) + Given slave surface elements, master surface elements, nodal coordinates and normal direction on nodes of slave surface elements, calculate contact segments. -Returns a dictionary, where key is slave element id and values is a list of -master elements giving contribution to that slave elements and xi-coordinates -of slave side element. +Return type is a dictionary, where key is slave element id and values is a +list of master elements giving contribution to that slave elements and +xi-coordinates of slave side element. + +# Example +```jldoctest +elements = Dict(1 => [1, 2], 2 => [3, 4]) +element_types = Dict(1 => :Seg2, 2 => :Seg2) +coords = Dict( + 1 => [1.0, 2.0], + 2 => [3.0, 2.0], + 3 => [0.0, 2.0], + 4 => [2.0, 2.0]) +normals = Dict( + 1 => [0.0, -1.0], + 2 => [0.0, -1.0]) +slave_ids = [1] +master_ids = [2] +segments = calculate_segments(slave_ids, master_ids, elements, + element_types, coords, normals) + +# output + +Dict{Int64,Array{Tuple{Int64,Array{Float64,1}},1}} with 1 entry: + 1 => Tuple{Int64,Array{Float64,1}}[(2, [-1.0, -0.0])] + +``` + +Here, output result means that slave element #1 has segment with master +element(s) #2 with dimensionless slave element coordinate xi = [-1, 0]. +That is, the start and end point of projection in physical coordinate +system is: + + x_start = 1/2*(1-xi[1])*xs1 + 1/2*(1+xi[1])*xs2 + x_stop = 1/2*(1-xi[2])*xs1 + 1/2*(1+xi[2])*xs2 + """ -function calculate_segments(slave_ids, master_ids, elements, element_types, coords, normals) - S = Dict() - for sid in slave_ids +function calculate_segments(slave_element_ids::Vector{Int}, + master_element_ids::Vector{Int}, + elements::Dict{Int, Vector{Int}}, + element_types::Dict{Int, Symbol}, + coords::Dict{Int, Vector{Float64}}, + normals::Dict{Int, Vector{Float64}}) + S = MortarSegmentation() + for sid in slave_element_ids @assert element_types[sid] == :Seg2 haskey(S, sid) || (S[sid] = []) scon = elements[sid] @@ -20,7 +67,7 @@ function calculate_segments(slave_ids, master_ids, elements, element_types, coor xs2 = coords[scon[2]] ns1 = normals[scon[1]] ns2 = normals[scon[2]] - for mid in master_ids + for mid in master_element_ids @assert element_types[mid] == :Seg2 mcon = elements[mid] xm1 = coords[mcon[1]] diff --git a/test/runtests.jl b/test/runtests.jl index a49b03b..c7ae227 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,7 @@ push!(test_files, "test_calculate_normals.jl") push!(test_files, "test_calculate_projections.jl") push!(test_files, "test_calculate_segments.jl") push!(test_files, "test_calculate_mortar_matrices.jl") -push!(test_files, "test_calculate_assembly.jl") +push!(test_files, "test_calculate_mortar_assembly.jl") @testset "Mortar2D.jl" begin for fn in test_files diff --git a/test/test_calculate_assembly.jl b/test/test_calculate_mortar_assembly.jl similarity index 92% rename from test/test_calculate_assembly.jl rename to test/test_calculate_mortar_assembly.jl index f034d37..dd2ef97 100644 --- a/test/test_calculate_assembly.jl +++ b/test/test_calculate_mortar_assembly.jl @@ -14,8 +14,8 @@ using Mortar2D: calculate_mortar_assembly Em = Dict(i => [i, i+1]+1 for i=ns:ns+nm-2) element_types = Dict(i => :Seg2 for i=1:(ns+nm)) - slave_element_ids = keys(Es) - master_element_ids = keys(Em) + slave_element_ids = collect(keys(Es)) + master_element_ids = collect(keys(Em)) elements = merge(Es, Em) coords = merge(Xs, Xm) diff --git a/test/test_calculate_mortar_matrices.jl b/test/test_calculate_mortar_matrices.jl index 030a2dd..454d26f 100644 --- a/test/test_calculate_mortar_matrices.jl +++ b/test/test_calculate_mortar_matrices.jl @@ -20,7 +20,7 @@ using Mortar2D: calculate_mortar_matrices 1 => [0.0, -1.0], 2 => [0.0, -1.0]) segmentation = Dict(1 => [(2, [-1.0, 0.0])]) - De, Me = calculate_mortar_matrices(1, [2], elements, element_types, + De, Me = calculate_mortar_matrices(1, elements, element_types, coords, normals, segmentation) De_expected = 1/24*[14 4; 4 2] Me_expected = 1/24*[13 5; 5 1]