inputs/dev/Kernel.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

\section{The {\classtitle Kernel} class}
The {\class Kernel} class is intended to deal with kernel, say function $G(x,y)$ where $x,y$ are two points. Basically, it manages the function itself using {\class Function} object that is able to deal with such function but it manages other functions such as $\nabla_xG$, $\nabla_yG$, $\nabla_{xy}G$, $\frac{\partial G}{\partial n_x}$, $\frac{\partial G}{\partial n_y}$:
\begin{lstlisting}[]{}
public:
Function kernel;
Function gradx;
Function grady;
Function gradxy;
Function ndotgradx;
Function ndotgrady;
Parameters userData;
\end{lstlisting}
All functions defined in {\class Kernel} class share the same parameters given by \textbf{userData}. Besides, the class provides some additional properties:
\begin{lstlisting}[]{}
SingularityType singularType;   // singularity type
real_t singularOrder;           // singularity order
complex_t singularCoefficient;  // singularity coefficient
SymType symmetry;               // kernel symmetry
string_t name;                  // kernel name
dimen_t dimPoint;               // dim of points x,y
//temporary data
mutable bool checkType_;       // to activate checking mode
mutable bool conjugate_;       // flag for conjugate operation
mutable bool transpose_;       // flag for transpose operation
mutable bool xpar;             // true if x is consider as a parameter
mutable Point xory;            // x/y as parameter of the kernel
 //to manage data requirements
bool requireNx;              // true if requires normal or x-normal vector
bool requireNy;              // true if requires y-normal vector
bool requireElt;             // true if requires element
bool requireDom;             // true if requires domain
bool requireDof;             // true if requires dof
\end{lstlisting}
The following singularity types are currently available:
\begin{lstlisting}[]{}
enum SingularityType {_notsingular, _r, _logr,_loglogr};
\end{lstlisting}
\vspace{.2cm}
It is also possible to specify the singular and regular part using {\class Kernel} pointer (0 by default):
\begin{lstlisting}[]{}
Kernel* singPart;
Kernel* regPart;
\end{lstlisting}
There are simple constructors from {\class Function} objects or from explicit C++ functions having a kernel form (see {\class Function} class):
\begin{lstlisting}[]{}
Kernel();                            // default constructor
Kernel(const Kernel&);               // copy constructor
Kernel& operator=(const Kernel&);    // assign operator
//from Function object
Kernel(const Function& ker, SingularityType st=_notsingular,
       real_t so=0, SymType sy=_noSymmetry);
Kernel(const Function& ker, const string_t& na, SingularityType st=_notsingular,
       real_t so=0, SymType sy=_noSymmetry);
//from explicit function having a kernel form
template <typename T>
Kernel(T(fun)(const Point&, const Point&, Parameters&),
       SingularityType st=_notsingular, real_t so=0, SymType sy=_noSymmetry);
template <typename T>
Kernel(T(fun)(const Point&, const Point&, Parameters&), const string_t& na,
       SingularityType st=_notsingular, real_t so=0, SymType sy=_noSymmetry);
...
virtual void initParameters();  //to share user parameters
\end{lstlisting}
Note that these simple constructors do not allow to set other functions (gradx, grady, ...).\\

This class provides some accessors:
\begin{lstlisting}[]{}
bool conjugate() const;
bool& conjugate();
void conjugate(bool v) const;
bool transpose() const;
bool& transpose();
void transpose(bool v) const;
virtual bool isSymbolic() const;
virtual ValueType valueType() const;
virtual StrucType structType() const;
virtual dimPair dims() const;
\end{lstlisting}
\vspace{.2cm}
It provides also some tools to manage parameters:
\begin{lstlisting}[]{}
void setParameters(const Parameters& par);
void setParameter(Parameter& par);
template<typename T>
 void setParameter(const string_t& name, T value);
template<typename T>
 T getParameter(const string_t& name, T value) const;
void setX(const Point& P) const;
void setY(const Point& P) const;
\end{lstlisting}
\vspace{.2cm}
and some functions to manage the kernel data required during computation:
\begin{lstlisting}[]{}
bool normalRequired() const;
bool xnormalRequired() const;
bool ynormalRequired() const;
void require(UnitaryVector uv, bool tf=true);
void require(const string_t& s, bool tf=true);
\end{lstlisting}
\vspace{.2cm}
In some computations, a {\class Kernel} object has to be viewed as a {\class Function} object, x or y being considered as a parameter. the following member functions do the job:
\begin{lstlisting}[]{}
//move to y->k(x,y) or  x->k(x,y)
const Function& operator()(const Point& x, variableName vn ) const;
const Function& operator()(variableName vn, const Point& y) const;
const Function& operator()(variableName vn) const;
\end{lstlisting}
\vspace{.1cm}
\textbf{variableName} is defined by the following enumeration:
\begin{lstlisting}[]{}
enum VariableName {_x, _x1=_x, _y, _x2=_y, _z, _x3=_z, _t};
\end{lstlisting}
\vspace{.2cm}
The class provides function to evaluate kernel:
\begin{lstlisting}[]{}
template<typename T>
 T& operator()(const Point&, const Point&, T&) const;//compute at (x,y)
