This document is aimed at explaining different available strategies to perform finite element computations on a periodic domain. To make our life a little simpler, this document is based on a two-dimensional problem where the unknowns at the nodes, $u$, are scalars. In this setting one can think of Poisson's equation
Using the finite element strategy of solving this in a weak form and an discretization using finite elements, one ends up with the following nodal equilibrium equations
\begin{equation}
\underline{f} ( \underline{u} ) = \underline{q}
\end{equation}
This equation is then usually linearized to obtain
For the simple finite element mesh shown in Figure~\ref{fig:mesh}(b), the system matrix, $\underline{\underline{A}}$, and the right-hand-side, $\underline{b}$, are specified in Appendix~\ref{sec:system}. The notation follows the typical assembly procedure of a finite element code (see Appendix~\ref{sec:code}), whereby the contributions of the different element matrices/vectors are included as $\left( . \right)^\mathrm{(element\, number)}_\mathrm{local\, node\, number(s)}$. Note that the definition of the local element numbers for this mesh are shown in Figure~\ref{fig:mesh}(a).
In assuming periodicity, $u$ is decomposed in a global contribution $\bar{u}$, plus a fluctuation $u^\star$:
\begin{equation}
u = \bar{u} + u^\star
\end{equation}
The global contribution is some finite affine field for which
\caption{Finite element used to exemplify the linear system: (a) the local node numbering inside each element, (b) the mesh including \textit{\textsf{element}} and (global) node numbers, and (c) the mesh including periodic repetitions complete with the \textit{\textsf{element}} numbers, the (global) node numbers of the independent nodes that remain, and all periodic repetitions.}
\label{fig:mesh}
\end{figure}
For this example, assuming periodicity implies in assuming the nodes on the right side of the mesh are the periodic repetitions of those on the left side, and likewise for the nodes on the top and bottom; whereby one side is shifted by the global affine field compared to the other. Let us ignore the latter for a moment, and just focus on the periodic fluctuations. In that case
\begin{equation}
\begin{aligned}
u^\star_{ 4} &= u^\star_{ 1} \\
u^\star_{13} &= u^\star_{ 1} \\
u^\star_{16} &= u^\star_{ 1} \\
u^\star_{ 8} &= u^\star_{ 5} \\
u^\star_{18} &= u^\star_{ 9} \\
u^\star_{14} &= u^\star_{ 2} \\
u^\star_{15} &= u^\star_{ 3}
\label{eq:tyings}
\end{aligned}
\end{equation}
We will start by dividing the system in independent, $i$, and dependent, $d$, degrees-of-freedom (DOFs). By this it is meant that the dependent DOFs are periodic repetitions, while independent DOFs are all the remaining ones. For the mesh from Fig.~\ref{fig:mesh}(b):
\begin{align}
i &= \big[ 1, 2, 3, 5, 6, 7, 9, 10, 11 \big] \\
d &= \big[ 4, 8, 12, 13, 14, 15, 16 \big]
\end{align}
We continue by renumbering the DOFs such that
\begin{align}
i &= 1 .. 9 \\
d &= 10 .. 16
\end{align}
as is shown in Fig.~\ref{fig:mesh}(c). The programmatic consequences can be found in Appendix~\ref{sec:code:per}. Once we have done this we can write Eq.~\eqref{eq:tyings} as a matrix equation
We can use this to eliminate the dependent DOFs from the system of Eq.~\ref{eq:system}. First we partition the system in independent and dependent DOFs
Starting from the original system in Appendix~\ref{sec:system}, this system is detailed in Appendix~\ref{sec:system:static-condensation}.
Once this system has been solved, the dependent DOFs can be reconstructed as follows
\begin{equation}
\delta \underline{u}^\star_{d}
=
\underline{\underline{C}}_{di}
\delta \underline{u}^\star_{i}
\end{equation}
From an implementation perspective this can be done much simpler by computing the element arrays as usual, but directly assembling to the independent DOFs only. In this case the dependent DOFs get the number of their periodic counterpart, see Appendix~\ref{sec:code:per:assembly}. The resulting system is included in Appendix~\ref{sec:system:periodic-assembly}, which is indeed exactly the same as in Appendix~\ref{sec:system:static-condensation}.
The only question which remains is to introduce $\bar{u}$ to the system. We can do this before the first iteration:
\begin{equation}
\underline{u}_{(0)} = \underline{\bar{u}}
\end{equation}
The iterative update will not change the mean, hence
\begin{equation}
\underline{u}
=
\underline{\bar{u}}
+
\delta \underline{u}^\star_{(0)}
+
\delta \underline{u}^\star_{(1)}
+
\delta \underline{u}^\star_{(2)}
+
\ldots
\end{equation}
One way to construct $\underline{\bar{u}}$ is as follows
The final thing we can to is to realize that there is an antagonist of $\vec{\bm{F}}$, namely $\vec{\bm{B}}$, which we can now compute for post-processing, but not control. Should we want to, we have to introduce extra degrees of freedom for $\vec{\bm{F}}$ (in this case the number of extra DOFs is equal to the number of dimensions). These DOFs can be used to apply mixed conditions on $\vec{\bm{F}}$ and $\vec{\bm{B}}$.
We will perform the implied volume averaging using a boundary integral. Ultimately this results in the following tying relations
where $L_x$ and $L_y$ are the dimensions of the cell in Fig.~\ref{fig:mesh}(b,c). Starting from the original system extended with the extra DOFs, in Appendix~\ref{sec:system:sec:system:static-condensation-virtual}, the final system is given in Appendix~\ref{sec:system:sec:system:static-condensation-virtual:reduced}. Compared to the system in Appendix~\ref{sec:system:static-condensation} one observes that the systems that belong the normal DOFs are identical, while the extra DOFs are clearly tied to the mean.
\appendix
\newpage
\section{Finite element mesh from Figure~\ref{fig:mesh} -- basic program structure -- non-periodic}
\label{sec:code}
\begin{python}
ndim = 2
x = [
[ 0. , 0. ],
[ 1./3. , 0. ],
[ 2./3. , 0. ],
[ 1. , 0. ],
[ 0. , 1./3. ],
[ 1./3. , 1./3. ],
[ 2./3. , 1./3. ],
[ 1. , 1./3. ],
[ 0. , 2./3. ],
[ 1./3. , 2./3. ],
[ 2./3. , 2./3. ],
[ 1. , 2./3. ],
[ 0. , 1. ],
[ 1./3. , 1. ],
[ 2./3. , 1. ],
[ 1. , 1. ],
]
u = zeros( ( coor.shape ) )
connectivity = [
[ 0, 1, 5, 4],
[ 1, 2, 6, 5],
[ 2, 3, 7, 6],
[ 4, 5, 9, 8],
[ 5, 6, 10, 9],
[ 6, 7, 11, 10],
[ 8, 9, 13, 12],
[ 9, 10, 14, 13],
[10, 11, 15, 14],
]
for ielem,elem in enumerate(connectivity):
for idim in range(ndim):
xe[:,idim] = x[elem,idim]
...
for i,idof in enumerate(elem):
for j,jdof in enumerate(elem):
A[idof,jdof] += Ael[i,j]
for i,idof in enumerate(elem):
b[idof] += bel[i]
\end{python}
\newpage
\section{Finite element mesh from Figure~\ref{fig:mesh} -- basic program structure -- periodic}