Update documentation

- Add theory very briefly and examples
- Rename calculate_assembly to calculate_mortar_assembly

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"
                  ])

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

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
-```
-

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

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.

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

src/calculate_assembly.jl → 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])