Functions and kernels#

Caution

If you are not familiar with C++ syntax, we recommend you to see C++ syntax before reading this page.

The Function class#

In order to deal with functions with parameters of any kind it is necessary to use an object function (Function class) which is related to a Parameters object (a list of parameters, see The Parameter and Parameters classes). This approach allows passing friendly, at low level of the code, some user’s functions, say functions defined in the main program.

User function and object function#

When you want to deal with the integral term:

\[\int_{\Omega}e^{ikx}\;u(x,y,z)\;v(x,y,z)\,d\Omega,\]

where the loop of finite element computation requires the computation of the function \(f(x)=e^{ikx}\) on quadrature points, it is necessary to pass the function \(f\) to the low-level code where the finite element loop is implemented. Most functions are function of a point or a list of points. So, only four kinds of function are concerned:

  • Function of an n-dimensional point (see The Point class for the class Point);

  • function of two n-dimensional points, usually named kernel;

  • function of a vector of n-dimensional point;

  • function of two vectors of n-dimensional points.

The functions may also have a Parameters as input argument if necessary (default value); for example, the real \(k\) in the next example. The output argument may be of any type (real, complex, vector, matrix, …), but has to be compatible which the type required in computation. The way to define such a function is the following. First, define a standard C++ function, for instance:

Complex f(const Point & P, Parameters& pa = defaultParameters)
{
  Real k=pa(1);      // k is the first parameter of the parameter's list pa
  // Real k=pa("k"); is available if you have a parameter with name "k"
  Real x=P(1);       // x is the first coordinate of the point P
  return exp(i*k*x); // return a complex value result,
  //return exp(i*pa(1)*P(1)); is also possible
}

Then create in your main program a Parameters object, a Function object from C++ function and you Parameters object:

Parameters pars(1, "k");  // function parameters
Function F(f, pars);      // Function object

If you have to deal with the integral term involving a real value matrix (with no parameter involved in this example):

\[\int_{\Omega}A\nabla u\,\cdot\nabla v\,d\Omega,\]

you may define:

Matrix<Real> A(const Point & P, Parameters& pa = defaultParameters)
{
  Matrix<Real> vA(3, 3);
  ...
  return vA;
}

Warning

Note that even if your function does not involve some parameters, the second argument of type Parameters is mandatory in the function definition.

For a kernel type function, it is quite similar. You have to specify two points as input arguments:

Complex Green_Helmholtz_3D(const Point& M, const Point& P, Parameters& pa = defaultParameters)
{
  Real r=distance(M, P); // we assume a distance function exists
  Real k=pa(1);
  Real eps=pa(2);
  if (r>eps) return exp(i*k*r)/r;
  else ...
}

The vector form of a Function is a function working with a vector of points and returning a vector of results (values on each point). It may be useful when computing n values together is really faster than computing \(n\) times a single value. Such a function should be declared, for instance, as follows:

Vector<Complex> vf(const Vector<Point> & vP, Parameters& pa = defaultParameters)
{
  Real k=pa(1);    //k is the first parameter of the parameter's list pa
  uint n=vP.size();
  Vector<Complex> res(n);
  for(int j=1;j<=n;j++) res(j)=exp(i*k*vP(j)(1));
  return res;
}

Note that the function has to return a Vector object. For a function involving a couple of vectors of points, the syntax of the declaration of the function is:

Vector<Complex> vf(const Vector<Point> & vP, const Vector<Point> & vQ, Parameters& pa = defaultParameters)

The C++ functions defined by the user may be used directly as argument of some functions of the library. But in most of cases it is necessary to define explicitly the object Function associated to the user function. This is the way to do this:

Parameters par;
par<<2.<<.000001;              // k and eps values, inserted in a parameter list
Function funcf(f, par);        // define a scalar function object using parameters
Function funcA(A);             // define a matrix function object
Function funcG(Green_Hemholtz_3D, par); // define a scalar kernel
Function funcvf(vf, par);      // define a scalar function in its vector form

Warning