template<typename T>
 T& operator()(const Point&, T&) const;  //compute at (x,P) or (P,y) with P=xory
template<typename T>
 Vector<T>& operator()(const Vector<Point>&,const Vector<Point>&,
                           Vector<T>&) const; //compute at a list of points
\end{lstlisting}
As the {\class Kernel} class is the base class of {\class TensorKernel} class which deals with particular kernel, there are some virtual functions:
\begin{lstlisting}[]{}
virtual KernelType type() const;                 //_generalKernel, _tensorKernel
virtual const TensorKernel* tensorKernel() const;//downcast to tensorkernel
virtual TensorKernel* tensorKernel();            //downcast to tensorkernel*
\end{lstlisting}
\vspace{.2cm}
To deal with some data that are available during FE computation such as normal vectors, current element, current dof or current domain, there is two things to do : tell to \xlifepp that the Kernel will use such data
\begin{lstlisting}[]{}
Kernel K(ker);     //link a C++ function(kernel) to a Kernel object
K.require(_nx);    //tell that function ker will use normal vector
\end{lstlisting}
and get data in the C++ function related to the Kernel object:
\begin{lstlisting}[]{}
Real ker(const Point& x, const Point&y, Parameters& pars=theDefaultParameters)
{...
Vector<Real>& n = getN();      //get the normal vector
...
}
\end{lstlisting}
This machinery is based on the {\class ThreadData} class that manages one instance of data by thread.

\vspace{.2cm}
As an example, we give the sketch of the program constructing the Laplace kernel in 3d:
\begin{lstlisting}[]{}
Kernel Laplace3dKernel(Parameters& pars)
{
    Kernel K;
    K.name="Laplace 3D kernel";
    K.singularType =_r;
    K.singularOrder = -1;
    K.singularCoefficient = over4pi;
    K.symmetry=_symmetric;
    K.userData.push(pars);
    K.kernel = Function(Laplace3d, K.userData);
    K.gradx = Function(Laplace3dGradx, K.userData);
    K.grady = Function(Laplace3dGrady, K.userData);
    K.ndotgradx = Function(Laplace3dNxdotGradx, K.userData);
    K.ndotgrady = Function(Laplace3dNydotGrady, K.userData);
    return K;
}

real_t Laplace3d(const Point& x, const Point& y, Parameters& pars)
{
  real_t r = x.distance(y);
  return over4pi / r;
}

//derivatives
Vector<real_t> Laplace3dGradx(const Point& x, const Point& y,Parameters& pars)
{
  real_t r2 = x.squareDistance(y);
  real_t r = std::sqrt(r2);
  Vector<real_t> g1(3);
  scaledVectorTpl(over4pi / (r*r2), x.begin(), x.end(), y.begin(), g1.begin());
  return g1;
}

Vector<real_t> Laplace3dGrady(const Point& x, const Point& y,Parameters& pars)
{
  real_t r2 = x.squareDistance(y);
  real_t r = std::sqrt(r2);
  Vector<real_t> g1(3);
  scaledVectorTpl(-over4pi / (r*r2), x.begin(), x.end(), y.begin(), g1.begin());
  return g1;
}

real_t Laplace3dNxdotGradx(const Point& x, const Point& y, Parameters& pars)
{
 Vector<real_t>& nxp = getNx(); //get normal at x point
 std::vector<real_t>::const_iterator itx=x.begin(),ity=y.begin(), itn=nxp.begin();
 real_t d1=(*itx++ - *ity++), d2=(*itx++ - *ity++), d3=(*itx - *ity);
 real_t r = d1* *itn++; r+= d2* *itn++; r+=d3* *itn;
 real_t r3= d1*d1+d2*d2+d3*d3; r3*=std::sqrt(r3);
 return -r*over4pi_/r3;
}
...
\end{lstlisting}
See in the \textit{mathRessources} directory, the kernels that are avalaible.\\

\displayInfos{library=utils, header=Kernel.hpp, implementation=Kernel.cpp, test=test\_Function.cpp,
header dep={config.h, Parameters.hpp, Point.hpp, String.hpp, Vector.hpp, Function.hpp}}

\section{The {\classtitle TensorKernel} class}

