• Introduction
• Algebra
• Quadrature
  • Singular integrals
• Common functions
• Mesh
• Finite Elements
• Computations
• Boundary conditions
• Post-processing
• Solved equations
• Time schemes

Singular integrals in Montjoie

A basic class concerns the integration of singular integrals : SingularDoubleQuadratureGalerkin_Base.

Methods for SingularDoubleQuadratureGalerkin_Base :

SetOperator	Selects the operator to integrate
SetAnalyticalIntegration	chooses a semi-analytical integration
SetNumericalIntegration	chooses a numerical integration
SetFiniteElement	sets the finite element used over the 2-D curve
SetTriangularFiniteElement	sets the triangular finite element used over the 3-D surface
SetQuadrilateralFiniteElement	sets the quadrilateral finite element used over the 3-D surface
ComputeElemMatrix	computes the elementary matrix for two panels
ComputeMassMatrix	computes the elementary mass matrix (single integral)
ComputeStiffnessMatrix	computes the elementary stiffness matrix (single integral)

Methods for EuclidianDistanceClass_Base :

GetDistance	returns the distance between two points
GetDiff	subtracts two points

SetOperator

Syntax

  void SetOperator(int type, Real_wp beta) const;

This method sets which operator to use and the associated constant β. The choices are the following:

• OPERATOR_DS : this operator is associated with the following elementary matrix:
  
• OPERATOR_GRAD : this operator is associated with the following elementary matrix:
  
• OPERATOR_DIFF : this operator is associated with the following elementary matrix:
  

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the two edges I and J : 
VectR2 PtsAi(2), PtsAj(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);
PtsAj(0).Init(0.4, 0.6); PtsAj(1).Init(0.8, 1.0);

// if the edge is curved, curved = true
bool curved = false;

// nodal points on the two edges and associated normales
VectR2 PointsNodalI(r+1), PointsNodalJ(r+1), NormaleNodalI(r+1), NormaleNodalJ(r+1);

// quadrature points on the two edges and associated normales
VectR2 PointsQuadI(r+1), PointsQuadJ(r+1), NormaleQuadI(r+1), NormaleQuadJ(r+1);
// lineic element (ds) on each quadrature point:
VectReal_wp Dsi, Dsj;

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobNodalI, MatJacobNodalJ, MatJacobI, MatJacobJ

// class describing the distance || X - Y || 
EuclidianDistanceClass_Base<Dimension2> var_dist;

// elementary matrices must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1), mat_elemI(r+1, r+1), mat_elemJ(r+1, r+1);

// finally elementary matrices are computed
sing_int.ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj, PointsNodalI, PointsNodalJ,
                  PointsQuadI, PointsQuadJ, NormaleNodalI, NormaleNodalJ, 
                  NormaleQuadI, NormaleQuadJ, MatJacobNodalI, MatJacobNodalJ,
                  MatJacobI, MatJacobJ, var_dist, mat_elem, mat_elemI, mat_elemJ);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

SetAnalyticalIntegration

Syntax

  void SetAnalyticalIntegration(Real_wp eta, int rp, int rj) const;

This method informs the object to use a semi-analytical method to evaluate singular integral. eta is the parameter used to determine if two edges Ei and Ej are close or far. They are assumed to be close if the following criterion is satisfied :

This parameter eta is only used in 2-D. For 3-D faces, the faces are assumed to be far if they do not share a common vertex. rp is the order of integration to use for close panels (hence only meaningful in 2-D), and rj the order of integration to use for joint panels (panels sharing points but not identical).

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5;
sing_int.SetAnalyticalIntegration(eta, rp, rj);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the two edges I and J : 
VectR2 PtsAi(2), PtsAj(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);
PtsAj(0).Init(0.4, 0.6); PtsAj(1).Init(0.8, 1.0);

// if the edge is curved, curved = true
bool curved = false;

// nodal points on the two edges and associated normales
VectR2 PointsNodalI(r+1), PointsNodalJ(r+1), NormaleNodalI(r+1), NormaleNodalJ(r+1);

// quadrature points on the two edges and associated normales
VectR2 PointsQuadI(r+1), PointsQuadJ(r+1), NormaleQuadI(r+1), NormaleQuadJ(r+1);
// lineic element (ds) on each quadrature point:
VectReal_wp Dsi, Dsj;

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobNodalI, MatJacobNodalJ, MatJacobI, MatJacobJ