Do not confuse the vector form of a Function and a function which returns a vector! The vector form means a function which computes a quantity (scalar, vector, matrix) on a set of a points or bipoints, the result being a vector of scalars, vectors or matrices. For most applications, scalar form of Function are sufficient. Vector form is an extension allowing the user to compute the function more efficiently in the case of multiple evaluations.

When the user wants to associate some parameters to his function, it is mandatory to define the object Function because it stores the list of parameters. To understand the role of the object Function, note that if P is a Point, the two following instructions are allowed and give the same result:

r=f(P, par); // call directly the user function f
funcf(P, r); // call the user function f using the object function funcf

In other words, the object Function shadows the Parameters object. In this example, using an object function seems to be artificial. The object function has a real interest when internal computational routines requires user’s functions, because it is easier to send one type of object encapsulating various type of function rather than many objects.

When you have to pass an object Function to a function which requires such an object, it is possible to use the constructor syntax Function(f,param), avoiding the explicit creation of the object Function:

Parameters par;
par<<2.<<.000001;          // value of k and eps
compute(Function(f, par));  // compute requiring an object function

Tip

When using a function returning a vector or a matrix, because XLiFE++ checks the dimension(s) of the returned object by using a fake point or a normal vector, it may occur some consistency problems with the dimension of points or normals that is set to 3 by default.

You can change it by specifying explicitly the dimension of the function when building it:

Function f(dnuinc, pars);
f.dimPoint=2;              // change the point dimension
// or
Function f(dnuinc, 2, pars); // specifying the point dimension when building

It works in the same way for kernels.

Dealing with normal vectors#

Sometimes, user functions have to deal with some normal vectors and obviously the normal vectors will depend on the point where the function is evaluated. For instance, a function computing the normal derivative of a given incident field (e.g. a plane wave \(e^{ikx}\)), will look like:

Complex dnuinc(const Point& P,  Parameters& pa = defaultParameters)
{
  Real x=P(1), k=pa("k");    // get k from parameters
  Reals n=getN();           // get the normal vector at P
  return i_*k*exp(i_*k*x)*n(1);
}

When this function is passed to FE computation routines, the normal vector will be evaluated and transmitted to the function only if the Function object associated to the function, has declared to use the normal vector. It is done by the following instructions in the main program:

Parameters pars(1, "k");      // declare k in the parameters
Function f(dnuinc,pars);     // associate parameters to Function object
f.require(_n);               // tell the function will use the normal
TermVector B(intg(Sigma,f*v));  // use function f

The normal vector refers to the domain on which the linear or bilinear form (involving the function) acts. The orientation of normal vector is described in section Dealing with normals.

If you want to call a function involving a normal vector out of the context of FE computations, you have to “transmit” explicitly the normal vectors required by your function. The following code shows how to do this:

Real f(const Point& x, Parameters& pa = defaultParameters)
{
  Vector<Real>& n=getN(); // get the normal vector
  return n(1);
}
...
Function F(f); F.require(_n);       //say the function is using normal vector
Point P(1.,0.);                     //point belonging to gamma
Vector<Real> n1=gamma.getNormal(P); //get the normal vector to gamma at P
setNx(n1);                          //pass the normal vector to the context
Real r = F(P,r);                    //compute F at P using n1

Note that the normal vector returned by the getN() function, is not interpolated. As a consequence, normal vector on an edge of the mesh is the normal vector of one of the elements supporting the edge!

The setNx() is thread dependent, so it can be used in a multithreading context.

Hint

Other data are available in user function, depending on the computation context :

All this data are returned as references.

Dealing with spectral family#

A spectral family is defined either from a XLiFE++ function or from a collection of TermVector. When it is defined from a function, the function has to propose the computation for each basis index, say \(n\). This basis index is managed by the XLiFE++ computation algorithms and transmitted to a global data structure (thread safe). To retrieve it as a Number, users have to call the global function getBasisIndex():

Real cosny(const Point& P, Parameters& pa = defaultParameters)
{
  Real y=P(2);
  Real h=pa("h");                    // get the parameter h (user definition)
  Number n=getBasisIndex()-1;        // get the index using the global function
  if (n==0) return std::sqrt(1./h);
  else      return std::sqrt(2./h)*std::cos(n*pi_*(y+h*0.5)/h);
}

