Computation of right hand sides in Montjoie

VarSourceProblem

VirtualSourceField

IncidentWaveField

VirtualProjectorFEM

VirtualSourceFEM

VolumetricSource

DiracSource

DiffractedWaveSource

TotalWaveSource

Unsteady problems

VirtualTimeSource

TimeSourceHyperbolic

Right hand sides for time-harmonic problems

The following source is implemented in Montjoie (in the case of Helmholtz equation) :

For other equations, the implemented sources are similar. More precisely, the implemented source is equal to the following expression :

where f_G and g_G are functions that you can set to 0. The source f_V may be set equal to a Dirac :

where E₀ is called polarization. This polarization can also be used in other expressions of f_V. It is also possible to compute the projection of a function h in the finite element basis (by calling AddVolumeProjection) :

In order to declare a new source, the following example can be used (see UserSource.hxx)

template<class T, class Dimension>
class MySource
 : public VirtualSourceFEM <T, Dimension>
{
public :
  typedef typename Dimension::R_N R_N; //!< point in R2 or R3
  typedef typename Dimension::VectR_N VectR_N; //!< vector of points

  // constructor with a given problem
  template<class TypeEquation>
  MySource(const EllipticProblem<TypeEquation>& var)
    : VirtualSourceFEM <T, Dimension>(var)
  {
    // you can initialize internal variables if present
  }

  // true if f_V is non null
  bool IsNonNullVolumetricSource(const VectR_N& s) { return true; }
  
  // true if f_G is non null
  bool IsNonNullGradientSource(const VectR_N& s) { return true; }

  // true if g_S is non null (for faces of reference ref)
  bool IsNonNullSurfacicSource(int ref) { return true; }

  // true if g_G is non null (for faces of reference ref)
  bool IsNonNullSurfacicSourceGradient(int ref) { return true; }

  // value of h (for Dirichlet condition or projection of a function h)
  void EvaluateFunction(int i, int j, const R_N& x, Vector<T>& h)
  {
    // h is a vector, its size being equal to the number of unknowns
    // of the equation
    h(0) = 2.0*x(0);
    h(1) = -1.5*x(1) + 4.0;
  }
  
  // value of f_V
  void EvaluateVolumetricSource(int i, int j, const R_N& x, Vector<T>& fV)
  {
    fV.Fill(0);
  }

  // value of f_G
  void EvaluateGradientSource(int i, int j, const R_N& x, Vector<T>& fG)
  {
    fG.Fill(0);
  }

  // value of g_S
  // to retrieve the normale, you can call MatricesElem.GetNormaleQuadratureBoundary(k)
  void EvaluateSurfacicSource(int k, const SetPoints<Dimension>& PointsElem,
                              const SetMatrices<Dimension>& MatricesElem,
                              Vector<T>& gS)
  {
    gS.Fill(0);
  }

  // value of g_G
  template<class Vector1>
  void EvaluateSurfacicSourceGradient(int k, const SetPoints<Dimension>& PointsElem,
                                     const SetMatrices<Dimension>& MatricesElem,
                                     Vector<T>& gG)
  {
    gG.Fill(0);
  }

};

// construction of EllipticProblem object
EllipticProblem<TypeEquation>var;
// ...

// then you call the computation of the source
MySource<Complex_wp, Dimension2> source;

// multiple right hand sides are allowed
Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > f(1);
f(0) = &source;                                                    

// the sources are computed
Vector<Vector<Complexe> > F(1);
var.ComputeGenericSource(F, source);

The following main classes are available in Montjoie

• VolumetricSource : generic source in which you can specify fV, fG, gS and gG
• DiffractedWaveSource : source associated with the computation of a diffracted field
• TotalWaveSource : source associated with the computation of a total field
• ModalSourceBoundary : an eigenmode of a section is provided as a source for this section
• DiracSource : source with only fV non-null and equal to a Dirac
• UserDefinedSource : source defined by the user

The last class is located in files UserSource.hxx/UserSource.cxx and can be modified by anyone to fit its needs.

Right hand sides for unsteady problems

In the case of unsteady problems, there are two classes of sources :

• Source with separable variables, i.e.
  
  In that case, g(x) is treated as described above (for time-harmonic problems). After the discretization in space , we will have a term
  
  where the vector F is precomputed (prior to time iterations). The function h(t) is provided through an instance of the class VirtualTimeSource
• Source with non-separable variables. In that case, the right hand side is re-computed at each iteration (for the required time). The time where the source needs to be evaluated is provided by the method Init of the class VirtualSourceFEM

Methods of class VarSourceProblem (a base class of EllipticProblem)

OnlyOneSource	returns true if there is only one source to compute
GetNbRhs	returns the number of right hand sides
GetSourceType	returns the type of source (i.e. generic source, diffracted field, etc)
GetThresholdSource	returns the threshold used to neglect higher modes
GetIncidentFieldType	returns the type of the incident field
GetNbParameterSource	returns the number of parameters defining a source
GetParameterSource	returns the parameters defining a source
SetParameterSource	sets the parameters defining a source
InitGaussianParameter	initializes a gaussian field
InitRandomGaussianParameter	initializes randomly a gaussian field
InitFftComputation	initialisation of fft computations
ModifyVolumetricSource	volume source f is modified (depending on the chosen formulation)
AddIncidentWave	adds projection of incident field on degrees of freedom to a vector
InitIncidentField	constructs the incident field
ComputeRightHandSide	computation of the source (right hand side of the equation)
GetCoefAB_Infinity	fills the coefficient a and b at infinity (for layered media)
ReadIncidentWaveParam	fills the parameters associated with an incident field
GetNewIncidentField	constructs and returns a new incident field
GetIncidentField	returns the incident field associated with the right hand side
GetNewIncidentProjector	constructs and returns a new class used to project the incident field on degrees of freedom
GetIncidentWaveProjector	returns the incident wave projector associated with the right hand side
AddVolumeProjection	adds projection of a function on degrees of freedom to a vector
SetSurfaceProjection	sets projection of a function on degrees of freedom (of a subset of boundaries) to a vector
AddSurfaceSource	adds to a vector the surface integral of a function against basis functions
AddVolumeSource	adds to a vector the volume integral of a function against basis functions
AddDiracSource	adds to a vector the values of basis functions at a given point (equivalent to a Dirac source)
ComputeGenericSource	computes the right hand side from functions fV, fG, gS, etc
GetNewVolumeSourceFunction	returns an instance of VirtualSourceField corresponding to parameters given in the arguments
ConstructVolumeSourceFunctions	modifies an instance of VolumetricSource (parameters contained in lines containing SRC_VOLUME)
GetNewSurfaceSourceFunction	returns an instance of VirtualSourceField corresponding to parameters given in the arguments
ConstructSurfaceSourceFunctions	modifies an instance of VolumetricSource (parameters contained in lines containing SRC_SURFACE)
FillVariableSource	computes the projection of a function on quadrature points of the surface
GetNewModalSourceEquation	returns the object used for a modal source
GetNewSourceEquationObject	returns the object used to compute the n-th right hand side

Methods of VirtualSourceField