// class describing the distance || X - Y || 
EuclidianDistanceClass_Base<Dimension2> var_dist;

// elementary matrices must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1), mat_elemI(r+1, r+1), mat_elemJ(r+1, r+1);

// finally elementary matrices are computed
sing_int.ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj, PointsNodalI, PointsNodalJ,
                  PointsQuadI, PointsQuadJ, NormaleNodalI, NormaleNodalJ, 
                  NormaleQuadI, NormaleQuadJ, MatJacobNodalI, MatJacobNodalJ,
                  MatJacobI, MatJacobJ, var_dist, mat_elem, mat_elemI, mat_elemJ);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

SetNumericalIntegration

Syntax

  void SetNumericalIntegration(Real_wp eta, int rp, int rj, int rj_ext, int n_int, int n_ext) const;

This method sets a numerical integration to evaluate the elementary matrix. eta is the parameter used to determine if two edges Ei and Ej are close or far. They are assumed to be close if the following criterion is satisfied :

This parameter eta is only used in 2-D. For 3-D faces, the faces are assumed to be far if they do not share a common vertex. rp is the order of integration to use for close panels (hence only meaningful in 2-D), and rj the order of integration to use for joint panels (panels sharing points but not identical) and intern integral. rj_ext is the order of integration to use for joint panels and extern integral. n_int and n_ext are the orders of integration for identical panels for the intern and extern integral.

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the two edges I and J : 
VectR2 PtsAi(2), PtsAj(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);
PtsAj(0).Init(0.4, 0.6); PtsAj(1).Init(0.8, 1.0);

// if the edge is curved, curved = true
bool curved = false;

// nodal points on the two edges and associated normales
VectR2 PointsNodalI(r+1), PointsNodalJ(r+1), NormaleNodalI(r+1), NormaleNodalJ(r+1);

// quadrature points on the two edges and associated normales
VectR2 PointsQuadI(r+1), PointsQuadJ(r+1), NormaleQuadI(r+1), NormaleQuadJ(r+1);
// lineic element (ds) on each quadrature point:
VectReal_wp Dsi, Dsj;

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobNodalI, MatJacobNodalJ, MatJacobI, MatJacobJ

// class describing the distance || X - Y || 
EuclidianDistanceClass_Base<Dimension2> var_dist;

// elementary matrices must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1), mat_elemI(r+1, r+1), mat_elemJ(r+1, r+1);

// finally elementary matrices are computed
sing_int.ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj, PointsNodalI, PointsNodalJ,
                  PointsQuadI, PointsQuadJ, NormaleNodalI, NormaleNodalJ, 
                  NormaleQuadI, NormaleQuadJ, MatJacobNodalI, MatJacobNodalJ,
                  MatJacobI, MatJacobJ, var_dist, mat_elem, mat_elemI, mat_elemJ);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

SetFiniteElement, SetTriangularFiniteElement, SetQuadrilateralFiniteElement

Syntax

  void SetFiniteElement(const ElementReference<Dimension1, 1>& elt);

  void SetTriangularFiniteElement(const ElementReference<Dimension2, 1>& elt, const VectR2& points, const VectReal_wp& weights);
  void SetQuadrilateralFiniteElement(const ElementReference<Dimension2, 1>& elt, const VectR2& points, const VectReal_wp& weights);

The method SetFiniteElement is available in 2-D only and provides basis functions on the unit interval [0, 1]. These basis functions are used when the elementary matrices are computed. The methods SetTriangularFiniteElement/SetQuadrilateralFiniteElement are available in 3-D only and provide basis functions on the unit triangle and the unit square, they take as argument the integration points and weights to use for far panels.

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the two edges I and J : 
VectR2 PtsAi(2), PtsAj(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);
PtsAj(0).Init(0.4, 0.6); PtsAj(1).Init(0.8, 1.0);

// if the edge is curved, curved = true
bool curved = false;

// nodal points on the two edges and associated normales
VectR2 PointsNodalI(r+1), PointsNodalJ(r+1), NormaleNodalI(r+1), NormaleNodalJ(r+1);

