Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in / Register
Toggle navigation
Menu
Open sidebar
Nektar
Nektar
Commits
5d601d93
Commit
5d601d93
authored
Apr 21, 2015
by
Chris Cantwell
Browse files
Merge branch 'master' into fix/linearisedwavespace
parents
26b1fe27
5b67aa99
Changes
68
Expand all
Hide whitespace changes
Inline
Sidebyside
docs/developerguide/coreconcepts/coreconcepts.tex
View file @
5d601d93
...
...
@@ 156,8 +156,9 @@ used within the code. We then initialise it in \inlsh{MyConceptImpl1.cpp}
\begin{lstlisting}
[style=C++Style]
string MyConceptImpl1::className
= GetMyConceptFactory().RegisterCreatorFunction(
"Impl1", MyConceptImpl1::create, "First implementation of my
concept.");
"Impl1",
MyConceptImpl1::create,
"First implementation of my concept.");
\end{lstlisting}
The first parameter specifies the value of the key which should be used to
select this implementation. The second parameter is a function pointer to our
...
...
@@ 239,3 +240,115 @@ a significant performance penalty (such as in tight inner loops). If needed,
Arrays allow access to the Cstyle array through the
\texttt
{
Array::data
}
member
function.
\section
{
Threading
}
\begin{notebox}
Threading is not currently included in the main code distribution. However, this
hybrid MPI/pthread functionality should be available within the next few months.
\end{notebox}
We investigated adding threaded parallelism to the already MPI parallel
Nektar++. MPI parallelism has multiple processes that exchange data using
network or networklike communications. Each process retains its own memory
space and cannot affect any other process’s memory space except through the MPI
API. A thread, on the other hand, is a separately scheduled set of instructions
that still resides within a single process’s memory space. Therefore threads
can communicate with one another simply by directly altering the process’s
memory space. The project's goal was to attempt to utilise this difference to
speed up communications in parallel code.
A design decision was made to add threading in an implementation independent
fashion. This was achieved by using the standard factory methods which
instantiate an abstract thread manager, which is then implemented by a concrete
class. For the reference implementation it was decided to use the Boost library
rather than native pthreads because Nektar++ already depends on the Boost
libraries, and Boost implements threading in terms of pthreads anyway.
It was decided that the best approach would be to use a thread pool. This
resulted in the abstract classes ThreadManager and ThreadJob. ThreadManager is
a singleton class and provides an interface for the Nektar++ programmer to
start, control, and interact with threads. ThreadJob has only one method, the
virtual method run(). Subclasses of ThreadJob must override run() and provide a
suitable constructor. Instances of these subclasses are then handed to the
ThreadManager which dispatches them to the running threads. Many thousands of
ThreadJobs may be queued up with the ThreadManager and strategies may be
selected by which the running threads take jobs from the queue. Synchronisation
methods are also provided within the ThreadManager such as wait(), which waits
for the thread queue to become empty, and hold(), which pauses a thread that
calls it until all the threads have called hold(). The API was thoroughly
documented in Nektar++’s existing Javadoc style.
Classes were then written for a concrete implementation of ThreadManager using
the Boost library. Boost has the advantage of being available on all Nektar++’s
supported platforms. It would not be difficult, however, to implement
ThreadManager using some other functionality, such as native pthreads.
Two approaches to utilising these thread classes were then investigated. The
bottomup approach identifies likely regions of the code for parallelisation,
usually loops around a simple and independent operation. The topdown approach
seeks to run as much of the code as is possible within a threaded environment.
The former approach was investigated first due to its ease of implementation.
The operation chosen was the multiplication of a very large sparse block
diagonal matrix with a vector, where the matrix is stored as its many smaller
sub matrices. The original algorithm iterated over the sub matrices multiplying
each by the vector and accumulating the result. The new parallel algorithm
sends ThreadJobs consisting of batches of sub matrices to the thread pool. The
worker threads pick up the ThreadJobs and iterate over the sub matrices in the
job accumulating the result in a thread specific result vector. This latter
detail helps to avoid the problem of cache pingpong which is where multiple
threads try to write to the same memory location, repeatedly invalidating one
another's caches.
Clearly this approach will work best when the sub matrices are large and there
are many of them . However, even for test cases that would be considered large
it became clear that the code was still spending too much time in its scalar
regions.
This led to the investigation of the topdown approach. Here the intent is to
run as much of the code as possible in multiple threads. This is a much more
complicated approach as it requires that the overall problem can be partitioned
suitably, that a mechanism be available to exchange data between the threads,
and that any code using shared resources be thread safe. As Nektar++ already
has MPI parallelism the first two requirements (data partitioning and exchange)
are already largely met. However since MPI parallelism is implemented by having
multiple independent processes that do not share memory space, global data in
the Nektar++ code, such as class static members or singleton instances, are now
vulnerable to change by all the threads running in a process.
To Nektar++’s communication class, Comm, was added a new class, ThreadedComm.
This class encapsulates a Comm object and provides extra functionality without
altering the API of Comm (this is the Decorator pattern). To the rest of the
Nektar++ library this Comm object behaves the same whether it is a purely MPI
Comm object or a hybrid threading plus MPI object. The existing data
partitioning code can be used with very little modification and the parts of the
Nektar++ library that exchange data are unchanged. When a call is made to
exchange data with other workers ThreadedComm first has the master thread on
each process (i.e. the first thread) use the encapsulated Comm object (typically
an MPI object) to exchange the necessary data between the other processes, and
then exchanges data with the local threads using direct memory to memory copies.
As an example: take the situation where there are two processes A and B,
possibly running on different computers, each with two threads 1 and 2. A
typical data exchange in Nektar++ uses the Comm method AllToAll(...) in which
each worker sends data to each of the other workers. Thread A1 will send data
from itself and thread A2 via the embedded MPI Comm to thread B1, receiving in
turn data from threads B1 and B2. Each thread will then pick up the data it
needs from the master thread on its process using direct memory to memory
copies. Compared to the situation where there are four MPI processes the number
of communications that actually pass over the network is reduced. Even MPI
implementations that are clever enough to recognise when processes are on the
same host must make a system call to transfer data between processes.
The code was then audited for situations where threads would be attempting to
modify global data. Where possible such situations were refactored so that each
thread has a copy of the global data. Where the original design of Nektar++ did
not permit this access to global data was mediated through locking and
synchronisation. This latter approach is not favoured except for global data
that is used infrequently because locking reduces concurrency.
The code has been tested and Imperial College cluster cx1 and has shown good
scaling. However it is not yet clear that the threading approach outperforms
the MPI approach; it is possible that the speedups gained through avoiding
network operations are lost due to locking and synchronisation issues. These
losses could be mitigated through more indepth refactoring of Nektar++.
\ No newline at end of file
docs/developerguide/datastructuresalgorithms/connectivity.tex
View file @
5d601d93
...
...
@@ 126,7 +126,7 @@ Both the local coordinate axis along an intersecting edge will then point in the
same direction. Obviously, these conditions will not be fulfilled by default.
But in order to do so, the direction of the local coordinate axis should be
reversed in following situations:
\begin{lstlisting}
\begin{lstlisting}
[style=C++Style]
if ((LocalEdgeId == 0)(LocalEdgeId == 1))
{
if( EdgeOrientation == Backward )
{
change orientation of local coordinate axis
...
...
@@ 184,7 +184,7 @@ mode within element e. This index i in this map array corresponds to the index o
\item
globalID represents the ID of the corresponding global degree of freedom.
\end{itemize}
However, rather than this twodimensional structure of the mapping array,
However, rather than this twodimensional structure of the mapping array,
\\
\texttt
{
LocalToGlobalMap2D::m
\_
locToContMap
}
stores the mapping array as a
onedimensional array which is the concatenation of the different elemental
mapping arrays map[e]. This mapping array can then be used to assemble the
...
...
@@ 215,7 +215,7 @@ the vertex mode, we will make a clear distinction between them in the
The fillin of the mapping array can than be summarised by the following part of
(simplified) code:
\begin{lstlisting}
\begin{lstlisting}
[style=C++Style]
for(e = 0; e < Number
_
Of
_
2D
_
Elements; e++)
{
for(i = 0; i < Number
_
Of
_
Vertices
_
Of
_
Element
_
e; i++)
{
offsetValue = ...
...
...
docs/developerguide/datastructuresalgorithms/timeintegration.tex
View file @
5d601d93
...
...
@@ 215,8 +215,8 @@ switch between some (depending on the problem) of the following
timeintegration schemes:
\begin{center}
\
small
\begin{tabular}
{
ll
}
\
footnotesize
\begin{tabular}
{
p
{
4cm
}
p
{
10cm
}
}
\toprule
AdamsBashforthOrder1
&
AdamsBashforth Forward multistep scheme of order 1
\\
AdamsBashforthOrder2
&
AdamsBashforth Forward multistep scheme of order 2
\\
...
...
docs/developerguide/developerguide.bib
deleted
100644 → 0
View file @
26b1fe27
@misc
{
nektarwebsite
,
title
=
{Nektar++: Spectral/hp element framework}
,
url
=
{http://www.nektar.info}
,
year
=
{2014}
}
@book
{
KaSh05
,
title
=
{Spectral/hp Element Methods for Computational Fluid Dynamics}
,
author
=
{Karniadakis, G. E. and Sherwin, S. J.}
,
publisher
=
{Oxford Science Publications}
,
year
=
{2005}
}
@article
{
Bu06
,
title
=
{General linear methods}
,
author
=
{Butcher, J. C.}
,
journal
=
{Acta Numerica}
,
volume
=
{15}
,
pages
=
{157256}
}
@article
{
VoEsBoChKi11
,
title
=
{A generic framework for timestepping partial differential equations (PDEs):
general linear methods, objectoriented implementation and application to fluid
problems}
,
author
=
{Vos, P. E. J. and Eskilsson, C. and Bolis, A. and Chun, S. and Kirby, R. M.
and Sherwin, S. J.}
,
journal
=
{International Journal of Computational Fluid Dynamics}
,
volume
=
{25}
,
issue
=
{3}
,
pages
=
{107125}
,
year
=
{2011}
}
docs/developerguide/developerguide.tex
View file @
5d601d93
...
...
@@ 31,6 +31,7 @@ openany, % A chapter may start on either a recto or verso page.
\usepackage
{
makeidx
}
\usepackage
{
import
}
%%% PAGE LAYOUT
%%%
\setlrmarginsandblock
{
0.15
\paperwidth
}{
*
}{
1
}
% Left and right margin
...
...
@@ 207,7 +208,9 @@ openany, % A chapter may start on either a recto or verso page.
columns=fullflexible,
backgroundcolor=
\color
{
black!05
}
,
linewidth=0.9
\linewidth
,
xleftmargin=0.1
\linewidth
xleftmargin=0.1
\linewidth
,
showspaces=false,
showstringspaces=false
}
\usepackage
{
tikz
}
...
...
@@ 371,7 +374,42 @@ Scientific Computing and Imaging Institute, University of Utah, USA}
\clearpage
\chapter
{
Introduction
}
Welcome to the developer's guide for Nektar++
\cite
{
nektarwebsite
}
.
Nektar++
\cite
{
CaMoCoBoRo15
}
is a tensor product based finite element package
designed to allow one to construct efficient classical low polynomial order
$
h
$
type solvers (where
$
h
$
is the size of the finite element) as well as higher
$
p
$
order piecewise polynomial order solvers. The framework currently has the
following capabilities:
\begin{itemize}
\item
Representation of one, two and threedimensional fields as a collection of
piecewise continuous or discontinuous polynomial domains.
\item
Segment, plane and volume domains are permissible, as well as domains
representing curves and surfaces (dimensionallyembedded domains).
\item
Hybrid shaped elements, i.e triangles and quadrilaterals or tetrahedra,
prisms and hexahedra.
\item
Both hierarchical and nodal expansion bases.
\item
Continuous or discontinuous Galerkin operators.
\item
Cross platform support for Linux, Mac OS X and Windows.
\end{itemize}
The framework comes with a number of solvers and also allows one to construct a
variety of new solvers.
Our current goals are to develop:
\begin{itemize}
\item
Automatic autotuning of optimal operator implementations based upon not
only
$
h
$
and
$
p
$
but also hardware considerations and mesh connectivity.
\item
Temporal and spatial adaption.
\item
Features enabling evaluation of highorder meshing techniques.
\end{itemize}
This document provides implementation details for the design of the libraries,
Nektar++specific data structures and algorithms and other development
information.
\begin{warningbox}
This document is still under development and may be incomplete in parts.
\end{warningbox}
\mainmatter
...
...
@@ 392,7 +430,7 @@ Welcome to the developer's guide for Nektar++\cite{nektarwebsite}.
%%% 
\bibliographystyle
{
plain
}
\bibliography
{
developerguide
}
\bibliography
{
../refs
}
\printindex
...
...
docs/developerguide/librarydesign/img/architecture.eps
0 → 100644
View file @
5d601d93
This diff is collapsed.
Click to expand it.
docs/developerguide/librarydesign/librarydesign.tex
View file @
5d601d93
...
...
@@ 3,22 +3,23 @@
A major challenge which arises when one aims to develop a software package that
implements the spectral/hp element method is to implement the mathematical
structure of the method in a digestible and coherent matter. Obviously, there
are many ways to encapsulate the fundamental concepts related to the
spectral/hp
element method, depending on e.g. the intended goal of the
developer or the
chosen programming language. We will (without going in too much
detail) give a
an overview of how we have chosen to abstract and implement
spectral/hp elements
in the
\nekpp
library. However, we want to emphasise that
this is not the
only possible choice.
are many ways to encapsulate the fundamental concepts related to the
spectral/
$
hp
$
element method, depending on e.g. the intended goal of the
developer or the
chosen programming language. We will (without going in too much
detail) give a
an overview of how we have chosen to abstract and implement
spectral/hp elements
in the
\nekpp
library. However, we want to emphasise that
this is not the
only possible choice.
Five different sublibraries, employing this characteristic pattern, are provided
in the full
\nekpp
library:
\begin{itemize}
\item
the supporting utilities sublibrary (LibUtilities library)
\item
the standard elemental region sublibrary (StdRegions library)
\item
the parametric mapping sublibrary (SpatialDomains library)
\item
the local elemental region sublibrary (LocalRegions library)
\item
the global region sublibrary (MultiRegions library)
\item
the s
upporting utilities
sublibrary (
LibUtilitie
s library)
\item
the s
olver support
sublibrary (
SolverUtil
s library)
\end{itemize}
This structure can also be related to the formulation of a global spectral/hp
...
...
@@ 31,18 +32,31 @@ element expansion, i.e.
library
}}}
\hat
{
u
}^
e
_
n
}_{
\mbox
{
\scriptsize
{
StdRegions library
}}}
\end{align*}
A more detailed overview of the
\nekpp
structure, including an overview of
the most important classes per sublibrary, is depicted in the figure below.
A more detailed overview of the
\nekpp
structure is given in
Figure~
\ref
{
f:library:structure
}
. A diagram showing the most important classes
in the core sublibraries, is depicted in the Figure~
\ref
{
f:library:overview
}
.
\begin{figure}
\centering
\includegraphics
[width=0.6\textwidth]
{
img/architecture
}
\caption
{
Structural overview of the Nektar++ libraries
}
\label
{
f:library:structure
}
\end{figure}
\begin{figure}
\centering
\includegraphics
[width=\textwidth]
{
img/overview.png
}
\caption
{
Diagram of the important classes in each library.
}
\label
{
f:library:overview
}
\end{figure}
\section
{
LibUtilities
}
This contains the underlying building blocks for constructing a Spectral Element
formulation including linear algebra, polynomial routines and memory management.
Contents
This contains the underlying building blocks for constructing a spectral element
formulation including linear algebra, polynomial routines and memory management
\cite
{
Ga39,AbSt64,CaHuYoQu88,GhOs70,KaSh05
}
. This includes:
\begin{itemize}
\setlength
{
\itemsep
}{
0em
}
\item
Basic Constants
\item
Basic Utilities ( Nektar++ Arrays)
\item
Expression Templates
...
...
@@ 61,98 +75,117 @@ These are routines for orthogonal polynomial calculus and interpolation based on
codes by Einar Ronquist and Ron Henderson.
\subsubsection
{
Abbreviations
}
\begin{itemize}
\item
z  Set of collocation/quadrature points
\item
w  Set of quadrature weights
\item
D  Derivative matrix
\item
h  Lagrange Interpolant
\item
I  Interpolation matrix
\item
g  Gauss
\item
k  Kronrod
\item
gr  GaussRadau
\item
gl  GaussLobatto
\item
j  Jacobi
\item
m  point at minus 1 in Radau rules
\item
p  point at plus 1 in Radau rules
\end{itemize}
\begin{tabular}
{
ll
}
\toprule
Character(s)
&
Description
\\
\midrule
z
&
Set of collocation/quadrature points
\\
w
&
Set of quadrature weights
\\
D
&
Derivative matrix
\\
h
&
Lagrange Interpolant
\\
I
&
Interpolation matrix
\\
g
&
Gauss
\\
k
&
Kronrod
\\
gr
&
GaussRadau
\\
gl
&
GaussLobatto
\\
j
&
Jacobi
\\
m
&
point at minus 1 in Radau rules
\\
p
&
point at plus 1 in Radau rules
\\
\bottomrule
\end{tabular}
\subsubsection
{
Main routines
}
Points and Weights
\begin{itemize}
\item
zwgj Compute GaussJacobi points and weights
\item
zwgrjm Compute GaussRadauJacobi points and weights (z=1)
\item
zwgrjp Compute GaussRadauJacobi points and weights (z= 1)
\item
zwglj Compute GaussLobattoJacobi points and weights
\item
zwgk Compute GaussKronrodJacobi points and weights
\item
zwrk Compute RadauKronrod points and weights
\item
zwlk Compute LobattoKronrod points and weights
\item
JacZeros Compute GaussJacobi points and weights
\end{itemize}
Derivative Matrices
\begin{itemize}
\item
Dgj Compute GaussJacobi derivative matrix
\item
Dgrjm Compute GaussRadauJacobi derivative matrix (z=1)
\item
Dgrjp Compute GaussRadauJacobi derivative matrix (z= 1)
\item
Dglj Compute GaussLobattoJacobi derivative matrix
\end{itemize}
Lagrange Interpolants
\begin{itemize}
\item
hgj Compute GaussJacobi Lagrange interpolants
\item
hgrjm Compute GaussRadauJacobi Lagrange interpolants (z=1)
\item
hgrjp Compute GaussRadauJacobi Lagrange interpolants (z= 1)
\item
hglj Compute GaussLobattoJacobi Lagrange interpolants
\end{itemize}
Interpolation Operators
\begin{itemize}
\item
Imgj Compute interpolation operator gj>m
\item
Imgrjm Compute interpolation operator grj>m (z=1)
\item
Imgrjp Compute interpolation operator grj>m (z= 1)
\item
Imglj Compute interpolation operator glj>m
\end{itemize}
Polynomial Evaluation
\begin{itemize}
\item
jacobfd Returns value and derivative of Jacobi poly. at point z
\item
jacobd Returns derivative of Jacobi poly. at point z (valid at z=1,1)
\end{itemize}
\begin{tabular}
{
ll
}
\toprule
\multicolumn
{
2
}{
l
}{
\textbf
{
Points and Weights
}}
\\
Routine
&
Description
\\
\midrule
zwgj
&
Compute GaussJacobi points and weights
\\
zwgrjm
&
Compute GaussRadauJacobi points and weights (z=1)
\\
zwgrjp
&
Compute GaussRadauJacobi points and weights (z= 1)
\\
zwglj
&
Compute GaussLobattoJacobi points and weights
\\
zwgk
&
Compute GaussKronrodJacobi points and weights
\\
zwrk
&
Compute RadauKronrod points and weights
\\
zwlk
&
Compute LobattoKronrod points and weights
\\
JacZeros
&
Compute GaussJacobi points and weights
\\
\bottomrule
\end{tabular}
\begin{tabular}
{
ll
}
\toprule
\multicolumn
{
2
}{
l
}{
\textbf
{
Derivative Matrices
}}
\\
Routine
&
Description
\\
\midrule
Dgj
&
Compute GaussJacobi derivative matrix
\\
Dgrjm
&
Compute GaussRadauJacobi derivative matrix (z=1)
\\
Dgrjp
&
Compute GaussRadauJacobi derivative matrix (z= 1)
\\
Dglj
&
Compute GaussLobattoJacobi derivative matrix
\\
\bottomrule
\end{tabular}
\begin{tabular}
{
ll
}
\toprule
\multicolumn
{
2
}{
l
}{
\textbf
{
Lagrange Interpolants
}}
\\
Routine
&
Description
\\
\midrule
hgj
&
Compute GaussJacobi Lagrange interpolants
\\
hgrjm
&
Compute GaussRadauJacobi Lagrange interpolants (z=1)
\\
hgrjp
&
Compute GaussRadauJacobi Lagrange interpolants (z= 1)
\\
hglj
&
Compute GaussLobattoJacobi Lagrange interpolants
\\
\bottomrule
\end{tabular}
\begin{tabular}
{
ll
}
\toprule
\multicolumn
{
2
}{
l
}{
\textbf
{
Interpolation Operators
}}
\\
Routine
&
Description
\\
\midrule
Imgj
&
Compute interpolation operator gj>m
\\
Imgrjm
&
Compute interpolation operator grj>m (z=1)
\\
Imgrjp
&
Compute interpolation operator grj>m (z= 1)
\\
Imglj
&
Compute interpolation operator glj>m
\\
\bottomrule
\end{tabular}
\begin{tabular}
{
ll
}
\toprule
\multicolumn
{
2
}{
l
}{
\textbf
{
Polynomial Evaluation
}}
\\
Routine
&
Description
\\
\midrule
jacobfd
&
Returns value and derivative of Jacobi poly. at point z
\\
jacobd
&
Returns derivative of Jacobi poly. at point z (valid at z=1,1)
\\
\bottomrule
\end{tabular}
\subsubsection
{
Local routines
}
\begin{itemize}
\item
jacobz Returns Jacobi polynomial zeros
\item
gammaf Gamma function for integer values and halves
\item
RecCoeff Calculates the recurrence coefficients for orthogonal poly
\item
TriQL QL algorithm for symmetrix tridiagonal matrix
\item
JKMatrix Generates the JacobiKronrod matrix
\end{itemize}
\subsubsection
{
Useful references
}
\begin{itemize}
\item
`[`1`]` Gabor Szego: Orthogonal Polynomials, American Mathematical
Society, Providence, Rhode Island, 1939.
\item
`[`2`]` Abramowitz
\&
Stegun: Handbook of Mathematical Functions, Dover,
New York, 1972.
\item
`[`3`]` Canuto, Hussaini, Quarteroni
\&
Zang: Spectral Methods in Fluid
Dynamics, SpringerVerlag, 1988.
\item
`[`4`]` Ghizzetti
\&
Ossicini: Quadrature Formulae, Academic Press, 1970.
\item
`[`5`]` Karniadakis
\&
Sherwin: Spectral/hp element methods for CFD, 1999
\end{itemize}
\begin{tabular}
{
ll
}
\toprule
Routine
&
Description
\\
\midrule
jacobz
&
Returns Jacobi polynomial zeros
\\
gammaf
&
Gamma function for integer values and halves
\\
RecCoeff
&
Calculates the recurrence coefficients for orthogonal poly
\\
TriQL
&
QL algorithm for symmetrix tridiagonal matrix
\\
JKMatrix
&
Generates the JacobiKronrod matrix
\\
\bottomrule
\end{tabular}
\subsubsection
{
Notes
}
\begin{itemize}
\setlength
{
\itemsep
}{
0em
}
\item
Legendre polynomial
$
\alpha
=
\beta
=
0
$
\item
Chebychev polynomial
$
\alpha
=
\beta
=