The class VirtualSourceField is a base class to specify a field in Montjoie. When you want to provide your own function f_V in the class VolumetricSource (by calling the function SetVolumeSource, you need to provide a class derived from VirtualSourceField. Below we list the virtual methods of this class.

EvaluateFunction	evaluates the field at a given point
EvaluateGradient	evaluates the gradient of the field at a given point

Currently in Montjoie, two classes derive from VirtualSourceField

• GaussianSourceField : gaussian function
• UniformSourceField : constant function

The class GaussianSource implements a scalar Gaussian function, whereas the class GaussianSourceField implements the multiplication of the scalar Gaussian with a vector. Below, we list the methods of the class GaussianSourceField (which derives from GaussianSource).

GaussianSource	constructors for GaussianSource
GetRadius	returns the distribution radius of the Gaussian
GetCutOffRadius	returns the cut-off radius of the Gaussian
GetCenter	returns the center of the Gaussian
Init	initializes parameters of the Gaussian
GetAmplitude	evaluates the scalar Gaussian function
GetGradAmplitude	evaluates the gradient of the scalar Gaussian function
GetHessianAmplitude	evaluates the hessian matrix of the scalar Gaussian function

Methods of class IncidentWaveField

The class IncidentWaveField is the base class for incident fields implemented in Montjoie. Most of the methods are virtual and overloaded in derived classes. Below we list the methods of this class

GetTimeSource	returns the time source
SetTimeSource	sets the time source
InitElement	initialization before computation of the field on a given element of the mesh
EvaluateFunction	evaluates the incident field
EvaluateFunctionGradient	evaluates the incident field and its gradient
UpdateCoefAB	updates coefficient a and b (for layered media) for the considered element

Currently in Montjoie, the following classes derive from IncidentWaveField

• PlaneWaveIncidentField : incident plane wave, given as
  
• HankelIncidentField : Hankel's function given as
  
  (in 3-D, we consider the spherical Hankel's function)
• GaussianBeamIncidentField : gaussian beam given as
  
  The class can be initialized with the constructor or with the method Init.
• LayeredPlaneWaveIncidentField : plane wave adapted for a layered media given as
  
  on each layer. The incident field satisfies Helmholtz equation :
  
  on each layer. The coefficients a_i and b_i are given in the constructor or by calling Init method.

Methods of class VirtualProjectorFEM

The class VirtualProjectorFEM is the base class for projectors (in the finite element base) implemented in Montjoie. Most of the methods are virtual and overloaded in derived classes. Below we list the methods of this class

GetNbUnknowns	returns the number of unknowns
InitElement	inits computation for a given element
EvaluateFunction	evaluates the fields at a given point
ApplyFFT	Performs FFTs (when the source is decomposed onto modes)
ModifyPoints	points are translated or rotated (for cyclic/periodic domains)
ModifyEvaluationProjection	modification of evaluations of h for cyclic/periodic domains

The class IncidentWaveProjector derives from the class VirtualProjectorFEM. It is used to project an incident field onto the finite element basis.

Methods of class VirtualSourceFEM (class inherited from VirtualProjectorFEM)

The class VirtualSourceFEM is the base class for all sources implemented in Montjoie (VolumetricSource, DiffractedWaveSource, etc). Most of the methods are virtual and overloaded in derived classes.

GetOrigin	returns the center of the source
GetPolarization	returns polarization of the source
GetPolarizationGrad	returns polarization of the source (for gradient terms)
SetPolarization	sets polarization of the source
SetPolarizationGrad	sets polarization of the source (for gradient terms
IsGradientDirac	returns true if the Dirac source is applied to gradient of basis functions
Init	initialization for unsteady problems
InitElement	initialization of some variables for computation of source on an element of the mesh
EvaluateFunction	evaluation of the function used for inhomogeneous Dirichlet condition
EvaluateVolumetricSource	evaluation of the volume source fV
IsNonNullVolumetricSource	returns true if the function fV is non-null
EvaluateGradientSource	evaluation of the volume source fG
IsNonNullGradientSource	returns true if the function fG is non-null
InitSurface	initialization of some variables for computation of source on a face of the mesh
EvaluateSurfacicSource	evaluation of the surface source gV
IsNonNullSurfacicSource	returns true if the function gV is non-null
EvaluateSurfacicSourceGradient	evaluation of the surface source gG
IsNonNullSurfacicSourceGradient	returns true if the function gG is non-null
IsNonNullDirichletSource	returns true if inhomogeneous Dirichlet condition is set
PresenceDirichlet	returns true if Dirichlet condition is present in all the mesh
IsDiracSource	returns true if a Dirac function is used for fV
GetCoefficientVolume	returns the multiplicative coefficient for volume source
ModifyEvaluationVolume	modification of evaluations of fV and fG for cyclic/periodic domains
ModifyEvaluationSurface	modification of evaluations of gV and gG for cyclic/periodic domains

Methods of class VolumetricSource (class inherited from VirtualSourceFEM)

The class VolumetricSource is the class implementing generic sources in Montjoie. Below we list methods specific to this class (other methods are inherited from VirtualSourceFEM).

SetVolumeSource	sets the function fV for a set of references
SetVolumeSourceGrad	sets the function fG for a set of references
SetVolumeSourceFunction	sets a global function fV
NullifyVolumeSourceFunction	nullifies pointers storing the functions fV
SetSurfaceSource	sets the function gS for a set of references
SetVariableSource	the volume source is given directly on quadrature points
SetVariableSurfaceSource	the surface source is given directly on quadrature points
SetVariableGradientSource	the volume source is given directly on quadrature points (for gradient of basis functions)

Class DiracSource (class inherited from VirtualSourceFEM)

The class DiracSource is the class implementing Dirac sources in Montjoie. The origin of the Dirac can be retrieved by the method GetOrigin. This kind of source is computed with a line beginning with TypeSource = SRC_DIRAC in the data file.

Class DiffractedWaveSource (class inherited from VirtualSourceFEM)

The class DiffractedWaveSource is the class implementing the computation of scattered fields in Montjoie. This kind of source is computed with a line beginning with TypeSource = SRC_DIFFRACTED_FIELD in the data file. It is implemented for only some equations (such as Helmholtz or Maxwell's equations). By writing

we obtain an equation with the diffracted field and source terms. For example, we have an hetereogeneous Dirichlet condition

Class TotalWaveSource (class inherited from VirtualSourceFEM)

The class TotalWaveSource is the class implementing the computation of scattered fields in Montjoie. This kind of source is computed with a line beginning with TypeSource = SRC_TOTAL_FIELD in the data file. It is implemented for only some equations (such as Helmholtz or Maxwell's equations). Contratry to DiffractedWaveSource, in this class the total field is computed (and not the diffracted field). Since only the diffracted field satisfies the outgoing radiation condition, we write

for PML layers and/or absorbing conditions to obtain an equation with the total field and source terms. For example, we have an hetereogeneous radiation condition

Class ModalSourceBoundary (class inherited from VirtualSourceFEM)

The class ModalSourceBoundary is the class implementing modal sources in Montjoie. We compute an eigenmode of the section, and we use this eigenmode as a source. More details will be explained later.

Methods of class VirtualTimeSource

The class VirtualTimeSource is the base class implementing time-dependency of sources in Montjoie. Below we list methods specific to this class

Evaluate	evaluates the function h(t)
EvaluateDerivative	evaluates the derivative of the function h(t)

The following classes derive from this base class :

• TimeRickerSource : Ricker source
• DerivativeTimeRickerSource : derivative of a Ricker
• TimeModifiedRickerSource : modified Ricker source
• TimeGaussianSource : gaussian source
• TimeHarmonicSource : sinusoidal source
• TimeSinuGaussianSource : sinusoidal source modulated by a gaussian
• TimeFileSource : source with discrete values given on a file
• TimeUserSource : source defined by the user (file UserSource.cxx)

The expression of the functions are given in the description of the keyword TemporalSource.

Methods of class TimeSourceHyperbolic

The class TimeSourceHyperbolic encapsulates an instance of VirtualTimeSource in order to compute the n-th derivative of the time function h(t). The function h(t) is interpolated with high-order polynomial to obtain this n-th derivative. Below, we list the methods of the class TimeSourceHyperbolic

Init	sets the function h(t)
EvaluateDerivative	computes the interpolation of the function h(t)

OnlyOneSource

This method returns true if there is only one source to be computed.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // then you can factorize the finite element matrix
    GlobalGenericMatrix<Complex_wp> nat_mat;
    solver->PerformFactorizationStep(nat_mat);

  if (var.OnlyOneSource())
    {
      // only one right hand side => a vector
      VectComplex_wp rhs;
      var.ComputeRightHandSide(rhs);

     // and solve the linear system A x = b
     VectComplex_wp x;
     x = rhs; solver->ComputeSolution(x);
    }
  else
    {
      // multiple right hand sides => a matrix
      // each source is stored in a column
      Matrix<Complex_wp, General, ColMajor> rhs;
      var.ComputeRightHandSide(rhs);

     // and solve the linear systems A x = b (for each column)
     Matrix<Complex_wp, General, ColMajor> x;
     x = rhs; solver->ComputeSolution(x);

    }
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNbRhs

This method returns the number of right hand sides that will be computed.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // then you can factorize the finite element matrix
    GlobalGenericMatrix<Complex_wp> nat_mat;
    solver->PerformFactorizationStep(nat_mat);

  if (var.OnlyOneSource())
    {
      // only one right hand side => a vector
      VectComplex_wp rhs;
      var.ComputeRightHandSide(rhs);

     // and solve the linear system A x = b
     VectComplex_wp x;
     x = rhs; solver->ComputeSolution(x);
    }
  else
    {
      // multiple right hand sides => a matrix
      // each source is stored in a column
      Matrix<Complex_wp, General, ColMajor> rhs;
      var.ComputeRightHandSide(rhs);

     // number of columns should be equal to GetNbRhs
     int nb_rhs = var.GetNbRhs();

     // and solve the linear systems A x = b (for each column)
     Matrix<Complex_wp, General, ColMajor> x;
     x = rhs; solver->ComputeSolution(x);

    }
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetSourceType

This method returns the type of the n -th source. It is an integer that can be equal to :

• SRC_NULL : null source
• SRC_VOLUME : generic source
• SRC_TOTAL_FIELD : source corresponding to the computation of the total field
• SRC_DIFFRACTED_FIELD : source corresponding to the computation of the scattered field
• SRC_USER : user-defined source (file UserSource.cxx)
• SRC_DIRAC : Dirac source

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // then you can factorize the finite element matrix
    GlobalGenericMatrix<Complex_wp> nat_mat;
    solver->PerformFactorizationStep(nat_mat);

  if (var.OnlyOneSource())
    {
      // only one right hand side => a vector
      VectComplex_wp rhs;
      var.ComputeRightHandSide(rhs);

     // and solve the linear system A x = b
     VectComplex_wp x;
     x = rhs; solver->ComputeSolution(x);

      // type of source ?
      int type_source = var.GetSourceType(0);
    }
  else
    {
      // multiple right hand sides => a matrix
      // each source is stored in a column
      Matrix<Complex_wp, General, ColMajor> rhs;
      var.ComputeRightHandSide(rhs);

     // loop over sources
     int nb_rhs = var.GetNbRhs();
     for (int k = 0; k < nb_rhs; k++)
       {
         int type_source = var.GetSourceType(k);
         if (type_source == var.SRC_TOTAL_FIELD)
           cout << "Computation of a total field for source " << k << endl;
       }      

    }
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetThresholdSource

This method returns the threshold used to neglect modes. If the source for the n-th has a norm below this threshold, this mode is dropped.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // threshold for sources
    Real_wp eps = var.GetThresholdSource();
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetIncidentFieldType

This method returns the type of the incident field associated with the n -th source. It is an integer that can be equal to :

• INCIDENT_PLANE_WAVE : plane wave
• INCIDENT_PLANE_WAVE_CPLX : plane wave with a complex wave number
• INCIDENT_GAUSSIAN_BEAM : gaussian beam
• INCIDENT_HANKEL : Hankel function
• INCIDENT_LAYERED_PLANE_WAVE : plane wave for layered media
• INCIDENT_NONE : no incident field

It is relevant only if the source is equal to SRC_TOTAL_FIELD or SRC_DIFFRACTED_FIELD. The second syntax can be used to obtain the integer associated with the name of the incident field.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // then you can factorize the finite element matrix
    GlobalGenericMatrix<Complex_wp> nat_mat;
    solver->PerformFactorizationStep(nat_mat);

  if (var.OnlyOneSource())
    {
      // only one right hand side => a vector
      VectComplex_wp rhs;
      var.ComputeRightHandSide(rhs);

     // and solve the linear system A x = b
     VectComplex_wp x;
     x = rhs; solver->ComputeSolution(x);

      // type of source ?
      int type_source = var.GetSourceType(0);
    }
  else
    {
      // multiple right hand sides => a matrix
      // each source is stored in a column
      Matrix<Complex_wp, General, ColMajor> rhs;
      var.ComputeRightHandSide(rhs);

     // loop over sources
     int nb_rhs = var.GetNbRhs();
     for (int k = 0; k < nb_rhs; k++)
       {
         int type_source = var.GetSourceType(k);
         if (type_source == var.SRC_TOTAL_FIELD)
           {
           cout << "Computation of a total field for source " << k << endl;
           if (var.GetIncidentFieldType(k) == var.INCIDENT_PLANE_WAVE)
             cout << " Incident field is a plane wave" << endl;
       }      

    }

   // second syntax
   int type_incident = VarSourceProblem_Base::GetIncidentFieldType("PLANE_WAVE");
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNbParameterSource

This method returns the number of lines (in the datafile) defining the n -th source. It corresponds to all the lines TypeSource associated with the n -the source.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // number of lines for the first source (n = 0)
    int n = var.GetNbParameterSource);

    // parameters for each line
    for (int k = 0; k < n; k++)
      cout << "Parameters = " << var.GetParameterSource(0, k);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetParameterSource

This method returns the k -th line (in the datafile) defining the n -th source. This line is usually set by inserting a line TypeSource in the datafile.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // number of lines for the first source (n = 0)
    int n = var.GetNbParameterSource);

    // parameters for each line
    for (int k = 0; k < n; k++)
      cout << "Parameters = " << var.GetParameterSource(0, k);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

SetParameterSource

This method sets the lines defining the sources. It corresponds to all the lines TypeSource in the data file.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // if you want to set manually a single source
    Vector<VectString> param(1);
    param(0).Reallocate(5);
    param(0)(0) = "SRC_VOLUME"; param(0)(1) = "GAUSSIAN"; param(0)(2) = "0.0"; param(0)(3) = "0.0"; param(0)(4) = "1.0";
    var.SetParameterSource(param);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

InitGaussianParameter

This method initializes a gaussian field from parameters contained in param. nb is the starting index that will be incremented. (param(nb), param(nb+1)) is the center of the gaussian (in 3-D, param(nb+2) is also used). param(nb+dim) is the radius of the gaussian.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // parameters of a Gaussian
    VectString param(3);
    param(0) = "0.0"; param(1) = "0.0"; // center
    param(2) = "1.0"; // radius

    GaussianSource<Dimension2> f;
    int nb = 0;
    var.InitGaussianParameter(f, param, nb);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

InitRandomGaussianParameter

This method initializes a random gaussian field. nb is the starting index that will be incremented. [param(nb), param(nb+1)] is the interval for the x-coordinate of the center. [param(nb+2), param(nb+3)] is the interval for the y-coordinate of the center. parm(nb+4) in 2-D is the radius of the gaussian. Only the center of the gaussian is chosen randomly (in the provided intervals).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // intervals for a Gaussian
    VectString param(5);
    param(0) = "-2.0"; param(1) = "1.0"; // x chosen in [-2, 1]
    param(2) = "-5.0"; param(3) = "-3.0"; // y chosen in [-5, -3]
    param(4) = "1.0"; // radius

    GaussianSource<Dimension2> f;
    int nb = 0;
    var.InitRandomGaussianParameter(f, param, nb);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

InitFftComputation

This method initializes the fft object in the case where the source is decomposed in Fourier modes (such as for a cyclic domain).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    FftInterface<Complex_wp> fft;
    var.InitFftComputation(fft);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ModifyVolumetricSource

This method modifies the evaluation of the volumetric source at a given point x. A modification of this vector is needed if you have a specific formulation of a given equation. In regular use, this method should not be called.

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

AddIncidentWave

This method adds to the vector u the projection of the incident field on degrees of freedom (multiplied by a coefficient alpha).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);


    // construction of incident field (projected on degrees of freedom)
    VectComplex_wp u(var.GetNbDof()); u.Zero();
    var.InitIncidentField();
    var.AddIncidentWave(Complex_wp(1, 0), u);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

InitIncidentField

This method constructs the incident field. It is automatically called in the function ComputeRightHandSide. If you did not call ComputeRightHandSide previously, you will need to call InitIncidentField if you want to evaluate the incident field.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);


    // construction of incident field (projected on degrees of freedom)
    VectComplex_wp u(var.GetNbDof()); u.Zero();
    var.InitIncidentField();
    var.AddIncidentWave(Complex_wp(1, 0), u);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ComputeRightHandSide

This method computes the right hand side of the linear system to solve. The first argument is the right hand side. The second argument is optional and meaningful only for parallel computation. If this argument is true, the right hand side will be assembled (between processors). On regular use, the assembly phase is not needed since it is performed when the method ComputeSolution is called.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   GlobalGenericMatrix<Complex_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients

   // to factorize the matrix (or prepare the computation if an iterative solver is selected)
   solver->PerformFactorizationStep(nat_mat);

    // to compute the right hand side
    VectComplex_wp b(var.GetNbDof());
    var.ComputeRightHandSide(b);

    // and solve the linear system A x = b
    VectComplex_wp x = b; solver->ComputeSolution(x);
   
   // when you no longer need the solver, you can release the memory
   delete solver;

}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetCoefAB_Infinity

This method fills the coefficients a and b at infinty. The incident field of a layered media satisfies

where a_i and b_i are coefficients that depend on the layer.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // to know a_i and b_i at infinity
    Real_wp a_inf, b_inf;
    var.GetCoefAB_Infinity(a_inf, b_inf);

}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ReadIncidentWaveParam