Warning

Note that basis index, set by internal algorithms, starts from 1!

Caution

Manipulating basis indexes does not work currently when two functions are involved at the same time with different indexes.

The basis index is automatically set by internal algorithms in their computation loops. But, in user context, it may be set by calling the global function setBasisIndex():

Function cos(cosny, params);
TermVectors cos_int(N);
for s(Number n=1; n<=N; n++)
{
  setBasisIndex(n); // update the basis index
  cos_int(n)=TermVector(u, sigmaP, cos, "c"+tostring(n-1));
}

Advanced user#

Delaying computations#

It may occur that the function you plan to use is a very complex one, involving some heavy computations that you want to compute only once by storing some intermediate results somewhere. In order to allow flexibility to the user, it is advised to use the capabilities of Parameters object to store void pointers. For instance, the first time the function is called, you can compute some reusable quantities and store them in any structure with dynamic memory allocation and store the pointer of this structure in the Parameters object. The next time the function is called, as you have an access to this void pointer (do not forget to recast it), you can recover your data.

Calling a Function object#

If you have to compute the values of the function via the object Function, there are mainly two ways to do it:

  • using an alias to the pointer function (requires that you know the type of function and arguments)

    Point P=Point(0, 0), Q=Point(1, 1);
Complex c=funcf.funSC(P);      //a function returning a complex scalar (funSC)
Complex g=funcG.kerSC(P, Q);   //a kernel returning a complex scalar (kerSC)
Matrix<Real> m=funcA.funMR(P); //a function returning a real matrix scalar(funMR)

    As this method uses a recasting of a void pointer with no checking, it can cause segmentation errors if there is a misfit between the type of function required and the real function stored in the void pointer! It is possible to check the type of arguments by using the utility functions typeReturned, structReturned, typeFunction and typeArg. This direct method is offered to developers in order to have the best performance.

  • using the safe overloaded operator (), allowing to deal with point or vector of points

    Point P=Point(0, 0), Q=Point(1, 1);
Complex c;            // complex to store the result
funcf.checkTypeOn();  // activate checking mode
funcf(P, c);          // compute a complex scalar function at point P
                      // checking mode is disabled after the computation
Vector<Point> pts;    // a vector of points
pts(1)=P;pts(2)=Q;
vector<Complex> vc;   // vector to store the result
funcf.checkTypeOn();  // activate the checking mode
funcf(Pts, vc);       // compute function at a vector of points

    This method does not require the knowledge of the exact type of the function (the output argument must be compatible !). It allows scalar or vector form independently of the form of the user function. Note that, contrary to the first method, this method uses a reference to return the values, so you have to manage its memory allocation. When the function is called with a vector of points as input, the vector result is resized if it is too small. Using the checkTypeOn function, it is possible to activate the checking of the type of argument. After computation, the checkType variable is reset to false in order to avoid unnecessary rechecking. As the checking process invokes RTTI functions (expensive in time), activate wisely this option. So, if you have to evaluate many times the function, activate the checking only for the first evaluation. Note that when the checking process is deactivated, this method is still slightly more expensive than the first one.

The Kernel class#

The Function class allows defining kernel type function, say function of two points x and y. But, to deal with integral equation, more information are required. It is the purpose of the Kernel class. A Kernel object manages mainly:

Function kernel, gradx, grady, gradxy, ndotgradx, ndotgrady, curlx, curly,
         curlxy, divx, divy, divxy, dx1, dx2, dx3;
Kernel* singPart;             // singular part of kernel
Kernel* regPart;              // regular part of kernel

Dimen dimPoint;               // dimension of points
SingularityType singularType; // singularity (_notsingular, _r, _logr, _loglogr)
Real singularOrder;           // order of singularity
Complex singularCoefficient;  // coefficient of singularity
SymType symmetry;             // kernel symmetry: _noSymmetry, _symmetric, ...
String name;                  // kernel name
Parameters userData;          // to store some additional information

