1%%%%%%%%%%%%%%%%%%%
2% XLiFE++ is an extended library of finite elements written in C++
3%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4%
5%     This program is free software: you can redistribute it and/or modify
6%     it under the terms of the GNU General Public License as published by
7%     the Free Software Foundation, either version 3 of the License, or
8%     (at your option) any later version.
9%     This program is distributed in the hope that it will be useful,
10%     but WITHOUT ANY WARRANTY; without even the implied warranty of
11%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12%     GNU General Public License for more details.
13%     You should have received a copy of the GNU General Public License
14%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15%%%%%%%%%%%%%%%%%%%
16
17\section{Kernel}
18
19The {\class Function} class allows to define kernel type function, say function of two points. But, to deal with integral equation, more informations are required. It is the role of the {\class Kernel} class. A {\class Kernel} object manages mainly:
20\begin{lstlisting}[deletekeywords={[3] kernel, symmetry, name}]
21Function kernel;              // kernel function
22Function gradx;               // x derivative
23Function grady;               // y derivative
24Function gradxy;              // x,y derivative
25Function ndotgradx;           // nx.gradx if available
26Function ndotgrady;           // ny.grady if available
27Function curlx;               // curl_x if available
28Function curly;               // curl_y if available
29Function curlxy;              // curl_x curl_y if available
30Function divx;                // div_x if available
31Function divy;                // div_y if available
32Function divxy;               // div_x div_y if available
33Function dx1;                 // d_x1 if available
34Function dx2;                 // d_x2 if available
35Function dx3;                 // d_x3 if available
36Kernel* singPart;             // singular part of kernel
37Kernel* regPart;              // regular part of kernel
38
39Dimen dimPoint;               // dimension of points
40SingularityType singularType; // singularity (_notsingular, _r, _logr,_loglogr)
41Real singularOrder;           // order of singularity
42Complex singularCoefficient;  // coefficient of singularity
43SymType symmetry;             // kernel symmetry :_noSymmetry, _symmetric ...
44String name;                  // kernel name
45Parameters userData;          // to store some additional informations
46\end{lstlisting}
47\vspace{.1cm}
48When dealing with a matrix kernel it may be useful to give {\cmd d\_x1, d\_x2, d\_x3} because it is not possible to define the gradient of a matrix kernel as a {\class Function}.\\
49
50So when defining a new one, you have to provide such informations. To understand how it works, this is  the example of Helmholtz3d kernel.\\
51First define all the functions as ordinary c++ functions :
52\begin{lstlisting}[deletekeywords={[3] x, y}]
53// kernel G(k; x, y)=exp(i*k*r)/(4*pi*r)
54Complex Helmholtz3d(const Point& x, const Point& y, Parameters& pa)
55{
56  Real k = real(pa("k"));
57  Real r = x.distance(y);
58  Complex ikr = Complex(0., 1.) * k * r;
59  return over4pi * std::exp(ikr) / r;
60}
61
62Vector<Complex> Helmholtz3dGradx(const Point& x, const Point& y, Parameters& pa)
63{
64  Real k = real(pa("k"));
65  Real r2 = x.squareDistance(y);
66  Real r = std::sqrt(r2);
67  Complex ikr = Complex(0., 1.) * k * r;
68  Complex dr = (ikr - 1.) / r2;
69  Vector<Complex> g1(3);
70  scaledVectorTpl(over4pi * exp(ikr)*dr / r, x.begin(), x.end(), y.begin(), g1.begin());
71  return g1;
72}
73
74Vector<Complex> Helmholtz3dGrady(const Point& x, const Point& y,Parameters& pa)
75{
76  Real k = real(pa("k"));
77  Real r2 = x.squareDistance(y);
78  Real r = std::sqrt(r2);
79  Complex ikr = Complex(0., 1.) * k * r;
80  Complex dr = (ikr - 1.) / r2;
81  Vector<Complex> g1(3);
82  scaledVectorTpl(-over4pi * exp(ikr)*dr / r, x.begin(), x.end(), y.begin(), g1.begin());
83  return g1;
84}
85\end{lstlisting}
86\vspace{.3cm}
87Define regular part functions :
88\begin{lstlisting}[deletekeywords={[3] x, y}]
89// regular part: G_reg(k; x, y)=(exp(i*k*r)-1)/(4*pi*r)
90Complex Helmholtz3dReg(const Point& x, const Point& y, Parameters& pa)
91{
92  Complex g;
93  Real k = real(pa("k"));
94  Real kr = k * x.distance(y);
95  Complex ikr = Complex(0., kr);
96  if (std::abs(kr) < 1.e-04)
97  {
98    int n=4; // for abs(kr)<1.e-4 this is a good choice for n (checked)
99    g = 1 + ikr / n--;
100    while (n > 1) {g = 1 + g * ikr / n--;}
101    return g *= Complex(0., over4pi * k);
102  }
103  else return over4pi * k * (std::exp(ikr) - 1.) / kr;
104}
105
106Vector<Complex> Helmholtz3dGradxReg(const Point& x, const Point& y, Parameters& pa)
107{
108  Real k = real(pa("k"));
109  Real r = x.distance(y);
110  Complex ikr = Complex(0., k*r);
111  Complex t = over4pi * (1. + std::exp(ikr)*(ikr - 1.))/ r;
112  Vector<Complex> g(3);
113  scaledVectorTpl(t/ r, x.begin(), x.end(), y.begin(), g.begin());
114  return g;
115}
116
117Vector<Complex> Helmholtz3dGradyReg(const Point& x, const Point& y, Parameters& pa)
118{
119  Real k = real(pa("k"));
120  Real r = x.distance(y);
121  Complex ikr = Complex(0., k*r);
122  Complex t = - over4pi * (1. + std::exp(ikr)*(ikr - 1.))/ r;
123  Vector<Complex> g(3);
124  scaledVectorTpl(t/ r, x.begin(), x.end(), y.begin(), g.begin());
125  return g;
126}
127\end{lstlisting}
128\vspace{.3cm}
129Define singular part functions :
130\begin{lstlisting}[deletekeywords={[3] x, y}]
131//construct Helmholtz3d Kernel singular part: G_sing(k; x, y)=1/(4*pi*r)
132Complex Helmholtz3dSing(const Point& x, const Point& y, Parameters& pa)
133{
134  Real r = x.distance(y);
135  return over4pi/r;
136}
137
138Vector<Complex> Helmholtz3dGradxSing(const Point& x, const Point& y, Parameters& pa)
139{
140  Real r = x.distance(y);
141  return -over4pi / (r*r);
142  Complex t = -over4pi / (r*r);
143  Vector<Complex> g(3);
144  scaledVectorTpl(t, x.begin(), x.end(), y.begin(), g.begin());
145  return g;
146}
147
148Vector<Complex> Helmholtz3dGradySing(const Point& x, const Point& y, Parameters& pa)
149{
150  Real r = x.distance(y);
151  Complex t = over4pi / (r*r);
152  Vector<Complex> g(3);
153  scaledVectorTpl(t, x.begin(), x.end(), y.begin(), g.begin());
154  return g;
155}
156\end{lstlisting}
157\vspace{.5cm}
158
159Now construct {\class Kernel} objects :
160\begin{lstlisting}[deletekeywords={[3] kernel, symmetry, name}]
161Parameters pars;
162pars<<Parameter(1.,"k");
163
164Kernel H3Dreg;   //regular part
165H3Dreg.name="Helmholtz 3D kernel regular part";
166H3Dreg.singularType =_notsingular;
167H3Dreg.singularOrder = 0;
168H3Dreg.singularCoefficient = over4pi;
169H3Dreg.symmetry=_symmetric;
170H3Dreg.userData = pars;
171H3Dreg.dimPoint = 3;
172H3Dreg.kernel = Function(Helmholtz3dReg, pars);
173H3Dreg.gradx = Function(Helmholtz3dGradxReg, pars);
174H3Dreg.grady = Function(Helmholtz3dGradyReg, pars);
175
176Kernel H3Dsing;   //singular part
177H3Dsing.name="Helmholtz 3D kernel, singular part";
178H3Dsing.singularType =_r;
179H3Dsing.singularOrder = -1;
180H3Dsing.singularCoefficient = over4pi;
181H3Dsing.symmetry=_symmetric;
182H3Dsing.userData = pars;
183H3Dsing.dimPoint = 3;
184H3Dsing.kernel = Function(Helmholtz3dSing, pars);
185H3Dsing.gradx  = Function(Helmholtz3dGradxSing, pars);
186H3Dsing.grady  = Function(Helmholtz3dGradySing, pars);
187
188Kernel H3D;   //kernel
189H3D.name="Helmholtz 3D kernel";
190H3D.singularType =_r;
191H3D.singularOrder = -1;
192H3D.singularCoefficient = over4pi;
193H3D.symmetry=_symmetric;
194H3D.userData = pars;
195H3D.dimPoint = 3;
196H3D.kernel = Function(Helmholtz3d, pars);
197H3D.gradx  = Function(Helmholtz3dGradx, pars);
198H3D.grady  = Function(Helmholtz3dGrady, pars);
199H3D.regPart  = &H3Dreg;
200H3D.singPart = &H3Dsing;
201\end{lstlisting}
202\vspace{.3cm}
203\begin{remark}
204If you do not define singular and regular part kernels, some computations will not be available.
205\end{remark}
206\\
207
208In fact the Helmholtz kernels are defined in {\lib mathsResources} library of \xlifepp. To load it, use the following code:
209\vspace{.1cm}
210\begin{lstlisting}
211Parameters pars;
212pars<<Parameter(1.,"k");
213Kernel H3D = Helmholtz3dKernel(pars);
214\end{lstlisting}
215\vspace{.2cm}
216For {\class Kernel} objects, one can pass to some kernel functions either \_nx vector or \_ny vector or both using the same method as one used for {\class Function} :
217\vspace{.1cm}
218\begin{lstlisting}[]{}
219Complex G(const Point& P,const Point& Q, Parameters& pa = default_Parameters)
220{
221Reals nx= getNx(); // get the normal vector at P
222Reals ny= getNy(); // get the normal vector at Q
223...
224}
225...
226Parameters pars(1,"k"); // declare k in the parameters
227Kernel K(g,pars);       // associate parameters to Kernel object
228K.require(_nx);         // declare that kernel uses x-normal
229K.require(_ny);         // declare that kernel uses y-normal
230TermMatrix B(intg(Sigma,Sigma,u*K*v));  // use kernel K
231\end{lstlisting}
232\vspace{.1cm}
233Kernel functions defined in \xlifepp managed the normal vectors, so the user has not to deal with.\\
234
235\vspace{.2cm}
236\begin{remark}
237	By default, the dimension of points of a Kernel is 3. When you define a 2D kernel, it may happen some troubles with point dimensions, in particular if the kernel function involves some operations sensitive to the point dimension. You can cure this problem by testing point dimensions in the kernel function or by specifying the point dimension when building Kernel:
238	\vspace{.1cm}
239	\begin{lstlisting}
240	Real ker(const Point& x, const Point& x, Parameters& pa=defaultParameters)
241	{...}
242	...
243	Kernel K(ker); K.dimPoint=2;
244	//or
245	Kernel K(ker,2);
246	\end{lstlisting}
247\end{remark}\\
248
249\vspace{2mm}
250\begin{remark}
251If you develop a new kernel for your own use, contact the administrators. May be they will be happy to integrate your work in \xlifepp.
252\end{remark}
253