This method fills parameters defining the incident field by using values contained in param (or by taking values from the class).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // parameters of the incident field
    // if param is empty, we take the values stored in the class
   VectString param;
   R2 kwave, R2_Complex_wp kc; R2 origin; Real_wp omega, w;
   var.ReadIncidentWaveParam(param, kwave, kc, origin, omega, w);

}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNewIncidentField

This method constructs a new incident field from values stored in param (or values in the class) and returns the incident field constructed. The integer n can be equal to :

• INCIDENT_PLANE_WAVE : the incident field is a plane wave
• INCIDENT_PLANE_WAVE_CPLX : the incident field is a plane wave with a complex wave vector
• INCIDENT_GAUSSIAN_BEAM : the incident field is a gaussian beam
• INCIDENT_HANKEL : the incident field is a Hankel function
• INCIDENT_LAYERED_PLANE_WAVE : the incident field is a combination of plane waves adapted for a layered media

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // parameters of the incident field
    // if param is empty, we take the values stored in the class
   Vector<VectString> param(1); // only param(0) is used

   IncidentWaveField<Complex_wp, Dimension2>* u_inc;
   u_inc = var.GetNewIncidentField(VarSourceProblem_Base::INCIDENT_PLANE_WAVE, param, Complex_wp(0));

   // evaluation of the incident field on a point
   R2 point(0.5, 0.8); Complex_wp val_u;
   u_inc->EvaluateFunction(point, val_u);
 
   // when you no longer need u_inc, you can delete it
   delete u_inc;
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetIncidentField

This method returns the incident field associated with the n-th source. It works only if InitIncidentField or ComputeRightHandSide has been called previously.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   var.InitIncidentField(); // to compute incident fields for all right hand sides
   IncidentWaveField<Complex_wp, Dimension2>* u_inc;
   u_inc = var.GetIncidentField(0, Complex_wp(0)); // first incident field

   // evaluation of the incident field on a point
   R2 point(0.5, 0.8); Complex_wp val_u;
   u_inc->EvaluateFunction(point, val_u);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNewIncidentProjector

This method constructs a new projector for the incident field from values stored in param (or values in the class) and returns the incident field constructed. The integer n can be equal to :

• INCIDENT_PLANE_WAVE : the incident field is a plane wave
• INCIDENT_PLANE_WAVE_CPLX : the incident field is a plane wave with a complex wave vector
• INCIDENT_GAUSSIAN_BEAM : the incident field is a gaussian beam
• INCIDENT_HANKEL : the incident field is a Hankel function
• INCIDENT_LAYERED_PLANE_WAVE : the incident field is a combination of plane waves adapted for a layered media

The difference with GetNewIncidentField is that the latter implements a scalar incident field, whereas the former implements a vector incident field (adapted to the considered equation).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // parameters of the incident field
    // if param is empty, we take the values stored in the class
   Vector<VectString> param(1); // only param(0) is used

   IncidentWaveProjector<Complex_wp, Dimension2>* u_inc;
   u_inc = var.GetNewIncidentProjector(VarSourceProblem_Base::INCIDENT_PLANE_WAVE, param, Complex_wp(0));

   // evaluation of the incident field on a point
   R2 point(0.5, 0.8); VectComplex_wp vec_u;
   u_inc->EvaluateFunction(point, vec_u);
 
   // when you no longer need u_inc, you can delete it
   delete u_inc;
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetIncidentWaveProjector

This method returns the incident projector associated with the n-th source. It works only if InitIncidentField or ComputeRightHandSide has been called previously. The difference with GetNewIncidentProjector is that the latter implements a scalar incident field, whereas the former implements a vector incident field (adapted to the considered equation).

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   var.InitIncidentField(); // to compute incident fields for all right hand sides
   IncidentWaveProjector<Complex_wp, Dimension2>* u_inc;
   u_inc = var.GetIncidentProjector(0, Complex_wp(0)); // first incident field

   // evaluation of the incident field on a point
   R2 point(0.5, 0.8); VectComplex_wp vec_u;
   u_inc->EvaluateFunction(point, vec_u);
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

AddVolumeProjection

This method adds the projection of alpha*f on degrees of freedom to the vector b. This procedure can be completed for a list of sources, that's why b is a vector of vectors.

Example :

    // you need to derive a class from the base class VirtualProjectorFEM
    class MyFunction : public VirtualProjectorFEM<Complex_wp, Dimension2>
    {
    public :
      void InitElement(int i, const VectR2& s) {}
      void EvaluateFunction(int i, int j, const R2& x, VectComplex_wp& f)
      {
        // definition of your function
        f(0) = exp(-x(0)*x(0) - x(1)*x(1)); // here f = exp(-x^2 - y^2)
      }
    
    };
    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the projection of f on degrees of freedom
   Vector<VectComplex_wp> f_proj(1); // only one vector here
   f_proj(0).Reallocate(var.GetNbDof()); f_proj(0).Zero();
   
   MyFunction my_f; // instance of the class defining f
   Vector<VirtualProjectorFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we project f on degrees of freedom
  var.AddVolumeProjection(Complex_wp(1, 0), f_proj, vec_f);
  DISP(f_proj(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

SetSurfaceProjection

This method sets the projection of alpha*f on degrees of freedom to the vector b, only for Dirichlet dofs. This procedure can be completed for a list of sources, that's why b is a vector of vectors.

Example :

    // you need to derive a class from the base class VirtualSourceFEM
    class MyFunction : public VirtualSourceFEM<Complex_wp, Dimension2>
    {
    public :
      void EvaluateFunction(int i, int j, const R2& x, VectComplex_wp& f)
      {
        // definition of your function
        f(0) = exp(-x(0)*x(0) - x(1)*x(1)); // here f = exp(-x^2 - y^2)
      }
    
    };
    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the projection of f on Dirichlet degrees of freedom
   Vector<VectComplex_wp> f_proj(1); // only one vector here
   f_proj(0).Reallocate(var.GetNbDof()); f_proj(0).Zero();
   
   MyFunction my_f; // instance of the class defining f
   Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we project f on degrees of freedom (only for Dirichlet dofs)
  var.SetSurfaceProjection(Complex_wp(1, 0), f_proj, vec_f);
  DISP(f_proj(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

AddSurfaceSource

This method adds the surface integrals of alpha*f against basis functions to the vector b. This procedure can be completed for a list of sources, that's why b is a vector of vectors.

Example :

    // you need to derive a class from the base class VirtualSourceFEM
    class MyFunction : public VirtualSourceFEM<Complex_wp, Dimension2>
    {
    public :
      bool IsNonNullSurfacicSource(int ref) { return true; } // you can specify a boundary where the integral will be computed
      void EvaluateSurfacicSource(int k, const SetPoints<Dimension2>& pts, const SetMatrices<& mat, VectComplex_wp& f)
    {
        R2 x = pts.GetPointQuadratureBoundary(k);
        R2 n = mat.GetNormaleQuadratureBoundary(k); // if you need outgoing normale
        // definition of your function
        f(0) = exp(-x(0)*x(0) - x(1)*x(1)); // here f = exp(-x^2 - y^2)
      }


      bool IsNonNullSurfacicSourceGradient(int ref) { return true; } // you can also integrate against gradient of basis functions
      void EvaluateSurfacicSourceGradient(int k, const SetPoints<Dimension2>& pts, const SetMatrices<& mat, VectComplex_wp& f)
      {
        f.Zero();
      }

    };
    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the surface integral of f against basis functions
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();
   
   MyFunction my_f; // instance of the class defining f
   Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we compute the surface integral of f with basis functions
  var.AddSurfaceSource(Complex_wp(1, 0), rhs, vec_f);
  DISP(rhs(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

AddVolumeSource

This method adds the volume integrals of alpha*f against basis functions to the vector b. This procedure can be completed for a list of sources, that's why b is a vector of vectors.

Example :

    // you need to derive a class from the base class VirtualSourceFEM
    class MyFunction : public VirtualSourceFEM<Complex_wp, Dimension2>
    {
    public :
      bool IsNonNullVolumetricSource(const VectR2& s) { return true; }
      void EvaluateVolumetricSource(int i, int j, const R2& x, VectComplex_wp& f)
    {
        // definition of your function
        f(0) = exp(-x(0)*x(0) - x(1)*x(1)); // here f = exp(-x^2 - y^2)
      }

      bool IsNonNullGradientSource(const VectR2&) { return true; } // you can also integrate against gradient of basis functions
      void EvaluateGradientSource(int i, int j, const R2& x, VectComplex_wp& f)
      {
        f.Zero();
      }

    };
    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the volume integral of f against basis functions
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();
   
   MyFunction my_f; // instance of the class defining f
   Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we compute the integral of f with basis functions
  var.AddVolumeSource(Complex_wp(1, 0), rhs, vec_f);
  DISP(rhs(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

AddDiracSource

This method adds alpha*phi(x₀) to the vector b (equivalent to integrate against a Dirac function). x₀ is the origin of the Dirac. This procedure can be completed for a list of sources, that's why b is a vector of vectors.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the volume integral of the Dirac against basis functions
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   Vector<VectString> param(1);
   param(0).Reallocate(3);
   param(0) = "-2.0"; param(1) = "-3.0"; // coordinates of the orign x_0
   DiracSource<HelmholtzEquation<Dimension2> > my_f(var, param); // instance of the class defining Dirac's source
   Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we compute the integral of Dirac with basis functions
  var.AddDiracSource(Complex_wp(1, 0), rhs, vec_f);
  DISP(rhs(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ComputeGenericSource

This method computes multiple right hand sides in the first argument. In the second argument, the user has to provide properties about functions to integrate. The third argument is optional. If equal to true, the right hand sides are assembled (relevant only for parallel computations).

Example :

// example of a class defining the different source terms
class MySource
 : public VirtualSourceFEM <Complex_wp, Dimension2>
{
public :
  typedef typename Dimension::R_N R_N; //!< point in R2 or R3
  typedef typename Dimension::VectR_N VectR_N; //!< vector of points

  // constructor with a given problem
  template<class TypeEquation>
  MySource(const EllipticProblem<TypeEquation>& var)
    : VirtualSourceFEM <T, Dimension>(var)
  {
    // you can initialize internal variables if present
  }

  // true if f_V is non null
  bool IsNonNullVolumetricSource(const VectR2& s) { return true; }
  
  // true if f_G is non null
  bool IsNonNullGradientSource(const VectR2& s) { return true; }

  // true if g_S is non null (for faces of reference ref)
  bool IsNonNullSurfacicSource(int ref) { return true; }

  // true if g_G is non null (for faces of reference ref)
  bool IsNonNullSurfacicSourceGradient(int ref) { return true; }

  // value of h (for Dirichlet condition or projection of a function h)
  void EvaluateFunction(int i, int j, const R2& x, Vector<T>& h)
  {
    // h is a vector, its size being equal to the number of unknowns
    // of the equation
    h(0) = 2.0*x(0);
    h(1) = -1.5*x(1) + 4.0;
  }
  
  // value of f_V
  void EvaluateVolumetricSource(int i, int j, const R2& x, Vector<T>& fV)
  {
    fV.Fill(0);
  }

  // value of f_G
  void EvaluateGradientSource(int i, int j, const R2& x, Vector<T>& fG)
  {
    fG.Fill(0);
  }

  // value of g_S
  // to retrieve the normale, you can call MatricesElem.GetNormaleQuadratureBoundary(k)
  void EvaluateSurfacicSource(int k, const SetPoints<Dimension2>& PointsElem,
                              const SetMatrices<Dimension2>& MatricesElem,
                              Vector<T>& gS)
  {
    gS.Fill(0);
  }

  // value of g_G
  template<class Vector1>
  void EvaluateSurfacicSourceGradient(int k, const SetPoints<Dimension2>& PointsElem,
                                     const SetMatrices<Dimension2>& MatricesElem,
                                     Vector<T>& gG)
  {
    gG.Fill(0);
  }

};

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   MySource my_f(var; // instance of the class defining the source
   Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &my_f;

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
  DISP(rhs(0));
}

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNewVolumeSourceFunction

This method returns a pointer to a class defining the function (for a volume source) depending on the provided parameters. If the source is only known at quadrature points, the returned pointer is null and the values of the source on quadrature points are placed in the object var. The parameters are detailed in the description of TypeSource.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    IVect ref_vol(2); ref_vol(0) = 1; ref_vol(1) = 3; // source non-null for references 1 and 3
    VectString param(5); int nb = 0; VectComplex_wp polar;
    param(0) = "GAUSSIAN"; param(1) = "0.0"; param(2) = "0.0"; param(3) = "1.0";
    VolumetricSource_Base<Complex_wp, Dimension2> var;
    VirtualSourceField<Complex_wp, Dimension2>* f = var.GetNewVolumeSourceFunction(ref_vol, param, nb, polar, var);

 }

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ConstructVolumeSourceFunctions

This method will construct all volumes functions defined in the parameters given as argument. These functions are stored in the output argument var.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    Vector<VectString> all_param(1);
    VectString param(5); int nb = 0; VectComplex_wp polar;
    param(0) = "GAUSSIAN"; param(1) = "0.0"; param(2) = "0.0"; param(3) = "1.0";
    all_param(0) = param;
    VolumetricSource_Base<Complex_wp, Dimension2> var;
     var.ConstructVolumeSourceFunctions(all_param, var);

 }

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNewVolumeSourceFunction

This method returns a pointer to a class defining the function (for a surface source) depending on the provided parameters. If the source is only known at quadrature points, the returned pointer is null and the values of the source on quadrature points are placed in the object var. The parameters are detailed in the description of TypeSource.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    IVect ref_surf(2); ref_surf(0) = 1; ref_surf(1) = 3; // source non-null for references 1 and 3
    VectString param(5); int nb = 0; VectComplex_wp polar;
    param(0) = "GAUSSIAN"; param(1) = "0.0"; param(2) = "0.0"; param(3) = "1.0";
    VolumetricSource_Base<Complex_wp, Dimension2> var;
    VirtualSourceField<Complex_wp, Dimension2>* f = var.GetNewSurfacSourceFunction(ref_surf, param, nb, polar, var);

 }

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

ConstructSurfaceSourceFunctions

This method will construct all surface functions defined in the parameters given as argument. These functions are stored in the output argument var.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    Vector<VectString> all_param(1); VectString param(5);
    param(0) = "GAUSSIAN"; param(1) = "0.0"; param(2) = "0.0"; param(3) = "1.0";
    all_param(0) = param;
    VolumetricSource_Base<Complex_wp, Dimension2> var;
     var.ConstructSurfaceSourceFunctions(all_param, var);

 }

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

FillVariableSource

This method projects a field known for a given set of points to all quadrature points of faces. It is assumed that the field is known on a surface, the array Points contain the points where the field is known, and val contain the value of the field on these points. The method computes the interpolation of the field on the quadrature points of the faces of reference ref.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // defining a field on discrete points
    VectR2 Points(4);
    Points(0).Init(0.0, 0.0);
    Points(1).Init(0.0, 0.2);
    Points(2).Init(0.0, 0.6);
    Points(3).Init(0.0, 1.0);

    int nb_comp = 2; // number of components of the field
    Matrix<Complex_wp> val(Points.GetM(), nb_comp);
    val(0, 0) = 0.2; val(0, 1) = 0.5;
    // etc

    // only faces such that the reference is equal to ref(i) are considered
    IVect ref(2); ref(0) = 1; ref(2) = 3;

    // allocating the output arrays
    Vector<Vector<VectComplex_wp> > > eval(nb_comp);
    eval(0).Reallocate(var.mesh.GetNbBoundaryRef());
    eval(1).Reallocate(var.mesh.GetNbBoundaryRef());

    // then we compute the interpolation on quadrature points of faces of reference 1 and 3
    var.FillVariableSource(ref, Points, val, eval);

 }

Location :

Source/DefineSourceElliptic.hxx
Source/DefineSourceElliptic.cxx

GetNewSourceEquationObject

This method constructs the object used to compute the n-th right hand side. It is used in the method ComputeRightHandSide.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // object to compute first right hand side ?
    VirtualSourceFEM<Complex_wp, DImension2>* src = var.GetNewSourceEquationObject(0);

 }

Location :

Source/DefineSourceElliptic.hxx
Harmonic/VarHarmonic.cxx

GetNewModalSourceEquation

This method constructs the object used to compute the a modal source for the considered equation.

Example :

    
int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    // object to compute a modal source
    ModalSourceBoundary_Dim<Complex_wp, DImension2>* src = var.GetNewModalSourceEquation();

 }

Location :

Source/DefineSourceElliptic.hxx
Harmonic/VarHarmonicInline.cxx

EvaluateFunction

This method evaluates the field at a given point.

Example :

// to define your own function :
class MyFunction : public VirtualSourceField<Real_wp, Dimension2>
{
  public :
    virtual void EvaluateFunction(const R2& x, VectReal_wp& f) const
    { f(0) = cos(x(0)*x(1)); f(1) = sin(x(1)) - cos(x(0)); }
};

int main()
{
   MyFunction func;
   VectReal_wp evalF(2); R2 x(0.4, 0.6);
   func.EvaluateFunction(x, evalF);
   cout << "Value of f at x = " <lt; x << " is equal to " << evalF << endl;
   
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateGradient

This method evaluates the gradient of a field at a given point. The gradient of each component are placed successively in the output vector. This function is not mandatory, and not used for most of equations. There are some equations (with a specific formulation) that will need the gradient.

Example :

// to define your own function :
class MyFunction : public VirtualSourceField<Real_wp, Dimension2>
{
  public :
    virtual void EvaluateFunction(const R2& x, VectReal_wp& f) const
    { f(0) = cos(x(0)*x(1)); f(1) = sin(x(1)) - cos(x(0)); }

    virtual void EvaluateGradient(const R2& x, VectReal_wp& df) const
    { df(0) = -x(1)*sin(x(0)*x(1)); df(1) = -x(0)*sin(x(0)*x(1)); // gradient of first component
      df(2) = sin(x(0)); df(3) = cos(x(1)); // then second component
    }
};

int main()
{
   MyFunction func;
   VectReal_wp evalF(2); R2 x(0.4, 0.6);
   func.EvaluateFunction(x, evalF);
   cout << "Value of f at x = " <lt; x << " is equal to " << evalF << endl;

   VectReal_wp evalDF(2); 
   func.EvaluateGradient(x, evalDF);
   cout << "Gradient of f at x = " <lt; x << " is equal to " << evalDF << endl;

}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

constructors for GaussianSource

The second constructor initializes the Gaussian with a center and a distribution radius r₁. The scalar field is equal to

where r is the distance to the center. The parameter alpha is set such that

r₂ is a cut-off radius, initialized to r₁ in the constructor. The parameter beta is given as

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5;
  GaussianSource<Dimension2> f(center, radius);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetRadius

This methods returns the distribution radius of the gaussian.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5;
  GaussianSource<Dimension2> f(center, radius);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));

  // if you want to retrieve the radius
  radius = f.GetRadius(); // returns r1
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetCenter

This methods returns the center of the gaussian.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5;
  GaussianSource<Dimension2> f(center, radius);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));

  // if you want to retrieve the center
  center = f.GetCenter(); 
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetCutOffRadius

This methods returns the cut-off radius of the gaussian. For r greater than the cut-off radius, the gaussian is replaced by zero.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5; Real_wp r2 = 1.0;
  GaussianSource<Dimension2> f;
  f.Init(center, radius, r2);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));

  // if you want to retrieve the cut-off radius
  r2 = f.GetCutOffRadius(); // returns r2
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

Init

This method initializes the gaussian with the center, and the two radii. You can look at the constructor for the mathematical expressions.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5; Real_wp r2 = 1.0;
  GaussianSource<Dimension2> f;
  f.Init(center, radius, r2);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetAmplitude

This method evaluates the gaussian at a point x.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5; Real_wp r2 = 1.0;
  GaussianSource<Dimension2> f;
  f.Init(center, radius, r2);

  // then you can evaluate the gaussian at any point
  Real_wp eval = f.GetAmplitude(R2(-1, -3));
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetGradAmplitude

This method evaluates the gradient of the gaussian at a point x.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5; Real_wp r2 = 1.0;
  GaussianSource<Dimension2> f;
  f.Init(center, radius, r2);

  // then you can evaluate the gaussian at any point and its gradient
  Real_wp eval; R2 eval_grad;
  f.GetGradAmplitude(R2(-1, -3), eval, eval_grad);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetHessianAmplitude

This method evaluates the hessian matrix of the gaussian at a point x.

Example :