Setting the kernel data is mandatory, but derivatives of kernel are often required. These can be provided by directly defining the data gradx, grady, … (Function) or the dx1, dx2, dx1 to obtain all first derivative operators. Integral calculations involving singular kernels can use quadrature methods handling singular and regular parts of the kernel, which are provided by singPart, regPart (Function pointers).

Note

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:

Real ker(const Point& x, const Point& x, Parameters& pa=defaultParameters)
{...}
...
Kernel K(ker); K.dimPoint=2;
//or
Kernel K(ker, 2);

Some classical Green’s kernel as offered by XLiFE++ as Kernel objects, see Available Green’s kernels. Besides, users can define their own Kernel objects, see User kernels.

Tensor kernel#

Tensor kernels are kernels of the following 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, the 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 XLiFE++, 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 SpectralBasis object. To define a TensorKernel, you have to specify either one or two Function defining the spectral families, one or two SpectralBasis, one or two Unknown defined on a spectral Space or one or two TermVectors handling numerical spectral family. The matrix \(\mathbb{A}\) must be also given as a Matrix or a Vector if \(\mathbb{A}\) is diagonal.

Real cosny(const Point& P, Parameters& pa = defaultParameters) // spectral functions
{
    Real y=P(2);
    Real h=pa("h");                //get the parameter h (user definition)
    Int n=basis_index();           //get the index of function to compute
    if(n==0) return sqrt(1./h);
    else     return sqrt(2./h)*std::cos(n*pi_*(y+h*0.5)/h);
}
...
Real k=2.; Number N=10;
Parameters params(1.,"h");
Vector<Complex> A(N);
for(Number n=0; n<N; n++)
    A[n]=sqrt(Complex(k*k-n*n*pi_*pi_/(h*h)));   // matrix A
Function cosN(cosny, params);                    // spectral Function's
TensorKernel tkp(cosN, A);                       // define TensorKernel

In an alternate way, you can also construct the TensorKernel as follows:

Space Sp(_domain=sigma, _basis=cosN, _dim=N, _name="vect{cos(n*pi*y)}");
Unknown phi(Sp, "phi");     // spectral unknown
TensorKernel tkp(phi, A);

Doing this, a spectral unknown is available to deal with problem involving such unknown. All the following constructors are available:

template <typename T>
TensorKernel(SpectralBasis&, vector<T>&, bool conj = false);
TensorKernel(SpectralBasis&, Matrix<T>&, bool conj = false);
TensorKernel(SpectralBasis&, vector<T>&, SpectralBasis&);
TensorKernel(SpectralBasis&, Matrix<T>&, SpectralBasis&);
TensorKernel(Function, vector<T>&, bool conj = false);
TensorKernel(Function, Matrix<T>&, bool conj = false);
TensorKernel(Function, vector<T>&, Function);
TensorKernel(Function, Matrix<T>&, Function);
TensorKernel(Unknown&, vector<T>&, bool conj = false);
TensorKernel(Unknown&, Matrix<T>&, bool conj = false);
TensorKernel(Unknown&, vector<T>&, Unknown&);
TensorKernel(Unknown&, Matrix<T>&, Unknown&);
TensorKernel(vector<TermVector>&, vector<T>&, bool conj = false);
TensorKernel(vector<TermVector>&, Matrix<T>&, bool conj = false);
TensorKernel(vector<TermVector>&, vector<T>&, vector<TermVector>&);
TensorKernel(vector<TermVector>&, Matrix<T>&, vector<TermVector>&);

By default, the second family (may be same as the first one) is not conjugated in computation. By activating the option conj to true, the second family will be conjugated in computation.

It is possible to associate some geometric maps (\(f,g\)) to a tensor kernel:

\[K(x,y)=\sum_{1\le m\le M}\sum_{1\le n\le N} \psi_m(f(x))\mathbb{A}_{mn}\phi_n(g(y)).\]
Vector<Real> f(const Point& P, Parameters& pa = defaultParameters){...}
Vector<Real> g(const Point& P, Parameters& pa = defaultParameters){...}
...
TensorKernel tkp(cosN, A);
tkp.xmap=Function(f);
tkp.ymap=Function(g);

