Algebraic representation#

From the previous chapter (Problem definition), a variational form have been defined:

\[\DeclareMathOperator{\spanop}{span} \mathrm{Find\;} u \in V = \spanop (w_j)_{j=1,n}, \mathrm{\;such\;that\;} \forall v \in W = \spanop (\tau_j)_{j=1,m}, a(u,v) = l(v)\]

Algebraic representation of the formulation consists in going to a representation in terms of matrices and vectors. It leads generally to the solving of a linear system or the computation of eigenvalues and eigenvectors.

Build the linear system#

Representing unknown function \(u\) as a linear combination of basis functions :

\[\DeclareMathOperator{\spanop}{span} u \in V = \spanop (w_j)_{j=1,n} \implies u = \sum_{j=1}^n u_j w_j\]

en choosing \(v=\tau_i\), the variational formulation reads:

\[\begin{split}\mathbb{A} U = L \text{ with } \ \left\{ \begin{array}{l} U=(u_j)_{j=1,n},\ L=\left(l(\tau_i)\right)_{i=1,m} \\ \mathbb{A}=\left(a(w_j, \tau_i)\right)_{i=1,m,j=1,n}. \end{array}\right.\end{split}\]

XLiFE++ provides two fundamental classes to deal with such vectors and matrices:

  • TermVector class which handles vector and space stuff (linear form, unknowns or test functions, dof numbering, …)

  • TermMatrix class which handles matrix and space stuff (bilinear form, unknowns and test functions, dof numbering, …)

Note

The TestFunction class inherits from the Unknown class. The two classes can therefore be used independently. However, it is advisable to use them as usual. In particular, linear forms should be constructed from TestFunction object leading to TermVector object related to a TestFunction. In the same vein, the TermVector solution to the problem, is linked to an Unknown object. XLiFE++ does not check strictly the consistency of unknowns/test functions.

These two classes support either single unknowns or multiple unknown representation. Multiple unknowns vector or matrix are represented by single unknown blocks:

\[\begin{split}\mathbb{A} = \begin{pmatrix} \mathbb{A}_{v_1,u_1} & \mathbb{A}_{v_1,u_2} & \cdots & \mathbb{A}_{v_1,u_N} \\ \mathbb{A}_{v_2,u_1} & \mathbb{A}_{v_2,u_2} & \cdots & \mathbb{A}_{v_2,u_N} \\ \vdots & \vdots & \cdots & \vdots \\ \mathbb{A}_{v_M,u_1} & \mathbb{A}_{v_M,u_2} & \cdots & \mathbb{A}_{v_M,u_N} \\ \end{pmatrix}, \qquad L=\begin{pmatrix} L_{v_1} \\ \vdots \\ L_{v_M} \\ \end{pmatrix}\end{split}\]

where \(u_1, \ldots u_N\) stands for unknowns and \(v_1, \ldots v_M\) stands for test functions.

Hint

Unknowns correspond to matrix columns and test functions to matrix rows!

General purpose on TermMatrix and TermVector#

Deal with real / complex coeffs#

The value type (real or complex) is automatically managed by TermMatrix and TermVector classes.

Tip

When the result of a linear combination of terms involves complex values, the resulting term will be set to real if the imaginary part is null.

Sometimes, in an advanced usage, it may be useful to convert a real term to a complex one or access to real part or imaginary part of a complex term.

For this purpose, XLiFE++ offers complex(), real(), imag() and abs() routines.

Delay computation#

When defined, TermVector and TermMatrix are automatically computed, except if the option _notCompute is set in definition of TermMatrix or TermVector.

To perform the real computation, the compute() routine has to be called explicitely.

...
BilinearForm auv=intg(omega,grad(u)|grad(v)));
TermMatrix A(auv, _notCompute, _name="A");
...
compute(A);

Size of matrix and vector terms#

