The concept of ``collections'' is that there are many operations that we (conceptually) accomplish on an element-by-element basis that can be aggregated together into a single operation. The
basic idea is that whenever you see a piece of code that consists of a loop in which, in the body of the loop, you are accomplishing the same operation (e.g. function evaluation) based upon
an input variable tied to the loop index -- you may be able to make the operation a ``collection''. For instance, if you were given the following snippet of code in Matlab notation:
\vspace{-5mm}
\begin{verbatim}
do i = 1:N,
y(i) = dot(a,x(i,:))
end
\end{verbatim}
\noindent where $y$ is a vector of length $N$, $x$ is an $N \times M$ matrix and $a$ is a vector of length $M$.
As mentioned in our discussions of the fundamentals of StdRegions (i.e. Section \ref{sec:stdregions-fundamentals}), one of the most
fundamental tools from calculus that we regularly employ is the idea of mapping from general domains to canonical domains. General
domains are the regions in world space over which we want to solve engineering problems, and thus want to be able to take derivatives and
compute integrals. But as we learned in calculus, it is often non-trivial to accomplish differentiation or integration over these regions. We resort
to the mapping arbitrary domains back to canonical domains over which we can define various operations. We introduced our canonical domains,
which we call standard regions, in Chapter \ref{chap:stdregions}. We refer to our world space regions as local regions, which we will present
in Chapter \ref{chap:localregions}. How these two are connected are via SpatialDomains.
As will be further discussed in Section \ref{sec:spatialdomains-datastructures}, there are two fundamental purposes served by SpatialDomains: (1) holding
basic geometric information (e.g. vertex values and curvature information) and (2) holding geometric factors information. The former information relates to
the geometric way we map standard regions to local regions. The latter information relates to how we use this map to allow us to accomplish
differentiation and integration of functions in world space via operations on standard regions with associated map (geometric) information. In this section,
we will highlight the important mathematical principles that are relevant to this section. We will first discuss the mapping itself: vertex and curvature information
and how it is used. We will then discuss how geometric factor information is computed. We break this down into two subsections following the convention of the
code. We will first discussed what is labeled in the code as {\em Regular}, and denotes mappings between elements of the same dimension
(i.e. standard region triangles to 2D triangles in world space). Although our notation will be slightly different, we will use \cite{KaSh05} as our guide; we
refer the interested reader there for a more in-depth discussion of these topics.
We will then discuss what is labeled in the code as {\em Deformed}, and denotes the mappings
between standard region elements and their world space variants in a higher embedded dimension (i.e. standard region triangles to triangles lying on a
surface embedded in a 3D space representing a manifold). Since the manifold work within {\nek} was introduced as an area of research, we will
use \cite{CantwellYKPS14} as our guide. The notation therein is slightly different than that of \cite{KaSh05} because of the necessity to use
broader coordinate system transformation principles (e.g. covariance and contravariance of vectors, etc.). We will abbreviate the detailed
mathematical derivations here, but encourage the interested reader to review \cite{CantwellYKPS14} and references therein as needed. For those
unfamiliar with covariant and contravariant spaces, we encourage the reader to review \cite{BorisenkoTS}.
\subsection{Vertex and Curvature Mapping Information}
When we load in a mesh into {\nek}, elements are often described in world space based upon their vertex positions. In traditional FEM formats, this can be as
simple as a list of d-dimensional vertex coordinates, followed by a list of element definitions: each row holding four integers (in the case of tetrahedra) denoting
the four vertices in the vertex list that comprise an element. {\nek} uses an HTML-based geometry file with a more rich definition of the basic geometric
information that just described; we encourage developers and users to review our User Guide(s) for the organization and conventions used within our geometry
files. For the purposes of this section, the important pieces of information are as follows. Let us assume that for each element, we have through our MeshGraph
data structure (described in the next section) access to the vertex positions of an element. In general, each vertex $\vec{v}_j$ is a n-tuple of dimension $d$ denoting
the dimension in which the points are specified in world space. For example, when considering a quadrilateral on the 2D plane, our vertex points each contain
two coordinates denotes the x-- and y--coordinates. As shown in Figure \ref{spatialdomains:planarmap}, we can express the mapping between a world space
quadrilateral and the standard region quadrilateral on $[-1,1]\times[-1,1]$ via a bi-linear mapping function $\chi(\vec{\xi})$.
\begin{figure}[htb]
\centering
\includegraphics[width=6in]{img/planarmap.pdf}
\caption{Reference space to world space mapping of a 2D quadrilateral to a straight-sided (2D) quadrilateral in the plane via a bi-linear (i.e., $Q(1)$) mapping.}
\label{spatialdomains:planarmap}
\end{figure}
In the case of segments, this mapping function $\chi(\vec{\xi})$ is a function of one variable and is merely an affine mapping. In two dimensions, the
mapping of a straight-sided triangle is a linear mapping (i.e., P(1) in the language of traditional finite elements -- a total degree $k=1$ space) and
the mapping of a straight-sided quadrilateral is a bi-linear mapping (i.e, Q(1) in the language of traditional finite elements -- a degree $k=1$ map along each
coordinate direction combined through tensor-product). In three dimensions, the mapping of a planar-sided tetrahedron is also a linear mapping, the
mapping of a planar-sided hexahedron is a tri-linear mapping, and the prism and pyramid are mathematically somewhere in-between these two canonical
types as given in \cite{KaSh05}. The key point is that in the case of straight-sided (or planar-sided) elements, the mapping between reference space and
world space can be deduced solely based upon the vertex positions. Furthermore in these cases,
as denoted in Figure \ref{spatialdomains:planarmap}, the form of the mapping function is solely determined by type (shape) of the element. If
only planar-sided elements are used, pre-computation involving the mapping functions can be done so that when vertex value information is available,
all the data structures can be finalized.
As presented in \ref{KaSh05}, there are many components of {\nek} that capitalize on the geometric nature of the basis functions we use. We often
speak in terms of vertex modes, edge modes, face modes and internal modes -- i.e., the coefficients that provide the weighting of vertex basis functions,
edge basis functions, etc. It is beyond the scope of this developer guide to go into all the mathematical details of their definitions, etc. However, we do
want to point out a few common developer-level features that are important. In the case of straight-sided (planar-sided) elements, the aforementioned
mapping functions can be fully described by vertex basis functions. The real benefit of this approach (of connecting the mapping representation with a
geometric basis) is seen when moving to curved elements.
Consider Figure \ref{spatialdomains:curvemap} in which we modify the example given above to accommodate on curved edge. From the
mathematical perspective, we know that the inclusion of this (quadratic) edge will require our mapping function to now be in $Q(2)$. If we were
not to use the fact that our basis is geometric in nature, we would be forced to form a Vandermonde system for a set of coefficients used
to combine the tensor-product quadratic functions (nine basis functions in all), and use the five pieces of information available to us (the four
vertex values and the one point $\vec{c}_0$ that informs the curve on edge $1$. As shown in Figure \ref{spatialdomains:curvemap}, we would
expect that this updated (to accommodate a curved edge) mapping function to consist of the bi-linear mapping
function with an additional term $C_1(\xi_1,\xi_2)$ that encompasses the new curvature information.
\begin{figure}[htb]
\centering
\includegraphics[width=6in]{img/curvemap.pdf}
\caption{Reference space to world space mapping of a 2D quadrilateral to a curve (2D) quadrilateral in the plane via a bi-quadratic (i.e., $Q(2)$) mapping.}
\label{spatialdomains:curvemap}
\end{figure}
The basis we use, following \cite{KaSh05}, allows us to precisely specify $C_1(\xi_1,\xi_2)$ using the edge basis function associated with edge $1$, and
to use the point value $\vec{c}_0$ to specify the coefficient to be used. In the figure, we assume that the form of the function is collocating, but in practice
it need not be so.
In practice, edge (and face) information can be given either as a set of point positions in world space that correspond to a particular point distribution
in the reference element (i.e., evenly-spaced points or GLL points) or modal information corresponding to the geometric basis we use internally. Our
geometric file formats assume the former -- that curve information is provided to us as physical values at specified positions from which we infer (calculate)
the modal values and store these values within SpatialDomain data structures.
Assuming we now have, for each element, a way of specifying the mapping function $\chi(\vec{\xi})$, we can now move to how we compute the geometric
factors for regular as well as for deformed mappings.
\subsection{Regular Mappings: Geometric Factors}
\subsection{Regular}
$$
\vec{x}=(x_j)^T,\, j=1,\ldots,D
$$
$$
\vec{\xi}=(\xi_j)^T,\, j=1,\ldots,D
$$
$$
\vec{x}=\chi(\vec{xi})
$$
$$
{\bf J}=\frac{\partial x_i}{\partial\xi_j}
...
...
@@ -34,10 +115,12 @@ gmat and jac
note that gmat is a tensor defined at each point, and jac is the determinant of the jacobian.