When TensorKernel is defined from vector spectral families (i.e. \(\boldsymbol{\phi}_n\) is a vector function), the vector tensor kernel reads

\[\mathbb{K}(x,y)=\sum_{1\le m\le M}\sum_{1\le n\le N} \boldsymbol{\phi}_n(y)\mathbb{A}_{mn}\boldsymbol{\psi}_m^t(x)\]

with \(\mathbb{K}\) a matrix dim \((\boldsymbol{\phi}_n)\times\) dim \((\boldsymbol{\psi}_m)\). This vector kernel can be used in bilinear form as follows:

\[\int_{\Gamma}\int_{\Sigma}\mathbf{u}(y)|(K(x,y)*\mathbf{v}(x))=\sum_{1\le m\le M}\sum_{1\le n\le N}\mathbb{A}_{mn}\int_{\Gamma}(\mathbf{u}|\boldsymbol{\phi}_n) \int_{\Sigma}(\mathbf{v}|\boldsymbol{\psi}_m)\]

with dim( \(\boldsymbol{u})=\) dim \((\boldsymbol{\phi}_n)\) and dim( \(\boldsymbol{v})=\) dim \((\boldsymbol{\psi}_m)\) and \((.|.)\) the inner product. In terms of FE (\(\boldsymbol{w}_j\) and \(\boldsymbol{\tau}_i\) basis functions related respectively to \(u\) and \(v\)), XLiFE++ produces the matrix

\[\mathbb{T}_{ij} = \sum_{1\le m\le M}\sum_{1\le n\le N}\mathbb{A}_{mn}\int_{\Gamma}(\boldsymbol{w}_j|\boldsymbol{\phi}_n) \int_{\Sigma}(\boldsymbol{\tau}_i|\boldsymbol{\psi}_m)\]

Dealing with normal vectors#

For 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 Function:

Complex G(const Point& P, const Point& Q, Parameters& pa = defaultParameters)
{
    Reals nx= getNx(); // get the normal vector at P
    Reals ny= getNy(); // get the normal vector at Q
    ...
}
...
Parameters pars(1, "k"); // declare k in the parameters
Kernel K(g, pars);       // associate parameters to Kernel object
K.require(_nx);         // declare that kernel uses x-normal
K.require(_ny);         // declare that kernel uses y-normal
TermMatrix B(intg(Sigma, Sigma, u*K*v));  // use kernel K

Kernel functions defined in XLiFE++ managed the normal vectors, so the user has not to deal with.

Tabulated functions: the Tabular class#

The Tabular class allows managing tabulated values of a function. More precisely, it handles the values of a function on a uniform grid of any dimension, mainly providing some tools to create the table, to save the table to a file, to load a table from a file, and to evaluate the values of the function at some points by interpolating from the values on the grid.

Up to now, only scalar (T is Real or Complex) functions of 1,2 or 3 real coordinates, may be with a Parameters, are handled:

T f1 (Real);
T f1p (Real, Parameters&);
T f2 (Real, Real);
T f2p (Real, Real, Parameters&);
T f3 (Real, Real, Real);
T f3p (Real, Real, Real, Parameters&);

To create a Tabular define your function and give the grid parameters (the starting coordinate x0, the step dx, the number of steps nbx, an optional name) along each coordinate:

Real f1(Real x){return sin(x);}
...
Tabular<Real> t1(0., 0.01, 100, f1);

Now, you can print the Tabular, save the Tabular to a file and evaluate the tabulated function at any point inside the grid domain:

theCout<<t1;
t1.saveToFile("f1.tab", "my optional comment");
Real y=t1(0.5);

If you saved a Tabular to a file, you can load it in a new Tabular of the same type:

Tabular<Real> t1("f1.tab");

To deal with a 2D or 3D grid, it is quite similar:

Real f2(Real x, Real y){return sin(x*y);}
...
Tabular<Real> t2(0., 0.01, 100, 0., 0.005, 50, f2);
theCout<<t2;
t2.saveToFile("f2.tab", "my optional comment");
Real y=t2(0.5, 0.5);

The general form of Tabular constructor and evaluation are ([ ] meaning an optional argument):