The computation algorithms find the minimal representation of vectors/matrices. It means that the size of vector/matrix is equal to the number of unknown/test function dofs (dofs = degrees of freedom) involved in the computation. For instance, a mass matrix on a boundary involves only dofs supported by the boundary.

The TermVector class#

A TermVector represents either a vector of coefficients of an unknown on a subset of degrees of freedom of an interpolation space (defined on a mesh) or the vector representing a linear form on the interpolation space (Risz theorem) . Then, the definition of a TermVector requires an Unknown (related to the approximation space and then the degrees of freedom), a Domain providing which subset of degrees of freedom is considered, and the values themselves.

There are two ways to define a TermVector:

  • from a LinearForm,

  • by evaluating an analytical function on the considered degrees of freedom.

Representation of linear forms#

This is the general case when building the right-hand-side vector of the linear system associated to the variational formulation in an approximation space (FE space for instance).

So in this context, a LinearForm and a name (of TermVector) must be given to the constructor.

Caution

The name, mainly useful for printing purpose, is also used to name stored internal objects and for post-processing. As a result, each name must be unique.

...
LinearForm lf=intg(omega, f*v);
TermVector tv(lf, _name="F");

LinearForm lf2 = intg(omega1, f*v1) + intg(omega2, f*v2);
TermVector tv2(lf2, _name="F2");

Define interpolated values of an analytical function#

Interpolation \(\Pi\) of a function \(f\) on shape functions \((w_i)\) is defined as follows:

\[\Pi : f \mapsto \Pi (f) = \sum_{i=1}^N d_i(f) w_i\]

\((d_i)_{i=1,N}\) are the generalized degrees of freedom (generalized means viewed as distributions) related to shape functions \((w_i)_{i=1,N}\).

For Lagrange finite elements, this results in:

\[d_i=\delta_{M_i}\ \rightarrow d_i(f)=f(M_i)\]

where \((M_i)_{i=1,N}\) are the Lagrange nodes. In this case, interpolated values are named “nodal values”.

Important

In the following, it will be explain how to build the interpolated values of a function, as defined above. To non Lagrange dofs (not nodal), some coordinates are virtually associated to, so the “nodal interpolation” may be enforced by using the key _nodal

Define interpolated values of C++ functions#

A scalar/vector real/complex interpolated value TermVector can be defined from a C++ function given by users. Here is an example:

Real sinxcosy(const Point& p, Parameters& pa=defaultParameters)
{ Real x=p.x, y=p.y; return sin(x)*cos(y); }
...
Domain omega=...;
Unknown u(...);
TermVector tv(u, omega, sinxcosy);

See Functions and kernels for more details on how to define such functions and manipulating parameters through the Parameters object.

For Morley elements, degrees of freedom involve first order derivatives, say \((\partial_x f, \partial_y f )\). For Argyris elements, they involve second order derivatives, say \(( \partial_{xx} f, \partial_{yy} f, \partial_{xy} f)\). This is the reason why, in order to compute interpolated values of an analytical function when using these finite elements, the analytical function but also its first order derivatives and / or second order dérivatives must be given:

