\documentclass[preprint,10pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

\addtolength{\evensidemargin}{-2cm}
\addtolength{\oddsidemargin}{-2cm}
\addtolength{\textwidth}{6cm}

\usepackage{url}
\usepackage{graphicx}
\usepackage{caption}
%\usepackage{subcaption}
\usepackage{subfig}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage{hyperref}
\usepackage{amsmath,amssymb,amsfonts,amsthm,stmaryrd}
\usepackage{dsfont}
\usepackage{eucal}
\usepackage{mathrsfs}
\usepackage{listings}    % for insert source code such as C++ code
\usepackage{color}
\usepackage{float}
\usepackage{natbib}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage{setspace}
\usepackage{algorithmicx}
\usepackage{algorithm}
\usepackage{algpseudocode}

\input{pre_commands.tex}

\newcommand{\grbf}[1]{\mbox{\boldmath{$#1$}}}

\floatstyle{ruled}
\newfloat{Program}{thp}{lop}
\floatname{Program}{Program}

\floatstyle{ruled}
\newfloat{Fbox}{thp}{lop}
\floatname{Fbox}{Box}

\newcolumntype{C}{>{\centering\arraybackslash}X}

\lstloadlanguages{C++,make,sh,csh,Matlab}

\definecolor{darkgray}{rgb}{0.95,0.95,0.95}
\definecolor{mygreen}{rgb}{0,0.6,0}

\lstset{backgroundcolor=\color{darkgray}}

\lstdefinestyle{Matlab}
{
 basicstyle=\footnotesize, numbers=none, numberstyle=\tiny,%
 showstringspaces=false, language=Matlab, escapechar=|,frame=tb,%
 commentstyle=\color{mygreen}
}

\lstdefinestyle{Matlab-num}
{
 basicstyle=\footnotesize, numbers=left, numberstyle=\tiny,%
 showstringspaces=false, language=Matlab, escapechar=|,frame=tb,%
commentstyle=\color{mygreen}
}

%-----------------------------------------------------------------------
\newcommand{\tty}[1]{\textnormal{\texttt{#1}}}
\newcommand{\sym}[1]{\textnormal{\textit{#1}}}

\lstnewenvironment{snippet}[1][]
{
 \lstset{style=Matlab, xleftmargin=5mm, gobble=4, #1}
}
{}

\lstnewenvironment{snippet1}[1][]
{
 \lstset{style=Matlab-num, xleftmargin=5mm, gobble=4, #1}
}
{}


% Uncomment the following line to make the final paper.!!!

%\definecolor{red}{gray}{0} 


\def\bibsection{\section*{References}}

%%%%%%%%%%%%%%%%%%%%%
\begin{document}

% for listing package

\definecolor{MyDarkBlue}{rgb}{1, 0.9, 1}
\lstset{language=Matlab,
       basicstyle=\footnotesize,
%       keywordstyle=\color{red},
       commentstyle=\itshape,
       stringstyle=\ttfamily,
       showstringspaces=false,
       tabsize=2}
%\lstset{backgroundcolor=\color{MyDarkBlue}}
\lstdefinestyle{commentstyle}{color=\color{green}}

\theoremstyle{remark}
\newtheorem{thm}{Theorem}[section]
\newtheorem{rmk}[thm]{Remark}

%*************************!!!!! uncomment this for the final version

\definecolor{red}{gray}{0}
\definecolor{blue}{gray}{0}

% double spacing for revising
%\doublespacing

\begin{frontmatter}

\title{A Nitsche's method method for mixed dimensional analysis: volume and surface solid-structure 
coupling}       

\author[cardiff]{Vinh Phu Nguyen  \fnref{fn1}}
\author[cardiff]{Pierre Kerfriden  \fnref{fn2}}
\author[cardiff]{St\'{e}phane P.A. Bordas \corref{cor1}\fnref{fn3}}

\cortext[cor1]{Corresponding author}

\address[cardiff]{School of Engineering, Institute of Mechanics and Advanced
Materials, Cardiff University, Queen's Buildings, The Parade, Cardiff \\
CF24 3AA}
%\address[weimar]{Institute of Structural Mechanics, Bauhaus-Universit\"{a}t
%Weimar, Marienstra\ss{}e 15 99423 Weimar}

\fntext[fn1]{\url nguyenpv@cardiff.ac.uk, ORCID: 0000-0003-1212-8311}
\fntext[fn2]{\url pierre@cardiff.ac.uk}
\fntext[fn3]{\url stephane.bordas@alum.northwestern.edu, ORCID: 0000-0001-7622-2193}

\maketitle
\tableofcontents

\begin{abstract}
A Nitche's method is presented to couple different mechanical models and discretisations. They 
include coupling of non-conforming two and three dimensional NURBS (Non Uniform Rational B-Splines) patches, 
coupling of a solid and a beam and of a solid
and a plate. Applications can be found in a multi-patch NURBS based isogeometric analysis,
multi-dimensional analysis in which parts of a structure are modeled with solid elements
and others are modeled using a coarser model with beam and/or plate elements. We present
various examples including composite delamination analysis.
\end{abstract}

\begin{keyword} 
   Nitsche \sep mixed dimensional analysis (MDA) \sep isogeometric analysis (IGA) 
\end{keyword}

\end{frontmatter}

% ============================================================================

% -----------------------
\section{Introduction}
% -----------------------


Nitsche's method \cite{nitsche} was originally proposed to weakly enforce Dirichlet boundary conditions
as an alternative to equivalent pointwise constraints.
The idea behind a Nitsche based approach is 
to replace the Lagrange multipliers arising in a dual formulation through their physical representation, 
namely the normal flux at the interface. Nitsche also added an extra penalty like term to restore 
the coercivity of the bilinear form. The method can be seen to lie in between the Lagrange multiplier
method and the penalty method.
The method has seen a resurgence in recent years and was applied for interface problems 
\cite{Hansbo20025537,Dolbow2009a},  for connecting overlapping meshes 
\cite{MZA:8203296,MZA:8203286,Sanders2012a,Sanders2011a}, for imposing Dirichlet boundary conditions 
in meshfree methods \cite{FernándezMéndez20041257}, in immersed boundary methods 
\cite{ruess2013weakly,NME:NME3339,embar_imposing_2010}, in fluid mechanics 
\cite{Bazilevs200712}
and for contact mechanics \cite{nitsche-wriggers}. It has also been applied for stabilising constraints
in enriched finite elements \cite{Sanders2008a}.

In this paper, a Nitsche's method is presented to (1) couple non-conforming Lagrange
and isogeometric discretisations, (2) couple two dimensional (2D) continua and beams
and (3) three dimensional (3D) continua and plates.
The need for the research presented in this paper stems from the problem of progressive
failure analysis of composite laminates of which recent studies were carried out by the authors
\cite{Nguyen2013,nguyen_cohesive_2013,nguyen-offset}. It was shown that using NURBS (Non Uniform
Rational B-splines) as the finite element basis functions--the concept coined as
isogeometric analysis (IGA) by Hughes and his co-workers \cite{hughes_isogeometric_2005,cottrel_book_2009}
results in speed up in both pre-processing and processing step of delamination analysis of composite laminates.
However in order for IGA to be applied for real problems such as large composite panels utilized in aerospace
industry Mixed Dimensional Analysis (MDA) \cite{Cuilliere2010} must be employed. In MDA 
some portions of the object of interest can be modeled using reduced-dimensional elements 
(typically beams and shells), while other portions must be modelled using volume (solid) elements.
In this way, large structures are likely tractable. 
The question is how to model the coupling of different models in a flexible and efficient manner.

Broadly speaking the coupling can be either surface coupling (non-overlapping coupling) or
volume coupling (overlapping coupling) \cite{Guidault2007b}, see Fig.~\ref{fig:concepts}. 
Volume coupling indicates the existence of 
a region in which both model co-exist and is usually realized using the Arlequin method 
\cite{Dhia2005a} of which Abaqus implementation can be found in \cite{Qiao2011a}.
Arlequin method is best suited for coupling different physical models such as continuum-particle
modeling see Refs. \cite{Bauman2008a,Pfaller2011a,Wellmann2012a,Rousseau2009a} among others.
In Arlequin method care must be taken when choosing the space of Lagrange multipliers
and the numerical integration of coupling terms on unstructured meshes is a non-trivial task, particularly
for three dimensional problems. In surface coupling there is
no overlapping of models and the two models can be coupled using Lagrange multipliers as in the mortar
method \cite{mortar-1999}, or the penalty method. Both Lagrange multiplier method and the penalty method
have their own disadvantages--in the former it is the introduction of extra unknowns 
(thus destroys positive definiteness of the augmented system of equation) and the proper choice
of the Lagrange multiplier discrete space. In the latter it is the choice of the penalty parameter.
Other method is to use constraint equations \cite{Unger,NME:NME967} (the common multipoint constraint MPC method
in commercial FE packages). MPC based coupling is a strong coupling method whereas all other coupling methods
are weak couplings.
In \cite{NME:NME967} constraint equations are determined by equating the work done by the stresses in each part of the model at the interface between dimensions. The authors also provided error estimations for their
MDA. The theory of MPC method is described in details in \cite{felippa:note}. In \cite{Unger} a comparative study of three coupling methods: MPC, mortar method and Arlequin method
was presented for the concurrent multiscale modeling of concrete material where the macro domain is
discretised by a coarse mesh and the micro domain (with heterogeneities) is discretised by a very fine mesh.
The authors concluded that the Arlequin method is too expensive to consider.
A somehow related approach is the global-local technique, see \eg \cite{felippa:global-local}
and references therein. 
In a global-local approach, the whole system is first analyzed as a global entity, discarding details 
deemed not to affect its overall behavior. Local details (cutoffs, cracks \etc) are then analyzed 
using the results of the global analysis as boundary conditions. In \cite{allix-non-intrusive} a global-local
technique was adopted as a non-intrusive coupling method.
The authors in \cite{Gabbert} compared the performance, through a simple numerical example, of
the MPC method, the Arlequin method and the global-local method and the conclusion was that the Arlequin method
is a promising coupling method for it results in lower stress jumps at the coupling boundary.
In aforementioned coupling methods either a modified weak formulation has to used
or a system of linear equations has to be solved with a set of (linear) constraints (for the MPC method).
Transition elements are a yet another method in which standard weak forms can be reused 
to couple solid elements and beam/shell elements 
see \eg \cite{NME:NME1620150704,transition-gmur} or in an IGA context the bending strip method \cite{kiendl_bending_2010,Raknes2013127}. 
A mesh-geometry-based solution to mixed-dimensional coupling was presented in \cite{Cuilliere2010}.
It should be emphasized that solving mixed dimensional models is much faster than solving full 3D models, but the time required to process these mixed dimensional models might eliminate that advantage.


\begin{figure}
  \centering
  \includegraphics[width=0.6\textwidth]{coupling-concepts}
  \caption{Discretisation (a) and model (b) coupling which can be either a surface or volume coupling. The bold lines represent the coupling surfaces and the grey domains are the coupling volumes.}
  \label{fig:concepts}
\end{figure}

In this paper, we decided to use a surface coupling method and a Nitsche based weak coupling formulation.
Surface coupling was chosen because  we are not solving a multi-physics model and the Nitsche's method
is adopted due to a rich mathematics behind it and the related discontinuous Galerkin methods \cite{Arnold01unifiedanalysis}. The method is symmetric for symmetric problems and do not need additional degrees of freedom. 
However, a user-defined stabilization parameter is required.
Our target is modeling large scale composite laminates and the
long-term goal is to start the analysis from a laminate shell model, build goal-oriented error estimates 
to identify the \textit{hot-spots} where intra or inter-laminar failure is likely to take place, to refine 
these zones with a continuum description, and, once the delaminations and intra-laminar cracks are fully open, replace them by an homogenised shell model. 

The contribution of the paper are listed as follows

\begin{enumerate}
\item A detailed description of Nitsche's method for coupling 2D/2D and 3D/3D Lagrangian and NURBS meshes (Section \ref{continua-continua});
\item Formulation and implementation details for 2D solid and beam coupling (Section \ref{sec:beam-continuum});
\item Formulation and implementation details for 3D solid and plate coupling (Section \ref{sec:plate-continuum});
\item Study of both symmetric Nitsche's method and non symmetric Nitsche's method.
\end{enumerate}
The formulation of continuum-continuum coupling is general and will be specialised to solid-beam and
solid-plate coupling. Numerical examples are provided in Section \ref{sec:examples}.
The presented algorithm applies for both Lagrange basis functions and NURBS basis functions. The latter
have been extensively used in isogeometric analysis (IGA) \cite{hughes_isogeometric_2005,cottrel_book_2009} and
facilitates the implementation of rotation-free beam/plate/shell elements 
see \eg \cite{kiendl_isogeometric_2009,benson_large_2011} which is also adopted in this work. 
Additionally, IGA is also the numerical framework
for our recent works on failure analysis of composite laminates 
\cite{Nguyen2013,nguyen_cohesive_2013,nguyen-offset}. We confine to quasi-static small strain problems 
although Nitsche's method was successfully applied to finite deformation \cite{Sanders2012a}.


We denote $d_p$ and $d_s$ as the number of parametric directions and spatial directions respectively.
Both tensor and matrix notations are used.
In tensor notation, tensors of order one or greater are written in boldface. Lower case bold-face letters 
are used for first-order tensor whereas upper case bold-face letters indicate high-order tensors. 
The major exception to this rule are the physical second order stress tensor and the strain tensor
which are written in lower case. In matrix notation, the same symbols as for tensors are used to denote
the matrices but the connective operator symbols are skipped. 

%-----------------------------------------------------------
\section{Problem description}\label{problem}
%-----------------------------------------------------------

%--------------------------------
\subsection{Solid-beam coupling}

We define the domain $\Omega \subset \mathbb{R}^{d_s}$ with boundary $\Gamma \equiv \partial \Omega$. 
The boundary is partitioned as $\Gamma = \overline{\Gamma_u \cup \Gamma_t}$, $\Gamma_u \cap \Gamma_t = \emptyset$ where $\Gamma_u$ and $\Gamma_t$ denote the Dirichlet and Neumann boundaries respectively with an 
overline representing a closed set. For sake of simplicity, we assume 
there is only one internal boundary denoted by $\Gamma_*$ that divides the domain into
two non-overlapping domains $\Omega^i, i=1,2$ such that $\Omega=\Omega^1 \cup \Omega^2$. 
The governing equations are 

\begin{itemize}

\item For the solid part

\begin{subequations}
\begin{alignat}{2}
%  \begin{array}{lllll}
-\text{div} \;\bsym{\sigma} &=  \vm{b}         &\quad\text{on} \quad \Omega \\ 
             \vm{u}         &=  \bar{\vm{u}}   &\quad\text{on} \quad \Gamma_u \\
  \bsym{\sigma} \cdot \vm{n} &=  \bar{\vm{t}}  &\quad\text{on} \quad \Gamma_t \label{eq:Neumann}\\
  \vm{u}^1 &=  \vm{u}^2                        &\quad\text{on} \quad \Gamma_* \\
  \bsym{\sigma}^1 \cdot \vm{n}^1 &= -\bsym{\sigma}^2 \cdot \vm{n}^2  &\quad\text{on} \quad \Gamma_* \label{eq:tr}
%  \end{array}
\end{alignat}
\end{subequations}

\item For the beam part

\end{itemize}

where $\vm{u}$ is the displacement field and $\bsym{\sigma}$ denotes the stress field; 
the last two equations express the continuity of displacements and tractions across $\Gamma_*$.
The prescribed displacement and traction are denoted by $\bar{\vm{u}}$ and $\bar{\vm{t}}$, respectively.
The outward unit normals to $\Omega^1$ and $\Omega^2$ are $\vm{n}^1$, $\vm{n}^2$. A super script, $(i)$, 
is used to denote a quantity that is valid over region $\Omega^i$, with $i = 1, 2$.

Under the small strain condition, the infinitesimal strain tensor reads $\epsilon_{ij}=0.5(u_{i,j}+u_{j,i})$.
Constitutive equations are given by

\begin{equation}
\bsym{\sigma}^i = \vm{C}^i : \bsym{\epsilon}^i, \quad i=1,2
\end{equation}
where the constitutive tensors are denoted by $\vm{C}^1$ and $\vm{C}^2$. For linear 
isotropic elastic materials, the constitutive tensor is written as

\begin{equation}
C_{ijkl} = \lambda \delta_{ik}\delta_{kl} + \mu (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})
\end{equation}
where $\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)}$ and $\mu=\frac{E}{2(1+\nu)}$ are the Lam{\'e} constants;
$E$ and $\nu$ are the Young's modulus and Poisson's ratio, respectively and $\delta_{ij}$ is the Kronecker
delta tensor. 

\begin{figure}
  \centering
  \includegraphics[width=\textwidth]{domain1}
  \caption{Computational domain with an internal interface.}
  \label{fig:domain}
\end{figure}

%--------------------------------
\subsection{Solid-plate coupling}

%---------------------
\section{Weak forms}

%--------------------------------
\subsection{Solid-beam coupling}

The derivation of the weak form starts with multiplying the equation of linear momentum with
a test function and integrating over the domain

\begin{equation}\label{weakform1}
 \int_\Omega \delta u_i ( \sigma_{ij,j} + b_i ) \mathrm{d}\Omega = 0
\end{equation}

\noindent The first term of the above equation can be written as follows according to the product
rule

\begin{equation}\label{weak2}
 \int_\Omega \delta u_i  \sigma_{ij,j}  \mathrm{d}\Omega = \int_\Omega \frac{\partial }{\partial x_j} ( \delta u_i  \sigma_{ij} ) \mathrm{d}\Omega
- \int_\Omega \frac{\partial \delta u_i}{\partial x_j}   \sigma_{ij}  \mathrm{d}\Omega
\end{equation}

\noindent Using the Gauss theorem, the first term of the above becomes

\begin{equation}
 \int_\Omega \frac{\partial }{\partial x_j} ( \delta u_i  \sigma_{ij} ) \mathrm{d}\Omega
= \int_\Gamma  \delta u_i  \sigma_{ij} n_j  \mathrm{d}\Gamma + 
   \int_{\Gamma_*^1}  \delta u_i^1  \sigma_{ij}^1 n_{j}^1  \mathrm{d}\Gamma +
   \int_{\Gamma_*^2}  \delta u_i^2  \sigma_{ij}^2 n_{j}^2  \mathrm{d}\Gamma
\end{equation}

\noindent By defining the jump operator $\jump{\cdot}\equiv(\cdot)^1-(\cdot)^2$ on $\Gamma_*$ and making use
of the Neumann boundary condition (Equation~\eqref{eq:Neumann}), we can rewrite the above as

\begin{equation}\label{weak3}
 \int_\Omega \frac{\partial }{\partial x_j} ( \delta u_i  \sigma_{ij} ) \mathrm{d}\Omega
= \int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma + 
\int_{\Gamma_*}  \jump{\delta u_i \sigma_{ij}} n_j^1 \mathrm{d}\Gamma 
\end{equation}

\noindent By using the following identity
\begin{equation}\label{weak3}
\jump{\delta u_i \sigma_{ij}} = \jump{\delta u_i}\{\sigma_{ij}\} + \{\delta u_i\}\jump{\sigma_{ij}}
\end{equation}
in which the average operator is defined as $\{\cdot\}\equiv 1/2[(\cdot)^1+(\cdot)^2]$ on $\Gamma_*$, we can write

\begin{equation}\label{weak3}
 \int_\Omega \frac{\partial }{\partial x_j} ( \delta u_i  \sigma_{ij} ) \mathrm{d}\Omega
= \int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma + 
\int_{\Gamma_*}  \jump{\delta u_i} \{\sigma_{ij}\} n_j^1 \mathrm{d}\Gamma 
\end{equation}
where the traction continuity (Equation~\eqref{eq:tr}), was used.

Substitution of Equations (\ref{weak2}) and (\ref{weak3}) into Equation(\ref{weakform1}), we obtain 

\begin{equation}
\int_\Omega \frac{\partial \delta u_i}{\partial x_j}   \sigma_{ij}  \mathrm{d}\Omega 
=  \int_\Omega \delta u_i  b_i \mathrm{d}\Omega+ 
\int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma +
\int_{\Gamma_*}  \jump{\delta u_i} \{\sigma_{ij}\} n_j^1 \mathrm{d}\Gamma 
\end{equation}

In order to have a symmetric method one more term is added to the weak form as follows

\begin{equation}
\int_\Omega \frac{\partial \delta u_i}{\partial x_j}   \sigma_{ij}  \mathrm{d}\Omega -
\int_{\Gamma_*}  \jump{\delta u_i} \{\sigma_{ij}\} n_j^1 \mathrm{d}\Gamma -
\int_{\Gamma_*}  \jump{u_i} \{\delta \sigma_{ij}\} n_j^1 \mathrm{d}\Gamma 
=  \int_\Omega \delta u_i  b_i \mathrm{d}\Omega+ 
\int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma 
\end{equation}


According to Nitsche's method, a penalty-like stabilisation term is added to the weak form resulting in

\begin{multline}
\int_\Omega \frac{\partial \delta u_i}{\partial x_j}   \sigma_{ij}  \mathrm{d}\Omega -
\int_{\Gamma_*}  \jump{\delta u_i} \{\sigma_{ij}\} n_j^1 \mathrm{d}\Gamma -
\int_{\Gamma_*}  \jump{u_i} \{\delta \sigma_{ij}\} n_j^1 \mathrm{d}\Gamma +
\int_{\Gamma_*}  \alpha \jump{\delta u_i} \jump{u_i} \mathrm{d}\Gamma 
=  \int_\Omega \delta u_i  b_i \mathrm{d}\Omega+ 
\int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma 
\label{eq:dg-weakform}
\end{multline}
in which $\alpha$ is a penalty-like stabilisation term served to stabilise the method. 
Although two extra terms are added, they vanish as $\jump{\vm{u}}=\vm{0}$ the weak form is thus equivalent
to the strong form. Therefore, the method is consistent.
Basically, mathematicians showed that there exists a minimum regularisation parameter 
that guarantees stability, see \eg \cite{Griebel}. 

For discretisation we rewrite Equation~\eqref{eq:dg-weakform} in a matrix form as follows

\begin{equation}
\int_\Omega  \delta \bsym{\epsilon}\trans  \bsym{\sigma}  \mathrm{d}\Omega -
\int_{\Gamma_*}  \jump{\delta \vm{u}}\trans \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\int_{\Gamma_*}  \{\delta\bsym{\sigma}\}\trans \vm{n}\trans \jump{\vm{u}} \mathrm{d}\Gamma +
\int_{\Gamma_*}  \alpha \jump{\delta \vm{u}}\trans \jump{\vm{u}} \mathrm{d}\Gamma 
=  \int_{\Gamma_t}  \delta \vm{u}\trans \bar{\vm{t}}  \mathrm{d}\Gamma 
\label{eq:dg-weakform-matrix}
\end{equation}
in which for sake of simplicity body force was neglected. Superscript T denotes the transpose operator.
Second order tensors ($\sigma_{ij}$ and $\epsilon_{ij}$) are written using the Voigt notation
as column vectors; $\bsym{\sigma}=[\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{xz}]\trans$, $\bsym{\epsilon}=[\epsilon_{xx}, \epsilon_{yy}, \epsilon_{zz}, 2\epsilon_{xy}, 2\epsilon_{yz}, 2\epsilon_{xz}]\trans$,
and $\vm{n}$ (note that we removed the subscript 1 for subsequent derivations) 
is a matrix that reads

\begin{equation}
\vm{n}_{2D} = \begin{bmatrix}
n_x & 0 & n_y \\ 0 & n_y & n_x
\end{bmatrix}, \quad
\vm{n}_{3D} = \begin{bmatrix}
n_x & 0 & 0 & n_y & 0 & n_z \\
0 & n_y & 0 & n_x & n_z & 0\\
0 & 0 & n_z & 0 & n_y & n_x  
\end{bmatrix}\label{eq:n-matrix}
\end{equation}
The average operator for the stress field can be written generally as \cite{Sanders2012a}

\begin{equation}
\{\bsym{\sigma}\}=\gamma \bsym{\sigma}^1 + (1-\gamma)\bsym{\sigma}^2
\label{eq:general-average}
\end{equation}
where $0 \le \gamma \le 1$. The usual average operator is reproduced if $\gamma=0.5$ is used.
Equation~\eqref{eq:general-average} is often utilized to join a soft model and a stiff one \cite{Sanders2011a}.
Taking $\gamma=1$ (or $\gamma=0$) results in the one-sided mortaring method. 
In this paper, the standard average operator is used unless otherwise stated.

%--------------------------------
\subsection{Solid-plate coupling}

%---------------------------
\section{Discretisation}

The two subdomains $\Omega^k (k=1,2)$ are discretised independently using finite elements.
At the interface $\Gamma_*$ there is a mismatch between the two meshes, cf. Fig.~\ref{fig:domain-mesh}.
The displacement fields and its variations are given by

\begin{figure}
  \centering
  \includegraphics[width=0.4\textwidth]{domain2}
  \caption{Independent discretisations of the domains.}
  \label{fig:domain-mesh}
\end{figure}

\begin{equation}
  \begin{array}{llllllll}
  u_i^1 & = & N_I^1 a_{iI}^1 & \text{and} & \delta u_i^1 & = & N_I^1 \delta a_{iI}^1 & \text{on $\Omega^1$}\\  
  u_i^2 & = & N_I^2 a_{iI}^2 & \text{and} & \delta u_i^2 & = & N_I^2 \delta a_{iI}^2 & \text{on $\Omega^2$}\\  
  \end{array}
  \label{eq:approximation}
\end{equation}
where $N^1_I$ and $N^2_I$ are the finite element shape functions (which can be either standard Lagrange basis
polynomials or NURBS basis functions to be specified in concrete examples) 
associated to mesh 1 and 2, respectively.
$a_{iI}^{1/2}$ denotes the $i$ component ($i=1,\ldots,d_s$) of the displacements at node I 
that belongs to domain $\Omega^{1/2}$.

The stresses, strains and displacements are given by 

\begin{equation}
\bsym{\sigma}^k = \vm{C}^k \vm{B}^k \vm{a}^k, \quad
\bsym{\epsilon}^k = \vm{B}^k \vm{a}^k, \quad 
\vm{u}^k=\vm{N}^k \vm{a}^k \; (k=1,2)
\end{equation}
where $\vm{B}$ is the standard strain-displacement matrix and $\vm{N}$ represents the standard
shape function matrix. For two dimensional element $e$, they are given by

\begin{equation}
\vm{B}_e^k = \begin{bmatrix}
N_{1,x}^k & 0 & N_{2,x}^k & 0 & \ldots\\
0 & N_{1,y}^k & 0 & N_{2,y}^k & \ldots \\
N_{1,y}^k & N_{1,x}^k & N_{2,y}^k & N_{2,x}^k & \ldots 
\end{bmatrix},\quad \vm{N}_e^k = \begin{bmatrix}
N_1^k & 0 & N_2^k & 0 & \ldots\\
0 & N_1^k & 0 & N_2^k & \ldots       
\end{bmatrix}
\end{equation}
Expressions for three dimensional elements can be found in many FEM textbooks \eg
\cite{hughes-fem-book}. The notation $N_{I,x}$ denotes
the derivative of shape function $N_I$ with respect to $x$. This notation for partial derivatives
will be used in subsequent sections.

The jump operator is given by

\begin{equation}
\jump{\vm{u}} = \vm{N}^1\vm{a}^1 - \vm{N}^2 \vm{a}^2
\end{equation}
Using those finite element approximations, Equation~\eqref{eq:dg-weakform-matrix} becomes

\begin{multline}
\int_\Omega  \delta \vm{a}\trans \vm{B}\trans \bsym{\sigma}  \mathrm{d}\Omega -
\int_{\Gamma_*} (\vm{N}^1\delta\vm{a}^1-\vm{N}^2\delta\vm{a}^2)\trans \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\int_{\Gamma_*}  \frac{1}{2}(\vm{C}^1\vm{B}^1\delta\vm{a}^1+\vm{C}^2\vm{B}^2\delta\vm{a}^2)\trans \vm{n}\trans \jump{\vm{u}} \mathrm{d}\Gamma \\ 
+\int_{\Gamma_*}  \alpha  (\vm{N}^1\delta\vm{a}^1-\vm{N}^2\delta\vm{a}^2)\trans \jump{\vm{u}} \mathrm{d}\Gamma 
=  \int_{\Gamma_t}  \delta \vm{a}\trans \vm{N}\trans \bar{\vm{t}}  \mathrm{d}\Gamma 
\label{eq:dg-weakform1}
\end{multline}
where the subscripts 1 and 2 refer to quantities that are associated to mesh 1 
and 2, respectively. 

Since $\delta\vm{a}$ are arbitrary\footnote{For simplicity, Dirichlet boundary
conditions are enforced as in standard FEM.}, Equation~\eqref{eq:dg-weakform1} leads to

\begin{equation}
\begin{array}{lll}
\D\int_\Omega  \vm{B}\trans \bsym{\sigma}  \mathrm{d}\Omega & = & \D \int_{\Gamma_t}   \vm{N}\trans \bar{\vm{t}}  \mathrm{d}\Gamma \\
-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{1\text{T}}\vm{C}^{1\text{T}}\vm{n}\trans\jump{\vm{u}} +
\D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \jump{\vm{u}} \mathrm{d}\Gamma & = & \vm{0}\\
\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{2\text{T}}\vm{C}^{2\text{T}}\vm{n}\trans\jump{\vm{u}} -
\D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \jump{\vm{u}} \mathrm{d}\Gamma & = & \vm{0}\\
\end{array}
\end{equation}
where the first equation is a standard one and the second and third equations are due to the
interface $\Gamma_*$.

Using the standard notation in nonlinear finite element method, the above is written as
\begin{equation}
\label{eq:fint-fext}
\vm{f}_\text{int}(\vm{a}) \equiv \vm{f}_\text{int}^{c} + \vm{f}_\text{int}^{d} = \vm{f}_\text{ext} \equiv  \int_{\Gamma_t}   \vm{N}\trans \bar{\vm{t}}  \mathrm{d}\Gamma
\end{equation}
in which the standard volumetric internal force vector reads
\begin{equation}
\vm{f}_\text{int}^c =\int_\Omega  \vm{B}\trans \bsym{\sigma}  \mathrm{d}\Omega  
\end{equation}
and the interfacial internal force vector $\vm{f}_\text{int}^d=\vm{f}_\text{int}^{d,1}+\vm{f}_\text{int}^{d,2}$ 
of which components are given by

\begin{equation}
\begin{split}
\vm{f}_\text{int}^{d,1}(\vm{a}^1,\vm{a}^2) &=-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{1\text{T}}\vm{C}^{1\text{T}}\vm{n}\trans\jump{\vm{u}} +
\D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \jump{\vm{u}} \mathrm{d}\Gamma \\
\vm{f}_\text{int}^{d,2}(\vm{a}^1,\vm{a}^2)  &=\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \{\bsym{\sigma}\} \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{2\text{T}}\vm{C}^{2\text{T}}\vm{n}\trans\jump{\vm{u}} -
\D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \jump{\vm{u}} \mathrm{d}\Gamma \\
\end{split}
\label{eq:internal-force-dg}
\end{equation}
where explicit dependence of the internal force vectors on the nodal displacements were made to ease
the linearisation to be presented.

Equation~\eqref{eq:fint-fext} can be solved iteratively using the Newton-Raphson method.
The tangent stiffness matrix is composed of two parts: one part from the interior domain
and the second part from the interface. The former, which is standard, is given by

\begin{equation}
\vm{K}^c = \int_\Omega  \vm{B}\trans \vm{C} \vm{B}  \mathrm{d}\Omega 
\end{equation}
As preparation for the derivation of the interfacial tangent matrices, 
Equation~\eqref{eq:internal-force-dg} is elaborated as follows

\begin{align}
\begin{split}
\vm{f}_\text{int}^{d,1}&= - \int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} 
\frac{1}{2}(\vm{C}^1\vm{B}^1\vm{a}^1+\vm{C}^2\vm{B}^2\vm{a}^2) \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{1\text{T}}\vm{C}^{1\text{T}}\vm{n}\trans(\vm{N}^1\vm{a}^1-\vm{N}^2\vm{a}^2) \\&+
\D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} (\vm{N}^1\vm{a}^1-\vm{N}^2\vm{a}^2)  \mathrm{d}\Gamma
\end{split}\\
\begin{split}
\vm{f}_\text{int}^{d,2} &=\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}(\vm{C}^1\vm{B}^1\vm{a}^1+\vm{C}^2\vm{B}^2\vm{a}^2)  \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{2\text{T}}\vm{C}^{2\text{T}}\vm{n}\trans (\vm{N}^1\vm{a}^1-\vm{N}^2\vm{a}^2)\\ & -
\D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} (\vm{N}^1\vm{a}^1-\vm{N}^2\vm{a}^2)  \mathrm{d}\Gamma 
\end{split}
\end{align}
and the tangent matrices associated to these internal force vectors are given by

\begin{align}
\begin{split}
\vm{K}^{d,11}_{[3n_1\times3n_1]}&=-\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{1\text{T}}\vm{C}^{1\text{T}}\vm{n}\trans \vm{N}^1 + \D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \vm{N}^1 \mathrm{d}\Gamma\\
\vm{K}^{d,12}_{[3n_1\times3n_2]} &=-\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma +
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{1\text{T}}\vm{C}^{1\text{T}}\vm{n}\trans \vm{N}^2 - \D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \vm{N}^2 \mathrm{d}\Gamma
\end{split}\\
\begin{split}
\vm{K}^{d,21}_{[3n_2\times3n_1]}  &=\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma -
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{2\text{T}}\vm{C}^{2\text{T}}\vm{n}\trans \vm{N}^1 - \D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \vm{N}^1 \mathrm{d}\Gamma\\
\vm{K}^{d,22}_{[3n_2\times3n_2]}  &=\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma +
\D\int_{\Gamma_*} \frac{1}{2}\vm{B}^{2\text{T}}\vm{C}^{2\text{T}}\vm{n}\trans \vm{N}^2 + 
\D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \vm{N}^2 \mathrm{d}\Gamma
\end{split}
\end{align}
where the subscripts denote the dimension of the matrices; $n_1$ and $n_2$ are the number of nodes
per element of $\Omega^1$ and $\Omega^2$. 
For 2D problems, the dimension of $\vm{K}^{d,11}$ is $[2n_1\times2n_1]$.
In a more compact form, 
the tangent matrix associated with the interface are then given by

\begin{equation}
\vm{K}^d = \vm{K}^d_{dg} + (\vm{K}^d_{dg})\trans + \vm{K}^d_{pe}
\end{equation}
in which the components are given by

\begin{equation}
\vm{K}^d_{dg} = \begin{bmatrix}
-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma \\
\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma 
\end{bmatrix}\label{eq:nitsche-kdg}
\end{equation}
and by
\begin{equation}
\vm{K}^d_{pe} = \begin{bmatrix}
\D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \vm{N}^1 \mathrm{d}\Gamma &
 - \D\int_{\Gamma_*}  \alpha  \vm{N}^{1\text{T}} \vm{N}^2 \mathrm{d}\Gamma\\
 - \D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \vm{N}^1 \mathrm{d}\Gamma&
 \D\int_{\Gamma_*}  \alpha  \vm{N}^{2\text{T}} \vm{N}^2 \mathrm{d}\Gamma
\end{bmatrix}\label{eq:nitsche-kpe}
\end{equation}

If the average operator defined in Equation~\eqref{eq:general-average} is used, 
we have   

\begin{equation}
\vm{K}^d_{dg} = \begin{bmatrix}
\gamma\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
(1-\gamma)\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \vm{C}^2\vm{B}^2 \mathrm{d}\Gamma \\
-\gamma\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
-(1-\gamma)\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \vm{C}^2\vm{B}^2 \mathrm{d}\Gamma 
\end{bmatrix}\label{eq:nitsche-kdg-general}
\end{equation}

%---------------------------
\subsection{Implementation}

For linear elasticity, the discrete system of equations is given by

\begin{equation}
\left(\vm{K}^c + \vm{K}^d_{dg} + (\vm{K}^d_{dg})\trans + \vm{K}^d_{pe}\right) \vm{a} = \vm{f}_\text{ext}
\end{equation}
Domain  integrals are carried out in a standard manner by looping over the elements

\begin{equation}
\int_\Omega f d \Omega = \int_{\Omega^1} f d \Omega + \int_{\Omega^2} f d \Omega 
=  \bigcup_{e=1}^{nel^1}\int_{\Omega_e^1} f d \Omega + 
\bigcup_{e=1}^{nel^2}\int_{\Omega_e^2} f d \Omega
\end{equation}
where $nel^k$ denotes the number of volume elements of $\Omega^k$ and $\bigcup$ denotes the
standard assembly operator. The contour integral over the Neumman
boundary to compute the external force is performed in a standard manner. 
What remains is the integrals over the internal
interface $\Gamma_*$ of which implementation complexity depends on the triangulations.

\begin{figure}
  \centering
  \includegraphics[width=0.8\textwidth]{domain4}
  \caption{Independent discretisations of the domains: hierarchical meshes. The interface $\Gamma_*$ is
  discretised using the element edges of $\Omega^2$ that intersect $\Gamma_*$. For the grey element, the Gauss point is denoted by the red star which is mapped to the GP in element 1 (green star). }
  \label{fig:domain-mesh-hier}
\end{figure}

%-----------------------------------
\subsubsection{Hierarchical meshes}

First, we consider hierarchical meshes as shown in Fig.~\ref{fig:domain-mesh-hier}.
In this case, the interface integrals can be straightforwardly calculated as explained in what follows.
Let assume that a fine mesh is adopted for $\Omega^2$ and a coarse mesh for $\Omega^1$, cf.
Fig.~\ref{fig:domain-mesh-hier}. We use the fine elements on $\Gamma_*$ to evaluate the interfacial integral

\begin{equation}
\int_{\Gamma_*} f(N^1,N^2) d \Gamma = \bigcup_{e=1}^{nbe} 
\int_{\Gamma_*^e} f(N^1,N^2) d \Gamma 
\end{equation}
where $\Gamma_*^e=\Gamma_* \cap \Omega_e^{2,b}$ and $\{\Omega_e^{2,b}\}_{1}^{nbe}$ denotes elements in $\Omega^2$
that intersect with $\Gamma_*$.
What makes hierarchical meshes attractive is that for a fine element on $\Gamma_*$ one knows
the element in the coarse mesh that locates the other side of the interface.

For the elemental interface integral, a Gauss quadrature rule for line 
elements is adopted. For example, two GPs are used for bilinear elements. Let the
GPs denoted by  $\{\xi_i\}_{i=1}^{ngp}$. These GPs have to be mapped to two parent elements--
one associated with $\Omega_e^{2,b}$ and  one associated with $\Omega_e^{1,b}$.
That is given $\xi_i$, one has to solve for $\bsym{\xi}_i^2$ and $\bsym{\xi}_i^1$
($\bsym{\xi}_i^2=(\xi_i^2,\eta_i^2)$)

\begin{equation}
\begin{split}
\vm{x}_i &= \vm{M}(\xi_i)\vm{x}_l \\
\vm{x}_i &= \vm{N}^2(\bsym{\xi}_i^2)\vm{x}_e^2 \rightarrow \bsym{\xi}_i^2\\
\vm{x}_i &= \vm{N}^1(\bsym{\xi}_i^1)\vm{x}_e^1 \rightarrow \bsym{\xi}_i^1
\end{split}
\end{equation}
where the first equation is used to compute the global coordinates of the GP
and the second and third equations are used to compute the natural coordinates
of the GP in the parent element associated with $\Omega_e^{k,b}$. In the above,
$\vm{M}$ denotes the vector of shape functions of a two-noded line element; $\vm{x}_l$
are the nodal coordinates of two boundary nodes of $\Gamma_*^e$;
$\vm{x}_e^k$ ($k=1,2$) denotes the nodal coordinates of $\Omega_e^{k,b}$. 
For the example given in Fig.~\ref{fig:domain-mesh-hier}, $\vm{x}_e^2$ stores the
coordinates of nodes 5,7,9 and 6. And, $\vm{x}_e^1$ stores the coordinates of nodes 10,22,20 and 16.

It is now ready to evaluate the interfacial integral as

\begin{equation}
\int_{\Gamma_*^e} f(N^1,N^2) d \Gamma = \sum_{i=1}^{ngp} f(N^1(\bsym{\xi}_i^1),N^2(
\bsym{\xi}_i^2) ) w_i
\end{equation}
where $w_i$ equals the weight multiplied with the Jacobian of the transformation from
the line parent element $[-1,1]$ to $\Gamma_*^e$.

Finally the coupling terms are assembled to the global stiffness matrix in a standard manner.
For example $\vm{K}^{d,11}$ is assembled using the connectivity of $\Omega_e^{1,b}$
and $\vm{K}^{d,22}$ is assembled using the connectivity of $\Omega_e^{2,b}$.

%--------------------------------------------
\subsubsection{Non-matching structured meshes}

Non-matching structured meshes are plotted in Fig.~\ref{fig:domain-mesh5}.
In those cases, the evaluation of the interfacial integrals are more complicated
for there might be two elements in $\Omega^1$ interacting with one element in $\Omega^2$.
For this case and especially for unstructured meshes, an efficient implementation of Nitsche's method requires
(1) a list of interface elements $\Omega^{k,b}$ and (2) an algorithm that gives, for a given
$\Omega_e^{1,b}$, elements on $\Omega^{2,b}$ that it interacts with. The latter can be implemented using
the concept of bounding boxes. Since our framework is based on NURBS-based IGA which involve only 
structured meshes, (1) is straightforward.

\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{domain5}
  \caption{Independent discretisations of the domains: non-matching structured meshes.}
  \label{fig:domain-mesh5}
\end{figure}

%-------------------------------------------
\subsubsection{Three dimensional formulations}

\begin{figure}
  \centering
  \includegraphics[width=0.55\textwidth]{3d-3d-coupling}
  \caption{Coupling of two three dimensional continuum models. For evaluating the coupling terms,
  the coupling surface $\Omega_*$ is discretised by $\{\Omega_*^e\}_{e=1}^{nbe}$ which are the faces
  of solid elements that have non-zero intersection with $\Omega_*$.}
  \label{fig:3d-3d-coupling}
\end{figure}

This section presents the implementation for 3D, we refer to Fig.~\ref{fig:3d-3d-coupling}.
We confine to hierarchical meshes to simplify the implementation. The interface is now 
a surface, hence we denote it by $\Omega_*$ rather than $\Gamma_*$. Coupling terms 
are computed according to

\begin{equation}
\int_{\Omega_*} f(N^1,N^2) d \Omega = \bigcup_{e=1}^{nbe} 
\int_{\Omega_*^e} f(N^1,N^2) d \Omega =
\bigcup_{e=1}^{nbe}
\sum_{i=1}^{ngp} f(N^1(\bsym{\xi}_i^1),N^2( \bsym{\xi}_i^2) ) \bar{w}_i
\label{eq:3d-3d-gauss}
\end{equation}
where 2D quadrature was used in the second equality of which the Gauss points are
determined using the following equations provided $\{(\xi_i,\eta_i,w_i)\}_{i=1}^{ngp}$
are the Gauss points defined on $[-1,1]\times[1,1]$     

\begin{alignat}{2}
\label{eq:hang}
\vm{x}_i &= \vm{M}(\xi_i,\eta_i)\vm{x}_l \\
\vm{x}_i &= \vm{N}^2(\xi_i^2,\eta_i^2,\zeta_i^2)\vm{x}_e^2 \rightarrow (\xi_i^2,\eta_i^2,\zeta_i^2) \label{eq:hang1}\\
\vm{x}_i &= \vm{N}^1 (\xi_i^1,\eta_i^1,\zeta_i^1) \vm{x}_e^1 \rightarrow (\xi_i^1,\eta_i^1,\zeta_i^1) \label{eq:hang2}
\end{alignat}
where $\vm{x}_l$ denotes the nodal coordinates of face $\Omega_*^e$.
In Equation~\eqref{eq:3d-3d-gauss}, $\bar{w}_i=w_i\norm{\vm{a}_1\times\vm{a}_2}$ with $\vm{a}_\alpha$ denote the tangent vectors
of $\Omega_*$ defined by

\begin{equation}
\vm{a}_1 = \vm{M}_{,\xi} \vm{x}_l,\quad
\vm{a}_2 = \vm{M}_{,\eta} \vm{x}_l
\label{eq:tangents}
\end{equation}
Equations~\eqref{eq:hang1} and \eqref{eq:hang2} determine the local coordinates of $\vm{x}_i$
in the parent elements $[-1,-1,-1]\times[1,1,1]$. Usually a Newton-Raphson method is used for this.

% ------------------------------
\subsection{Stabilisation parameter}

The penalty-like term in the weak form restores coercivity (in the continuous problem), or positive definiteness 
(in the discrete problem) to the system of equations. A stabilisation parameter that is too small produces unstable results, whereas a stabilisation parameter  that is too large essentially forces the problem to behave
as a penalty method, which can produce the same locking problems as typical Lagrange multipliers.
There exists techniques to determine the optimal value of the stabilisation parameter $\alpha$.
The idea is based on numerical analysis of the discrete forms and lead to the global \cite{Griebel}
and local generalized eigenvalue approaches \cite{embar_imposing_2010}.

It can be shown that for linear elasticity $\alpha$ depends on the element size $h_e$ and the material parameters, see for example \cite{fritz-nitsche-mortar,Bazilevs200712}

\begin{equation}
\alpha = \frac{\lambda + \mu}{2} \frac{\theta(p)}{h_e}
\label{eq:alpha-estimation}
\end{equation}
where $\theta(p)$ is a positive number that depends only on the polynomial order $p$ of the finite element
approximation. %Equation~\eqref{eq:alpha-estimation}.

%--------------------------------------------
\subsection{Non symmetric Nitsche's method}

In \cite{Oden:1998:DHF:301787.301788} a non symmetric discontinuous Galerkin for elliptic problems
was proposed.
This method does not contain the stabilizing (penalty) term as the classical symmetric Nitsche's method.
Stability and error estimate of this non symmetric method was studied by \cite{Larson_analysisof}.
In a non symmetric Nitsche's method the diffusive term is skipped and the
sign of the term which was added to make the method symmetric is changed. The
resulting weak form is thus given by 

\begin{equation}
\int_\Omega \frac{\partial \delta u_i}{\partial x_j}   \sigma_{ij}  \mathrm{d}\Omega -
\int_{\Gamma_*}  \jump{\delta u_i} \{\sigma_{ij}\} n_j^1 \mathrm{d}\Gamma +
\int_{\Gamma_*}  \jump{u_i} \{\delta \sigma_{ij}\} n_j^1 \mathrm{d}\Gamma
=  \int_\Omega \delta u_i  b_i \mathrm{d}\Omega+ 
\int_{\Gamma_t}  \delta u_i \bar{t}_i  \mathrm{d}\Gamma 
\label{eq:dg-weakform-non-symmetry}
\end{equation}
which leads to the interfacial stiffness matrix that reads

\begin{equation}
\vm{K}^d = \vm{K}^d_{dg} - (\vm{K}^d_{dg})\trans 
\end{equation}
in which the components are given by
\begin{equation}
\vm{K}^d_{dg} = \begin{bmatrix}
-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
-\D\int_{\Gamma_*} \vm{N}^{1\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma \\
\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^1\vm{B}^1 \mathrm{d}\Gamma &
\D\int_{\Gamma_*} \vm{N}^{2\text{T}} \vm{n} \frac{1}{2}\vm{C}^2\vm{B}^2 \mathrm{d}\Gamma 
\end{bmatrix}\label{eq:nitsche-kdg-non}
\end{equation}

%-----------------------------------------
\subsection{Extension to NURBS elements}

NURBS are a commonly used basis function in isogeometric analysis 
\cite{hughes_isogeometric_2005,cottrel_book_2009}.
In general an object is represented by multiple NURBS patches. Connecting those patches of which
parametrisation is not compatible at the interface is not a trivial task. 
The Nitche's method as previously described can be used to join non-conforming NURBS patches.
In this case the shape functions $N_I^k(\bar{\bsym{\xi}})$ are the NURBS basis functions of which
details can be found elsewhere \cite{hughes_isogeometric_2005}.
Note that NURBS basis are defined in a parameter space which is different from the parent space
where Lagrange polynomials are defined. Therefore, for IGA there is a slight modification to
Equations~\eqref{eq:hang1} and \eqref{eq:hang2}. 




%-------------------------------------
\section{Beam-continuum coupling}\label{sec:beam-continuum}
%-------------------------------------

We present a Nitsche based method to couple a beam and a 2D continuum, see Fig.~\ref{fig:beam-continuum}.
To make it concrete, we assume that the beam is modeled with the Timoshenko beam theory.
Formulation for the Euler-Bernoulli beam can be obtained using the presented ideas. Only treatment of the coupling
matrices is presented since beam stiffness matrices are standard and can be found for example 
in \cite{reddy_book}. The idea is  to couple the continuum and the real beam $\Omega^\text{beam}$. 
Therefore the formulation presented in Section \ref{continua-continua} can be used.
The coupling interface (line) $\Gamma_*=\Omega^\text{beam}\cap\Omega^\text{solid}_b$.
Firstly, the formulation is presented for the case in which the beam coordinate system is
identical to the global one. Generalization is straightforward and not presented here.
We emphasize that,  cf. Fig.~\ref{fig:beam-continuum}, the unknowns at beam nodes and continuum nodes are 
three and two, respectively.

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{beam-2d-coupling}
  \caption{Coupling a beam and a continuum using a Nitsche based method.}
  \label{fig:beam-continuum}
\end{figure}

The displacement field of the beam is given by (we will use the superscripts $s$ and $b$
      for solid and beam respectively)

\begin{equation}
\begin{split}
u_x^b &= u(x) -y\theta(x)\\
u_y^b &= w(x)
\end{split}
\end{equation}
where $u,w$ and $\theta$ are the axial, transverse displacement and the rotation of the beam mid-line, 
respectively.
The strain field is therefore given by

\begin{equation}
\epsilon_{xx}^b  = u_{,x}-y \theta_{,x}, \quad
\epsilon_{yy}^b  = 0, \quad
2\epsilon_{xy}^b = -\theta + w_{,x}
\end{equation}
The stresses are defined as
\begin{equation}
\begin{bmatrix}
\sigma_{xx}^b\\
\sigma_{yy}^b\\
\sigma_{xy}^b
\end{bmatrix}=
\begin{bmatrix}
E &  0 & 0\\
0 &  0 & 0\\
0 & 0 & kG  
\end{bmatrix}
\begin{bmatrix}
\epsilon_{xx}^b\\
\epsilon_{yy}^b\\
2\epsilon_{xy}^b
\end{bmatrix} \equiv \vm{C}^\text{b}\bsym{\epsilon}^\text{b}
\end{equation}
where $k$ denotes the shear correction factor (a value of $5/6$ is usually used) and $G$ is the shear modulus.


For sake of simplicity, assume that the beam is discretised by two-noded
linear elements. In this case, the displacements of any beam element $e$ are given by

\begin{equation}
\begin{bmatrix}
u_{x}^b\\
u_{y}^b
\end{bmatrix}=\underbrace{
\begin{bmatrix}
N_1 & 0 & -y N_{1} & N_2 & 0 & - y N_{2}\\
0 & N_{1} & 0 & 0 & N_{2} & 0     
\end{bmatrix}}_{\vm{N}^b}
\begin{bmatrix}
u_1\\w_1\\\theta_1\\u_2\\w_2\\ \theta_2
\end{bmatrix}\label{eq:beam-Nb}
\end{equation}
where $(u_I,w_I,\theta_I)$ are the nodal unknowns at node $I$ of a beam element; $N_I(\xi)$ are the
univariate shape functions, and the strains are given by

\begin{equation}
\bsym{\epsilon}^b=
\begin{bmatrix}
\epsilon_{xx}^b\\
\epsilon_{yy}^b\\
2\epsilon_{xy}^b
\end{bmatrix}=\underbrace{
\begin{bmatrix}
N_{1,x} & 0 & -y N_{1,x} & N_{2,x} & 0 & - y N_{2,x}\\
   0 & 0 & 0 & 0 & 0 & 0\\
0 & N_{1,x} & -N_1 & 0 & N_{2,x} & - N_2     
\end{bmatrix}}_{\vm{B}^b}
\begin{bmatrix}
u_1\\w_1\\\theta_1\\u_2\\w_2\\ \theta_2
\end{bmatrix}
\end{equation}

The coupling matrices given in Equations ~\eqref{eq:nitsche-kdg} and \eqref{eq:nitsche-kpe}
are now given by

\begin{equation}
\vm{K}^d_{dg} = \begin{bmatrix}
-\D\int_{\Gamma_*} (\vm{N}^{s})\trans \vm{n} \frac{1}{2}\vm{C}^s\vm{B}^s \mathrm{d}\Gamma &
-\D\int_{\Gamma_*} (\vm{N}^{s})\trans \vm{n} \frac{1}{2}\vm{C}^b\vm{B}^b \mathrm{d}\Gamma \\
\D\int_{\Gamma_*}  (\vm{N}^{b})\trans  \vm{n} \frac{1}{2}\vm{C}^s\vm{B}^s \mathrm{d}\Gamma &
\D\int_{\Gamma_*}  (\vm{N}^{b})\trans   \vm{n} \frac{1}{2}\vm{C}^{b}\vm{B}^{b} 
\mathrm{d}\Gamma 
\end{bmatrix}
\label{eq:2D-beam-stiffness1}
\end{equation}
and by
\begin{equation}
\vm{K}^d_{pe} = \begin{bmatrix}
   \D\int_{\Gamma_*}  \alpha  (\vm{N}^{s})\trans \vm{N}^s \mathrm{d}\Gamma &
 - \D\int_{\Gamma_*}  \alpha  (\vm{N}^{s})\trans \vm{N}^{b} \mathrm{d}\Gamma\\
 - \D\int_{\Gamma_*}  \alpha  (\vm{N}^{b})\trans \vm{N}^s \mathrm{d}\Gamma&
   \D\int_{\Gamma_*}  \alpha  (\vm{N}^{b})\trans \vm{N}^{b} \mathrm{d}\Gamma
\end{bmatrix}
\label{eq:2D-beam-stiffness2}
\end{equation}
where $\vm{C}^s$, $\vm{B}^s$ and $\vm{N}^s$ denote the constitutive matrix, the strain-displacement
matrix and the shape function matrix of the plane continuum elements, respectively.

%-----------------------------------
\subsection{Stabilisation parameter}

\begin{equation}
a(\vm{u},\vm{u})=\tilde{a}(\vm{u},\vm{u}) + \alpha \int_{\Gamma_*} \jump{\vm{u}}\jump{\vm{u}}d\Gamma -  
\int_{\Gamma_*} \jump{\vm{u}}\left[ \vm{D}^c \bsym{\epsilon}^c \vm{n}^c + \vm{D}^b \bsym{\epsilon}^b
\vm{n}^c \right]d\Gamma
\label{eq:eq1}
\end{equation}

For any scalar product and any parameter $\epsilon$, one has

\begin{equation}
\norm{\vm{u}-\epsilon\vm{v}}^2=\norm{\vm{u}}^2 + \epsilon^2 \norm{\vm{v}}^2 - 2\epsilon (\vm{u},\vm{v})
\end{equation}
where $\norm{\cdot}$ is the norm with $(\cdot,\cdot)$ \ie $\norm{\cdot}=\sqrt{(\cdot,\cdot)}$. From this equation,one can write
\begin{equation}
(\vm{u},\vm{v}) = \frac{1}{2\epsilon}\norm{\vm{u}}^2 + \frac{\epsilon}{2}\norm{\vm{v}}^2 - \norm{\vm{u}-\epsilon\vm{v}}^2 \le \frac{1}{2\epsilon}\norm{\vm{u}}^2 + \frac{\epsilon}{2}\norm{\vm{v}}^2
\label{eq:bdt1}
\end{equation}
If we apply Equation~\eqref{eq:bdt1} to $(\jump{\vm{u}},\bar{\vm{t}})_{\Gamma_*}$, we have

\begin{equation}
 (\jump{\vm{u}},\bar{\vm{t}})_{\Gamma_*} \le
\frac{1}{2\epsilon} \norm{\jump{\vm{u}}}^2 + \frac{\epsilon}{2} \norm{\bar{\vm{t}}}^2
\label{eq:eq2}
\end{equation}

Substituting Equation~\eqref{eq:eq2} into Equation~\eqref{eq:eq1} leads to the following

\begin{equation}
a(\vm{u},\vm{u}) \ge \tilde{a}(\vm{u},\vm{u}) + \left(\alpha -\frac{1}{2\epsilon}\right) 
    \norm{\jump{\vm{u}}}^2 -  \frac{\epsilon}{2}\norm{\bar{\vm{t}}}^2
\label{eq:eq3}
\end{equation}
The main idea of the method is to suppose the existence of a constant $C$ such that
\begin{equation}
\norm{\bar{\vm{t}}}^2 \le C^2 \tilde{a}(\vm{u},\vm{u})
\label{eq:eq4}
\end{equation}
which means that the flux are bounded by the gradient inside the domain.

Substituting Equation~\eqref{eq:eq4} into Equation~\eqref{eq:eq3} leads to the following

\begin{equation}
a(\vm{u},\vm{u}) \ge \left( 1 - \frac{C^2 \epsilon}{2}\right) \tilde{a}(\vm{u},\vm{u}) + 
\left(\alpha -\frac{1}{2\epsilon}\right) \norm{\jump{\vm{u}}}^2
\label{eq:eq5}
\end{equation}

\begin{equation}
 1 - \frac{C^2 \epsilon}{2} \ge 0, \quad \alpha -\frac{1}{2\epsilon} \ge 0
\label{eq:eq6}
\end{equation}
If we take $\epsilon=1/C^2$, the first inequality in Equation~\eqref{eq:eq6} will be satisfied.
Therefore the bilinear form $a$ is positive definite if $\alpha > C^2/2$. In what follows we present
a way to determine the constant C satisfying Equation~\eqref{eq:eq4}. This is where one has to solve
an eigenvalue problem.

\subsubsection{Finding C}

\begin{equation}
\int_{\Gamma_*} \left[ \vm{D}^c \bsym{\epsilon}^c \vm{n}^c + \vm{D}^b \bsym{\epsilon}^b \vm{n}^c 
\right]^2 d\Gamma \le C^2 
\int_\Omega \bsym{\epsilon} \vm{D} \bsym{\epsilon} d \Omega
\end{equation}
Equation~\eqref{eq:eq4} becomes
\begin{equation}
\frac{\vm{u}\trans \vm{H} \vm{u}}{\vm{u}\trans \tilde{\vm{K}} \vm{u}} \le C^2 
\end{equation}
So we need to find the maximum of the Rayleigh quotient $R$ defined as $R=\frac{\vm{u}\trans \vm{H} \vm{u}}{\vm{u}\trans \tilde{\vm{K}} \vm{u}}$. By theorem this maximum is $\lambda_1$, the largest eigenvalue of 
matrix $\tilde{\vm{K}}^{-1}\vm{H}$. Finally 

\begin{equation}
\alpha = \frac{\lambda_1}{2}
\end{equation}

%------------------------------------------------------------
\section{Plate-continuum coupling}\label{sec:plate-continuum}
%--------------------------------------------------------------

In this section, we present the coupling of a plate and a 3D continuum, see Fig.~\ref{fig:plate-continuum}.
We present the formulation for the Reissner-Mindlin plate element with three unknowns at 
node--the deflection and the two rotations $(w_I,\beta_{1I},\beta_{2I})$ and
the Kirchoff plate elements with only one unknowns per node.
In the same manner with beam-continuum coupling, we couple
the 3D plate with the solid although only the mid-surface of the plate is discretised.
In what follows, we will use the superscripts $s$ and $p$ for solid and plate respectively.


\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{plate-3D-coupling}
  \caption{Coupling a plate and a 3D continuum using a Nitsche based method.}
  \label{fig:plate-continuum}
\end{figure}

%------------------------------------
\subsection{Reissner-Mindlin plate}

The displacement field of the plate is given by \cite{reddy_book}

\begin{equation}
\begin{split}
u_1(x_1,x_2,x_3) &= -x_3 \beta_1\\
u_2(x_1,x_2,x_3) &= -x_3 \beta_2\\
u_3(x_1,x_2,x_3) &= w(x_1,x_2)
\end{split}\label{eq:plate-displacement}
\end{equation}
where $w,\beta_1,\beta_2$ are the transverse displacement and rotations of the mid-plane
of the plate.

The strain field is then given by

\begin{equation}
\begin{split}
\epsilon_{11}  &= -x_3 \beta_{1,1}\\
\epsilon_{22}  &= -x_3 \beta_{2,2}\\
2\epsilon_{12} &= -x_3 (\beta_{1,2}+\beta_{2,1})\\
2\epsilon_{13} &= - \beta_{1} + w_{,1}\\
2\epsilon_{23} &= - \beta_{2} + w_{,2}\\
\end{split}\label{eq:plate-strain}
\end{equation}

The deflection $w$ and the rotation are approximated as follows

\begin{equation}
w = N_I(\xi,\eta) w_I,\quad
\beta_\alpha = N_I(\xi,\eta) \beta_{\alpha I}, \quad \alpha=1,2
\end{equation}

\noindent where $N_I(\xi,\eta)$ is the shape function associated with node I
and $(w_I,\beta_{1I},\beta_{2I})$ denote the nodal unknowns that include the nodal deflection and two rotations.

Hence, the displacement field given in Equation~\eqref{eq:plate-displacement} can be written as

\begin{equation}
\vm{u}^p = \begin{bmatrix}
0 & -x_3N_I & 0\\
0 & 0 & -x_3 N_I\\
N_I & 0 & 0 
\end{bmatrix}
\begin{bmatrix}
w\\ \beta_1 \\ \beta_2
\end{bmatrix}_I \equiv \vm{N}_I^p\vm{a}^p_I
\end{equation}
and the strain field given in Equation~\eqref{eq:plate-strain} as

\begin{equation}
\begin{bmatrix}
\epsilon_{11}\\
\epsilon_{22}\\
2\epsilon_{12}\\
2\epsilon_{23}\\
2\epsilon_{13}\\
\end{bmatrix}=
\begin{bmatrix}
0 & -x_3 N_{I,1} & 0 \\
0 & 0 & -x_3 N_{I,2}\\
0 & -x_3N_{I,2} & -x_3 N_{I,1}\\  
N_{I,2} & 0 & -N_I \\
N_{I,1} & -N_I & 0      
\end{bmatrix}
\begin{bmatrix}
w_I\\ \beta_{1I}\\ \beta_{2I}
\end{bmatrix}\equiv\vm{B}_{I}^p  \vm{a}_I^p
\label{plate-curvature1}
\end{equation}
And finally the stress field is given by
\begin{equation}
\begin{bmatrix}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{12}\\
\sigma_{23}\\
\sigma_{13}\\
\end{bmatrix}=
\begin{bmatrix}
\frac{E}{(1-\nu^2)}\begin{bmatrix}
1 & \nu & 0\\
\nu & 1 & 0\\
0 & 0 & 0.5(1-\nu)
\end{bmatrix} & 
\begin{bmatrix}
0 & 0 \\
0 & 0 \\
0 & 0 
\end{bmatrix}\\
\begin{bmatrix}
0 & 0 & 0\\
0 & 0 & 0
\end{bmatrix} &
\frac{k E}{2(1+\nu)}\begin{bmatrix}
1 & 0\\
0 & 1
\end{bmatrix}
\end{bmatrix}
\begin{bmatrix}
\epsilon_{11}\\
\epsilon_{22}\\
2\epsilon_{12}\\
2\epsilon_{23}\\
2\epsilon_{13}\\
\end{bmatrix} \equiv \vm{C}^p\bsym{\epsilon}^p
\label{plate-curvature1}
\end{equation}
where $k$ is the shear correction factor.

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.65\textwidth]{plate-3D-Gauss}
  \caption{Coupling a plate and a 3D continuum using a Nitsche based method: determination of 
           Gauss points for evaluating the coupling terms.}
  \label{fig:plate-continuum-gauss}
\end{figure}

The coupling matrices are identical with Equations~\eqref{eq:2D-beam-stiffness1}
and \eqref{eq:2D-beam-stiffness2} by replacing $\vm{N}^b,\vm{B}^b,\vm{C}^b$ by  $\vm{N}^p,\vm{B}^p,\vm{C}^p$, 
respectively. Note that the elasticity matrix of the solid $\vm{C}^s$ is now a $5\times6$ matrix in which
the third row that corresponds with $\sigma_{33}$ is removed. Similarly the normal vector $\vm{n}_{3D}$ given
in Equation~\eqref{eq:n-matrix} is reduced to a $3\times5$ matrix where the third column was removed.
Coupling terms are computed as follows given $(\xi_i,\eta_i,w_i)$--the Gauss points

\begin{equation}
\int_{\Omega_*} f(N^s,N^p) d \Omega = \bigcup_{e=1}^{nbe} 
\int_{\Omega_*^e} f(N^s,N^p) d \Omega =
\bigcup_{e=1}^{nbe}
\sum_{i=1}^{ngp} f(N^s(\bsym{\xi}_i^s),N^p( \bsym{\xi}_i^p) ) \bar{w}_i
\end{equation}
where $\bsym{\xi}_i^s=(\xi,\eta,\zeta)_i^s$ and $\bsym{\xi}_i^p=(\xi,\eta)_i^p$ are
the Gauss point defined in the unit cube and unit square that are associated with the
solid and plate element, cf. Fig.~\ref{fig:plate-continuum-gauss}.
They are computed using the following equations
\begin{equation}
\begin{split}
\vm{x}_i &= \vm{M}(\xi_i,\eta_i)\vm{x}_l \\
\vm{x}_i &= \vm{N}^s(\xi_i^s,\eta_i^s,\zeta_i^s)\vm{x}_e^s \rightarrow (\xi_i^s,\eta_i^s,\zeta_i^s) \\
\vm{x}_i^* &= \vm{N}^p (\xi_i^p,\eta_i^p) \vm{x}_e^p \rightarrow (\xi_i^p,\eta_i^p)
\end{split}
\end{equation}
with $\bar{w}_i = w_i \norm{\vm{a}_1 \times \vm{a}_2}$ where $\vm{a}_\alpha$ are the tangent
vectors defined in Equation~\eqref{eq:tangents}. Again $\vm{x}_i$ denotes the global coordinates
of the Gauss point $i$ and $\vm{x}_i^*$ is its projection on the plate mid-plane.

%--------------------------
\subsection{Kirchhoff plate}

According to the Kirchhoff plate theory, the displacement field is uniquely described in terms
of the transverse displacement of the mid-plane

\begin{equation}
\begin{split}
u_1(x_1,x_2,x_3) &= -x_3 \pderiv{w}{x_1}=-x_3w_{,1}\\
u_2(x_1,x_2,x_3) &= -x_3 \pderiv{w}{x_2}=-x_3w_{,2}\\
u_3(x_1,x_2,x_3) &= w(x_1,x_2)
\end{split}
\label{eq:Kirchhoff-disp}
\end{equation}
The strain field is then given by

\begin{equation}
\begin{split}
\epsilon_{11}  &= -x_3 \frac{\partial^2 w}{\partial x_1^2} =-x_3w_{,11}\\
\epsilon_{22}  &= -x_3 \frac{\partial^2 w}{\partial x_2^2} =-x_3w_{,22}\\
2\epsilon_{12} &= -2x_3 \frac{\partial^2 w}{\partial x_1 x_2} =-2x_3w_{,12}\\
\end{split}
\label{eq:Kirchhoff-strain}
\end{equation}
The deflection $w$ is approximated as follows

\begin{equation}
w = N_I(\xi,\eta) w_I
\end{equation}

\noindent where $N_I(\xi,\eta)$ is the NURBS basis function associated with node I
and $w_I$ is the nodal deflection--the only unknown per node/control point. We note that any
smooth basis functions of which $C^1$ continuity can be straightforwardly constructed can be used.
Besides NURBS, subdivision surfaces can be adopted as well see \eg \cite{Cirak_2000}.

The displacement field in Equation~\eqref{eq:Kirchhoff-disp}
is then given by

\begin{equation}
\vm{u}^p = \begin{bmatrix}
-x_3 N_{I,1}\\
-x_3 N_{I,2}\\
N_I 
\end{bmatrix}
\begin{bmatrix}
w_I
\end{bmatrix}\equiv\vm{N}_I^p\vm{a}^p_I
\end{equation}
and the strain field given in Equation~\eqref{eq:Kirchhoff-strain} as

\begin{equation}
\begin{bmatrix}
\epsilon_{11}\\
\epsilon_{22}\\
2\epsilon_{12}\\
\end{bmatrix}=
\begin{bmatrix}
-x_3 N_{I,11}\\
-x_3 N_{I,22}\\
-2x_3 N_{I,12}  
\end{bmatrix}
\begin{bmatrix}
w_I
\end{bmatrix}\equiv\vm{B}_{I}^p  \vm{a}_I^p
\end{equation}
And finally the stress field is given by
\begin{equation}
\begin{bmatrix}
\sigma_{11}\\
\sigma_{22}\\
\sigma_{12}\\
\end{bmatrix}=
\frac{E}{(1-\nu^2)}\begin{bmatrix}
1 & \nu & 0\\
\nu & 1 & 0\\
0 & 0 & 0.5(1-\nu)
\end{bmatrix}
\begin{bmatrix}
\epsilon_{11}\\
\epsilon_{22}\\
2\epsilon_{12}\\
\end{bmatrix}=\vm{C}^p\bsym{\epsilon}^p
\end{equation}
Note that the elasticity matrix of the solid $\vm{C}^s$ is now a $3\times6$ matrix in which
the rows that correspond with $\sigma_{33},\sigma_{13}, \sigma_{23}$ are removed.



%------------------------------------------------
\section{Numerical examples}\label{sec:examples}
%------------------------------------------------

In accordance to the three topics presented previously, in this section, numerical examples
on discretisation coupling, beam/continuum and plate/continuum coupling are presented to
demonstrate the performance of the proposed formulation. Specially, the following examples
are provided together with their raison d'etre.

\begin{enumerate}
\item Discretisation coupling
   \begin{enumerate}
   \item Timoshenko beam (2D/2D coupling)
   \item Cantilever beam (3D/3D coupling)
   \end{enumerate}
\item Beam-continuum coupling
   \begin{enumerate}
   \item Cantilever beam (2D/beam coupling): convergence analysis, stabilisation, non-symmetric Nitsche...
   \end{enumerate}
\item Plate-continuum coupling
   \begin{enumerate}
   \item Cantilever beam (3D/plate coupling)
   \item Composite laminate
   \end{enumerate}
\end{enumerate}

Unless otherwise stated, we use 
MIGFEM--an open source Matlab IGA code\footnote{which is available at \url{https://sourceforge.net/projects/cmcodes/}} for our computations and  the Timoshenko beam theory and the Reissner-Mindlin plate theory.
The visualisation was performed in Paraview \cite{para}.
In examples of mixed dimensional analyses domains where reduced models such as beams and plates
are utilized are assumed to be known \textit{a priori}. For all considered examples, 
the symmetric and non symmetric versions of the Nitsche's method are both utilized.
%We use a C++ code \cite{jemjive} for computations.

%-------------------------------------
\subsection{Discretisation coupling}

% ------------------------------
\subsubsection{Timoshenko beam}

Consider a beam of dimensions $L \times D$ (unit thickness), subjected to a
parabolic traction at the free end as shown in Fig.~\ref{fig:beam-geo}.
The beam is considered to be of unit depth and is in plane stress
state. The parabolic traction is given by

\begin{equation}
    t_y(y) = -\frac{P}{2I} \biggl ( \frac{D^2}{4} - y^2 \biggr)\label{eq:ty}
\end{equation}

\noindent where $I = D^3/12$ is the moment of inertia. The exact displacement
field of this problem is \cite{elasticity_book}

\begin{equation}
\begin{split}
    u_x(x,y) &=  \frac{Py}{6EI} \biggl [ (6L-3x)x + (2+\nu)\biggl(y^2-\frac{D^2}{4}\biggr) \biggr] \\
    u_y(x,y) &=  - \frac{P}{6EI} \biggl [ 3\nu y^2(L-x) + (4+5\nu)\frac{D^2x}{4} +(3L-x)x^2 \biggr] \\
\end{split}
\label{eq:tBeamExactDisp}
\end{equation}

\noindent and the exact stresses are

\begin{equation}
	\sigma_{xx}(x,y) =  \frac{P(L-x)y}{I}; 
	\quad \sigma_{yy}(x,y) =  0, \quad
\sigma_{xy}(x,y) =  -\frac{P}{2I} \biggl ( \frac{D^2}{4}-y^2\biggr)
\end{equation}

\noindent In the computations, material properties are taken as $E=
3.0 \times 10^7$, $\nu = 0.3$ and the beam dimensions are $D=6$ and
$L=48$. The shear force is $P = 1000$.  In order to model the clamping condition, 
the displacement defined by Equation~\eqref{eq:tBeamExactDisp} is prescribed as essential boundary 
conditions at $x=0, -D/2 \le y \le D/2$. This problem is solved with bilinear Lagrange elements
and high order B-splines elements. The former helps to verify the implementation in addition to
the ease of enforcement of Dirichlet boundary conditions (BCs). For the latter, care must be taken
in enforcing the Dirichlet BCs given in Equation~\eqref{eq:tBeamExactDisp} since the B-splines
are not interpolatory. The beam is divided into two domains by a vertical line at $x=L/2$.\\

\begin{figure}[htbp]
  \centering 
  \psfrag{p}{P}\psfrag{l}{$L$}\psfrag{d}{$D$}\psfrag{x}{$x$}\psfrag{y}{$y$}
  \includegraphics[width=0.5\textwidth]{beam}
  \caption{Timoshenko beam.}
  \label{fig:beam-geo}
\end{figure}

\noindent \textbf{Lagrange elements} Firstly, a conforming mesh is considered and each domain is
discretised by a mesh of $20\times4$ elements as given in Fig.~\ref{fig:beam-meshes}a. Then, a non-conforming
mesh where the left domain is discretised by $20\times8$ elements and
the right domain is meshed by $20\times4$ is considered, cf. Fig.~\ref{fig:beam-meshes}b. 
The vertical displacements are plotted in Fig.~\ref{fig:beam-disp} together with the exact solution
and the stresses are plotted in Fig.~\ref{fig:beam-stress}.

\begin{figure}[htbp]
  \centering 
   \subfloat[Conforming mesh]{\includegraphics[width=0.45\textwidth]{tbeam-mesh1}}\\
   \subfloat[Non conforming mesh]{\includegraphics[width=0.45\textwidth]{tbeam-mesh2}}
  \caption{Timoshenko beam: conforming and non-conforming meshes. Note that even with the conforming mesh,
  there are double nodes at the coupling interface $x=L/2, -D/2 \le y \le D/2$.}
  \label{fig:beam-meshes}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \subfloat[Conforming mesh]{\includegraphics[width=0.45\textwidth]{2D-2D-conforming-disp}}
   \subfloat[Conforming mesh]{\includegraphics[width=0.45\textwidth]{2D-2D-nonconforming-disp}}
  \caption{Timoshenko beam.}
  \label{fig:beam-disp}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \subfloat[Stresses along the beam length]{\includegraphics[width=0.45\textwidth]{2D-2D-stress-length}}
   \subfloat[Stresses over the beam height]{\includegraphics[width=0.45\textwidth]{2D-2D-stress-width}}
  \caption{Timoshenko beam: stresses conforming mesh.}
  \label{fig:beam-stress}
\end{figure}

\noindent \textbf{B-splines elements} Next, we study the performance of the B-splines elements of
which one mesh is given in Fig.~\ref{fig:beam-bspline-mesh}. 
Dirichlet BCs are enforced using the least square projection method see \eg \cite{nguyen_iga_review}.
As detailed in \cite{hughes-fem-book} for Lagrangian basis functions,
a rule of $(p+1)\times(q+1)$ Gaussian quadrature can be applied for
two-dimensional elements in which $p$ and $q$ denote the orders of
the chosen basis functions in the $\xi$ and $\eta$ direction. The same
procedure is also used for NURBS basis functions in the present work,
although it should be emphasised that Gaussian quadrature is not optimal for IGA
\cite{hughes_efficient_2010,Auricchio201215}.

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{2D-2D-Bspline-mesh}
  \caption{Timoshenko beam: B-spline cubic mesh with $4\times4$ elements for the left
  domain and $2\times2$ elements for the right one. The filled circles denote the control points.}
  \label{fig:beam-bspline-mesh}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.45\textwidth]{2D-2D-Bspline-stress-length}
   \includegraphics[width=0.45\textwidth]{2D-2D-Bspline-stress}
  \caption{Timoshenko beam: stresses with B-splines elements. The left domain
  is meshed by $8\times8$ cubic elements and the right domain with $2\times2$ cubic elements.}
  \label{fig:beam-geo}
\end{figure}

Finally we present results obtained with a non-hierarchical B-spline mesh as given in 
Fig.~\ref{fig:beam-bspline-general}. When evaluating the coupling terms, we use the \textit{FindSpan}
algorithm--which is a commonly used algorithm in CAD.


\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{tbeam-general}
  \caption{Timoshenko beam: non-hierarchical B-spline mesh ($8\times6$ cubic elements
        for the left domain and $4\times3$ cubic elements for the right domain.}
  \label{fig:beam-bspline-general}
\end{figure}

% -------------------------------
\subsubsection{3D-3D coupling}
% -------------------------------

In order to test the implementation for 3D problems, we consider the 3D cantilever beam shown in 
Fig.~\ref{fig:beam3D-geo}.
The data are: Young modulus $E=1000$, Poisson ratio $\nu=0.3$, $L=10$, $W=H=1$ and the imposed
displacement in the $z$-direction is $1$. The non-conforming 
B-splines discretisation is given in Fig.~\ref{fig:beam3D-mesh} where the beam is divided into two
equal parts.
In Fig.~\ref{fig:beam3D-res} the contour plot of $\sigma_{xx}$ is given and a comparison was made with
a standard Galerkin discretisation of $32\times4\times4$ tri-cubic B-splines elements.
The visualisation was performed in Paraview \cite{para} and a good agreement was obtained.

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.65\textwidth]{3d-3d-geo}
  \caption{A 3D cantilever beam subjected to an imposed vertical displacement.}
  \label{fig:beam3D-geo}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.65\textwidth]{nitsche-3D3D-spline-mesh}
  \caption{A 3D cantilever beam subjected to an imposed vertical displacement: $16\times4\times4$ tri-cubic
     B-splines elements for the left domain and  $16\times1\times2$ tri-cubic
     B-splines elements for the right domain.}
  \label{fig:beam3D-mesh}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \subfloat[Nitsche]{\includegraphics[width=0.55\textwidth]{nitsche-3D3D-spline-stress}}\\
   \subfloat[Galerkin]{\includegraphics[width=0.55\textwidth]{iga-3D3D-spline-stress}}
  \caption{Timoshenko beam.}
  \label{fig:beam3D-res}
\end{figure}

%--------------------------------------
\subsection{Continuum-beam coupling}

In this example, the Timoshenko beam problem as shown in Fig.~\ref{fig:beam-geo} is solved
using the continuum-beam coupling presented in Section \ref{sec:beam-continuum}. The mixed
continuum-beam model is given in Fig.~\ref{fig:beam-continuum-example}.
Geometry, material and boundary conditions are kept unaltered. The end shear force given in
Equation~\eqref{eq:ty} is now written as (for it is applied to a point)

\begin{equation}
    F = \int_{-D/2}^{D/2}t_y(y) dy = - \int_{-D/2}^{D/2} \frac{P}{2I} \biggl ( \frac{D^2}{4} - y^2 \biggr) dy
    = -\frac{PD^3}{12I}
\end{equation}

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{beam-coupling-example}
  \caption{Timoshenko beam:  mixed continuum-beam model.}
  \label{fig:beam-continuum-example}
\end{figure}

\noindent \textbf{Lagrange elements} In the first calculation we take $l_c=12$ and a mesh of $40\times10$ Q4 elements (40 elements in the length direction)  was used for the continuum part
and 29 two-noded elements for the beam part. 
Fig.~\ref{fig:beam-continuum-res} plots the transverse displacement (taken as nodal values) along the beam length at $y=0$ together with the exact solution given in Equation~\eqref{eq:tBeamExactDisp}. An excellent agreement 
with the exact solution can be observed.\\

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{beam-2d-result}
  \caption{Mixed dimensional analysis of the Timoshenko beam: comparison of numerical transverse displacements
     and exact solution.}
  \label{fig:beam-continuum-res}
\end{figure}


\begin{figure}[htbp]
  \centering 
  \subfloat[]{\includegraphics[width=0.5\textwidth]{beam-2D-stress1}}
  \subfloat[]{\includegraphics[width=0.5\textwidth]{beam-2D-stress2}}
  \caption{Cantilever beam subjects to an end shear force}
  \label{fig:beam-continuum-example-res}
\end{figure}

\noindent \textbf{B-splines elements} are used to discretise the continuum part (with bi-variate B-splines
      elements) and the beam part (with uni-variate B-splines elements).

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{beam-2D-spline-mesh}
  \caption{Mixed dimensional analysis of the Timoshenko beam: discretisation with B-splines elements.
  The continuum part is meshed by $2\times2$ cubic B-splines and the beam part is with 2 cubic elements.}
  \label{fig:beam-2D-spline-mesh}
\end{figure}

%--------------------------------------
\subsection{Continuum-plate coupling}

For verification of the continuum-plate coupling, we consider the 3D cantilever beam given in 
Fig.~\ref{fig:plate-continuum-geo}. The material properties are $E=1000$ N/mm$^2$, $\nu=0.3$.
The end shear traction is $\bar{t}=10$ N/mm in case of continuum-beam model and is
$\bar{t}=10/20$ N/mm$^2$ in case of continuum model which is referred to as the reference model.
We use B-splines elements to solve both the MDA and the reference model.
A mesh of $64\times4\times5$ tri-cubic elements is utilized for the reference model
and a mesh of $32\times4\times5$/ $16\times2$ cubic elements is utilized for the mixed dimensional models,
cf. Fig.~\ref{fig:plate-continuum-meshes}. The plate part of the mixed dimensional model is discretised using
the Reissner-Mindlin plate theory with three unknowns per node and the Kirchhoff plate theory with only one
unknown per node. Fig.~\ref{fig:plate-continuum-deform}.
  Fig.~\ref{fig:plate-continuum-stresses}.

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{plate-3D-example1}
  \caption{Cantilever beam subjects to an end shear force.}
  \label{fig:plate-continuum-geo}
\end{figure}


\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{plate-3D-bspline-mesh}
  \caption{Cantilever beam subjects to an end shear force.}
  \label{fig:plate-continuum-meshes}
\end{figure}

\begin{figure}[htbp]
  \centering 
   \includegraphics[width=0.55\textwidth]{plate-3D-meshes}
  \caption{Cantilever beam subjects to an end shear force.}
  \label{fig:plate-continuum-deform}
\end{figure}

\begin{figure}[htbp]
  \centering 
  \subfloat[reference model]{\includegraphics[width=0.5\textwidth]{plate-3D-iga-stress}}\\
  \subfloat[mixed dimensional model, Mindlin plate]{\includegraphics[width=0.5\textwidth]{plate-3D-nitsche-stress-mindlin}}
  \subfloat[mixed dimensional model, Kirchhoff plate]{\includegraphics[width=0.5\textwidth]{plate-3D-nitsche-stress-kirchhoff}}
  \caption{Cantilever beam subjects to an end shear force}
  \label{fig:plate-continuum-stresses}
\end{figure}

%-------------------------------------------
\section{Conclusions}\label{sec:conclusions}

We presented a Nitsche's method to couple (1) two and three dimensional non-matching
discretisatins, (2) two dimensional continua and beams and (3) three dimensional continua
and plates. A unified implementation of those coupling methods was given. 
Numerical examples using low order Lagrange finite elements and high order
NURBS isogeometric finite elements were provided. The first application is to glue non-compatible
NURBS patches in isogeometric analysis. The second application was in the context of mixed dimensional
analysis.

What have not been addressed in this manuscript  are (1) dynamics problems (2) finite deformation
problems (3) continua and shell coupling and (4) error estimation. 

% +++++++++++++++++++++++++++
\section*{Acknowledgements}

The authors would like to acknowledge the partial financial support of the
Framework Programme 7 Initial Training Network Funding under grant number
289361 ``Integrating Numerical Simulation and Geometric Design Technology".
St\'{e}phane Bordas also thanks partial funding for his time provided by
1) the EPSRC under grant EP/G042705/1 Increased Reliability for Industrially
Relevant Automatic Crack Growth Simulation with the eXtended Finite Element
Method and 2) the European Research Council Starting Independent Research
Grant (ERC Stg grant agreement No. 279578) entitled ``Towards real time multiscale simulation of cutting in
non-linear materials with applications to surgical simulation and computer
guided surgery''.
The authors would like to express the gratitude towards Drs. Erik Jan Lingen
and Martijn Stroeven at the Dynaflow Research Group, Houtsingel 95, 2719 EB Zoetermeer, The Netherlands 
for providing us the numerical toolkit jem/jive.


\bibliography{isogeometric}
%\bibliographystyle{elsarticle-num}
\bibliographystyle{unsrt}
%\bibliographystyle{plainnat}


% ============================================================================

\end{document}