// quadrature points on the two edges and associated normales
VectR2 PointsQuadI(r+1), PointsQuadJ(r+1), NormaleQuadI(r+1), NormaleQuadJ(r+1);
// lineic element (ds) on each quadrature point:
VectReal_wp Dsi, Dsj;

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobNodalI, MatJacobNodalJ, MatJacobI, MatJacobJ;

// class describing the distance || X - Y || 
EuclidianDistanceClass_Base<Dimension2> var_dist;

// elementary matrices must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1), mat_elemI(r+1, r+1), mat_elemJ(r+1, r+1);

// finally elementary matrices are computed
sing_int.ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj, PointsNodalI, PointsNodalJ,
                  PointsQuadI, PointsQuadJ, NormaleNodalI, NormaleNodalJ, 
                  NormaleQuadI, NormaleQuadJ, MatJacobNodalI, MatJacobNodalJ,
                  MatJacobI, MatJacobJ, var_dist, mat_elem, mat_elemI, mat_elemJ);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

ComputeElemMatrix

Syntax

  void ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj,
                         PointsNodalI, PointsNodalJ, PointsQuadI, PointsQuadJ,
                         NormaleNodalI, NormaleNodalJ, NormaleQuadI, NormaleQuadJ,
                         MatJacobNodalI, MatJacobNodalJ, MatJacobQuadI, MatJacobQuadJ,
                         var_dist, mat_elem, mat_elemI, mat_elemJ);