Real fx(const Point& P, Parameters& pa = defaultParameters)
{
  return P(1);
}
Vector<Real> gx(const Point& P, Parameters& pa = defaultParameters)
{
  Vector<Real> g(2, 0.);
  g(1)=1.;
  return g;
}
Vector<Real> g2x(const Point& P, Parameters& pa = defaultParameters)
{
  return {Vector<Real> g(3,0.);
}
...
TermVector ux(u, omega, Function(fx)); //order 0 dofs
TermVector vx(u, omega, Function(fx), Function(gx)); //order 0/1 dofs involved
TermVector wx(u, omega, Function(fx)  Function(gx), Function(hx)); // order 0/1/2 dofs involved

It is also working for vector unknowns; in that case derivatives must be given as a unique vector collecting all the derivatives for each component \(( \partial_{x} f_1, \partial_{y} f_1, \partial_{x} f_2, \partial_{y} f_2, \cdots )\).

Define nodal values of a symbolic function#

For classical functions (trigonometrical functions, hyperbolic functions, exponential functions, …), there is an easier way to define nodal values: with symbolic functions. Let’s see the previous example rewritten with a symbolic function:

Domain omega=...;
Unknown u(...);
TermVector tv(u, omega, sin(x_1)*cos(x_2));

A symbolic function, handled by the SymbolicFunction class, can use the cartesian coordinates (x_1, x_2, x_3) or the spherical coordinates (r_ theta_ phi_), and can be linear combination of elementary functions.

As for analytical functions, a TermVector can be defined from a symbolic function on TermVector (up to 4). The purpose here is to build a TermVector by a symbolic function applied directly on TermVector values. This supposes that all considered TermVector are defined on the same approximation space but with possible different unknowns.

Warning

This works only for single unknown TermVector defined on the same approximation space (but possibly on different unknowns).

Combine nodal values#

In a same way, it is possible to combine values of 2 TermVector using a C++ function taking up to 2 values (real or complex):

Real f(const Real& r1, const Real& r2)
{ ... }
...
Domain omega=...;
TermVector tv1(...), tv2(...);
TermVector tv(tv1, tv2, f, _name="F"); // tv=f(tv1, tv2)

Warning

This works only for single unknown TermVector defined on the same unknown.

Operations on TermVector#

Algebraic operations#

Algebraic operations (\(+=, -=, *=, /=, +, -, *, /\)) are available on TermVector to build a new one.

Important

In order to be more efficient, the linear combination of terms is delayed up to the assign (=) operation or a constructor operation. It means that some expression may not be evaluated and produce warning/error message related to LcTerm class.

...
TermVector W=U+V; // OK
theCout << U+V;   // NOT EVALUATED

Inner products and norms#

The general routine for inner/hermitian product is the | operator between 2 TermVector. If both are real, | operator will compute the inner product. If both are complex, | operator will compute the hermitian product. Whatever values of TermVector, | operator always returns a complex value.

Hint

The specific routine for inner product is innerProduct(), and the specific routine for hermitian product is hermitianProduct(). To compute an inner product between a real TermVector and a complex TermVector, use the routine dotRC().

To compute norms, XLiFE++ provides norm1(), norm2() and norminfty().

Restriction#

In some circumstances it may usefull to restrict a TermVector to dofs related to a smaller domain (subdomain) than its domain definition. Use the | operator:

Domain gamma=...;
TermVector U=...;
TermVector UonGamma=U|gamma;

In case of a multiple unknowns TermVector, an unknown block may be extracted as follows:

Space V(_domain=omega, _interpolation=P1, _name="V"), H(_domain=omega, _interpolation=P0, _name="V");
Unknown u(V, _name="u", _dim=2), p(W, _name="p");
LinearForm fv=intg(omega,f|u)+intg(omega,p);  // multiple unknowns form
TermVector B(fv);                             // multiple unknowns TermVector
TermVector B_u=B(u);                          // extract u part

The extraction process creates a copy of the extracted part, which is a single unknown TermVector.

Deal with collections of TermVector#

The class TermVectors allows to deal with a collection of TermVector. One of the first objective was to write transient problems the easiest way possible (seee Solving wave equation with Leap-frog scheme), but there are a plenty of situations where it is friendly to use TermVectors.

TermVectors tvs(20);
tvs(1)=...;
tvs(2)=...;
...

TermVectors can be constructed by giving its size and then, each TermVector of the collection can be accessed with the parentheses operator (index numbering starting from 1).

There are lot of constructors and operations on TermVector. See the Language reference documentation.

The TermMatrix class#

A TermMatrix represents a matrix of values where column index is related to unknown and row index to test function, each on a subset of degrees of freedom of an interpolation space (defined on a domain). Then, a TermMatrix is mainly defined from an Unknown (also providing the related approximation space and then the related degrees of freedom), a TestFunction (also providing the related approximation space and then the related degrees of freedom), one or two Domain providing which subsets of degrees of freedom are considered, and the values themselves stored in various ways (dense, skyline, row or column compressed sparse storage regarding the matrix structure).

There are two ways to define a TermMatrix, according to the type of object represented:

Representation of bilinear forms#

This is the general case when building the matrix of the linear system associated to a PDE variational formulation. So in this context, a BilinearForm and a name has to be given.

Caution

The name, mainly useful for printing purpose, is also used to name stored internal objects and for post-processing. As a result, each name must be unique.

...
BilinearForm auv=intg(omega, grad(u)|grad(v));
TermMatrix A(auv, _name="A");

BilinearForm auv2=intg(omega1, grad(u1)|grad(v1))+intg(omega2, grad(u2)|grad(v2));
TermMatrix A2(auv2, _name="A2");

Deal with essential conditions#

As previously mentioned, essential conditions are attached directly to TermMatrix. They must be specified when constructing a TermMatrix from a bilinear form:

BilinearForm auv=intg(omega, grad(u)|grad(v));
EssentialConditions ecs= (u|sigmaM = 1) & (u|sigmaP = 1);
TermMatrix A(auv, ecs, _name="A");

Hint

Essential conditions are never attached to a TermVector! When solving a system involving essential conditions, the TermVector representing the right-hand side of the system is automatically corrected to take into account essential conditions effects.

In few words, the process to take into account essential conditions is the following:

  • the essential conditions are transformed in a constaints system \(\mathbb{C}U=D\) where \(\mathbb{C}\) is a \(p\times n\) rectangular matrix (\(n\) the dimension of \(U\), \(p < n\))

  • the constraints system is then reduced to an upper triangular system using a QR algorithm : \(\mathbb{C}=[\mathbb{E}\ \mathbb{R}]\) with \(\mathbb{E}\) a \(q\times q\) invertible matrix (\(q\le p\)) and \(\mathbb{R}\) an upper triangular \(q\times n\) matrix

  • the constraints system, now reading \(U_e=\mathbb{E}^{-1}D-\mathbb{E}^{-1}\mathbb{R}U_r\), variables \(U_e\) may be substituted in terms of \(U_r\) variables to change the matrix and may be the right hand side.

Note

The QR factorization process can detect redundant or conflicting constraints; in any case only one constraint is kept and a warning message is raised. It is the responsability of the user to check if the result is correct when such situation arises.

The QR process stops when the pivot coefficient is below a tolerance factor (by default \(100\times \varepsilon_{mach}\approx 2.10^{-14}\) in double precision). If the tolerance factor is too low, due to rounding errors, some quasi-redundant conditions are not reduced, leading to quasi-ill-posed linear system. To overcome this problem, the tolerance factor may be increased:

EssentialConditions ecs = (ncross(ud)-ncross(ug)|Sigma=0);
ecs.redundancyTolerance(1.E-10);

Manage reduction methods#

Essential conditions, corresponding to an elimination of some dofs and a reduction of the matrix and may be of right hand side vector. There are several ways to deal with this reduction:

  • pseudo-elimination, changing some matrix coefficients (row or column index corresponding to eliminated dofs) without changing the size of the system but may be the storage (default).

  • real elimination, eliminating row and column corresponding to the eliminated dofs, decreasing the size of the system and changing the storage as a consequence.

  • penalization, adding to the matrix, the matrix constraints affected with a large coefficient, may be changing the storage and having a significant impact on matrix conditionning.

  • dualization, using a Lagrange multiplier related to the constraints system, thus increasing the size of the system by adding matrix blocks.

Up to now, only the pseudo-reduction method is available for any set of essential conditions and the real elimination, penalization methods are only available for Dirichlet condition. The dualization method is not yet implemented in XLiFE++ in a automatic way but can be handled by using an augmented variational formulation.

If necessary, it is possible to change the diagonal coefficient (by default 1) of the pseudo eliminated block matrix by invoking the ReductionMethod object:

BilinearForm auv=intg(omega, grad(u)|grad(v));
EssentialConditions ecs= (u|sigmaM = 0);
TermMatrix A(auv, ecs, _reduction=ReductionMethod(_pseudoReduction,10.), _name="A");

The penalization method adds \(\alpha\) to the diagonal coefficients related to dofs where the Dirichlet conditions act. To use the penalization method, change the method in the ReductionMethod object:

BilinearForm auv=intg(omega, grad(u)|grad(v));
EssentialConditions ecs= (u|sigmaM = 0);
TermMatrix A(auv, ecs, _reduction=ReductionMethod(_penalizationReduction,10000.), _name="A");

Pseudo-reduction and eigen problems#

When dealing with eigen problems be very cautious when there are some essential conditions. Up to now, only the pseudo-reduction method is available to take into account some essential conditions. Just recall that essential conditions are expressed in terms of a linear constraints’ system which is reduced to a minimal form using a QR factorization, say

\[U_e = g_e + \mathbb{S} \, U_r\]

where \(e\) stands for \(n_e\) eliminated dofs, \(r\) for \(n_r\) reduced dofs and \(\mathbb{S}\) is a \(n_e \times n_r\) matrix. The pseudo-reduction of a matrix, say \(\mathbb{A}\), consists in using the previous constraints’ system for all rows corresponding to reduced dofs and replacing rows corresponding to eliminated dofs by the constraints’ system, possibly scaled by a factor \(\alpha\). After reduction, the matrix formally looks like:

\[\begin{split}\begin{pmatrix} \alpha \mathbb{I}_{ee} & \alpha \mathbb{C}_{er} \\ \mathbb{O}_{re} & \mathbb{A}_{rr} \end{pmatrix}\end{split}\]

When \(\mathbb{C}_{er}\) is a null matrix (Dirichlet condition for instance), the eigen vectors are well computed but may be spoiled by the artificial part added to recover constraints. For instance, if the smallest eigen values are required, these artificial eigen values can be shifted by choosing the scale factor \(\alpha\) large enough:

...
BilinearForm auv=intg(omega, grad(u)|grad(v));
BilinearForm uv=intg(omega, u*v);
EssentialConditions ecs = (u|Gamma=0);
TermMatrix A(auv, ecs, _reduction=ReductionMethod(_pseudoReduction,1000.));
TermMatrix M(uv, ecs); // not scaled !
EigenElements eigs=eigenSolve(A, M, _nev=10, _which="SM");

When \(\mathbb{C}_{er}\) is not a null matrix (case of transmission or periodic condition for instance), the previous reduction is no longer working. Another scaling parameter \(\beta\) is available in the ReductionMethod object. It corresponds to the following generalized scaling:

\[\begin{split}\begin{pmatrix} (\alpha + \beta)\mathbb{I}_{ee} & \alpha \mathbb{C}_{er} \\ \mathbb{O}_{re} & \mathbb{A}_{rr} \end{pmatrix}\end{split}\]

Now choosing \(\alpha=0\), the matrix will be reduced to a correct matrix for searching eigen values

\[\begin{split}\begin{pmatrix} \beta \mathbb{I} & \mathbb{O} \\ \mathbb{O} & \mathbb{A}_{rr} \end{pmatrix}\end{split}\]
...
BilinearForm auv=intg(omega, grad(u)|grad(v));
BilinearForm uv=intg(omega, u*v);
EssentialConditions ecs = (u|Gamma) -(u|Sigma)=0;
TermMatrix A(auv, ecs, _reduction=ReductionMethod(_pseudoReduction,0.,1000.));
TermMatrix M(uv, ecs, _reduction=ReductionMethod(_pseudoReduction,0.,1));
EigenElements eigs(A, M, _nev=10, _which="SM");

Although the eigen system has been “reduced” by “eliminating” some dofs, the eigen vectors well satisfy the essential conditions because the applyEssentialConditions() function is automatically called to recover the essential conditions.

Define nodal values of an analytical kernel#

A scalar/vector real/complex nodal value TermMatrix can be constructed form a two variables C++ function (i.e. a kernel \(K\)):

\[\mathbb{G}_{ij}=K(M_i,M_j)\quad M_i,\ M_j\text{ nodes of domains } \Sigma,\ \Gamma.\]

Here is an example:

Kernel K = Helmholtz3dKernel(...);

...
Domain sigma=..., gamma...;
Unknown us(...), ug(...);
...
TermMatrix G(ug, gamma, us, sigma, K, _name="G");
TermMatrix dnG(ug, gamma, us, sigma, grad_y(G)|_ny, _name="dnG");

This construction works also with an operator on a kernel.

About storages#

XLiFE++ manages several matrix storages in several access types. Let’s begin with the latter:

  • row: matrix coefficients are stored row by row one after the other;

  • column: matrix coefficients are stored column by column one after the other;

  • symmetric: matrix coefficients of the lower triangular part are stored row by row, one after the other;

  • dual: diagonal matrix coefficients are stored first, then matrix coefficients of the strict lower triangular part are stored row by row, one after the other, then matrix coefficients of the strict upper triangular part are stored column by column, one after the other.

Let’s now talk about available storages:

  • dense: all matrix coefficients are stored. It is available in all 4 access types.

  • compressed sparse: Only non zero coefficients are stored using compressed sparse format. It is available in all 4 access types. This is the default storage, well suited to to iterative solvers.

  • skyline: Leading and trailing zero coeeficients on row or column are not stored. It is only available in symmetric and dual access types. This is the dedicated storage for matrix factorizations and direct solvers.

  • sparse: Only non zero coefficients are stored using coordinates access. Access type is not concerned with this storage. It is only for export purpose, in particular to plug computations with other softwares.

See Matrix storages for more details.

By default, XLiFE++ chooses an appropriate storage for TermMatrix, according to the type of computation and the the symmetry of the associated BilinearForm:

  • for FEM matrix, the dual or column compressed sparse storage is used,

  • for BEM or spectral matrix, the dense storage is used,

  • matrix in compressed sparse storage move to skyline storage when a direct method is called,

  • when combining matrices, the dual or column compressed sparse storage is used except if the storage size is greater than the half of the dense storage size, in which case it will be used.

To set the storage use the dedicated _storage key. Possible values are csRow, csCol, csDual, csSym, denseRow, denseCol, denseDual, skylineSym, and skylineDual.

BilinearForm auv=intg(omega, grad(u)|grad(v));
TermMatrix A(auv, _storage=skylineDual, _name="A");

Note

TermMatrix related to BEM computation may also use a hierarchical storage based on a dof clustering, see Hmatrix.

Standard operations#

Algebraic operations#

Algebraic operations (\(+=, -=, *=, /=, +, -, *, /\)) are available on TermMatrix and TermVector to build a new one.

...
Space V(_domain=omega, _interpolation=P1, _name="V");
Unknown u(V, _name="u"); TestFunction v(u, _name="v");
TermMatrix M(intg(omega,u*v));               // mass TermMatrix
TermMatrix K(intg(omega,grad(u)|grad(v)));   // rigidity TermMatrix
TermMatrix H = K - k2*M;                     // Helmholtz matrix
theCout << K+M;         // NOT EVALUATED

Important

In order to be more efficient, the linear combination of terms is delayed up to the assign (=) operation or a constructor operation. It means that some expression may not be evaluated and produce warning/error message related to LcTerm class.

Hint

When combining more than two terms, especially for TermMatrix, it is better to write the summation in one step rather than in several steps.

Caution

The product beetween two TermMatrix is always generated in dense format.

The operator * can also be used to do the product TermMatrix x TermVector or TermVector x TermMatrix:

Space V(_domain=omega, _interpolation=P1, _name="V");
Unknown u(V, _name="u"); TestFunction v(u, _name="v");
TermMatrix M(intg(omega,u*v));  // mass TermMatrix
TermVevtor U(u,omega,1.);       // TermVector (1,1,...,1)
theCout<<sqrt(M*U|U);           // L2 norm of U as element of space V
theCout<<norm(U);               // l2 norm of U as an element of Rn

XLiFE++ provides the \(l^2\) norm (Frobenius norm \(\sqrt{\sum|A_{ij}|^2}\)) and the \(l^\infty\) norm (\(\max_{i,j}|A_{i,j}|\)):

...
TermMatrix M(intg(omega,u*v));
theCout<<norm2(M)<<" "<<norminfty(M);

Extraction#

In case of a multiple unknowns TermMatrix, an unknown block may be extracted as follows:

Space V(_domain=omega, _interpolation=P1, _name="V");
Unknown u1(V, _name="u1"), u2(V, _name="u2");
TestFunction v1(u1, _name="v1"), v2(u2, _name="v2");
BilinearForm auv=intg(omega1,grad(u1)|grad(v1))+intg(omega2,grad(u2)|grad(v2));
TermMatrix A(auv);            // multiple unknowns TermMatrix
TermMatrix A11=A(u1,v1);      // extract (u1, v1) part, single unknown TermMatrix

The returned TermMatrix is a copy of the extracted block.

An example with DtN map as transparent boundary condition#

The following example illustrates multiple uses of algebraic operators on TermMatrix.

...
// define spectral space to deal with DtN
Number N=10;
Space Sp(_domain=sigmaP, _basis=Function(cosny, params), _dim=N, _name="cos(n∗pi∗y)");
Unknown phiP(Sp, _name="phiP");
Complexes lambda(N);
for (Number n=0; n<N; n++) lambda[n]=sqrt(Complex(k*k-n*n*pi*pi/(h*h)));
// define TermMatrix with no DtN
BilinearForm auv=intg(omega, grad(u)|grad(v)) - k*k*intg(omega, u*v);
TermMatrix A(auv, _storage=csDual, _name="A");
// contruct DtN TermMatrix using matrix product
BilinearForm euv=intg(sigmaP, phiP*v), fuv=intg(sigmaP, u*phiP);
TermMatrix E(euv, _name="E"), F(fuv, _name="F");
TermMatrix L(phiP, sigmaP, lambda, _name="L"); // diagonal matrix
TermMatrix ELF=E*L*F;
TermMatrix A2=A-i_*ELF;

This DtN approach using product of matrices was the Melina approach. In XLiFE++, it is better to use TensorKernel approach. See Tensor kernels for more details.

There are lot of constructions or operations on TermMatrix, see the Language reference documentation.

Projectors#

When a Projector is constructed (see Spaces and projections), the matrix \(\mathbb A\) and \(\mathbb B\) are computed and the matrix \(\mathbb A\) is factorized. So the projection of TermVector (consistent with the projector) can be computed:

TermVector B2(u2, omega, fx2, _name="B2"); // fx2 user function
TermVector B1=P2toP1(B2, u1);

By specifying the unknown u1 as second argument, the result B1 will have u1 as unknown. It is also possible to omit the unknown argument or to pass the TermVector result as argument:

TermVector B1=P2toP1(B2);
TermVector B1b;
P2toP1(B2, B1b);

When no unknown is given, XLiFE++ will choose the first unknown that has been defined on result space.

By using the projection() function, the standard projections can be done without managing explicitely a Projector object:

TermVector B1=projection(B2, V1);

In that case the Projector object is constructed in the back and kept in memory.

Hint

Only single unknown TermVector can be projected. In case of multiple unknowns TermVector, extract the part related to the unknown of interest before projection (see Restriction).