The {\class TensorKernel} class, inherited from the {\class Kernel} class, deals with kernels of the form:
$$K(x,y)=\sum_{1\le m\le M}\sum_{1\le n\le N} \psi_m(x)\mathbb{A}_{mn}\phi_n(y)$$
with $\mathbb{A}$ a $M\times N$ matrix. Such kernel is involved, for instance, in Dirichlet to Neumann method using spectral expansion. In many cases, this matrix is a diagonal one:
$$K(x,y)=\sum_{1\le n\le N} \psi_n(x)\lambda_n\phi_n(y).$$
In the meaning of \xlifepp, the families $(\psi_m)_{1\le m\le M}$ and  $(\psi_n)_{1\le n\le N}$ may be considered as some basis of spectral space, say {\class SpectralBasis} object. So, the  {\class TensorKernel} class is defined as follow:
\begin{lstlisting}[]{}
protected :
 const SpectralBasis* phi_p;
 const SpectralBasis* psi_p;
 VectorEntry* matrix_p;
public :
 bool isDiag;       // true if matrix is a diagonal one
 Function xmap;     // geometric transformation applied to x
 Function ymap;     // geometric transformation applied to y
 bool isConjugate;  // conjugate phi_p in computation (default=false)
\end{lstlisting}
\vspace{.1cm}
The \textbf{xmap} and \textbf{ymap} functions allow to map the basis functions from a reference space to an other one. For instance, suppose  the basis $\cos(n\pi x)$ on the segment $]0,1[$ is given:
\begin{lstlisting}[]{}
Real cosn(const Point& P, Parameters& pars)
{ Real x=P(1);
  Int n=pa("basis index")-1;  //get the index of function to compute
  return cos(n*pi*x);
}
\end{lstlisting}
and has to be used on the boundary $\Sigma=\{(0,y), 0<y<h\}$. The basis will be map onto using the function:
\begin{lstlisting}[]{}
Point sigmaTo(const Point& P, Parameters& pars)
{
  return Point(P(2)/h);
}
\end{lstlisting}
\vspace{.1cm}
The {\class SpectralBasis} class manages either a basis given by an explicit C++ function (analytic form) or a basis given by a set of \textbf{TermVector} objects (numerical form). So the constructors provided by the {\class TensorKernel} class manages these different cases too:
\begin{lstlisting}[]{}
TensorKernel();    //default constructor
template<typename T>
 //from one or two SpectralBasis objects
 TensorKernel(const SpectralBasis&, const vector<T>&, bool conj=false);
 TensorKernel(const SpectralBasis&, const vector<T>&, const SpectralBasis&);
 //from one or two Unknown objects belonging to spectral spaces
 TensorKernel(const Unknown&, const vector<T>&, bool conj=false);
 TensorKernel(const Unknown&, const vector<T>&, const Unknown&);
 //from one or two Function objects
 TensorKernel(Function, const vector<T>&, bool conj=false);
 TensorKernel(Function, const vector<T>&, Function);
 //from one or two set of TermVector objects
 TensorKernel(const vector<TermVector>&, const vector<T>&, bool conj =false);
 TensorKernel(const vector<TermVector>&, const vector<T>&,
              const vector<TermVector>&);
\end{lstlisting}
Note that these constructors does not allow to set maps! They have to be set after construction. When the user defines {\class TensorKernel} from an explicit function, it has to get the index from the parameters related to the function. This index is automatically set by the computation functions and its name in parameters is "basis index".\\

The {\class TensorKernel} class provides some accessors ans properties:
\begin{lstlisting}[]{}
const VectorEntry& vectorEntry() const;
const SpectralBasis* phip() const;
const SpectralBasis* psip() const;
const SpectralBasis& phi()  const;
const SpectralBasis& psi()  const;
dimen_t dimOfPhi() const;
dimen_t dimOfPsi() const;
number_t nbOfPhi() cons;
number_t nbOfPsi() cons;
bool phiIsAnalytic() const;
bool psiIsAnalytic() const;
virtual KernelType type() const;
virtual bool isSymbolic() const;
virtual ValueType valueType() const;;
virtual StrucType structType() const;;
bool sameBasis() const;
virtual const TensorKernel* tensorKernel() const;
virtual TensorKernel* tensorKernel();
\end{lstlisting}
\vspace{.1cm}
It implements the computation of kernel:
\begin{lstlisting}[]{}
template<typename T>
 T kernelProduct(number_t, const T&) const;
template<typename T>
 T kernelProduct(number_t, number_t, const T&) const;
\end{lstlisting}
\displayInfos{library=terms, header=TensorKernel.hpp, implementation=TensorKernel.cpp, test=test\_Function.cpp,
header dep={config.h, Parameters.hpp, Point.hpp, String.hpp, Vector.hpp, Function.hpp, TermVector.hpp, SpectralBasis.hpp}}