Tabular<T> tab(x0, dx, nx, [y0, dy, ny], [z0, dz, nz],
               fun, ["name_x"], ["name_y"], ["name_z"]);
T r = tab(x, [y], [z]);

Tabular class inherits from std::vector<Real>. On a grid \(\small x_i=x_0+i*dx,\ [y_j=y_0+j*dy,\ [z_k=z_0+k*dz]]\), \(\small i=0,n_x,\ [j=0,n_y,\ [k=0,n_z]]\), tabulated values \(\small f(x_i, [y_j, [z_k]])\) are stored in vector at place \(\small [[k*nx*ny]+j*nx]+i\).

Tip

If you need a grid of dimension greater than 3, contact the developers. But be care, a 4D grid with 100 steps in each direction requires to store \(10^8\) values, almost 1 Go of memory for double values.

Tabulated Function#

In order to use tabulated function in FE computation, it is possible to associate a Tabular to a Function:

Real f(Real x, Real y){return sin(x*y);}
...
Tabular<Real> tab(0., 0.01, 100, 0., 0.005, 50, f);  // create tabular
Function F(tab, 2); //associate tabular to f

To be consistent, the dimension of the input Point of the function has to be specified to the dimension of the Tabular ! Up to now, Tabular dimension is limited to 3.

It may be of interest to use a Tabular with a dimension less than the Point dimension; for instance a Tabular defined only for the radial coordinate. In that case, a Function relating Point to radial coordinates must be defined:

Real f(Real x, Real y){return sin(x*y);}
Real radial(const Point& X, Parameters& pars=defaultParameters) {return norm(X);}
...
Tabular<Real> tab(0., 0.01, 100, 0., 0.005, 50, f); // create tabular from f
Function R(radial);                // create R function from radial function
Function F(tab, 2, R);             // associate R map

Note that it is also possible to associate a Tabular to the function itself:

Real f(const Point& X, Parameters& pars=defaultParameters){return sin(x*y);}
...
Tabular<Real> tab(0., 0.01, 100, 0., 0.005, 50, f); // create tabular from f
Function F(f, 2);                                      // create Function
F.createTabular(0., 0.01, 100, 0., 0.005, 50);         // F is now tabulated

Once a Tabular is handled by a Function, it will be used in any circumstances!

Warning

Tabulated Function works only for scalar functions !

The SymbolicFunction class#

Some finite element constructions require passing as argument a SymbolicFunction object, that is a symbolic expression of a function. Such object is built in by writting any expression involving

  • some constants (real or complex)

  • some variables : x_1, x_2, x_3

  • some algebraic operators : +, -, *, /, ^ (power)

  • some boolean operators : &&, ||, <, <=, >, >=, ==, !=, ! (not)

  • some mathematical functions : abs, realPart, imagPart, sqrt, squared, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, exp, log, log10

The result of a boolean expression is either 0 (false) or 1 (true).

Tip

the functions asinh, acosh, atanh are only available when compiling with C++ 2011 standard.

To deal with a symbolic function is very easy. For instance consider the function

\[f(x_1, x_2)=[1-x_1+cos(x_2)]^3\]

To define it as a symbolic function, write:

SymbolicFunction fs = (1-x_1+cos(x_2))^3;

To evaluate it, write:

Real r = fs(-1, 0);
Complex c = fs(i_, 0);

You can mix algebraic and boolean expressions. Consider the following function:

\[\begin{split}f(x_1)=\left\{ \begin{array}{cc} e^{x_1} & x_1 \le 0 \\ 1 & x_1 > 0 \end{array} \right.\end{split}\]

To define it, write:

SymbolicFunction fs = (x_1<=0)*exp(x_1)+(x_1>0);

Symbolic function can be derived:

SymbolicFunction fs = power(sin(x_1),cos(x_2));
SymbolicFunction dx1fs = derivative(fs,_x1);
SymbolicFunction dx2fs = derivative(fs,_x2);
/* gives :
     ((cos(x2)*cos(x1))/sin(x1))*pow(sin(x1),cos(x2))
    -((sin(x2)*log(sin(x1)))*pow(sin(x1),cos(x2))) */