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