0
.
5
$
\item
All routines are double precision.
\item
All array subscripts start from zero, i.e. vector
$
0
..N

1
$
\end{itemize}
\section
{
StdRegions
}
The StdRegions library, see the figure below, bundles all classes that mimic a
spectral/hp element expansion on a standard region. Such an expansion, i.e.
The StdRegions library, a summary of which is shown in
Figure~
\ref
{
f:library:stdregions
}
, bundles all classes that mimic a
spectral/
$
hp
$
element expansion on a standard region. Such an expansion, i.e.
\begin{align*}
u(
\boldsymbol
{
\xi
}_
i) =
\sum
_{
n
\in\mathcal
{
N
}}
\phi
_
n(
\boldsymbol
{
\xi
}_
i)
\hat
{
u
}_
n,
...
...
@@ 169,30 +202,34 @@ at the quadrature points $\boldsymbol{\xi}_i$.
All standard expansions, independent of the dimensionality or shape of the
standard region, can be abstracted in a similar way. Therefore, it is possible
to define these data structures in an
\emph
{
abstract
}
base class, i.e. the
class
!StdExpansion. This base class can also contain the implementation of methods
to define these data structures in an
\emph
{
abstract
}
base class, i.e. the
class
!StdExpansion. This base class can also contain the implementation of methods
that are identical across all shapes. Derived from this base class is another
level of abstraction, i.e. the abstract classes StdExpansion1D, StdExpansion2D
and StdExpansion3D. All other shapespecific classes (such as e.g. StdSegExp or
StdQuadExp) are inherited from these abstract base classes. These
shapespecific
classes are the classes from which objects will be instantiated.
StdQuadExp) are inherited from these abstract base classes. These
shapespecific
classes are the classes from which objects will be instantiated.
They also contain the shapespecific implementation for operations such as
integration or differentiation.
\begin{center}
\includegraphics
{
img/StdRegions.png
}
\end{center}
\begin{figure}
\centering
\includegraphics
[width=\textwidth]
{
img/StdRegions.png
}
\caption
{
Main classes in the StdRegions library.
}
\label
{
f:library:stdregions
}
\end{figure}
\section
{
SpatialDomains
}
The most important family of classes in the SpatialDomains library is the
Geometry family, as can also be seen in the figure below. These classes are the
representation of a (geometric) element in
\emph
{
physical space
}
. It has been
indicated before that every local element can be considered as an image of the
standard element where the corresponding onetoone mapping can be represented
as an elemental standard spectral/hp expansion. As such, a proper encapsulation
should at least contain data structures that represent such an expansion in
order to completely define the geometry of the element. Therefore, we have
equipped the classes in the Geometry family with the following data structures:
The most important family of classes in the SpatialDomains library is the
Geometry family, as can also be seen in Figure~
\ref
{
f:library:spatialdomains
}
.
These classes are the representation of a (geometric) element in
\emph
{
physical
space
}
. It has been indicated before that every local element can be considered
as an image of the standard element where the corresponding onetoone mapping
can be represented as an elemental standard spectral/hp expansion. As such, a
proper encapsulation should at least contain data structures that represent such
an expansion in order to completely define the geometry of the element.
Therefore, we have equipped the classes in the Geometry family with the
following data structures:
\begin{itemize}
\item
an object of !StdExpansion class, and
...
...
@@ 208,9 +245,12 @@ However, for every shapespecific geometry class, it needs to be initialised
according to the corresponding StdRegions class (e.g. for the QuadGeom class,
it needs to be initialised as an StdQuadExp object).
\begin{center}
\begin{figure}
\centering
\includegraphics
{
img/SpatialDomains.png
}
\end{center}
\caption
{
Main classes in the SpatialDomains library.
}
\label
{
f:library:spatialdomains
}
\end{figure}
\section
{
LocalRegions
}
...
...
@@ 247,8 +287,8 @@ x = \chi(\xi) = \frac{1\xi}{2}x_{e1} + \frac{1+\xi}{2}x_e \quad \xi
In two dimensions, for a quadrilateral, each coordinate is given by
\begin{align*}
x
_
i =
\chi
(
\xi
_
1,
\xi
_
2) = x
_
i
^
A
\frac
{
1
\xi
_
1
}{
2
}
\frac
{
1
\xi
_
2
}{
2
}
+
x
_
i
^
B
\frac
{
1+
\xi
_
1
}{
2
}
\frac
{
1
\xi
_
2
}{
2
}
+
x
_
i =
\chi
(
\xi
_
1,
\xi
_
2)
&
= x
_
i
^
A
\frac
{
1
\xi
_
1
}{
2
}
\frac
{
1
\xi
_
2
}{
2
}
+
x
_
i
^
B
\frac
{
1+
\xi
_
1
}{
2
}
\frac
{
1
\xi
_
2
}{
2
}
\\
&
\qquad
+
x
_
i
^
D
\frac
{
1
\xi
_
1
}{
2
}
\frac
{
1+
\xi
_
2
}{
2
}
+
x
_
i
^
C
\frac
{
1+
\xi
_
1
}{
2
}
\frac
{
1+
\xi
_
2
}{
2
}
,
\quad
i=1,2
\end{align*}
...
...
@@ 263,10 +303,11 @@ shape functions, $f^A(\xi_1), f^B(\xi_2), f^C(\xi_1)$ and
$
f
^
D