This method computes elementary matrices corresponding to the following integrals (it depends on the operator specified when calling SetOperator:

• OPERATOR_DS : this operator is associated with the following elementary matrix:
  
• OPERATOR_GRAD : this operator is associated with the following elementary matrix:
  
• OPERATOR_DIFF : this operator is associated with the following elementary matrix:
  

This method takes the following arguments:

• PtsAi, PtsAj: Vertices of the two panels E_i and E_j.
• curved : true if one (or both) of these two panels is curved
• Dsi, Dsj: surface elements on quadrature points.
• PointsNodalI, PointsNodalJ, PointsQuadI and PointsQuadJ : nodal points and quadrature points of these two panels (corresponding to far panels).
• NormaleNodalI, NormaleNodalJ, NormaleQuadI, NormaleQuadJ : associated normales
• MatJacobNodalI, MatJacobNodalJ, MatJacobQuadI, MatJacobQuadJ : jacobian matrices DF_i^-1 which are used to transform a 2-D gradient on the reference element to a 3-D surface gradient
• var_dist : object providing the euclidian norm || X- Y || (of type EuclidianDistanceClass_Base)
• mat_elem : elementary matrix between panel I and panel J
• mat_elemI : elementary matrix of panel I
• mat_elemJ : elementary matrix of panel J

mat_elemI and mat_elemJ are filled only for OPERATOR_DIFF.

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the two edges I and J : 
VectR2 PtsAi(2), PtsAj(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);
PtsAj(0).Init(0.4, 0.6); PtsAj(1).Init(0.8, 1.0);

// if the edge is curved, curved = true
bool curved = false;

// nodal points on the two edges and associated normales
VectR2 PointsNodalI(r+1), PointsNodalJ(r+1), NormaleNodalI(r+1), NormaleNodalJ(r+1);

// quadrature points on the two edges and associated normales
VectR2 PointsQuadI(r+1), PointsQuadJ(r+1), NormaleQuadI(r+1), NormaleQuadJ(r+1);
// lineic element (ds) on each quadrature point:
VectReal_wp Dsi, Dsj;

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobNodalI, MatJacobNodalJ, MatJacobI, MatJacobJ

// class describing the distance || X - Y || 
EuclidianDistanceClass_Base<Dimension2> var_dist;

// elementary matrices must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1), mat_elemI(r+1, r+1), mat_elemJ(r+1, r+1);

// finally elementary matrices are computed
sing_int.ComputeElemMatrix(PtsAi, PtsAj, curved, Dsi, Dsj, PointsNodalI, PointsNodalJ,
                  PointsQuadI, PointsQuadJ, NormaleNodalI, NormaleNodalJ, 
                  NormaleQuadI, NormaleQuadJ, MatJacobNodalI, MatJacobNodalJ,
                  MatJacobI, MatJacobJ, var_dist, mat_elem, mat_elemI, mat_elemJ);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

ComputeMassMatrix

Syntax

  void ComputeMassMatrix(PtsAi, curved, Dsi, mat_elem);

This method computes the elementary mass matrix corresponding to the following integral :

This method fills the elementary matrix mat_elem and asks the following input arguments:

• PtsAi : Vertices of the panel E_i .
• curved : true if the panel E_i is curved
• Dsi: surface elements on quadrature points.

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the edge I
VectR2 PtsAi(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);

// if the edge is curved, curved = true
bool curved = false;

// lineic element (ds) on each quadrature point:
VectReal_wp Dsi;

// mass elementary matrix must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1);

// and computed
sing_int.ComputeMassMatrix(PtsAi, curved, Dsi, mat_elem);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

ComputeStiffnessMatrix

Syntax

  void ComputeStiffnessMatrix(PtsAi, curved, Dsi, NormaleQuadI, MatJacobQuadI, kinf, epsilon, mat_elem);

This method computes the elementary stiffness matrix corresponding to the following integral :

This method fills the elementary matrix mat_elem and asks the following input arguments:

• PtsAi : Vertices of the panel E_i .
• curved : true if the panel E_i is curved
• Dsi : surface elements on quadrature points.
• NormaleQuadI : associated normales.
• MatJacobQuadI : jacobian matrices DF_i^-1 which are used to transform a 2-D gradient on the reference element to a 3-D surface gradient
• kinf : k_ε for ε = 0
• mat_elem : elementary stiffness matrix

Example :

// for 1-D integrals over 2-D curves :
SingularDoubleQuadratureGalerkin_Base<Real_wp, Dimension1> sing_int;

// you select which operator you want to compute with the exponent beta
sing_int.SetOperator(sing_int.OPERATOR_DS, Real_wp(0.5));

// order of approximation
int r = 4;

// which quadrature strategy to use
Real_wp eta = 2.0;
int rp = r+2, rj = r+5, rj_ext = 2*r+6, n_int = 2*r, n_ext = r+4;
sing_int.SetNumericalIntegration(eta, rp, rj, rj_ext, n_int, n_ext);

// 1-D finite element describing basis functions
EdgeLobattoReference elt;
elt.ConstructFiniteElement(r);
sing_int.SetFiniteElement(elt);

// points of the edge I
VectR2 PtsAi(2);
PtsAi(0).Init(-0.5, 0.2); PtsAi(1).Init(-0.2, 0.3);

// if the edge is curved, curved = true
bool curved = false;

// normale on quadrature points
VectR2 NormaleQuadI(r+1);

// matrices DF_i^{-1} (not used in 2-D) on nodal points and quadrature points
Vector<TinyMatrix<Real_wp, General, 3, 2> > MatJacobQuad;

// lineic element (ds) on each quadrature point:
VectReal_wp Dsi;

// mass elementary matrix must be allocated
Matrix<Real_wp> mat_elem(r+1, r+1);

// wave number for epsilon = 0
Real_wp kinf = 2.0*pi_wp;

// damping parameter epsilon (k_epsilon = k_inf + i epsilon for
epsilon > 0)
// (k_epsilon = k_inf + 0.4 i k_inf^{1/3} \Kappa^{2/3} for epsilon < 0)
Real_wp epsilon  = -1.0;

// elementary matrix is computed
sing_int.ComputeStiffnessMatrix(PtsAi, curved, Dsi, NormaleQuad, MatJacobQuad, kinf, epsilon, mat_elem);

Location :

SingularIntegration2D.cxx
SingularIntegration3D.cxx

GetDistance

Syntax

  Real_wp GetDistance(R_N x, R_N y);

This method computes the distance between two points X and Y.

Example :

EuclidiantDistanceClass_Base<Dimension2> var_dist;

R2 X, Y;
// distance between X and Y :
Real_wp dist = var_dist.GetDistance(X, Y);

Location :

SingularIntegration2D.cxx

GetDiff

Syntax

  R_N GetDiff(R_N x, R_N y);

This method computes the difference between points X and Y.

Example :

EuclidiantDistanceClass_Base<Dimension2> var_dist;

R2 X, Y;
// X-Y
R2 diff = var_dist.GetDiff(X, Y);

Location :

SingularIntegration2D.cxx