int main()
{
  R2 center(1, 2); Real_wp radius = 0.5; Real_wp r2 = 1.0;
  GaussianSource<Dimension2> f;
  f.Init(center, radius, r2);

  // then you can evaluate the gaussian at any point and its gradient (and hessian)
  Real_wp eval; R2 eval_grad; Matrix2_2sym eval_hess;
  f.GetHessianAmplitude(R2(-1, -3), eval, eval_grad, eval_hess);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetTimeSource, SetTimeSource

These methods can be used to get/set the time function used for plane waves. In time-domain, the plane wave, need a time source, i.e. the incident field is given as

where h is the time function. This function is usually specified in the data file (keyword TemporalSource).

Example :

int main()
{
  R2 origin, k(3.0, 1.4);
  PlaneWaveIncidentField<Real_wp, Dimension2> u_inc(origin, k);
  Real_wp freq = 1.0;
  TimeRickerSource ricker(freq);
  u_inc.SetTimeSource(&ricker);
  
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

InitElement

This function is called in Montjoie to specify the element where the incident field will be evaluated. The function EvaluateFunction is called with points x that belong to the element.

Example :

int main()
{
  R2 origin, k(3.0, 1.4);
  PlaneWaveIncidentField<Real_wp, Dimension2> u_inc(origin, k);

  int i = 2; VectR2 s; mesh.GetVerticesElement(i, s);
  u_inc.InitElement(i, s);

  // then EvaluateFunction will be called with points inside element i
  Complex_wp f; u_inc.EvaluateFunction(R2(1, 1), f);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateFunction

This function evaluates the incident field at a given point.

Example :

int main()
{
  R2 origin, k(3.0, 1.4);
  PlaneWaveIncidentField<Real_wp, Dimension2> u_inc(origin, k);

  int i = 2; VectR2 s; mesh.GetVerticesElement(i, s);
  u_inc.InitElement(i, s);

  // then EvaluateFunction will be called with points inside element i
  Complex_wp f; u_inc.EvaluateFunction(R2(1, 1), f);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateFunctionGradient

This function evaluates the incident field and its gradient at a given point.

Example :

int main()
{
  R2 origin, k(3.0, 1.4);
  PlaneWaveIncidentField<Real_wp, Dimension2> u_inc(origin, k);

  int i = 2; VectR2 s; mesh.GetVerticesElement(i, s);
  u_inc.InitElement(i, s);

  // then EvaluateFunctionGradient will be called with points inside element i
  Complex_wp f; R2_Complex_wp df;
  u_inc.EvaluateFunctionGradient(R2(1, 1), f, df);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

UpdateCoefAB

This function retrieves the coefficients a and b inside the considered element. It is meaningful for an incident plane wave for a layered media.

Example :

int main()
{
  R2 origin, k(3.0, 1.4);
  PlaneWaveIncidentField<Real_wp, Dimension2> u_inc(origin, k);

  int i = 2; VectR2 s; mesh.GetVerticesElement(i, s);
  u_inc.InitElement(i, s);

  // then EvaluateFunctionGradient will be called with points inside element i
  Complex_wp f; R2_Complex_wp df;
  u_inc.EvaluateFunctionGradient(R2(1, 1), f, df);

  // coefficients a and b for this element ?
  Complex_wp a, b;
  u_inc.UpdateCoefAB(a, b);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

Init

This function initializes a gaussian beam.

Example :

int main()
{
  R2 origin, k(3.0, 1.4); Real_wp w = 0.2;
  // you can use the constructor
  GaussianBeamIncidentField<Real_wp, Dimension2> u_inc(origin, k, w);

  // evaluate the field
  Complex_wp f;
  u_inc.EvaluateFunction(R2(1, 1), f);

  // if you want to change parameters of the gaussian beam, use Init
  origin.Init(0.2, 0.8); k.Init(2.0, -0.9); w = 0.3;
  u_inc.Init(origin, k, w);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

Init

This function initializes an incident plane wave adapted for layered media. The incident field solves Helmholtz equation on each layer :

The position of the interfaces are located at y = d_i. The coefficients a₀ and b₀ correspond to the bottom layer y < d₀. The coefficients a_i and b_i correspond to the layer d_i-1 < y < d_i. The layer y > d_N has coefficients a_infty and b_infty.

parameters contains the coefficients in the following order :

• parameters(2) : number of intermediate layers N
• parameters(3) : coefficient b₀
• parameters(4) : coefficient a₀
• parameters(5) : position d₀
• parameters(6) : coefficient b₁
• parameters(7) : coefficient a₁
• parameters(8) : position d₁
• ...
• parameters(end) : position d_N

Example :

int main()
{
  R2 origin, k(3.0, 1.4); Real_wp omega = 0.2; Real_wp a_inf(1), b_inf(1);
  int N = 1; // one intermediate layer
  VectString parameters(3*N + 6);
  parameters(2) = "1"; parameters(3) = "2.0"; parameters(4) = "0.8"; parameters(5) = "-1.0";
  parameters(6) = "3.0"; parameters(7) = "1.4"; parameters(8) = "1.0";
  // you can use the constructor
  LayeredPlaneWaveIncidentField<Complex_wp, Dimension2> u_inc(origin, k, omega, a_inf, b_inf, param);

  // evaluate the field
  Complex_wp f;
  u_inc.EvaluateFunction(R2(1, 1), f);

  // if you want to change parameters, use Init
  param(3) = "2.5";
  u_inc.Init(origin, k, omega, a_inf, b_inf, param);
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetNbUnknowns

This function returns the number of unknowns to project. In some cases, we do not want to project all the unknowns of the equation, but only the first ones. In that case, we specify the number of unknowns that are projected (starting from the first unknown).

Example :

int main()
{
    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

  // projector for the incident field
   Vector<VectString> param(1); // only param(0) is used

   IncidentWaveProjector<Complex_wp, Dimension2>* u_inc;
   u_inc = var.GetNewIncidentProjector(VarSourceProblem_Base::INCIDENT_PLANE_WAVE, param, Complex_wp(0));

  int nb_u = u_inc->GetNbUnknowns();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

InitElement

This function is called when the projection is perfomed. The code calls InitElement for the considered element and then calls EvaluateFunction for all the needed points of the element.

Example :

class MyProjector : public VirtualProjectorFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things 
    }

   // evaluates the field at x (a point belonging to the element i)
   void EvaluateFunction(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateFunction

This function evaluates the field at a given point x.

Example :

class MyProjector : public VirtualProjectorFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things 
    }

   // evaluates the field at x (a point belonging to the element i)
   void EvaluateFunction(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

ApplyFFT

This function applies Fourier transform to the vector f in order to handle the case where the solution is decomposed into Fourier modes (cyclic or periodic domains). This method does not need to be called on regular use, it is already called in ComputeRightHandSide.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

ModifyPoints

This function modifies the points where the solution has to be evaluated in order to handle the case of periodic/cyclic domains. On regular use, this method does not need to be called.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

ModifyPoints

This function modifies the values of a function in order to handle the case of periodic/cyclic domains. On regular use, this method does not need to be called.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetOrigin

This function returns the center of the source, it will correspond for instance to the origin of the Dirac if the source is a Dirac. For a diffracted/total field, the center corresponds to the origin of the phase.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    Vector<VectString> param(1);
    param(1).Reallocate(3);
    param(0)(0) = "SRC_DIRAC"; param(0)(1) = "0.0"; param(0)(2) ="2.0"; // origin is (0, 2)

    DiracSource<HelmholtzEquation<Dimension2> > src(var, param);
    R2 x0 = src.GetOrigin(); // to retrieve the origin of the Dirac
      
}
  

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetPolarization, SetPolarization

The function GetPolarization returns the polarization of the source. The function SetPolarization sets the polarization.

Example :

int main()
{

    EllipticProblem<HarmonicElasticEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    Vector<VectString> param(1);
    param(1).Reallocate(6);
    param(0)(0) = "SRC_DIRAC"; param(0)(1) = "0.0"; param(0)(2) ="2.0"; // origin is (0, 2)
    param(0)(3) = "Polarization"; param(0)(4) = "-1.0"; param(0)(5) = "2.5";

    DiracSource<HarmonicElasticEquation<Dimension2> > src(var, param);
    R2 x0 = src.GetOrigin(); // to retrieve the origin of the Dirac
    R2_Complex_wp polar = src.GetPolarization(); // direction

    // you can also change the polarization
    polar.Init(Complex_wp(3, -2), Complex_wp(-2, 0.5));  
    src.SetPolarization(polar);
}
  

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetPolarizationGrad, SetPolarizationGrad

The function GetPolarizationGrad returns the polarization of the source (used for a gradient of a Dirac). The function SetPolarizationGrad sets this polarization.

Example :

int main()
{

    EllipticProblem<HarmonicElasticEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

    Vector<VectString> param(1);
    param(1).Reallocate(8);
    param(0)(0) = "SRC_DIRAC"; param(0)(1) = "0.0"; param(0)(2) ="2.0"; // origin is (0, 2)
    param(0)(3) = "PolarizationGrad"; param(0)(4) = "-1.0"; param(0)(5) = "2.5";
    param(0)(6) = "3.0"; param(0)(7) = "0.8";

    DiracSource<HarmonicElasticEquation<Dimension2> > src(var, param);
    R2 x0 = src.GetOrigin(); // to retrieve the origin of the Dirac
    R2_Complex_wp polar_grad = src.GetPolarizationGrad(); // direction

    // you can also change the polarization
    polar(1) = Complex_wp(3, -2);
    src.SetPolarizationGrad(polar_grad);
}
  

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsGradientDirac

This function returns true if the source is a gradient of a Dirac.

Example :

int main()
{

    EllipticProblem<HarmonicElasticEquation<Dimension2> > var;;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // presence of PolarizationGrad => gradient of a Dirac
    Vector<VectString> param(1);
    param(1).Reallocate(8);
    param(0)(0) = "SRC_DIRAC"; param(0)(1) = "0.0"; param(0)(2) ="2.0"; // origin is (0, 2)
    param(0)(3) = "PolarizationGrad"; param(0)(4) = "-1.0"; param(0)(5) = "2.5";
    param(0)(6) = "3.0"; param(0)(7) = "0.8";

    DiracSource<HarmonicElasticEquation<Dimension2> > src(var, param);
    R2 x0 = src.GetOrigin(); // to retrieve the origin of the Dirac
    R2_Complex_wp polar_grad = src.GetPolarizationGrad(); // direction

    bool grad = src.IsGradientDirac();
}
  

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

Init

This function is called in Montjoie for unsteady simulations when the source does depend on time. It provides mainly to the class defining the time-dependent source, the time for which the source needs to be evaluated.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateVolumetricSource

This function evaluates the volume source (function f_V) at a given point x.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things 
    }

   // evaluates the volume source at x (a point belonging to the element i)
   void EvaluateVolumetricSource(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsNonNullVolumetricSource

This function returns true if the volume source (function f_V) is non-null for the considered element. The argument s contains the list of vertices of the element.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  bool source_null;

  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things
      source_null = true;
      // the source will be non-null only for elements of reference 3
      if (this->var_problem.mesh.Element(i).GetReference == 3)
        source_null = false;
    }

   // to inform if the source is null or not
   bool IsNonNullVolumetricSource(const VectR2& s)
   {
     return source_null;
   }
   
   // evaluates the volume source at x (a point belonging to the element i)
   void EvaluateVolumetricSource(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateVolumetricSource

This function evaluates the volume source (function f_G) at a given point x.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things 
    }

   // evaluates the volume source at x (a point belonging to the element i)
   // this source is multiplied by gradient of basis functions
   void EvaluateGradientSource(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
     f(1) = sin(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsNonNullGradientSource

This function returns true if the volume source (function f_G) is non-null for the considered element. The argument s contains the list of vertices of the element.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  bool source_null;

  public :
    // i : element number, s : vertices of the element
    void InitElement(int i, const VectR2& s)
    {
      // if you need to precompute things
      source_null = true;
      // the source will be non-null only for elements of reference 3
      if (this->var_problem.mesh.Element(i).GetReference == 3)
        source_null = false;
    }

   // to inform if the source is null or not
   bool IsNonNullGradientSource(const VectR2& s)
   {
     return source_null;
   }
   
   // evaluates the volume source at x (a point belonging to the element i)
   // this source will be multiplied by gradient of basis functions
   void EvaluateGradientSource(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = cos(x(0))*sin(x(1));
     f(1) = sin(x(0))*sin(x(1));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsNonNullSurfacicSource

This function returns true if the surface source (function g_S) is non-null for the considered element. The argument ref contains the reference of the face.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if the surfacic source is null or not
   bool IsNonNullSurfacicSource(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the surface source at x
   void EvaluateSurfacicSource(int k, const SetPoints<Dimension2>& pts,
                               const SetMatrices<Dimension2>& mat, VectComplex_wp& f)
   {
     R2 x = pts.GetPointQuadratureBoundary(k); // to retrieve point x
     R2 n = mat.GetNormaleQuadratureBoundary(k); // outgoing normale n
     f(0) = cos(x(0))*n(0) + sin(x(1))*n(1);
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateSurfacicSource

This function evaluates the surface source (function g_S) at a given point. .

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if the surfacic source is null or not
   bool IsNonNullSurfacicSource(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the surface source at x
   void EvaluateSurfacicSource(int k, const SetPoints<Dimension2>& pts,
                               const SetMatrices<Dimension2>& mat, VectComplex_wp& f)
   {
     R2 x = pts.GetPointQuadratureBoundary(k); // to retrieve point x
     R2 n = mat.GetNormaleQuadratureBoundary(k); // outgoing normale n
     f(0) = cos(x(0))*n(0) + sin(x(1))*n(1);
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

InitSurface

This function is called before evaluating the surface source (function g_S or g_G) on quadrature points of a given face.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if the surfacic source is null or not
   bool IsNonNullSurfacicSource(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the surface source at x
   void EvaluateSurfacicSource(int k, const SetPoints<Dimension2>& pts,
                               const SetMatrices<Dimension2>& mat, VectComplex_wp& f)
   {
     R2 x = pts.GetPointQuadratureBoundary(k); // to retrieve point x
     R2 n = mat.GetNormaleQuadratureBoundary(k); // outgoing normale n
     f(0) = cos(x(0))*n(0) + sin(x(1))*n(1);
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsNonNullSurfacicSourceGradient

This function returns true if the surface source (function g_G) is non-null for the considered element. The argument ref contains the reference of the face.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if the surfacic source is null or not (with gradient of basis functions)
   bool IsNonNullSurfacicSourceGradient(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the surface source at x
   void EvaluateSurfacicSourceGradient(int k, const SetPoints<Dimension2>& pts,
                                       const SetMatrices<Dimension2>& mat, VectComplex_wp& f)
   {
     R2 x = pts.GetPointQuadratureBoundary(k); // to retrieve point x
     R2 n = mat.GetNormaleQuadratureBoundary(k); // outgoing normale n
     f(0) = cos(x(0))*n(0);
     f(1) = sin(x(1))*n(1);
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

EvaluateSurfacicSourceGradient

This function evaluates the surface source (function g_G) at a given point. .

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if the surfacic source is null or not
   bool IsNonNullSurfacicSourceGradient(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the surface source at x (integral against gradient of basis functions)
   void EvaluateSurfacicSourceGradient(int k, const SetPoints<Dimension2>& pts,
                                       const SetMatrices<Dimension2>& mat, VectComplex_wp& f)
   {
     R2 x = pts.GetPointQuadratureBoundary(k); // to retrieve point x
     R2 n = mat.GetNormaleQuadratureBoundary(k); // outgoing normale n
     f(0) = cos(x(0))*n(0);
     f(1) = sin(x(1))*n(1);
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsNonNullDirichletSource

This function returns true if the Dirichlet condition (function h) is hetereogeneous for the considered face. The argument ref contains the reference of the face.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
    // i : element number, s : vertices of the element
    void InitSurface(int i, int num_face, int num_elem, int num_loc)
    {
      // if you need to precompute things
    }

   // to inform if there is an heteregeneous Dirichlet condition
   bool IsNonNullDirichletSource(int ref)
   {
     // source only for faces of reference 3
     if (ref == 3)
       return true;

     return false;
   }
   
   // evaluates the Dirichlet function h
   void EvaluateFunction(int i, int j, const R2& x, VectComplex_wp& f)
   {
     f(0) = sin(x(1))*cos(x(0));
   }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

PresenceDirichlet

This function returns true if some hetereogeneous Dirichlet conditions (function h) may be present in the source term.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
   // no hetereogeneous Dirichlet => return false
   bool PresenceDirichlet() const { return false; }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

IsDiracSource

This function returns true if a Dirac term is present in the current source. In this case, GetOrigin returns the origin of the Dirac, and GetPolarization its direction.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
   // no Dirac => return false
   bool IsDiracSource() const { return false; }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

GetCoefficientVolume

This function returns a coefficient used to multiply the volume source f_V.

Example :

class MySource : public VirtualSourceFEM<Complex_wp, Dimension2>
{
  public :
   // multiplicative coefficient
   Complex_wp GetCoefficientVolume() const { return Complex_wp(2, 0); }

};

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

ModifyEvaluationVolume

This function is used to compute the Fourier component of the source (when the source is decomposed onto modes). On regular use, it should not be called.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

ModifyEvaluationSurface

This function is used to compute the Fourier component of the source (when the source is decomposed onto modes). On regular use, it should not be called.

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVolumeSource

This function can be used to set directly the function f_V of the right hand side. The source term is equal to

The first argument is used to set f_V only on a part of mesh. You provide the references for which the source is set.

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_V can be defined as a derived class of VirtualSourceField
    GaussianSourceField<Complex_wp, Dimension2> fV;    
    fV.Init(R2(0.1, 0.3), 1.0, 2.0);
    VectComplex_wp coef(1); coef(0) = Complex_wp(1, 0);
    fV.SetPolarization(coef);

    // SetVolumeSource to provide fV on a set of references
    IVect ref(1); ref(0) = 1; // here only for elements of reference 1
    f_obj.SetVolumeSource(ref, &fV);    

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

    // to avoid that fV is destructed twice :
    f_obj.NullifyVolumeSourceFunction();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVolumeSourceGrad

This function can be used to set directly the function f_G of the right hand side. The source term is equal to

The first argument is used to set f_G only on a part of mesh. You provide the references for which the source is set.

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_G can be defined as a derived class of VirtualSourceField
    GaussianSourceField<Complex_wp, Dimension2> fG;
    fG.Init(R2(0.1, 0.3), 1.0, 2.0);
    VectComplex_wp coef(2); coef(0) = Complex_wp(1, 0); coef(1) = Complex_wp(2, 0);
    fG.SetPolarization(coef);

    // SetVolumeSource to provide fG on a set of references
    IVect ref(1); ref(0) = 1; // here only for elements of reference 1
    f_obj.SetVolumeSourceGrad(ref, &fG);    

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

    // to avoid that fG is destructed twice :
    f_obj.NullifyVolumeSourceFunction();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVolumeSourceFunction

This function can be used to set directly the function f_V of the right hand side. The source term is equal to

The difference between this function and SetVolumeSource is that SetVolumeSource modifies f_V for only a set of references whereas SetVolumeSourceFunction modifies f_V for all the references.

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_V can be defined as a derived class of VirtualSourceField
    GaussianSourceField<Complex_wp, Dimension2> fV;    
    fV.Init(R2(0.1, 0.3), 1.0, 2.0);
    VectComplex_wp coef(1); coef(0) = Complex_wp(1, 0);
    fV.SetPolarization(coef);

    // SetVolumeSourceFunction to provide fV on all the mesh
    f_obj.SetVolumeSourceFunction(fV);    

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

    // to avoid that fV is destructed twice :
    f_obj.NullifyVolumeSourceFunction();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

NullifyVolumeSourceFunction

This function nullifies the pointers storing f_V and f_G. This method has to be called if you modified these pointers by calling SetVolumeSource, SetVolumeSourceGrad or SetVolumeSourceFunction.

Example :

int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_V can be defined as a derived class of VirtualSourceField
    GaussianSourceField<Complex_wp, Dimension2> fV;    
    fV.Init(R2(0.1, 0.3), 1.0, 2.0);
    VectComplex_wp coef(1); coef(0) = Complex_wp(1, 0);
    fV.SetPolarization(coef);

    // SetVolumeSourceFunction to provide fV on all the mesh
    f_obj.SetVolumeSourceFunction(fV);    

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

    // to avoid that fV is destructed twice :
    f_obj.NullifyVolumeSourceFunction();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetSurfaceSource

This function can be used to set directly the function g_S of the right hand side. The source term is equal to

The first argument is used to set G_S only on a part of mesh. You provide the references for which the source is set.

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function g_S can be defined as a derived class of VirtualSourceField
    GaussianSourceField<Complex_wp, Dimension2> gS;    
    gS.Init(R2(0.1, 0.3), 1.0, 2.0);
    VectComplex_wp coef(1); coef(0) = Complex_wp(1, 0);
    gS.SetPolarization(coef);

    // SetSurfaceSource to provide gS on a set of references
    IVect ref(1); ref(0) = 1; // here only for faces of reference 1
    f_obj.SetSurfaceSource(ref, &gS);    

  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

    // to avoid that gS is destructed twice :
    f_obj.NullifyVolumeSourceFunction();
}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVariableSurfaceSource

This function is used to specify the function g_S of the right hand side by giving directly the values of g_S on quadrature points. The source term is equal to

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function g_S can be defined directly on quadrature points
    IVect ref(1); ref(0) = 1; // here non-null source only for faces of reference 1
    f_obj.SetVariableSurfaceSource(ref);    
  
    // then you have to fill f_obj.evalSurf for the quadrature points of these faces
    R2 kwave(2, 1); VectR2 s;
    SetPoints<Dimension2> pts; SetMatrices<Dimension2> mat;
    for (int i = 0; i < var.mesh.GetNbBoundaryRef(); i++)
      if (var.mesh.BoundaryRef(i).GetReference() == 1)
        {
          int num_elem = var.mesh.BoundaryRef(i).numElement(0);
          const ElementReference_Dim<Dimension2>& Fb = var.GetReferenceElement(num_elem);
          var.mesh.GetVerticesElement(num_elem, s);
          Fb.FjElemQuadrature(s, pts, var.mesh, num_elem);
          Fb.DFjElemQuadrature(s, pts, mat, var.mesh, num_elem);

          int num_loc = var.mesh.Element(num_elem).FindPositionBoundary(i);
          Fb.FjSurfaceElem(s, pts, var.mesh, num_elem, num_loc);
          Fb.DFjSurfaceElem(s, pts, mat, var.mesh, num_elem, num_loc);

         int nb_pts_quad = pts.GetNbPointsQuadratureBoundary();
         f_obj.evalSurf(0)(i).Reallocate(nb_pts_quad);
         for (int j = 0; j < nb_pts_quad; j++)
           {
             R2 x = pts.GetPointQuadratureBoundary(j);
             Complex_wp vloc = exp(Iwp*DotProd(x, kwave));
             f_obj.evalSurf(0)(i)(j) = vloc;
           }
       }
    
  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
  DISP(rhs(0));

}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVariableSource

This function is used to specify the function f_V of the right hand side by giving directly the values of f_V on quadrature points. The source term is equal to

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_V can be defined directly on quadrature points
    IVect ref(1); ref(0) = 1; // here non-null source only for elements of reference 1
    f_obj.SetVariableSource(ref);    
  
    // then you have to fill f_obj.evalS for the quadrature points of these elements
    R2 kwave(2, 1); VectR2 s;
    SetPoints<Dimension2> pts; SetMatrices<Dimension2> mat;
    for (int num_elem = 0; num_elem < var.mesh.GetNbElt(); num_elem++)
      if (var.mesh.Element(num_elem).GetReference() == 1)
        {
          const ElementReference_Dim<Dimension2>& Fb = var.GetReferenceElement(num_elem);
          var.mesh.GetVerticesElement(num_elem, s);
          Fb.FjElemQuadrature(s, pts, var.mesh, num_elem);
          Fb.DFjElemQuadrature(s, pts, mat, var.mesh, num_elem);

          int nb_pts_quad = Fb.GetNbPointsQuadratureInside();
          f_obj.evalS(0)(num_elem).Reallocate(nb_pts_quad);
          for (int j = 0; j < nb_pts_quad; j++)
           {
             R2 x = pts.GetPointQuadrature(j);
             Complex_wp vloc = exp(Iwp*DotProd(x, kwave));
             f_obj.evalS(0)(num_elem)(j) = vloc;
           }
       }
    
  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

SetVariableGradientSource

This function is used to specify the function f_G of the right hand side by giving directly the values of f_G on quadrature points. The source term is equal to

Example :

    int main()
{

    EllipticProblem<HelmholtzEquation<Dimension2> > var;

    // constructing the problem (mesh, finite element)
    All_LinearSolver* solver;
    var.ConstructAll("test.ini", "QUADRANGLE_LOBATTO", "HELMHOLTZ", solver);

   // vector containing the right hand sides
   Vector<VectComplex_wp> rhs(1); // only one vector here
   rhs(0).Reallocate(var.GetNbDof()); rhs(0).Zero();

   // generic source => VolumetricSource
    VolumetricSource<HelmholtzEquation<Dimension2> >  f_obj(var); // instance of the class defining the source
    Vector<VirtualSourceFEM<Complex_wp, Dimension2>* > vec_f(1);
   vec_f(0) = &f_obj;

    // then function f_G can be defined directly on quadrature points
    IVect ref(1); ref(0) = 1; // here non-null source only for elements of reference 1
    f_obj.SetVariableGradientSource(ref);    
  
    // then you have to fill f_obj.evalG for the quadrature points of these elements
    R2 kwave(2, 1); VectR2 s;
    SetPoints<Dimension2> pts; SetMatrices<Dimension2> mat;
    for (int num_elem = 0; num_elem < var.mesh.GetNbElt(); num_elem++)
      if (var.mesh.Element(num_elem).GetReference() == 1)
        {
          const ElementReference_Dim<Dimension2>& Fb = var.GetReferenceElement(num_elem);
          var.mesh.GetVerticesElement(num_elem, s);
          Fb.FjElemQuadrature(s, pts, var.mesh, num_elem);
          Fb.DFjElemQuadrature(s, pts, mat, var.mesh, num_elem);

          int nb_pts_quad = Fb.GetNbPointsQuadratureInside();
          f_obj.evalG(0)(num_elem).Reallocate(nb_pts_quad);
          for (int j = 0; j < nb_pts_quad; j++)
           {
             R2 x = pts.GetPointQuadrature(j);
             Complex_wp vloc_x = exp(Iwp*DotProd(x, kwave));
             Complex_wp vloc_y = vloc_x*kwave(1); vloc_x *= kwave(0);
             f_obj.evalG(0)(num_elem)(j) = vloc_x;
             f_obj.evalG(1)(num_elem)(j) = vloc_y;
           }
       }
    
  // then we compute the source
  var.ComputeGenericSource(rhs, vec_f);
    DISP(rhs(0));

}

Location :

Source/SourceSpatiale.hxx
Source/SourceSpatiale.cxx

Evaluate

This function returns the function h(t)

Example :

    int main()
{
    VirtualTimeSource<Real_wp>* fsrc;

    // you choose the time source you want
    fsrc = new TimeRickerSource(1.0);

    // then you can call the virtual method Evaluate
    Real_wp t = 2.5;
    Real_wp h = fsrc->Evaluate(t);

}

Location :

Source/UserSource.hxx
Source/UserSource.cxx

EvaluateDerivative

This function returns h'(t)

Example :

    int main()
{
    VirtualTimeSource<Real_wp>* fsrc;

    // you choose the time source you want
    fsrc = new TimeRickerSource(1.0);

    // then you can call the virtual method Evaluate
    Real_wp t = 2.5;
    Real_wp h = fsrc->Evaluate(t);

    // and derivative :
    Real_wp hp = fsrc->EvaluateDerivative(t);

    delete fsrc;
}

Location :

Source/UserSource.hxx
Source/UserSource.cxx

Init

By calling this method, you provide the function h(t) to the class.

Example :

    int main()
{
    VirtualTimeSource<Real_wp>* fsrc;

    // you choose the time source you want
    Real_wp freq = 1.0;
    fsrc = new TimeRickerSource(freq);

    // you encapsulate this source in TimeSourceHyperbolic
    TimeSourceHyperbolic var;
    var.Init(fsrc, freq, 1e-12, 0.0);

    // and you can evaluate the second derivative of h for instance
    Real_wp hpp; Real_wp t = 0.2;
    var.EvaluateDerivative(t, 2, hpp);

    // no need to delete fsrc, it will be deleted by the destructor of TimeSourceHyperbolic
}

Location :

Source/UserSource.hxx
Source/UserSource.cxx

EvaluateDerivative

This method computes the n-th derivative of the function at time t.

Example :

    int main()
{
    VirtualTimeSource<Real_wp>* fsrc;

    // you choose the time source you want
    Real_wp freq = 1.0;
    fsrc = new TimeRickerSource(freq);

    // you encapsulate this source in TimeSourceHyperbolic
    TimeSourceHyperbolic var;
    var.Init(fsrc, freq, 1e-12, 0.0);

    // and you can evaluate the second derivative of h for instance
    Real_wp hpp; Real_wp t = 0.2;
    var.EvaluateDerivative(t, 2, hpp);

    // no need to delete fsrc, it will be deleted by the destructor of TimeSourceHyperbolic
}

Location :

Source/UserSource.hxx
Source/UserSource.cxx