Finite element computations in Montjoie

Computation of finite element matrix

The general variational formulation for continuous elements is equal to (see the section devoted to the description of equations if you want more details about that formulation)

The tensors A, B, C, D and E are provided by the class defining the solved equation. These tensors will be of the form :

where

The variational formulation leads to the following linear system

where the finite element matrix Ah is equal to

The matrices Mh, Sh and Kh are respectively the mass matrix, damping matrix and stiffness matrix. The finite element matrix is computed by calling method AddMatrixWithBC :

// The definition of the problem is constructed via EllipticProblem class
EllipticProblem<LaplaceEquation<Dimension2> > var;
var.InitIndices(100);
var.SetTypeEquation("LAPLACE"); // name of the equation, it can be used to use an equivalent formulation of the same equation

var.ComputeMeshAndFiniteElement("QUADRANGLE_LOBATTO"); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can call AddMatrixWithBC
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
Matrix<Real_wp, Symmetric, ArrayRowSymSparse> A;

// but you can change them
Real_wp alpha = 2.0, beta = 0.5, gamma = 0.25;
nat_mat.SetCoefMass(alpha);
nat_mat.SetCoefDamping(beta);
nat_mat.SetCoefStiffness(gamma);
// matrix is added, so you need to clear it if you do not want to keep
// previous non-zero entries
A.Clear();

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Real_wp(0));
delete Ah;


The method AddMatrixWithBC calls method AddBoundaryConditionsTerms for terms due to boundary conditions and other extra terms then, the method AddMatrixFEM is called. In this last method, the boundary integrals (appearing only in discontinuous Galerkin formulation) are added to the matrix by the function AddElementaryFluxesDG, then the volume integrals are added by the function AssembleMatrix (which calls the function ComputeElementaryMatrix for each element of the mesh).

Functions related to computation of finite element matrix

 AssembleMatrix assembles elementary matrices in order to form the finite element matrix ComputeElementaryMatrix computes elementary matrix for continuous or discontinuous elements AddElementaryFluxesDG adds boundary integrals that appear in DG formulation to a given matrix

Public attributes of class VarProblem_Base

The class VarProblem_Base is the base class for steady (or time-harmonic) problems (class EllipticProblem). This class is abstract and non-template. It does not depend on the dimension or the type of equation. Below, we detail the public attributes of this class. Methods are listed in the next table.

 print_level verbosity level nb_unknowns_scal number of scalar unknowns nb_unknowns_vec number of vectorial unknowns nb_unknowns number of unknowns nb_components_en number of components for the trace of the solution nb_components_hn number of components for the trace of the gradient of solution nb_unknowns_hdg number of surface unknowns (for HDG formulation) type_element type of the finite element used for the main unknown other_type_element type of the finite elements used for other unknowns finite_element_name name of the finite element used for the main unknown name_other_elements name of the finite elements used for other unknowns mesh_num_unknown list of mesh numberings to use for each unknown offset_dof_unknown incremental number of degrees of freedom for each unknown offset_dof_condensed incremental number of degrees of freedom for each unknown (after static condensation) dg_formulation Formulation used (continuous, discontinous or HDG) sipg_formulation true for Interior Penalty Discontinuous Galerkin compute_dfjm1 true if the jacobian matrices DFi need to be computed and stored alpha_penalization penalty parameter for discontinuous Galerkin formulations delta_penalization penalty parameter for discontinuous Galerkin formulations upwind_fluxes true if upwind fluxes will be used for a discontinuous Galerkin formulation automatic_choice_penalization true if penalty parameters are set automatically Glob_CoefPenalDG Penalty parameters for each face (IPDG) mesh_data Parameters for constructing the mesh exit_if_no_boundary_condition if true, the simulation is stopped if there is a boundary face without boundary condition var_chrono Object used to compute timings

Methods of class VarComputationProblem_Base

The class VarComputationProblem_Base is an abstract class that is used for the assembly of matrices (see details in the description of AssembleMatrix). Below we list the methods of this class, most of them are virtual and are overloaded in the derived classes.

 GetThresholdMatrix returns the threshold used to drop negligible entries of the matrix SetThresholdMatrix sets the threshold used to drop negligible entries of the matrix GetNbElt returns the number of elements GetNbRows returns the number of rows of the assembled matrix ComputeElementaryMatrix computes the elementary matrix of a given element GetInternalNodesElement retrieves degrees of freedom that can be eliminated for a given element GetNewCondensationSolver constructs a new object handling static condensation

Methods of class VarComputationProblem (class inherited from VarComputationProblem_Base)

The class VarComputationProblem implements the computation of finite element matrices.

Public attributes of class VarGeometryProblem (class inherited from VarProblem_Base)

 mesh mesh used for the simulations Glob_PointsQuadrature quadrature points for all the elements (if stored) write_quadrature_points if true, quadrature points are written in a file write_quad_points_pml if true, quadrature points (of PMLs) are written in a file Glob_jacobian Determinant of jacobian matrices on quadrature points Glob_decomp_jacobian Decomposition of the determinant of jacobian matrices on a polynomial basis Glob_normale Outgoing normales on quadrature points of the boundaries of the mesh Glob_dsj Surface integration elements on quadrature points of the boundaries of the mesh OrthogonalElement type of orthogonality for each element Glob_DFjm1 inverse of jacobian matrices (multiplied by the determinant) IsNewFace true if the face j of element i has not been considered

Methods of class VarGeometryProblem (class inherited from VarProblem_Base)

 ConstructFiniteElement constructs finite element classes ConstructExactIntegrationFiniteElement constructs finite element classes with exact integration CopyFiniteElement copies finite element classes from another object EllipticProblem ClearFiniteElement clears finite element classes ComputeNumberOfDofs computes the number of degrees of freedom for only the main unknowns CheckContinuity checks continuity of basis functions (or continuity of tangential traces for edge elements) UseNumericalIntegration if true, quadrature points are needed to compute integrals because of varying coefficients FaceHasToBeConsideredForBoundaryIntegral if true the i-th is involved in a boundary integral of the variational formulation GetDofNumberOnElement retrieves dof numbers of an element of the mesh GetLocalUnknownVector extracts the local components of u on an element of the mesh AddLocalUnknownVector adds a local (with only values on degrees of freedom of an element of the mesh) vector to the global vector ModifyLocalComponentVector intermediary function used by GetLocalUnknownVector ModifyLocalUnknownVector intermediary function used by AddLocalUnknownVector GetGlobalUnknownVector modification of local projection on reference element to obtain a global projection (respecting signs and orientations of global dofs) GetFaceBasis returns finite element class of order r for unit triangle/tetrahedron GetFaceBasis2 returns finite element class of order r for unit square/pyramid GetFaceBasis3 returns finite element class of order r for unit prism GetFaceBasis4 returns finite element class of order r for unit cube HybridMesh returns true if the mesh is hybrid ComputeValueNodalBoundary computes the values of shape functions on the boundary of an element ComputeGaussIntegralSurface computes integral on a face of the mesh against basis functions with exact integration ComputeValuesPhi computes the values of basis functions of an element on a local point of this element ComputeValuesGradientPhi computes the gradients of basis functions of an element on a local point of this element

Public attributes of class DistributedProblem

 comm_group_mode MPI communicator for nodes sharing the computation of the same mode comm_group_master MPI communicator with master nodes associated with the computation of different modes splitting_algorithm returns the splitting algorithm used to part the mesh (Metis or Scotch) OverlapDofNumber_Subdomain array storing the dof numbers which are already handled by another processor of lower rank OverlapProcNumber_Subdomain array storing the processor ranks associated with dofs which are already handled by another processor of lower rank OverlapNumber_Subdomain dof numbers already handled by another processor PeriodicDofComm dof numbers for periodic dofs on the interface between two processors OriginalPeriodicDofComm original dof numbers for periodic dofs on the interface between two processors nodl_mesh total number of dofs associated with the mesh nodl_all total number of dofs (size of the global matrix)

Public attributes of class VarHarmonic

 Glob_jacobian jacobian for affine elements and weighted jacobian for non-affine elements Glob_decomp_jacobian decomposition of jacobian on a polynomial basis for straight elements ElementRho local number of elements among elements of the same reference Glob_normale unit normale on quadrature points of the faces Glob_dsj element of surfacic integration on quadrature points of the faces Glob_CoefPenalDG coefficients for penalty terms in DG formulation Glob_DFjm1 Ji DF-1 on quadrature points of each element IsNewFace if IsNewFace(i)(j) is true the face j of the element i has not been seen before storage_matrix_asked if true, the finite element matrix will be stored in a file file_name_matrix_stored name of the file where the finite element matrix may be stored threshold_matrix entries of the matrix below that value are discarded static_condensation if true, the internal dofs are eliminated through a static condensation process var_transmission the object containing arrays needed to compute terms due to transmission conditions

TO DO

Methods of class VarHarmonic

 GetWaveVector returns the wave vector SetWaveVector sets the wave vector GetPolarization returns the polarization vector GetPhaseOrigin sets the origin of the phase (for plane waves) IsComplexProblem returns true if the finite element matrix will contain complex non-zero entries StaticCondensation returns true if static condensation is performed to eliminate interior dofs PutAdditionalDofsAndOtherInitializations method to handle additional degrees of freedom and other initializations for various models RunAll runs a complete simulation (from reading the input file until writing the results in output files) AddMatrixFEM adds volume integrals in finite element matrix AddMatrixWithBC adds all terms for the computation of finite element matrix ComputeDiagonalMatrix computes only the diagonal of the finite element matrix GetProfilBlockDiagonal returns the sparsity pattern for a block-diagonal approximation of the finite element matrix ComputeBlockDiagonalMatrix computes only a block-diagonal submatrix of the finite element matrix ComputeMassMatrix computes geometric quantities Ji, DFi, normale on quadrature points, etc needed for the computation of the finite element matrix AllocateMassMatrices allocation for arrays storing geometric quantities needed for the computation of the finite element matrix ComputePhysicalCoefficients computes the physical coefficients (rho, epsilon, mu, etc) on quadrature points of the mesh FinalizeComputationVaryingIndices method called after the computation of physical coefficients GetInverseSquareRootMassMatrix computes the inverse of the square root of mass matrix (if diagonal) GetMassMatrix computes the mass matrix (if diagonal) ClearMassMatrix clears arrays containing geometric quantities Ji, DFi, etc

Matrix-vector product with finite element matrix

The finite element matrix is represented by the class FemMatrixFreeClass. This class may store the matrix or not, depending on the finite element, the order of approximation, and the solved equation. The method MltVector is overloaded for this class.

// The definition of the problem is constructed via EllipticProblem class
EllipticProblem<TypeEquation> var;
var.InitIndices(100);
var.SetTypeEquation("HELMHOLTZ");

var.PerformOtherInitializations();
var.ComputeMassMatrix();
var.ComputeQuasiPeriodicPhase();

// once var is constructed, you can call AddMatrixWithBC
GlobalGenericMatrix<Complex_wp> nat_mat;
Matrix<Complex_wp, Symmetric, ArrayRowSymSparse> A;

// By default, alpha = beta = gamma = 1
// but you can change them
Complex_wp alpha = 2.0, beta = 0.5, gamma = 0.25;
nat_mat.SetCoefMass(alpha);
nat_mat.SetCoefDamping(beta);
nat_mat.SetCoefStiffness(gamma);
// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass&_Baselt;Real_wp>* Ah = var.GetNewIterativeMatrix(Real_wp(0));

// then you can compute the matrix-vector product
Vector<Complex_wp> x(Ah->GetM()), y(Ah->GetM()); x.FillRand();
Ah->MltVector(x, y);

delete Ah;


The function MltAddFree is overloaded in the leaf classes (depending on the equation).

Public attributes of FemMatrixFreeClass

 mat_boundary* additional sparse matrix coming from boundary conditions or other models mat_iterative* main sparse matrix if the finite element is stored matCSR_boundary* mat_boundary in CSR (Compressed Sparse Row) format matCSR_iterative* mat_iterative in CSR (Compressed Sparse Row) format var problem associated with the finite element matrix

Methods of FemMatrixFreeClass (inherited from VirtualMatrix)

 constructor of FemMatrixFreeClass SetCoefficientDirichlet changes coefficient on diagonal entry for Dirichlet dofs SetCoefficientMatrix changes mass, stifness and damping coefficients GetCoefMass returns the coefficient associated with the mass matrix IsSymmetric returns true if the current matrix is symmetric FormulationDG returns the type of discontinuous Galerkin formulation SetCondensedSolver sets the condensed solver used to compute the condensed matrix DirichletDofIgnored returns true if Dirichlet dofs will be ignored in the matrix vector product IgnoreDirichletDof informs that Dirichlet dofs should be ignored SetScaling sets scalings to be used for the rows and columns of the matrix SucceedInAffectingPointer returns true if the method succeeded in addressing the pointer to the current matrix InitSymmetricMatrix initializes the matrix as a symmetric matrix InitUnsymmetricMatrix initializes the matrix as a non-symmetric matrix ApplyLeftScaling applies row scaling on a vector ApplyRightScaling applies column scaling on a vector CompressMatrix converts sparse matrices to CSR matrices to reduce memory usage AddExtraBoundaryTerms adds terms due to boundary conditions SetNbDirichletCondition sets the number of right hand sides (for Dirichlet conditions) ApplyDirichletCondition modifies the right hand side such that it vanishes for Dirichlet dofs MltAddHetereogeneousDirichlet multiplies the matrix with only columns associated with Dirichlet dofs SetDirichletCondition modifies the matrix due to Dirichlet condition InitDirichletCondition initializes Dirichlet condition ImposeDirichletCondition imposes a null Dirichlet condition to the vector given on input MltAddFree performs the matrix-vector product (matrix-free implementation) GetExtrapolVariables returns the intermediary object used to perform the matrix vector product

Methods of MatrixVectorProductLevel

 SetLevel specifies which elements are to be considered for the matrix-vector product GetLevelArray returns the element numbers for each level SetLevelArray sets the element numbers for each level GetNbElt returns the number of elements for a given level GetElementNumber returns the element number of the selected level GetLocalElementNumber returns the local element number of the selected level TreatElement returns true if the element i should be considered in the matrix-vector product SetNbElt sets the number of elements in the mesh GetMemorySize returns the memory used by the object in bytes

Methods of CondensationBlockSolver_Base

 ModifyElementaryMatrix applies the static condensation to an elementary matrix SetElementNumber sets the element number (global and condensed) GetCondensedElementNumber returns the condensed number of the selected element GetGlobalElementNumber returns the global number of the selected element SetNbCondensedElt sets the number of condensed elements GetNbCondensedElt returns the number of condensed elements GetMemorySize returns the memory used to store the object (in bytes)

Methods of GlobalGenericMatrix

 GetCoefMass returns the mass coefficient GetCoefStiffness returns the stiffness coefficient GetCoefDamping returns the damping coefficient SetCoefMass sets the mass coefficient SetCoefStiffness sets the stiffness coefficient SetCoefDamping sets the damping coefficient

AssembleMatrix

Syntax

 void AssembleMatrix( Matrix& A, Matrix& mat_elem, GlobalGenericMatrix nat_mat VarComputationProblem_Base& var, CondensationBlockSolver_Base& solver, int offset_row, int offset_col)

Parameters

A (inout)
matrix to modify
mat_elem(inout)
elementary matrix
nat_mat (in)
coefficients
var (inout)
class defining the computation of elementary matrices
solver (inout)
solver (for static condensation)
offset_row (in)
offset for row numbers
offset_col (in)
offset for column numbers

This function adds elementary matrices to the global matrix A. The elementary matrices are computed by AssembleMatrix by calling the method ComputeElementaryMatrix of the object var given as a parameter. Below we reproduce an example located in the file src/Program/Unit/Computation/assemble_matrix.cc. If the global matrix A is not allocated, the function will allocate it and fill it.

Example :


class MyExample : public VarComputationProblem_Base
{
int Nx, Ny; Real_wp dx, dy;

public :
MyExample(double Lx, double Ly, int Nx_, int Ny_)
{
Nx = Nx_; Ny = Ny_;
dx = Lx / Nx; dy = Ly / Ny;
}

// number of elementary matrices to compute
int GetNbElt() const { return Nx*Ny; }

// number of rows of the global matrix
int GetNbRows() const { return (Nx+1)*(Ny+1); }

// verbosity level
int GetPrintLevel() const { return 0; }

// elementary matrix for real coefficients
// i : element number, num_row : row numbers of the element i
// mat_elem : elementary matrix, solver : static condensation solver, coef : coefficients
void ComputeElementaryMatrix(int i, IVect& num_row,
VirtualMatrix<Real_wp>& mat_elem,
CondensationBlockSolver_Base<Real_wp>& solver,
const GlobalGenericMatrix<Real_wp>& coef)
{
// basic exemple alpha M + beta K
// where M is a mass matrix and K stiffness matrix

mat_elem.Reallocate(4, 4);
mat_elem.Zero();

int ie = i %Nx, je = i / Nx;
// row numbers of the considered element
num_row.Reallocate(4);
num_row(0) = je*(Nx+1) + ie;
num_row(1) = num_row(0) + 1;
num_row(2) = num_row(1) + Nx;
num_row(3) = num_row(2) + 1;

Real_wp K00 = 0.5*dy/dx + 0.5*dx/dy;
Real_wp K01 = -0.5*dy/dx, K10 = -0.5*dx/dy;
Real_wp M00 = 0.25*dx*dy;
Real_wp alpha = coef.GetCoefMass();
Real_wp beta = coef.GetCoefStiffness();
mat_elem.SetEntry(0, 0, beta*K00 + alpha*M00);
mat_elem.SetEntry(0, 1, beta*K01);
mat_elem.SetEntry(0, 2, beta*K10);

mat_elem.SetEntry(1, 1, beta*K00 + alpha*M00);
mat_elem.SetEntry(1, 0, beta*K01);
mat_elem.SetEntry(1, 3, beta*K10);

mat_elem.SetEntry(2, 2, beta*K00 + alpha*M00);
mat_elem.SetEntry(2, 3, beta*K01);
mat_elem.SetEntry(2, 0, beta*K10);

mat_elem.SetEntry(3, 3, beta*K00 + alpha*M00);
mat_elem.SetEntry(3, 2, beta*K01);
mat_elem.SetEntry(3, 1, beta*K10);
}

// elementary matrix for complex matrices
void ComputeElementaryMatrix(int i, IVect& num_row,
VirtualMatrix<Complex_wp>& mat_elem,
CondensationBlockSolver_Base<Complex_wp>& solver,
const GlobalGenericMatrix<Complex_wp>& coef)
{
cout << "Not implemented" << endl;
abort();
}

};

int main(int argc, char** argv)
{
InitMontjoie(argc, argv);

MyExample var(2.0, 2.0, 5, 5);

// global matrix
Matrix<Real_wp, Symmetric, ArrayRowSymSparse> Aref;
Matrix<Real_wp, Symmetric, RowSymPacked> mat_elem; // elementary matrix
CondensationBlockSolver_Base<Real_wp> cond_solver; // object to handle static condensation
GlobalGenericMatrix<Real_wp> nat_mat; // coefficients

// setting coefficients (if needed)
nat_mat.SetCoefMass(0.4); nat_mat.SetCoefStiffness(1.3);

// we assemble the matrix A
AssembleMatrix(A, mat_elem, nat_mat, var, cond_solver, 0, 0);
}



ModifyElementaryMatrix

Syntax

 void ModifyElementaryMatrix( int i, IVect& num_dof, Matrix& mat_elem, GlobalGenericMatrix nat_mat)

Parameters

i (in)
element number
num_dof (inout)
row numbers that are kept
mat_elem (inout)
elementary matrix that can be modified due to static condensation
nat_mat (in)
coefficients

This methods modifies the elementary matrix previously computed by the method ComputeElementaryMatrix. If static condensation is applied, it should produce a smaller matrix with only rows that cannot be eliminated. The parameter num_dof contains the rows that are kept.

Example :


class MySolver<Real_wp> : public CondensationBlockSolver_Base<Real_wp>
{
VectReal_wp schur_coef;

public :
void ModifyElementaryMatrix(int i, IVect& num_ddl, VirtualMatrix<Real_wp> mat_elem, const GlobalGenericMatrix<Real_wp>& nat_mat)
{
int j = mat_elem.GetM()-1;
// for instance the last dof is eliminated
VectReal_wp last_row, last_col
mat_elem.GetDenseRow(j, last_row);
mat_elem.GetDenseCol(j, last_col);
Real_wp invA22 = Real_wp(1) / last_row(j);

// updating Schur complement
mat_elem.Resize(j, j);
for (int i = 0; i < j; i++)
if (last_row(i) != Real_wp(0))
for (int k = 0; k < j; k++)
if (last_col(k) != Real_wp(0))

// last dof is removed
num_ddl.Resize(j);

// we can store invA22
schur_coef(i) = invA22;
}

};


SetTreatmentStiffnessInside

Syntax

 void SetTreatmentStiffnessInside( bool t)

This method informs if there are stiffness terms to handle for the static condensation. It is used by local implicit schemes.

SetElementNumber

Syntax

 void SetElementNumber(int local_num , int global_num )

This method sets the local element number (among condensed elements) and the global element number. Usually these two numbers are equal (if all the elements contribute to the condensed matrix).

GetCondensedElementNumber

Syntax

 int GetCondensedElementNumber() const

This method returns the current condensed element number.

GetGlobalElementNumber

Syntax

 int GetGlobalElementNumber() const

This method returns the current global element number.

GetNbCondensedElt

Syntax

 int GetNbCondensedElt() const

This method returns the number of condensed elements (elements that contribute to the condensed matrix).

SetNbCondensedElt

Syntax

 void SetNbCondensedElt(int n )

This method sets the number of condensed elements (elements that contribute to the condensed matrix).

GetMemorySize

Syntax

 size_t GetMemorySize() const

This method returns the memory used to store the object (in bytes).

GetCoefMass

Syntax

 T GetCoefMass() const

This method returns the coefficient associated with the mass matrix. If the equation is used in time-domain, it corresponds to the coefficient for the second derivatives in time. If the equation only involves first derivatives in time, the coefficient applies to first derivatives. This coefficient is the coefficient α detailed in the equation below.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


GetCoefDamping

Syntax

 T GetCoefDamping() const

This method returns the coefficient associated with the damped matrix. If the equation is used in time-domain, it corresponds to the coefficient for the first derivatives in time if the equation is a second-order formulation in time. If the equation is a first-order formulation, it corresponds to coefficients associated with damping terms. It correspdonds to the coefficient β in the equations below.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


GetCoefStiffness

Syntax

 T GetCoefStiffness() const

This method returns the coefficient associated with the stiffness matrix. If the equation is used in time-domain, it corresponds to the coefficient for terms without time-derivatives. It correspdonds to the coefficient γ in the equations below.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


SetCoefMass

Syntax

 void SetCoefMass(T coef) const

This method sets the coefficient associated with the mass matrix. This coefficient is the coefficient α detailed in the description of GetCoefMass.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


SetCoefDamping

Syntax

 void SetCoefDamping(T coef) const

This method sets the coefficient associated with the damped matrix. This coefficient is the coefficient β detailed in the description of GetCoefDamping.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


SetCoefStiffness

Syntax

 void SetCoefStiffness(T coef) const

This method sets the coefficient associated with the stiffness matrix. This coefficient is the coefficient γ detailed in the description of GetCoefStiffness.

Example :

   GlobalGenericMatrix<Real_wp> nat_mat;

// you can set coefficients
nat_mat.SetCoefMass(Real_wp(0.4));
nat_mat.SetCoefDamping(Real_wp(0.8));
nat_mat.SetCoefStiffness(Real_wp(1.5));

// and retrieve them
Real_wp m = nat_mat.GetCoefMass();
Real_wp sig = nat_mat.GetCoefDamping();
Real_wp s = nat_mat.GetCoefStiffness();


SetLevel

Syntax

 void SetLevel(int lvl)

This method sets the level, such that the matrix-vector product will be performed only for elements of the selected level. There are predefined levels (which correspond to negative numbers):

• ALL_LEVELS : all the elements of the mesh are considered
• LVL_PML : all the elements in PMLs are considered
• LVL_NOPML : all the elements outside PMLs are considered

Example :

    // class for solving acoustic equation
HyperbolicProblem<AcousticEquation<Dimension2> > var;

// assuming that var is correctly constructed
// you can retrieve the different levels
// each level is usually associated with an area with a given time step
MatrixVectorProductLevel& list_level = var.GetTimeLevelDistribution();

// you can select a level with SetLevel
list_level.SetLevel(1);

// and loop over selected elements
for (int i0 = 0; i0 < list_level.GetNbElt(); i0++)
{
int i = list_level.GetElementNumber(i0);
// and perform the computation for the element i
}

// to select a predefined level
list_level.SetLevel(MatrixVectorProductLevel::ALL_LEVELS);


GetLevelArray

Syntax

 Vector GetLevelArray()

This method returns the list of elements for each level.

Example :

    // class for solving acoustic equation
HyperbolicProblem<AcousticEquation<Dimension2> > var;

// assuming that var is correctly constructed
// you can retrieve the different levels
// each level is usually associated with an area with a given time step
MatrixVectorProductLevel& list_level = var.GetTimeLevelDistribution();

// you can select a level with SetLevel
list_level.SetLevel(1);

// you can display the elements of each level
Vector<IVect>& lvl_array = list_level.GetLevelArray();
cout << "Elements for level 0 = " << lvl_array(0) << endl;
cout << "Elements for level 1 = " << lvl_array(1) << endl;


SetLevelArray

Syntax

 void SetLevelArray(Vector& liste )

This method sets the list of elements for each level.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);


GetNbElt

Syntax

 int GetNbElt() int GetNbElt(int level )

This method returns the number of elements of a given level. If no argument is given, it returns the number of elements of the selected level (with method SetLevel).

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);

// number of elements for level 1 ?
int n1 = prod_level.GetNbElt(1);

// setting a level
prod_level.SetLevel(2);
// and number of elements on the selected level (2 here)
int n2 = prod_level.GetNbElt();


SetNbElt

Syntax

 void SetNbElt(int num_elt , int lnum_elt_pml )

This method sets the number of elements (in total) and the number of elements inside PML.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);

// you can select all elements inside the PML
prod_level.SetLevel(MatrixVectorProductLevel::LVL_PML);

// and loop over them
int n0 = prod_level.GetNbElt()
for (int i0 = 0; i0 < n0; i0++)
{
int i = prod_level.GetElementNumber(i0);
}


GetElementNumber

Syntax

 int GetElementNumber(int i )

This method returns the global element number of the i-th element of the selected level.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);

// you can select all elements inside the PML
prod_level.SetLevel(MatrixVectorProductLevel::LVL_PML);

// and loop over them
int n0 = prod_level.GetNbElt()
for (int i0 = 0; i0 < n0; i0++)
{
int i = prod_level.GetElementNumber(i0); // global element number
}


GetLocalElementNumber

Syntax

 int GetLocalElementNumber()

This method returns the local element number. This method is used in an alternative way to browse the elements of a given level, as detailed in the example below.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);

// you can select all elements inside the PML
prod_level.SetLevel(MatrixVectorProductLevel::LVL_PML);

// and loop over them
int n0 = prod_level.GetNbElt()
// first method : loop over local numbers
for (int i0 = 0; i0 < n0; i0++)
{
int i = prod_level.GetElementNumber(i0); // global element number (in the mesh)
}

// second method : loop over global numbers
// in that case, you need to browse these elements in ascending order
prod_level.SetLevel(0);

for (int i = 0; i < var.mesh.GetNbElt(); i++)
if (var.TreatElement(i))
{
// you can retrieve the local number i0
int i0 = prod_level.GetLocalElementNumber();
}


TreatElement

Syntax

 bool TreatElement(int i )

This method returns true is the element i belongs to the selected level. This method is used in an alternative way to browse the elements of a given level, as detailed in the example below.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

MatrixVectorProductLevel prod_level;
// setting the number of elements (total number and number in PMLs)
prod_level.SetNbElt(var.mesh.GetNbElt(), var.GetNbEltPML());

// setting the element numbers for differents levels
Vector<IVect> list_level(3);
for (int i = 0; i < var.mesh.GetNbElt(); i++)
list_level(i%3).PushBack(i)

prod_level.SetLevelArray(list_level);

// you can select all elements inside the PML
prod_level.SetLevel(MatrixVectorProductLevel::LVL_PML);

// and loop over them
int n0 = prod_level.GetNbElt()
// first method : loop over local numbers
for (int i0 = 0; i0 < n0; i0++)
{
int i = prod_level.GetElementNumber(i0); // global element number (in the mesh)
}

// second method : loop over global numbers
// in that case, you need to browse these elements in ascending order
prod_level.SetLevel(0);

for (int i = 0; i < var.mesh.GetNbElt(); i++)
if (var.TreatElement(i))
{
// you can retrieve the local number i0
int i0 = prod_level.GetLocalElementNumber();
}

// TreatElement can only be used in a loop over all the elements (in ascending order)
// it cannot be used solely to know if a single elements belongs to the selected level


print_level

This attribute is the verbosity level. It is usually modified through data file (by inserting a line PrintLevel = lvl).

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

// you can modify the attribute manually
var.print_level = 2;


nb_unknowns

This attribute is the number of unknowns of the solved equation. It cannot be modified, this attribute is actually set in the class defining the equation. It is equal to

• The number of continuous unknowns (that need to be numbered) for a continuous formulation. For instance, for the Helmholtz equation, the number of unknowns is equal to 1 (even though a mixed formulation is used with a discontinuous unknown v).
• The number of discontinuous unknowns for a discontinuous Galerkin formulation. For instance, for the Helmholtz equation (solved with LDG method), the number of unknowns is equal to d+1 (where d is the dimension).
• The number of volume unknowns for HDG. For instance for the Helmholtz equation (solved with HDG method), the number of unknowns is equal to d+1, the surface unknown λ is not counted.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

int n = var.nb_unknowns;


nb_unknowns_scal

This attribute is the number of scalar unknowns of the solved equation. It cannot be modified, this attribute is actually set in the class defining the equation. It is equal to

• The number of continuous unknowns (that need to be numbered) for a continuous formulation. For a continuous formulation, nb_unknowns and nb_unknowns_scal are equal. For instance, for the Helmholtz equation, the number of scalar unknowns is equal to 1.
• The number of scalar discontinuous unknowns for a discontinuous Galerkin formulation. For instance, for the Helmholtz equation (solved with LDG method), the number of scalar unknowns is equal to 1 (only unknown u is assumed to be scalar).
• The number of scalar volume unknowns for HDG. For instance for the Helmholtz equation (solved with HDG method), the number of scalar unknowns is equal to 1.

For discontinuous Galerkin formulation, this number is used in order to use staggered time-schemes (with a scalar unknown u and a vectorial unknown v).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

int n = var.nb_unknowns_scal;


nb_unknowns_vec

This attribute is the number of vectorial unknowns of the solved equation. It cannot be modified, this attribute is actually set in the class defining the equation. It is equal to

• The number of discontinuous unknowns (for a mixed formulation) for a continuous formulation. For instance, for the Helmholtz equation, the number of vectorial unknowns is equal to d where d is the dimension.
• The number of vectorial discontinuous unknowns for a discontinuous Galerkin formulation. For instance, for the Helmholtz equation (solved with LDG method), the number of vectorial unknowns is equal to d (only unknown v is assumed to be vectorial).
• The number of vectorial volume unknowns for HDG. For instance for the Helmholtz equation (solved with HDG method), the number of vectorial unknowns is equal to d.

For discontinuous Galerkin formulation, this number is used in order to use staggered time-schemes (with a scalar unknown u and a vectorial unknown v).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

int n = var.nb_unknowns_vec;


nb_unknowns_hdg

This attribute is the number of surface unknowns of the solved equation. It cannot be modified, this attribute is actually set in the class defining the equation. It is significant only for HDG formulation. If an equation is not solved with HDG, it is usually set to 0.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

int n = var.nb_unknowns_hdg;


nb_components_en

This attribute is the number of components of the trace of the solution. It cannot be modified, this attribute is actually set in the class defining the equation. It is mainly used for the transparent condition, wich needs to compute the solution on a closed surface : the trace of the solution and derivatives. For instance, for Helmholtz equation, it will need to compute u and on the surface. For Helmholtz equation, it will be equal to 1, whereas for Maxwell's equations, it is equal to three (in 3-D) because the electric field has three components.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

// you can access to the number of components to store E \times n
int n = var.nb_components_en;


nb_components_hn

This attribute is the number of components of the trace of the derivative of the solution. It cannot be modified, this attribute is actually set in the class defining the equation. It is mainly used for the transparent condition, wich needs to compute the solution on a closed surface : the trace of the solution and derivatives. For instance, for Helmholtz equation, it will need to compute u and on the surface. For Helmholtz equation, it will be equal to 1 (du/dn is scalar), whereas for Maxwell's equations, it is equal to three (in 3-D) because the magnetic field has three components.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

// you can access to the number of components to store H \times n
int n = var.nb_components_hn;


type_element

This attribute is the type of the finite element to use for the main unknown. This type is an integer that can be equal to 1, 2 or 3 :

• 1 : use of scalar finite element for the main unknown (discretization of H1 or L2 space)
• 2 : use of edge element (discretization of H(curl) space)
• 3 : use of facet element (discretization of H(div) space)

The finite elements classes are described in the section devoted to finite elements. They depend on a template parameter which corresponds to the type of the finite element. The attribute type_element cannot be modified and is defined in the class defining the equation. The main unknown is the unknown associated with the first mesh numbering. Most of solved equations rely on a single mesh numbering.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

// you can access to the type of finite element used to solve the equation
int type_elt = var.type_element;


other_type_element

This attribute is the type of the finite element to use for other unknowns. This type is an integer that can be equal to 1, 2 or 3 :

• 1 : use of scalar finite element for another unknown (discretization of H1 or L2 space)
• 2 : use of edge element (discretization of H(curl) space)
• 3 : use of facet element (discretization of H(div) space)

The finite elements classes are described in the section devoted to finite elements. They depend on a template parameter which corresponds to the type of the finite element. The attribute other_type_element cannot be modified and is defined in the class defining the equation. The main unknown is the unknown associated with the first mesh numbering, while other unknowns are associated with the next mesh numberings. Most of solved equations rely on a single mesh numbering. In that case, the array other_type_element is void. If there are at least two mesh numberings, other_type_element(i) will be the type of the finite element for the i+1-th mesh numbering.

Example :

    EllipticProblem<HarmonicMaxwellEquation_HcurlAxi> var;

// you can access to the type of finite element used for the main unknown
int type_elt = var.type_element;

int type_elt2 = var.other_type_element(0);


finite_element_name

This attribute stores the finite element name (used for the main unknown). It corresponds to the parameter given when filling TypeElement in the data file.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;

// after constructing the mesh

string name_elt = var.finite_element_name ; // should be QUADRANGLE_LOBATTO


name_other_elements

This attribute stores the finite element name used for other unknowns. It corresponds to the parameter given when filling TypeElement in the data file.

Example :

    EllipticProblem<HarmonicMaxwellEquation_HcurlAxi> var;;

// after setting finite element

string name_elt = var.finite_element_name ; // for the main unknown
string name_elt2 = var.name_other_elements(0); // and second unknown


mesh_num_unknown

This attribute stores the mesh numberings to use for each unknown. For most equations with only one mesh numbering, this array is filled with zeros.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;
var.InitIndices(20);
var.SetTypeEquation("LAPLACE");

// mesh numbering to use for the unknown 1 ?
int n = var.mesh_num_unknown(1);


offset_dof_unknown

This attribute stores the cumulated number of degrees of freedom for each unknown. For instance, let us consider three unknowns and N1, N2, N3 degrees of freedom for each of them. offset_dof_unknown will be equal to (0, N1, N1+N2, N1+N2+N3).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// starting row for the first unknown
int n0 = var.offset_dof_unknown(0); // should be 0

// starting row for the second unknown
int n1 = var.offset_dof_unknown(1); // should be equal to N1

// starting row for the third unknown
int off2 = var.offset_dof_unknown(2); // should be equal to N1+N2


offset_dof_condensed

This attribute stores the cumulated number of degrees of freedom for each unknown (with static condensation). For instance, let us consider three unknowns and N1, N2, N3 condensed degrees of freedom for each of them. offset_dof_condensed will be equal to (0, N1, N1+N2, N1+N2+N3). The condensed number of degrees of freedom is the number of degrees of freedom after static condensation (internal dofs of elements are removed).

Example :

    EllipticProblem<ElasticEquation<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// starting row for the first unknown
int n0 = var.offset_dof_condensed(0); // should be 0

// starting row for the second unknown (condensed system)
int n1 = var.offset_dof_condensed(1); // should be equal to N1

// starting row for the third unknown (condensed system)
int off2 = var.offset_dof_condensed(2); // should be equal to N1+N2


dg_formulation

This attribute stores the type of the formulation used to solve the equation. It can be equal to

• ElementReference_Base::CONTINUOUS :the equation is solved with continuous finite elements.
• ElementReference_Base::DISCONTINUOUS :the equation is solved with discontinuous Galerkin formulation.
• ElementReference_Base::HDG :the equation is solved with Hybridizable Discontinuous Galerkin formulation.

The formulation is specified in the definition of the equation.

Example :

    EllipticProblem<ElasticEquation<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// which formulation is used ?
int dg_form = var.dg_formulation;


sipg_formulation

This attribute is equal to true if Interior Penalty Discontinuous Galerkin is used. If it is false, Local Discontinuous Galerkin is used

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// dg_form should be equal to ElementReference_Base::DISCONTINUOUS
int dg_form = var.dg_formulation;

// sipg should be true (LAPLACE_SIPG given in ConstructAll)
bool sipg = var.sipg_formulation;


compute_dfjm1

This attribute is equal to true if the jacobian matrices will be computed (and stored). For most of equations, they are computed.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;;

// if you want to force that jacobian matrices are stored
var.compute_dfjm1 = true;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;


alpha_penalization

This attribute is a penalty parameter used in discontinuous Galerkin formulations. It can be modified through PenalizationDG in the data file.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to know parameter alpha for penalization
Real_wp alpha = var.alpha_penalization;


delta_penalization

This attribute is a penalty parameter used in discontinuous Galerkin formulations. It can be modified through PenalizationDG in the data file.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to know parameter alpha for penalization
Real_wp alpha = var.alpha_penalization;

// to know parameter delta for penalization
Real_wp delta = var.delta_penalization;


upwind_fluxes

If true, upwind fluxes are used in the discontinuous Galerkin formulation. It can be specified through PenalizationDG in the data file. If false, centered fluxes are used with penalty terms. Depending on the penalty coefficients (alpha_penalization, delta_penalization) and the considered equation, it can coincide with upwind fluxes. It is the case for Helmholtz equation, if alpha and delta are set to -1. If alpha and delta are set to zero, and upwind_fluxes is false, it correspond to pure centered fluxes (without penalty terms).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to know if upwind fluxes are used
Real_wp upwind = var.upwind_fluxes;


automatic_choice_penalization

If this attribute is true, the penalty coefficient is automatically chosen (depending on the order of approximation and the mesh size). It is usually set through CoefficientPenalization in the data file. It is significant only for Interior Penalty Discontinuous formulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to know if the penalty coefficients are automatically chosen
bool auto_param = var.automatic_choice_penalization;


Glob_CoefPenalDG

This array stores the penalty coefficients for each boundaries of the mesh (edges in 2-D, faces in 3-D). It is significant only for Interior Penalty Discontinuous formulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to display the penalty coefficients
cout << "Penalty coefficients " << var.Glob_CoefPenalDG << endl;


Glob_CoefPenalDG

This array stores the penalty coefficients for each boundaries of the mesh (edges in 2-D, faces in 3-D). It is significant only for Interior Penalty Discontinuous formulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to display the penalty coefficients
cout << "Penalty coefficients " << var.Glob_CoefPenalDG << endl;


mesh_data

This array stores the parameters used to construct the mesh. It is usually a vector of length since only one mesh is used in the simulations. It corresponds to parameters given in the field FileMesh of the data file.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// to display the mesh parameters
cout << "Mesh parameters " << var.mesh_data(0) << endl;


exit_if_no_boundary_condition

If true, the computation will be stopped is there are faces (or edges in 2-D) on a boundary without boundary conditions. By default, this attributed is true, in order to detect if the user forgot lines ConditionReference in the data file. It can be modified by setting Exit_IfNo_BoundaryCondition in the data file.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// if you do not want to stop the computation because of isolated faces without boundary conditions
var.exit_if_no_boundary_condition = false;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;


var_chrono

It is the stopwatch used in Montjoie to compute the time elapsed during the different stages of the simulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve a computation time
Real_wp t_mesh = var.var_chrono.GetSeconds("MeshGeneration");

// or display all stopwatches
var.var_chrono.DisplayAll();


GetNbDof, SetNbDof

Syntax

 int GetNbDof() const void SetNbDof(int N)

The method GetNbDof returns the total number of degrees of freedom of the solved problem. This number is computed by adding the degrees of freedom for each unknown and optional degrees of freedom due to boundary conditions (or other models). The method SetNbDof changes this number.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the total number of degrees of freedom
int N = var.GetNbDof();

// you want to add a row to the linear system
var.SetNbDof(N+1);


GetNbMeshDof

Syntax

 int GetNbMeshDof() const int GetNbMeshDof(int n) const

This method returns the number of degrees of freedom for the n-th mesh numbering. If no argument is given, it returns the number of degrees for the first mesh numbering.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the number of dofs for the main unknown
int Nvol = var.GetNbMeshDof();


GetNbMeshNumberings

Syntax

 int GetNbMeshNumberings() const

This method returns the number of mesh numberings.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the number of mesh numberings
int nb_mesh_num = var.GetNbMeshNumberings(); // should be one


GetMeshNumberingBase

Syntax

 const MeshNumbering_Base& GetMeshNumberingBase(int n = 0) const

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the mesh numbering for the main unknown
const MeshNumbering_Base<Real_wp>& mesh_num = var.GetMeshNumberingBase();


GetOffsetDofUnknown

Syntax

 int GetOffsetDofUnknown(int n) const

This method returns offset_dof_unknown(n) (see details in the description of offset_dof_unknown).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the first row for unknown 2
int N = var.GetOffsetDofUnknown(2);


GetOffsetDofCondensed

Syntax

 int GetOffsetDofCondensed(int n) const

This method returns offset_dof_condensed(n) (see details in the description of offset_dof_condensed).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the first row for unknown 2
int N = var.GetOffsetDofCondensed(2);


GetDefaultOrder

Syntax

 int GetDefaultOrder(int r) const

This method returns the order that will be used for the mesh. It corresponds to the order given when filling OrderGeometry in the data file. If this field is not present, the default order is initialized with the order given in OrderDiscretization.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the order used for the geometry
int r = var.GetDefaultOrder();


FormulationDG

Syntax

 int FormulationDG() const

This method returns the attribute dg_formulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// if you want to retrieve the formulation used to solve the equation
int dg_form = var.FormulationDG();


ComputeDFjm1

Syntax

 bool ComputeDFjm1() const

This method returns the attribute compute_dfjm1.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// jacobian matrices are computed ?
bool dfj_computed = var.ComputeDFjm1();


FirstOrderFormulation

Syntax

 bool FirstOrderFormulation() const

This method returns true if a mixed formulation is used. For example, the Helmholtz equation :

is transformed into

which is a mixed formulation of Helmholtz equation. This method is relevant only for continuous formulations. Only some equations have an implementation of a mixed formulation.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// mixed formulation is used ?
bool mix_form = var.FirstOrderFormulation();


FirstOrderFormulationDG

Syntax

 bool FirstOrderFormulationDG() const

This method returns true if only first-order derivatives (in space) appear in the discontinuous Galerkin formulation.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;

// only first-order derivatives in the formulation ?
bool first_order_deriv = var.FirstOrderFormulationDG();


SetFirstOrderFormulation

Syntax

 bool SetFirstOrderFormulation(bool mix = true)

This method enables (or disables) the use of a mixed formulation (as detailed in FirstOrderFormulation).

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;;

// if you want to force a first-order formulation :
var.SetFirstOrderFormulation(true);

// constructing the problem (mesh, finite element)
All_LinearSolver* solver;


UseExactIntegrationElement

Syntax

 bool UseExactIntegrationElement() const

This method returns true if we consider that the used finite element will ensure an exact integration. It is used for some equations, in order to use more efficients formulations. This feature can be used by inserting a line with ExactIntegration in the data file. It is up to the user to be sure that the used finite element ensures an exact integration.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// if the user put a line ExactIntegration = YES, it should return true :
bool exact_integration = var.UseExactIntegrationElement();


GetOverIntegration

Syntax

 int GetOverIntegration() const

This method returns the over-integration used in the simulations. This integer can be modified by inserting a line with ExactIntegration in the data file. Otherwise, it is equal to 0 (no over-integration).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// is there over-integration ?
// integration or order r+p (where r is the order of the finite element)
int p = var.GetOverIntegration();


GetXmin

Syntax

 Real_wp GetXmin() const

This method returns the minimum of x-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp xmin_d = var.mesh.GetXmin();


GetYmin

Syntax

 Real_wp GetYmin() const

This method returns the minimum of y-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp ymin_d = var.mesh.GetYmin();


GetZmin

Syntax

 Real_wp GetZmin() const

This method returns the minimum of z-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension3> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp zmin_d = var.mesh.GetZmin();


GetXmax

Syntax

 Real_wp GetXmax() const

This method returns the maximum of x-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp xmax_d = var.mesh.GetXmax();


GetYmax

Syntax

 Real_wp GetYmax() const

This method returns the maximum of y-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension2> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp ymax_d = var.mesh.GetYmax();


GetZmax

Syntax

 Real_wp GetZmax() const

This method returns the maximum of z-coordinates of the physical domain. The physical domain does not include PML layers.

Example :

    EllipticProblem<LaplaceEquationDG<Dimension3> > var;;

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

// bounds of the physical domain
Real_wp xmin = var.GetXmin();
Real_wp xmax = var.GetXmax();
Real_wp ymin = var.GetYmin();
Real_wp ymax = var.GetYmax();
Real_wp zmin = var.GetZmin();
Real_wp zmax = var.GetZmax();

// if you want to know bounds of the computational domain
// (physical domain + PMLs), use mesh
Real_wp zmax_d = var.mesh.GetZmax();


SetComputationalDomain

Syntax

 void SetComputationalDomain(Real_wp xmin, Real_wp xmax, Real_wp ymin, Real_wp ymax, Real_wp zmin, Real_wp zmax)

This method sets the bounds of the physical domain. The physical domain does not include PML layers. The name of the method is misleading since the methods changes the bounds of physical domain (and not the computational domain).

Example :

    EllipticProblem<LaplaceEquationDG<Dimension3> > var;;

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

// if you want to change the bounds of the physical domain
var.SetComputationalDomain(-2.0, 2.0, -2.0, 2.0, -1.0, 1.0);


GetOmega, SetOmega

Syntax

 Real_wp GetOmega() const void SetOmega(Real_wp omega)

The method GetOmega returns the pulsation (which is equal to the frequency multiplied by 2 π) . The method SetOmega sets the pulsation.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// pulsation omega ?
Real_wp omega = var.GetOmega();


GetSquareOmega

Syntax

 Real_wp GetSquareOmega() const

This method returns the square of the pulsation .

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// square of pulsation omega^2 ?
Real_wp omega2 = var.GetSquareOmega();


GetMiomega

Syntax

 void GetMiomega(Real_wp& z) const void GetMiomega(Complex_wp& z) const

This method sets z to - i ω where ω is the pulsation. If z is a real number, it is set to 1.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// -i omega ?
Complex_wp m_iomega; var.GetMiomega(m_iomega);


GetMomega2

Syntax

 void GetMomega2(Real_wp& z) const void GetMomega2(Complex_wp& z) const

This method sets z to - ω2 where ω is the pulsation. If z is a real number, it is set to 1.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// -omega^2 ?
Complex_wp m_omega2; var.GetMomega2(m_omega2);


GetFrequency, SetFrequency

Syntax

 Real_wp GetFrequency() const void SetFrequency(Real_wp f)

The method GetFrequency returns the frequency. The method SetFrequency sets the frequency.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// frequency ?
Real_wp f = var.GetFrequency();


Syntax

This method returns the characteristical length. It represents the length of 1 in the mesh. For example if this characteristical length is equal to 1e-9, it means that the units in the mesh are in nanometers (instead of meters). It is modified by adding a line with WavelengthAdim. It is used only if Adimensionalization has been set to YES, and for Helmholtz equations or Maxwell's equations.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// characteristical length


GetDimension

Syntax

 int GetDimension() const

This method returns the dimension of the solved problem (2 or 3).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// dimension should be three here
int d = var.GetDimension();


GetNbComponentsUnkown

Syntax

 int GetNbComponentsUnkown(int n ) const

This method returns the number of components for an unknown associated with the mesh numbering n . It can be 1 (for a scalar unknown), 2 or 3 (vectorial unknown).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of components for the first mesh numbering ?
int p = var.GetNbComponentsUnknown(0);


Syntax

This method returns the number of components for the gradient of an unknown associated with the mesh numbering n . It can be 1 (for a scalar unknown), 2 or 3 (vectorial unknown). For edge elements, we consider the curl (instead of the gradient) and for facet elements, the divergence is considered.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of "gradient" components for the first mesh numbering ?


GetNbLocalDof

Syntax

 int GetNbLocalDof(int i ) const int GetNbLocalDof(int i , int n ) const

This method returns the number of degrees of freedom for the element i . The second argument n is the number of the considered mesh numbering. If not provided, n is equal to zero.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of degrees of freedom for a given element
int i = 3;
int nb_dof = var.GetNbLocalDof(i);


GetNbSurfaceDof

Syntax

 int GetNbSurfaceDof(int i ) const int GetNbSurfaceDof(int i , int n ) const

This method returns the number of degrees of freedom for the surfacic element i . The second argument n is the number of the considered mesh numbering. If not provided, n is equal to zero. This method is relevant for HDG formulation where the surface unknowns λ are discretized with surfacic finite elements.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of degrees of freedom for a given face
int i = 3;
int nb_dof = var.GetNbSurfaceDof(i);


GetNbDofBoundaries

Syntax

 int GetNbDofBoundaries(int i ) const int GetNbDofBoundaries(int i , int n ) const

This method returns the number of degrees of freedom for the element i , only degrees of freedom associated with the vertices/edges/faces are counted. We do not count internal degrees of freedom (associated with the inside of the element). The second argument n is the number of the considered mesh numbering. If not provided, n is equal to zero.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of degrees of freedom on the boundaries of an element
int i = 3;
int nb_dof_boundaries = var.GetNbDofBoundaries(i);


Syntax

This method returns the number of quadrature points inside the element i .

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of quadrature points for volume integrals of element i
int i = 3;


WeightsND

Syntax

 const VectReal_wp& WeightsND(int i ) const

This method returns the weights associated with quadrature points inside the element i .

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// number of quadrature points for volume integrals of element i
int i = 3;

// associated weights
const VectReal_wp& weights = var.WeightsND(i);


ElementInsidePML

Syntax

 bool ElementInsidePML(int i ) const

This method returns returns true if the element i is located in PMLs.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// is element i a PML element ?
int i = 3;
bool inside_pml = var.ElementInsidePML(i);


GetReferenceElementBase

Syntax

 const ElementReference_Base& GetReferenceElementBase(int i ) const const ElementReference_Base& GetReferenceElementBase(int i , int n ) const

This method returns the finite element associated with the element i . It returns a reference to the base class of the finite element. If you already know the dimension, you can call GetReferenceElement instead. The second argument is optional (if you want to select another mesh numbering).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// reference element associated with an element
int i = 3;
const ElementReference_Base& Fb = var.GetReferenceElementBase(i);


GetSurfaceElementBase

Syntax

 const ElementReference_Base& GetSurfaceElementBase(int i ) const const ElementReference_Base& GetSurfaceElementBase(int i , int n ) const

This method returns the finite element associated with the surfacic element i . It returns a reference to the base class of the finite element. If you already know the dimension, you can call GetReferenceElement instead. The second argument is optional (if you want to select another mesh numbering). This method is relevant for HDG formulation (the surfacic elements are associated with surfacic unknowns λ).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// reference element associated with a face
int i = 3;
const ElementReference_Base& Fb = var.GetSurfaceElementBase(i);


WriteMesh

Syntax

 void WriteMesh(string file_name) const

This method writes the mesh in a file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

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

// if you want to write the mesh
var.WriteMesh("test.mesh");


Syntax

This method forces to use same numbers for periodic dofs. This can also be achieved by inserting a line UseSameDofsForPeriodicCondition = YES in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// if you want to force same numbers for periodic dofs

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


InitIndices

Syntax

 void InitIndices(int N)

This method allocates the arrays storing physical indexes (such as ϵ, μ for Maxwell's equations). The size of the arrays should be equal to N+1 such that you can use a reference equal to N.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
// we usually put the maximal reference number (Physical Volume in 3-D)
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file
// the lines with MateriauDielec or PhysicalMedia must have a number between 1 and N (included)

// you can continue the computations


GetNbPhysicalIndices

Syntax

 int GetNbPhysicalIndices()

This method returns the number of physical indices stored in the class. It corresponds to N+1 where N is the argument given when InitIndices has been called.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
// we usually put the maximal reference number (Physical Volume in 3-D)
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file
// the lines with MateriauDielec or PhysicalMedia must have a number between 1 and N

// you can continue the computations

// if you need to retrieve N+1
int np1 = var.GetNbPhysicalIndices();


SetIndices

Syntax

 void SetIndices(int ref, const VectString& param)

This method sets the physical indexes associated with the reference ref. It is the method called when a line MateriauDielec is present in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
// we usually put the maximal reference number (Physical Volume in 3-D)
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file
// the lines with MateriauDielec or PhysicalMedia must have a number between 1 and N (included)

// if you want to add other indexes (not given in the data file)
VectString param(4);
// for Helmholtz equation, we have isotropy rho mu sigma
param(0) = "ISOTROPE"; param(1) = "2.4"; param(2) = "1.5"; param(3) = "0.1";
// it is equivalent to place : MateriauDielec = 4 ISOTROPE 2.4 1.5 0.1
var.SetIndices(4, param);

// you can continue the computations


SetPhysicalIndex

Syntax

 void SetPhysicalIndex(string index_name, int ref, const VectString& param)

This method sets a single physical index associated with the reference ref. It is the method called when a line PhysicalMedia is present in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
// we usually put the maximal reference number (Physical Volume in 3-D)
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file
// the lines with MateriauDielec or PhysicalMedia must have a number between 1 and N (included)

// if you want to add/modify a single index(not given in the data file)
VectString param(4);
// for Helmholtz equation, we have rho mu sigma
param(0) = "ISOTROPE"; param(1) = "2.4";
// it is equivalent to place : PhysicalMedia = rho 4 ISOTROPE 2.4
var.SetPhysicalIndex("rho", 4, param);

// you can continue the computations


IsVaryingMedia

Syntax

 bool IsVaryingMedia(int ref) const bool IsVaryingMedia(int m, int ref) const

This method returns true if one of the physical index of reference ref is variable. In the second syntax, (which is used only for unsteady simulations with discontinuous Galerkin), the method returns true if the physical index associated with the unknown m is variable (mass and/or damping terms) for the reference ref.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file

// you can continue the computations

// variable index for reference 5 ?
bool variable = var.IsVaryingMedia(5);


GetPhysicalIndexName

Syntax

 string GetPhysicalIndexName(int m) const

This method returns the name of the m-th physical index. It depends on the solved equation. For example, it will return "rho" for Helmholtz equation and m = 0.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ");

// then you can read the data file

// name of the physical index 1 ?
string name = var.GetPhysicalIndexName(1);


GetCoefficientPenaltyStiffness

Syntax

 Real_wp GetCoefficientPenaltyStiffness(int ref) const

This method returns the physical coefficient associated with reference ref. This coefficient will be used for the penalty terms of the Interior Penalty Discontinuous Galerkin method. For Helmholtz equation, it will correspond to the amplitude of μ.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ_SIPG");

// then you can read the data file

// coefficient for penalty terms
Real_wp mu_coef = var.GetCoefficientPenaltyStiffness(5);


GetVelocityOnElements

Syntax

 void GetVelocityOnElements(VectReal_wp& velocity, Mesh& mesh)

This method fills the array velocity with the velocities for each element of the mesh.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ_SIPG");

// then you can read the data file

// the mesh is constructed

// velocities for all elements
VectReal_wp velocity;
var.GetVelocityOnElements(velocity, var.mesh);


GetVelocityOfInfinity

Syntax

 void GetVelocityOfInfinity()

This method returns the velocity associated with the infinity. This quantity is used for absorbing boundary conditions. It is correctly computed if the user has inserted a line ReferenceInfinity in the data file, such that the physical indexes are known at the infinity.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ_SIPG");

// then you can read the data file

// the mesh is constructed

// velocity for reference 5 ?
Real_wp c = var.GetVelocityOfMedia(5);

// and for the media at infinity ?
Real_wp c_inf = var.GetVelocityOfInfinity();


CopyInputData

Syntax

 void CopyInputData(const VarProblem_Base& var)

This method copies input parameters from another similar object.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var, var_bis;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("HELMHOLTZ_SIPG");

// then you can read the data file (input parameters are read from a file)

// you can copy input parameters to object var_bis
var_bis.CopyInputData(var);


IsSymmetricProblem

Syntax

 bool IsSymmetricProblem(bool eigen=false) const

This method returns true if the finite element matrix is symmetric. The second argument is optional. If eigen is equal to true, we ask if the eigenproblem is symmetric (with a positive definite mass matrix).

Example :

    EllipticProblem<LaplaceEquation<Dimension3> > var;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("LAPLACE");

// then you can read the data file (input parameters are read from a file)

// symmetric finite element matrix ?
bool sym = var.IsSymmetricProblem();


IsSymmetricMassMatrix

Syntax

 bool IsSymmetricMassMatrix() const

This method returns true if the mass matrix is symmetric.

Example :

    EllipticProblem<LaplaceEquation<Dimension3> > var;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("LAPLACE");

// then you can read the data file (input parameters are read from a file)

// symmetric mass matrix ?
bool sym = var.IsSymmetricMassMatrix();


IsComplexProblem

Syntax

 bool IsComplexProblem() const

This method returns true if the solved equation has to be solved with complex numbers.

Example :

    EllipticProblem<LaplaceEquation<Dimension3> > var;

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation("LAPLACE");

// then you can read the data file (input parameters are read from a file)

// complex numbers ? LaplaceEquation is solved with real numbers => we should get false
bool cplx = var.IsComplexProblem();


ComputeMeshAndFiniteElement

Syntax

 void ComputeMeshAndFiniteElement(string name_finite_element) void ComputeMeshAndFiniteElement(string name_finite_element, bool split_mesh)

This method computes the mesh and constructs the finite elements. The name of the main finite element (to use for the main unknown) is given as argument. It corresponds to the line TypeElement of the data file. The second argument is optional and is meaningful in parallel. If true, the mesh is distributed between the different processors.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);


PerformOtherInitializations

Syntax

 void PerformOtherInitializations()

This method performs other initializations (if there are specific models to initialize or specific boundary conditions). This method is called after constructing the mesh and finite elements.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();


SetTypeEquation

Syntax

 void SetTypeEquation(string name_finite_element)

This method sets the type of equation to solve. The name of the equation gives also the formulation used to solve it (discontinuous, hdg, etc). This method is called just after InitIndices such that the array is correctly filled.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();


GetThresholdMatrix

Syntax

 void GetThresholdMatrix() const

This method returns the threshold used to drop values in the finite element matrix. All values below this threshold (in magnitude) are dropped. It corresponds to the value given by the field ThresholdMatrix in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can call AddMatrixWithBC
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
Matrix<Real_wp, Symmetric, ArrayRowSymSparse> A;
// which threshold is set in the data file ?
cout << "Threshold for matrix = " << var.GetThresholdMatrix() << endl;



SetThresholdMatrix

Syntax

 void SetThresholdMatrix(Real_wp eps)

This method sets the threshold used to drop values in the finite element matrix. All values below this threshold (in magnitude) are dropped. It corresponds to the value given by the field ThresholdMatrix in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can call AddMatrixWithBC
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
Matrix<Real_wp, Symmetric, ArrayRowSymSparse> A;
// if you want a sparser matrix, you can drop small values
// (if the threshold is too large, the solution will be less accurate
var.SetThresholdMatrix(1e-12);



GetNbElt

Syntax

 int GetNbElt() const

This method returns the number of elements of the mesh (number of faces in 2-D, number of volumes in 3-D).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

// number of elements
int nb_elt = var.GetNbElt();



GetNbRows

Syntax

 int GetNbRows() const

This method returns the number of rows of the finite element matrix.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

// size of matrix
int nb_rows = var.GetNbRows();



GetPrintLevel

Syntax

 int GetPrintLevel() const /td>

This method returns the verbosity level used in Montjoie. Higher values induces more informations to be displayed during the simulation. A print level negative or null implies that Montjoie is silencious (no messages are displayed). It corresponds to the field PrintLevel in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// print level ?
int print_level = var.GetPrintLevel();


ComputeElementaryMatrix

Syntax

 void ComputeElementaryMatrix(int i, IVect& num_ddl, VirtualMatrix& A, CondensationBlockSolver_Base& solver, const GlobalGenericMatrix& nat_mat) const

This method computes the elementary matrix of a given element. Actually, this method is overloaded for each equation solved in Montjoie. If you want to compute the global finite element matrix, AddMatrixWithBC should be called.

Parameters

i (in)
element number
num_ddl (out)
global row numbers
A (out)
elementary matrix
solver (out)
solver handling static condensation
nat_mat (in)
mass, damping and stiffness coefficients

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Complex_wp> nat_mat;

// if you need only an elementary matrix
// (for the global matrix, call AddMatrixWithBC)
int num_elem = 5;
IVect num_row; Matrix<Complex_wp> mat_elem;
CondensationBlockSolver_Base<Complex_wp> cond_solver;
var.ComputeElementaryMatrix(num_elem, num_row, mat_elem, cond_solver, nat_mat);



GetInternalNodesElement

Syntax

 void GetInternalNodesElement(int i, int nb_dof_loc, int& nb_dof_edges, int& nb_dof_int, IVect& intern_node) const

This method is used to know which rows of the elementary matrix will be condensed.

Parameters

i (in)
element number
nb_dof_loc (in)
number of degrees of freedom of element i
nb_dof_edges (out)
number of degrees of freedom that cannot be eliminated
nb_dof_int (out)
number of degrees of freedom that will be condensed
intern_node (in)
If intern_node(i) is nonnegative, it is the local number among rows that are conserved after condensation. If intern_node(i) is negative, -intern_node(i)-1 is the local number among rows that are eliminated after condensation.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Complex_wp> nat_mat;

// if you need only an elementary matrix
// (for the global matrix, call AddMatrixWithBC)
int num_elem = 5;
IVect num_row; Matrix<Complex_wp> mat_elem;
CondensationBlockSolver_Base<Complex_wp> cond_solver;
var.ComputeElementaryMatrix(num_elem, num_row, mat_elem, cond_solver, nat_mat);

// if you want to know which rows can be eliminated
int nb_dof_edges, nb_dof_inside; Vector<int> intern_node;
var.GetInternalNodesElement(i, mat_elem.GetM(), nb_dof_edges, nb_dof_inside, intern_node);


GetNewCondensationSolver

Syntax

 CondensationBlockSolver_Base* GetNewCondensationSolver(Real_wp) CondensationBlockSolver_Base* GetNewCondensationSolver(Complex_wp)

This method constructs a new object handling static condensation. This new object is suited for the solved equation.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Complex_wp> nat_mat;

// if you want an adapted condensation solver
CondensationBlockSolver_Base<Complex_wp>* cond_solver;
cond_solver = var.GetNewCondensationSolver(Complex_wp(0));

// if you need only an elementary matrix
// (for the global matrix, call AddMatrixWithBC)
int num_elem = 5;
IVect num_row; Matrix<Complex_wp> mat_elem;
var.ComputeElementaryMatrix(num_elem, num_row, mat_elem, *cond_solver, nat_mat);

// the created object must be deleted after use
delete cond_solver;


UseMatrixFreeAlgorithm

Syntax

 bool UseMatrixFreeAlgorithm() const

This method returns true if the finite element matrix will not be stored (if an iterative matrix is constructed).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

bool matrix_stored = var.UseMatrixFreeAlgorithm(); // matrix will be stored ?

// computation of the iterative matrix
FemMatrixFreeClass_Base<Complex_wp>* A = var.GetNewIterativeMatrix(Complex_wp(0));
GlobalGenericMatrix<Complex_wp> nat_mat;

delete A;


IsSymmetricGlobalMatrix

Syntax

 bool IsSymmetricGlobalMatrix() const

This method returns true if the global finite element matrix is symmetric.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

bool sym = var.IsSymmetricGlobalMatrix(); // matrix is symmetric ?

GlobalGenericMatrix<Complex_wp> nat_mat;
if (sym)
{
DistributedMatrix<Complex_wp, Symmetric, ArrayRowSymSparse> A;
}
else
{
DistributedMatrix<Complex_wp, General, ArrayRowSparse> A;
}



GetStorageFiniteElementMatrix

Syntax

 int GetStorageFiniteElementMatrix() const

This method returns the storage used for the computation of the iterative matrix. The following storages are possible :

• MATRIX_AUTO_STORAGE : Montjoie selects automatically if the matrix is stored or not
• MATRIX_STORED : the iterative matrix is stored
• MATRIX_FREE : the iterative matrix is not stored

This storage is chosen by inserting a line with the field ExplicitMatrixFEM in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// storage for the iterative matrix
int storage = var.GetStorageFiniteElementMatrix();



SetStorageFiniteElementMatrix

Syntax

 void SetStorageFiniteElementMatrix(int type)

This method sets the storage used for the computation of the iterative matrix. The following storages are possible :

• MATRIX_AUTO_STORAGE : Montjoie selects automatically if the matrix is stored or not
• MATRIX_STORED : the iterative matrix is stored
• MATRIX_FREE : the iterative matrix is not stored

This storage can also be chosen by inserting a line with the field ExplicitMatrixFEM in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// if you want to change the storage of the iterative matrix
var.SetStorageFiniteElementMatrix(var.MATRIX_STORED);

// computation of the iterative matrix
FemMatrixFreeClass_Base<Complex_wp>* A = var.GetNewIterativeMatrix(Complex_wp(0));
GlobalGenericMatrix<Complex_wp> nat_mat;

delete A;



SetSymmetricElementaryMatrix

Syntax

 void SetSymmetricElementaryMatrix(bool sym=true)

This method can be used to modify the type of elementary matrix.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// if you want to change the type of symmetry for the elementary matrix
// it is not needed to set it, since the correct type is usually automatically selected
var.SetSymmetricElementaryMatrix(false); // unsymmetric elementary matrix



SetLeafStaticCondensation

Syntax

 void SetLeafStaticCondensation(bool condensed=true)

This method is used to tell Montjoie that the static condensation has to been performed. If condensed is true (and if static condensation has been enabled), the condensed matrix will be computed if AddMatrixWithBC is called, otherwise the non-condensed matrix is computed. This method does not enable or disable static condensation, it is done by inserting a line StaticCondensation in the data file. Thanks to the method SetLeafStaticCondensation, the user can compute the non-condensed matrix even though the static condensation has been enabled (to compute efficiently the solution).

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// the mesh and finite elements are constructed
var.ComputeMeshAndFiniteElement(type_element);

// then other initializations (if needed)
var.PerformOtherInitializations();

var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// to tell that we want to compute the non-condensed matrix
var.SetLeafStaticCondensation(false);

// computation of the matrix
GlobalGenericMatrix<Complex_wp> nat_mat;
DistributedMatrix<Complex_wp, General, ArrayRowSparse> A;



GetLeafStaticCondensation

Syntax

 bool GetLeafStaticCondensation() const

This method returns true if the static condensation has to been effectively performed. More details are given in the description of SetLeafStaticCondensation.

LightStaticCondensation

Syntax

 bool LightStaticCondensation() const

This method returns true if a light static condensation is enabled. A light static condensation consists of removing discontinuous unknowns (labelled as vectorial unknowns) whereas continuous unknowns (labelled as scalar unknowns) are conserved. This method makes sense only for a continuous formulation of the equation. Internal degrees of freedom are not eliminated as done by a regular static condensation. A light static condensation is enabled if the user inserts the following line in the data file :

   TypeCondensation = Light


Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// light static condensation ?
bool light = var.LightStaticCondensation();


GetSymmetrizationUse

Syntax

 bool GetSymmetrizationUse() const

This method returns true if a symmetrization is used (if possible). For example, Helmholtz equation (in the mixed formulation) given as

can be symmetrized by multiplying the second equation by -1. With this modification, the finite element matrix is symmetric. This symmetrization is possible for specific equations.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// use of a symmetrization ?
bool use_sym = var.GetSymmetrizationUse();


SetSymmetrizationUse

Syntax

 void SetSymmetrizationUse(bool sym=true) const

This method can be used if the user wants to use a symmetrization (if possible). For example, Helmholtz equation (in the mixed formulation) given as

can be symmetrized by multiplying the second equation by -1. With this modification, the finite element matrix is symmetric. This symmetrization is possible for specific equations. It is advised to check that the finite element matrix will be symmetric by calling IsSymmetricProblem.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// we want to symmetrize the matrix if possible
var.SetSymmetrizationUse();

// checking that the matrix will be symmetric
bool sym = var.IsSymmetricProblem();


SetHomogeneousDirichlet

Syntax

 void SetHomogeneousDirichlet(bool hg_dir) const

This method informs Montjoie that all the Dirichlet conditions are homogeneous (i.e. u = 0) if hg_dir is true. If a Dirichlet condition is hetereogeneous (i.e. u = f), hg_dir should be set to false.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// if you are sure that all Dirichlet conditions are homogeneous (u = 0)
// then you can say it :
var.SetHomogeneousDirichlet(true);

// it will save time for the computation of the solution


IsHomogeneousDirichlet

Syntax

 bool IsHomogeneousDirichlet() const

This method returns true if only Homogeneous Dirichlet conditions (i.e. u = 0) are present.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// if you are sure that all Dirichlet conditions are homogeneous (u = 0)
// then you can say it :
var.SetHomogeneousDirichlet(true);

// it will save time for the computation of the solution
// then IsHomogeneousDirichlet should return true
bool hg = var.IsHomogeneousDirichlet();


SetPrintLevel

Syntax

 void SetPrintLevel(int level) const /td>

This method sets the verbosity level used in Montjoie. Higher values induces more informations to be displayed during the simulation. A print level negative or null implies that Montjoie is silencious (no messages are displayed). It corresponds to the field PrintLevel in the data file.

Example :

    EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

// print level ?
int print_level = var.GetPrintLevel();

// if you want to change it (for example to keep Montjoie silent for a specific computation
var.SetPrintLevel(0);

// and you can get back to the previous print level
var.SetPrintLevel(print_level);


Syntax

 void AddMatrixWithBC( Matrix& A, const GlobalGenericMatrix& nat_mat, int offset_row = 0, int offset_col = 0, CondensationBlockSolver_Fem* solver = NULL, bool diag_matrix = false) const

This method computes the finite element matrix (with boundary terms). It works if A is an iterative matrix or if A is a distributed matrix. In this last case, the finite element matrix is added to the previously matrix given as argument. The object nat_mat stores the mass coefficient α, damping coefficient β, stiffness coefficient γ such that the following matrix is computed :

For time-harmonic equations, these coefficients are usually set to one. For unsteady equations, these coefficients depend on the used time scheme.

Parameters

A (inout)
finite element matrix
nat_mat (in)
mass, damping and stiffness coefficients
offset_row (optional)
offset for row numbers
offset_col (optional)
offset for column numbers
solver (optional)
solver handling static condensation
diag_matrix (optional)
true if only the diagonal of the matrix has to be computed

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can call AddMatrixWithBC
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
Matrix<Complex_wp, Symmetric, ArrayRowSymSparse> A;
var.AddMatrixWithBC(A, nat_mat); // A is a sparse matrix => it is stored

// if you only need an iterative matrix (the matrix is not necessarily stored)
// ie only the matrix-vector product is needed
FemMatrixFreeClass_Base<Complex_wp>* Aiter = var.GetNewIterativeMatrix(Complex_wp(0));

delete Aiter;


Syntax

 void AddMatrixFEM( Matrix& A, const GlobalGenericMatrix& nat_mat, int offset_row = 0, int offset_col = 0, CondensationBlockSolver_Fem* solver = NULL, bool diag_matrix = false) const

This method computes the finite element matrix (without boundary terms). It works if A is an iterative matrix or if A is a distributed matrix. In this last case, the finite element matrix is added to the previously matrix given as argument. The object nat_mat stores the mass coefficient α, damping coefficient β, stiffness coefficient γ such that the following matrix is computed :

For time-harmonic equations, these coefficients are usually set to one. For unsteady equations, these coefficients depend on the used time scheme.

Parameters

A (inout)
finite element matrix
nat_mat (in)
mass, damping and stiffness coefficients
offset_row (optional)
offset for row numbers
offset_col (optional)
offset for column numbers
solver (optional)
solver handling static condensation
diag_matrix (optional)
true if only the diagonal of the matrix has to be computed

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can call AddMatrixFEM
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
Matrix<Complex_wp, Symmetric, ArrayRowSymSparse> A;
var.AddMatrixFEM(A, nat_mat); // A is a sparse matrix => it is stored
// here only volume integrals are computed (no boundary conditions)

// if you only need an iterative matrix (the matrix is not necessarily stored)
// ie only the matrix-vector product is needed
FemMatrixFreeClass_Base<Complex_wp>* Aiter = var.GetNewIterativeMatrix(Complex_wp(0));

delete Aiter;


ComputeDiagonalMatrix

Syntax

 void ComputeDiagonalMatrix( Vector& diagonal, const GlobalGenericMatrix& nat_mat, bool assemble = true); void ComputeDiagonalMatrix( Vector& diagonal, Matrix& A, const GlobalGenericMatrix& nat_mat, bool assemble = true);

This method computes the diagonal of the finite element matrix. In the first syntax, only mass, damping and stiffness coefficients are provided. In the second syntax, the computed matrix is provided (iterative or distributed).

Parameters

diagonal (out)
diagonal of the matrix
A (optional)
finite element matrix previously computed
nat_mat (in)
mass, damping and stiffness coefficients
assemble (optional)
if true, the diagonal is assembled (between processors)

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// once var is constructed, you can compute the diagonal of the matrix
GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to set coefficients alpha, beta and gamma
// By default, alpha = beta = gamma = 1
VectComplex_wp diagonal;
var.ComputeDiagonalMatrix(diagonal, nat_mat);


IsSymmetricElementaryMatrix

Syntax

 bool IsSymmetricElementaryMatrix( const GlobalGenericMatrix& nat_mat) const;

This method returns true if the elementary matrix (with coefficients given in nat_mat) is symmetric. This method is called to select the correct type of elementary matrix when the matrix is computed with AddMatrixWithBC.

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// is the elementary matrix symmetric ?
GlobalGenericMatrix<Complex_wp> nat_mat;
nat_mat.SetCoefDamping(Complex_wp(0)); // you can change some coefficients
bool sym = var.IsSymmetricElementaryMatrix(nat_mat);



IsDiagonalElementaryMatrix

Syntax

 bool IsDiagonalElementaryMatrix( const GlobalGenericMatrix& nat_mat) const;

This method returns true if the elementary matrix (with coefficients given in nat_mat) is symmetric. This method is called to select the correct type of elementary matrix when the matrix is computed with AddMatrixWithBC. For example, it will return true if you ask to compute the mass matrix with mass-lumped elements (for most of equations).

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// is the elementary matrix diagonal ?
GlobalGenericMatrix<Complex_wp> nat_mat;
nat_mat.SetCoefDamping(Complex_wp(0)); // you can change some coefficients
nat_mat.SetCoefStiffness(Complex_wp(0));
// alpha = 1, beta = gamma = 0 => mass matrix is asked
bool diag = var.IsDiagonalElementaryMatrix(nat_mat);



IsSparseElementaryMatrix

Syntax

 bool IsSparseElementaryMatrix( const GlobalGenericMatrix& nat_mat) const;

This method returns true if the elementary matrix (with coefficients given in nat_mat) is sparse. This method is called to select the correct type of elementary matrix when the matrix is computed with AddMatrixWithBC.

Example :

   EllipticProblem<HelmholtzEquationDG<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// is the elementary matrix sparse ?
GlobalGenericMatrix<Complex_wp> nat_mat;
nat_mat.SetCoefDamping(Complex_wp(0)); // you can change some coefficients
nat_mat.SetCoefStiffness(Complex_wp(0));
// alpha = 1, beta = gamma = 0 => mass matrix is asked
bool sparse = var.IsSparseElementaryMatrix(nat_mat);



GetStaticCondensedRows

Syntax

 void GetStaticCondensedRows( IVect& IndexCondensedRows, IVect& global_row, IVect& overlap_row, IVect& overlap_proc, int& nb_scalar_dof, int& nb_global_dof, IVect& sharing_procs, Vector& sharing_rows) const

This method computes the row numbers of degrees of freedom after a static condensation. It is used to construct a distributed matrix with only condensed degrees of freedom (hence, the matrix has a small size).

Parameters

IndexCondensedRows (inout)
IndexCondensedRows(i) is equal to -1 if the dof is eliminated, and to index if the dof is kept. The condensed matrix will be a smaller matrix using "index" as row/column numbers.
global_row (out)
global row numbers of the condensed matrix
overlap_row (out)
dofs that are owned by another processor
overlap_proc (out)
processor that owns the overlapped dofs
nb_scalar_dof (out)
number of dofs per unknown for the condensed matrix
nb_global_dof (out)
global size of the condensed matrix
sharing_procs (out)
numbers of processors that share dofs with current processor
sharing_rows (out)
row numbers of shared dofs

Example :

   EllipticProblem<HelmholtzEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// to construct a distributed matrix with condensed dofs :
IVect IndexCondensedRows, global_row, overlap_row, overlap_proc, sharing_procs;
int nb_scalar_dof, nb_global_dof; Vector<IVect> sharing_rows;
var.GetStaticCondensedRows(IndexCondensedRows, global_row, overlap_row, overlap_proc,
nb_scalar_dof, nb_global_dof, sharing_procs, sharing_rows);

int nb_u = 1;
if (var.GetNbMeshNumberings() == 1)
nb_u = var.nb_unknowns_scal;

DistributedMatrix<Complex_wp, General, ArrayRowSparse> A;
A.Init(nb_global_dof, global_row, overlap_row, overlap_proc, nb_scalar_dof, nb_u,
sharing_procs, sharing_rows, var.comm_group_mode);


Syntax

 void AddElementaryFluxesDG( Matrix& A, const GlobalGenericMatrix& nat_mat, int offset_row = 0, int offset_col = 0)

This method adds to the given matrix A the terms due to numerical fluxes (for a Discontinuous Galerkin formulation). It is called by AddMatrixWithBC.

Parameters

A (inout)
matrix to modify
nat_mat (in)
mass, damping and stiffness coefficients
offset_row (optional)
offset for row numbers
offset_col (optional)
offset for column numbers

Example :

   EllipticProblem<HelmholtzEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// matrix with only numerical fluxes
GlobalGenericMatrix<Complex_wp> nat_mat;
DistributedMatrix<Complex_wp, General, ArrayRowSparse> A;
A.Reallocate(var.GetNbDof(), var.GetNbDof());


Syntax

This method updates the shifts (for the computation of eigenvalues) because of adimensionalization. The method modifies the shifts such that they are adimensionalized. This methods modifies the shifts only if there is a line Adimensionalization = YES in the data file.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// physical shift (square of pulsation) for second-order formulation
Real_wp freq = 4000; // frequency in Hz
Real_wp omega = 2*pi*freq;
Real_wp sr = square(omega), si = 0;

// we want to know adimensionalized shifts (which correspond to the computed finite element matrix which is adimensionalized)



Syntax

 void UpdateEigenvaluesAdimensionalization( Vector& Lr, Vector& Li, Matrix& V)

This method updates the eigenvalues (stored in Lr and Li) and the eigenvectors (stored in V) . The obtained eigenvalues are the physical ones. This methods modifies the eigenvalues only if there is a line Adimensionalization = YES in the data file.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// physical shift (square of pulsation) for second-order formulation
Real_wp freq = 4000; // frequency in Hz
Real_wp omega = 2*pi*freq;
Real_wp sr = square(omega), si = 0;

// we want to know adimensionalized shifts (which correspond to the computed finite element matrix which is adimensionalized)

DistributedMatrix<Real_wp, Symmetric, ArrayRowSymSparse> Mh, Kh;
GlobalGenericMatrix<Real_wp> nat_mat;
nat_mat.SetCoefDamping(0.0); nat_mat.SetCoefStiffness(0.0);
nat_mat.SetCoefMass(0.0); nat_mat.SetCoefStiffness(1.0);

// we compute eigenvalues and eigenvectors
SparseEigenProblem<Real_wp, DistributedMatrix<Real_wp, Symmetric, ArrayRowSymSparse>,
DistributedMatrix<Real_wp, Symmetric, ArrayRowSymSparse> > eigen_solver;

eigen_solver.InitMatrix(Kh, Mh);
eigen_solver.SetTypeSpectrum(eigen_solver.CENTERED_EIGENVALUES, sr, eigen_solver.SORTED_MODULUS);

VectReal_wp Lr, Li; Matrix<Real_wp, General, ColMajor> V;
GetEigenvaluesEigenvectors(eigen_solver, Lr, Li, V);

// then we get back to physical eigenvalues


FindIntervalDofSignSymmetry

Syntax

 void FindIntervalDofSignSymmetry( int& i0, int& i1, int& j0, int& j1) const

This method is relevant only if a symmetrization is enabled (see SetSymmetrizationUse). It provides the rows that are multiplied by -1 (to obtain a symmetric matrix). The rows between i0 and i1 (i1 being exlcuded) and the rows between j0 and j1 are multiplied by -1.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// rows modified because of symmetrization ? (mixed formulation)
int i0=0, i1=0, j0=0, j1=0;
var.FindIntervalDofSignSymmetry(i0, i1, j0, j1);


ModifySourceSymmetry

Syntax

 void ModifySourceSymmetry( Vector& rhs) const void ModifySourceSymmetry( Matrix& rhs) const

This method is relevant only if a symmetrization is enabled (see SetSymmetrizationUse). It modifies the right hand side by multiplying some rows by -1 (some rows of the finite element matrix have been multiplied by -1 to obtain a symmetric matrix). A matrix can be provided in the case of multiplie right hand sides.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

// filling a right hand side
VectComplex_wp rhs(var.GetNbDof()); rhs.FillRand();

// rows affected by the symmetrization are multiplied by -1
var.ModifySourceSymmetry(rhs);


GetNewIterativeMatrix

Syntax

 FemMatrixFreeClass_Base* GetNewIterativeMatrix( T& s) const

This method constructs a new object (representing the iterative matrix) and returns the address of the created object. The object is usually an instance of FemMatrixFreeClass (that depends on the solved equation). The method GetNewIterativeMatrix returns a pointer of the base class FemMatrixFreeClass_Base such that it can be used in generic functions.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// do not forget to release the memory when Ah is no longer needed
delete Ah;


GetNewLinearSolver

Syntax

 All_LinearSolver* GetNewLinearSolver() const

This method constructs a new object (representing the linear solver) and returns the address of the created object. This solver can be used to compute the solution of a linear system with the finite element matrix.

Example :

   EllipticProblem<HarmonicElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// creation of the linear solver
All_LinearSolver* solver = var.GetNewLinearSolver();

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

// and solve the linear system A x = b
VectComplex_wp b(var.GetNbDof()), x; b.FillRand();
x = b; solver->ComputeSolution(x);

// when you no longer need the solver, you can release the memory
delete solver;


GetNewPreconditioning

Syntax

 All_Preconditioner_Base* GetNewPreconditioning(T)

This method constructs a new object (representing the preconditioning) and returns the address of the created object. This preconditioning can be used to compute efficiently the solution of a linear system with the finite element matrix.

Example :

   EllipticProblem<HarmonicElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

FemMatrixFreeClass_Base<Complex_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// creation of the preconditioning
All_Preconditioner_Base<Complex_wp> prec = var.GetNewPreconditioning();

// to apply the preconditioning to a vector
VectComplex_wp r(var.GetNbDof()), z(var.GetNbDof()); r.FillRand();
prec->Solve(*Ah, r, z); // z = M^{-1} r

// when you no longer need the preconditioning, you can release the memory
delete prec;


ComputeElementaryMatrix

Syntax

 void ComputeElementaryMatrix( int i, IVect& num_ddl, VirtualMatrix& A, const GlobalGenericMatrix& nat_mat, const EllipticProblem& var, const ElementReference& Fb) const

This generic function computes the elementary matrix of a given element. It is implemented if all finite elements are scalar. It is also implemented if all finite elements are 2-D edge elements. This function has to be distinguished from the method ComputeElementaryMatrix of the class EllipticProblem. Most of the times, the method ComputeElementaryMatrix will call the generic function ComputeElementaryMatrix (except for specific equations).

Parameters

i (in)
element number
num_ddl (out)
global row numbers
A (out)
elementary matrix
nat_mat (in)
mass, damping and stiffness coefficients
var (in)
object describing the problem to solve
Fb (in)
finite element associated with the element i

Example :

   // to define your own equation
class MyEquation : GenericEquation_Base<Complex_wp>
{
public:
typedef Dimension3 Dimension;

static const bool FirstOrderFormulation = true;

enum {nb_unknowns = 3, nb_unknowns_hdg=0,
nb_components_en = 1, nb_components_hn = 1,
nb_unknowns_scal = 3, nb_unknowns_vec = 0};

static inline bool SymmetricGlobalMatrix() { return false; }
static inline bool SymmetricElementaryMatrix() { return false; }

// other functions ( such as GetTensorMass, GetGradPhiTensor, etc)
};

// in the specialization of EllipticProblem
// you can call the generic function ComputeElementaryMatrix
template<>
class EllipticProblem<MyEquation>
: public VarHarmonic <MyEquation>
{
public :
void ComputeElementaryMatrix(int i, IVect& num_dof, VirtualMatrix<Complex_wp>& mat_elem,
CondensationBlockSolver_Base<Complex_wp>&,
const GlobalGenericMatrix<Complex_wp>& nat_mat)
{
// we call generic function
Montjoie::ComputeElementaryMatrix(i, num_dof, mat_elem, nat_mat, *this, this->GetReferenceElementH1(i));
}

};



Syntax

 void AddElementaryFluxesDG(VirtualMatrix& A, const GlobalGenericMatrix& nat_mat, const EllipticProblem& var, int offset_row, int offset_col) const

This generic function computes the boundary terms associated with numerical fluxes (for a Discontinuous Galerkin formulation) and adds them to a given matrix. It is implemented if all finite elements are scalar. This function has to be distinguished from the method AddElementaryFluxesDG of the class EllipticProblem. Most of the times, the method AddElementaryFluxesDG will call the generic function AddElementaryFluxesDG (except for specific equations).

Parameters

A (out)
elementary matrix
nat_mat (in)
mass, damping and stiffness coefficients
var (in)
object describing the problem to solve
offset_row (in)
offset for row numbers
offset_col (in)
offset for column numbers

Example :

    // to define your own equation
class MyEquation : GenericEquation_Base<Complex_wp>
{
public:
typedef Dimension3 Dimension;

static const bool FirstOrderFormulation = true;

enum {nb_unknowns = 3, nb_unknowns_hdg=0,
nb_components_en = 1, nb_components_hn = 1,
nb_unknowns_scal = 3, nb_unknowns_vec = 0};

static inline bool SymmetricGlobalMatrix() { return false; }
static inline bool SymmetricElementaryMatrix() { return false; }

// other functions ( such as GetTensorMass, GetGradPhiTensor, etc)
};

// in the specialization of EllipticProblem
// you can call the generic function ComputeElementaryMatrix and ElementaryFluxesDG
template<>
class EllipticProblem<MyEquation>
: public VarHarmonic <MyEquation>
{
public :
void ComputeElementaryMatrix(int i, IVect& num_dof, VirtualMatrix<Complex_wp>& mat_elem,
CondensationBlockSolver_Base<Complex_wp>&,
const GlobalGenericMatrix<Complex_wp>& nat_mat)
{
// we call generic function
Montjoie::ComputeElementaryMatrix(i, num_dof, mat_elem, nat_mat, *this, this->GetReferenceElementH1(i));
}

const GlobalGenericMatrix<Complex_wp>& nat_mat,
int offset_row, int offset_col)
{
}

};



mat_boundary_sym, mat_boundary_unsym

The attribute mat_boundary_sym stores the terms due to boundary conditions (or other models) in the case where the finite element matrix is symmetric. If the matrix is unsymmetric, the attribute mat_boundary_unsym is used.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// you can write terms due to boundary conditions
Ah->mat_boundary_sym.WriteText("mat_boundary.dat");

// do not forget to release the memory when Ah is no longer needed
delete Ah;


matCSR_boundary_sym, matCSR_boundary_unsym

The attribute matCSR_boundary_sym stores the terms due to boundary conditions (or other models) in the case where the finite element matrix is symmetric. If the matrix is unsymmetric, the attribute matCSR_boundary_unsym is used. These matrices are constructed if the method CompressMatrix has been called previously.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// matrices are compressed (use of CSR storage)
Ah->CompressMatrix();

// you can write terms due to boundary conditions
Ah->matCSR_boundary_sym.WriteText("mat_boundary.dat");

// do not forget to release the memory when Ah is no longer needed
delete Ah;


mat_iterative_sym, mat_iterative_unsym

The attribute mat_iterative_sym stores the finite element matrix if it is stored (and symmetric). If the matrix is unsymmetric, the attribute mat_iterative_unsym is used. To force the storage of the matrix, you can insert a line ExplicitMatrixFEM = YES in the data file.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// you can write the matrix (if stored)
Ah->mat_iterative_sym.WriteText("mat_iterative.dat");

// do not forget to release the memory when Ah is no longer needed
delete Ah;


matCSR_iterative_sym, matCSR_iterative_unsym

The attribute matCSR_iterative_sym stores the finite element matrix if it is stored (and symmetric). If the matrix is unsymmetric, the attribute matCSR_iterative_unsym is used. To force the storage of the matrix, you can insert a line ExplicitMatrixFEM = YES in the data file. These matrices are constructed if the method CompressMatrix has been called previously.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// using CSR storage (to save memory)
Ah->CompressMatrix();

// you can write the matrix (if stored)
Ah->matCSR_iterative_sym.WriteText("mat_iterative.dat");

// do not forget to release the memory when Ah is no longer needed
delete Ah;


var

This attribute stores a reference to the instance EllipticProblem.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> gt; Ah(var);

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah.MltVector(x, y); // y = A x

// var is present in the object Ah (here var and Ah.var refer to the same object)
Ah.var.SetPrintLevel(3);


Constructor of FemMatrixFreeClass

The only constructor takes an object EllipticProblem as argument. If you only have a base class of EllipticProblem (with a generic function that works for different equations), you can use the method GetNewIterativeMatrix.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> gt; Ah(var);

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah.MltVector(x, y); // y = A x

// if you are in a function with a base class of EllipticProblem :
FemMatrixFreeClass_Base<Real_wp>* Ah0 = var.GetNewIterativeMatrix(Real_wp(0));

delete Ah0;


SetCoefficientDirichlet

Syntax

 void SetCoefficientDirichlet(const Real_wp& coef) const

This method changes the coefficient for the diagonal of Dirichlet dofs.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// computation of the finite element matrix
Matrix<Real_wp, General, ArrayRowSparse> A;

// we want void columns for Dirichlet => coef = 0
// default value is 1.0
var.SetCoefficientDirichlet(Real_wp(0));



SetCoefficientMatrix

Syntax

 void SetCoefficientMatrix(const GlobalGenericMatrix& coefs) const

This method changes the mass, damping and stiffness coefficient for the iterative matrix. Usually it is not needed to call this method since AddMatrixWithBC already initializes these coefficients.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// if you want to change the coefficient afterwards (not recommended because boundary terms do not change)
nat_mat.SetCoefficientStiffness(2.0);
Ah.SetCoefficientMatrix(nat_mat);


GetCoefMass

Syntax

 T GetCoefMass() const

This method returns the mass coefficient for the iterative matrix.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// should return 2.0
Real_wp mass = Ah.GetCoefMass();


IsSymmetric

Syntax

 bool IsSymmetric() const

This method returns true if the matrix is symmetric.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// symmetric matrix ?
bool sym = Ah.IsSymmetric();


FormulationDG

Syntax

 int FormulationDG() const

This method returns the type of the formulation used. It can be equal to ElementReference_Base::CONTINUOUS, ElementReference_Base::DISCONTINUOUS or ElementReference_Base::HDG as explained in the description of dg_formulation.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// DG formulation ?
int dg_form = Ah.FormulationDG();


SetCondensedSolver

Syntax

 void SetCondensedSolver(CondensationBlockSolver_Fem*)

This method sets the condensed solver used to assemble the matrix. This method should not be used in regular use (since it is already called when AddMatrixBC is called).

DirichletDofIgnored

Syntax

 bool DirichletDofIgnored() const

This method returns true if Dirichlet dofs are ignored. If Dirichlet dofs are present (= not ignored), we set yi = 0 for rows associated with Dirichlet condition (where y = A x is the result of the matrix-vector product).

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// Dirichlet dofs treated ?
bool dirichlet = Ah.DirichletDofIgnored();


IgnoreDirichletDof

Syntax

 void IgnoreDirichletDof()

This method forces Dirichlet dofs to be ignored. If Dirichlet dofs are present (= not ignored), we set yi = 0 for rows associated with Dirichlet condition (where y = A x is the result of the matrix-vector product.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// if you want to ignore Dirichlet dofs
// (default : Dirichlet dofs are treated)
Ah.IgnoreDirichletDof();


SetScaling

Syntax

 void SetScaling(VectReal_wp& scale_left, VectReal_wp& scale_right)

This method sets the scaling factors to use if the matrix is not stored. This method is called when ScaleMatrix is used. ScaleMatrix handles both cases (matrix stored or not).

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// if you want to replace A by L A R (L and R are diagonal matrices)
// where L is the left scaling and R the right scaling :
VectReal_wp L(Ah.GetM()), R(Ah.GetM());
L.FillRand(); R.FillRand();

// best option : use ScaleMatrix (it will call SetScaling)
ScaleMatrix(Ah, L, R);


SucceedInAffectingPointer

Syntax

 bool SucceedInAffectingPointer(Matrix*& ptr_A, Matrix*& ptr_Acsr)

This method tries to initialize ptr_A (or ptr_Acsr if the matrix has been compressed) with the matrix stored in the object. If the method returns false, it means that no matrix has been found (pointers have not been initialized). In that case, the matrix is probably not stored, or the symmetry is not correct.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// for a symmetric matrix
Matrix<Real_wp, Symmetric, RowSymSparse>* A_csr;
Matrix<Real_wp, Symmetric, ArrayRowSymSparse>* A_stored;
bool success = Ah.SucceedInAffectingPointer(A_stored, A_csr);

// if success is false, it means that the matrix is either unsymmetric or not stored


InitSymmetricMatrix

Syntax

 void InitSymmetricMatrix()

This method inits the object to a symmetric matrix. This method should not be called in regular use, since it is already called by AddMatrixWithBC.

InitUnsymmetricMatrix

Syntax

 void InitUnsymmetricMatrix()

This method inits the object to an unsymmetric matrix. This method should not be called in regular use, since it is already called by AddMatrixWithBC.

ApplyRightScaling

Syntax

 void ApplyRightScaling(const Vector& B0, Vector& C0, Vector& B, Vector& C)

This method applies the right scaling to vector B0 (B = R B0 where R is the right scaling) and inits C = C0. If no right scaling is required, B and C are initialized with B0 and C0. This method is used in the method MltAddFree to handle scaling of the matrix.

Example :

     // if you specialize the class with your equation
class FemMatrixFreeClass<Real_wp, MyEquation> : public FemMatrixFreeClass_Eq<Real_wp, MyEquation>
{
public :
FemMatrixFreeClass(const EllipticProblem<MyEquation>& var_)
: FemMatrixFreeClass_Eq<Real_wp, MyEquation>(var_) {}

// product C0 = (L A R) B0
const SeldonTranspose&, int level,
const VectReal_wp& B0, VectReal_wp& C0) const
{
VectReal_wp B, C;
// right scaling for B = R B 0
this->ApplyRightScaling(B0, C0, B, C);

// you do the matrix vector product without scaling C = A B
// ... (code to insert)

// finally left scaling for C0 = L C
this->ApplyLeftScaling(B0, C0, B, C);
}

};


ApplyLeftScaling

Syntax

 void ApplyLeftScaling(const Vector& B0, Vector& C0, Vector& B, Vector& C)

This method applies the left scaling to vector C (C0 = L C where L is the left scaling). This method is used in the method MltAddFree to handle scaling of the matrix.

Example :

     // if you specialize the class with your equation
class FemMatrixFreeClass<Real_wp, MyEquation> : public FemMatrixFreeClass_Eq<Real_wp, MyEquation>
{
public :
FemMatrixFreeClass(const EllipticProblem<MyEquation>& var_)
: FemMatrixFreeClass_Eq<Real_wp, MyEquation>(var_) {}

// product C0 = (L A R) B0
const SeldonTranspose&, int level,
const VectReal_wp& B0, VectReal_wp& C0) const
{
VectReal_wp B, C;
// right scaling for B
this->ApplyRightScaling(B0, C0, B, C);

// you do the matrix vector product without scaling C = A B
// ... (code to insert)

// finally left scaling for C0 = L C
this->ApplyLeftScaling(B0, C0, B, C);
}

};


CompressMatrix

Syntax

 void CompressMatrix()

This method compresses the matrices stored in the object. It consists of converting them into CSR storage (e.g. RowSparse instead of ArrayRowSparse) which requires less memory. The matrix vector product is also more efficient.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

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

// for an iterative matrix (the matrix is not necessary stored, use FemMatrixFreeClass)
FemMatrixFreeClass_Base<Real_wp>* Ah = var.GetNewIterativeMatrix(Complex_wp(0));

// with Ah, you can only compute matrix-vector products
VectReal_wp x(var.GetNbDof()), y(var.GetNbDof());
x.FillRand();
Ah->MltVector(x, y); // y = A x

// matrices are compressed (use of CSR storage)
Ah->CompressMatrix();

// do not forget to release the memory when Ah is no longer needed
delete Ah;


Syntax

 void AddExtraBoundaryTerms(const T& alpha, const Vector& B, Vector& C)

This method adds the terms due to boundary conditions (or other models) contained in mat_boundary* to the vector C. It can be used to write a specialization of the class for your own equation.

Example :

     // if you specialize the class with your equation
class FemMatrixFreeClass<Real_wp, MyEquation> : public FemMatrixFreeClass_Eq<Real_wp, MyEquation>
{
public :
FemMatrixFreeClass(const EllipticProblem<MyEquation>& var_)
: FemMatrixFreeClass_Eq<Real_wp, MyEquation>(var_) {}

// product C0 = (L A R) B0
const SeldonTranspose&, int level,
const VectReal_wp& B0, VectReal_wp& C0) const
{
VectReal_wp B, C;
// right scaling for B
this->ApplyRightScaling(B0, C0, B, C);

// you do the matrix vector product without scaling C = A B
// ... (code to insert)

// contributions of boundary terms

// finally left scaling for C0 = L C
this->ApplyLeftScaling(B0, C0, B, C);
}

};


SetNbDirichletCondition

Syntax

 void SetNbDirichletCondition(int n)

This method sets the number of Dirichlet conditions, usually it corresponds to the number of right-hand-sides.

Example :

        EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// for a single right hand side
Ah.SetNbDirichletCondition(1);


ApplyDirichletCondition

Syntax

 void ApplyDirichletCondition(const SeldonTranspose& trans, Vector& b, int k=0)

This method modifies the right hand side if hetereogenous Dirichlet condition is present. The new right hand side is equal to

With this transformation, the hetereogeneous Dirichlet is converted into an homogeneous Dirichlet condition (with a symmetric finite element matrix). The values bk where k is a Dirichlet dof are stored in the object and used when ImposeDirichletCondition is called. If trans is equal to SeldonTranspose, we assume that the solution with the transpose of A is searched.

Example :

    EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// for a single right hand side
Ah.SetNbDirichletCondition(1);
VectReal_wp rhs(Ah.GetM()); rhs.FillRand();

// we modify the right hand side to convert Dirichlet conditions to homogeneous ones
Ah.ApplyDirichletCondition(SeldonNoTrans, rhs);

// then you have to solve A x = rhs
// to do
VectReal_wp x(rhs.GetM()); x.Zero(); // etc
Cg(Ah, x, rhs, id_precond, iter); // solving with Cg for instance

// to get back to the solution (with hetereogeneous Dirichlet)
Ah.ImposeDirichletCondition(SeldonNoTrans, x);


ImposeDirichletCondition

Syntax

 void ImposeDirichletCondition(const SeldonTranspose& trans, Vector& x, int k=0)

This method sets the values of x with Dirichlet conditions given when ApplyDirichletCondition has been called.

Example :

    EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// for two right hand sides
Ah.SetNbDirichletCondition(2);
VectReal_wp rhs(Ah.GetM()), rhs2(Ah.GetM()); rhs.FillRand(); rhs2.FillRand();

// we modify the right hand sides to convert Dirichlet conditions to homogeneous ones
Ah.ApplyDirichletCondition(SeldonNoTrans, rhs, 0);
Ah.ApplyDirichletCondition(SeldonNoTrans, rhs2, 1);

// then you have to solve A x = rhs
// to do
VectReal_wp x(rhs.GetM()), x2(rhs2.GetM()); x.Zero(); x2.Zero(); // etc
Cg(Ah, x, rhs, id_precond, iter); // solving with Cg for instance
Cg(Ah, x2, rhs2, id_precond, iter);

// to get back to the solution (with hetereogeneous Dirichlet)
Ah.ImposeDirichletCondition(SeldonNoTrans, x, 0);
Ah.ImposeDirichletCondition(SeldonNoTrans, x2, 1);


SetDirichletCondition

Syntax

 void SetDirichletCondition(Matrix& A, int offset_row=0, int offset_col=0)

This method modifies the matrix by removing rows (and/or columns) associated with Dirichlet dofs. For a symmetric matrix, the Dirichlet rows and columns are removed and replaced by a diagonal with a constant coefficient (usually one). This method does not need to be called, since it is already done in AddMatrixWithBC.

Syntax

 void MltAddHetereogeneousDirichlet(const T& alpha, const SeldonTranspose& trans, const Vector& B, Vector& C)

This method computes C = C + alpha A B with only columns a A associated with Dirichlet dofs.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

// product with only Dirichlet columns
VectReal_wp b(Ah.GetM()); b.FillRand();
VectReal_wp c(Ah.GetM()); c.Zero();


InitDirichletCondition

Syntax

 void InitDirichletCondition(const Vector& b, int k=0)

This method stores hetereogeneous Dirichlet condition contained in x. The values can then be retrieve by calling ImposeDirichletCondition.

Example :

   EllipticProblem<ElasticEquation<Dimension3> > var;

// we retrieve the type of finite element and equation in the data file
string type_element, type_equation;
getElement_Equation("test.ini", type_element, type_equation);

// allocating the arrays for physical indexes
var.InitIndices(100); // we choose N = 100
var.SetTypeEquation(type_equation);

// then you can read the data file

var.ComputeMeshAndFiniteElement(type_element); // mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations
var.ComputeMassMatrix(); // computation of geometric quantities (such as jacobian matrices)
var.ComputeQuasiPeriodicPhase(); // for quasi-periodic conditions

GlobalGenericMatrix<Real_wp> nat_mat; // this object is used to store mass, damping and stiffness coefficients
nat_mat.SetCoefMass(2.0); // chaning the mass coefficient

// computation of the finite element matrix
FemMatrixFreeClass<Real_wp, ElasticEquation<Dimension3> > Ah;

VectReal_wp x(Ah.GetM()); x.FillRand();

// if you want to store values contained in x on Dirichlet dofs :
Ah.InitDirichletCondition(x);

// and retrieve them
VectReal_wp y(Ah.GetM()); y.Zero();
Ah.ImposeDirichletCondition(SeldonNoTrans, y);



Syntax

 void MltAddFree(const GlobalGenericMatrix& nat_mat, const SeldonTranspose& trans, int level, const Vector& X, Vector& Y)

This method adds to Y the matrix-vector product A X where A is the matrix associated with the current object. If trans is equal to SeldonTrans, we add AT X. MltAddFree is usually not called directly, but rather MltVector (or MltAddVector). These two functions will call MltAddFree in the case where the matrix is not stored. You can specialize MltAddFree to perform the matrix-vector product with your own equation (matrix-free implementation).

Example :

     // if you specialize the class with your equation
class FemMatrixFreeClass<Real_wp, MyEquation> : public FemMatrixFreeClass_Eq<Real_wp, MyEquation>
{
public :
FemMatrixFreeClass(const EllipticProblem<MyEquation>& var_)
: FemMatrixFreeClass_Eq<Real_wp, MyEquation>(var_) {}

// product C0 = (L A R) B0
const SeldonTranspose&, int level,
const VectReal_wp& B0, VectReal_wp& C0) const
{
VectReal_wp B, C;
// right scaling for B = R B 0
this->ApplyRightScaling(B0, C0, B, C);

// you do the matrix vector product without scaling C = A B
// ... (code to insert)

// finally left scaling for C0 = L C
this->ApplyLeftScaling(B0, C0, B, C);
}

};


GetExtrapolVariables

Syntax

 ExtrapolVariablesProductFEM& GetExtrapolVariables()

This method returns the object storing additional variables needed to perform the matrix-vector product. You can also specialize the class ExtrapolVariablesProductFEM to store the variables needed to perform the matrix-vector product.

Example :

     // if you specialize the class with your equation
class FemMatrixFreeClass<Real_wp, MyEquation> : public FemMatrixFreeClass_Eq<Real_wp, MyEquation>
{
public :
FemMatrixFreeClass(const EllipticProblem<MyEquation>& var_)
: FemMatrixFreeClass_Eq<Real_wp, MyEquation>(var_) {}

// product C0 = (L A R) B0
const SeldonTranspose&, int level,
const VectReal_wp& B0, VectReal_wp& C0) const
{
VectReal_wp B, C;
// right scaling for B = R B 0
this->ApplyRightScaling(B0, C0, B, C);

// you do the matrix vector product without scaling C = A B
// ... (code to insert)
ExtrapolVariablesProductFEM<Real_wp, MyEquation>& var_extra = this->GetExtrapolVariables();

// finally left scaling for C0 = L C
this->ApplyLeftScaling(B0, C0, B, C);
}

};


mesh

This attribute is the mesh used to complete the simultation. It is usually modified through data file (by inserting a line FileMesh = file_name).

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed

// you can write the mesh on a file
var.mesh.Write("test.mesh");


Location :

This attribute stores the quadrature points on all the elements. These points are computed when the method ComputeMassMatrix is called.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
// (second argument is false to keep quadrature points)
// (if the second argument is true, the quadrature points are cleared)
var.ComputeMassMatrix(true, false);

// you can display quadrature points of a given element
int num_elem = 5;

// to retrieve one point :
int j = 7;


Location :

If the attribute write_quadrature_points is true, the quadrature points are written in the file quadrature_points.dat when the method ComputeMassMatrix is called. If the attribute write_quad_points_pml is true, the quadrature points of PMLs are written, otherwise only quadrature points of elements in the physical domain are written. These attributes are usually modified by inserting a line WriteQuadraturePoints = YES NO_PML_POINTS in the data file.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// if you want to force quadrature points to be written
var.write_quad_points_pml = true; // with PMLs

// quadrature points and jacobian matrices are computed (and written if asked)
var.ComputeMassMatrix(true, false);


Glob_jacobian

This attribute stores the following quantities :

An element is affine if the transformation Fi is affine. In 2-D, it corresponds to straight triangles and parallelograms.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

int num_elem = 5;
bool affine = var.mesh.IsElementAffine(num_elem);
if (affine)
{
// affine elements, the jacobian is stored (without weight)
Real_wp jacob = var.Glob_jacobian(num_elem)(0);
}
else
{
// weighted jacobian on quadrature points
int j = 4;
Real_wp omJi = var.Glob_jacobian(num_elem)(j);
// if you want to recover the jacobian on the quadrature point => divide by the weight
const ElementReference<Dimension2, 1>& Fb = var.GetReferenceElementH1(num_elem);
Real_wp om = Fb.WeightsND(j); Real_wp jacob = omJi / om;
}


Glob_decomp_jacobian

This attribute stores the polynomial coefficients of the jacobian. They are computed only for elements that have a sparse mass matrix for straight elements. It can be checked by calling the method LinearSparseMassMatrix of the finite element.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

// polynomial coefficients of the jacobian
int num_elem = 11;
const ElementReference<Dimension2, 1>& Fb = var.GetReferenceElementH1(num_elem);
if (Fb.LinearSparseMassMatrix())
cout << "coefficients = " << var.Glob_decomp_jacobian(num_elem) << endl;


Glob_normale

This attribute stores the normales on quadrature points of faces (edges in 2-D) of the mesh. They are computed when the method ComputeMassMatrix is called.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

// normales for a given face
int num_face = 11; int num_pt = 2;
R2 normale = var.Glob_normale(num_face)(num_pt);


Glob_dsj

This attribute stores the surface integration elements on quadrature points of faces (edges in 2-D) of the mesh. They are computed when the method ComputeMassMatrix is called.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

// surface jacobian (used to compute surface integrals)
int num_face = 11; int num_pt = 2;
Real_wp ds = var.Glob_dsj(num_face)(num_pt);


OrthogonalElement

This attribute stores the type of orthogonality of the considered element. It is computed only for elements inside PML elements. In 2-D, it can take the following values:

• 0 : the element has diagonal jacobian matrices (the quadrilateral is a rectangle)
• 1 : at least one jacobian matrix is not diagonal

In 3-D, it can take the following values:

• 0 : the element has diagonal jacobian matrices (the hexahedron is a rectangular parallelepiped)
• 1 : the element has jacobian matrices with zeros on last row and last column (except the diagonal coefficient). It corresponds to hexahedra that are obtained by extrusion in z-direction.
• 2 : other elements

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

// orthogonal element ?
int num_elem = 13;
int ortho = var.OrthogonalElement(num_elem);


Glob_DFjm1

This attribute stores the following quantities :

An element is affine if the transformation Fi is affine. In 2-D, it corresponds to straight triangles and parallelograms.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

int num_elem = 5;
bool affine = var.mesh.IsElementAffine(num_elem);
if (affine)
{
// affine elements, only one matrix is stored
Matrix2_2 jacob_dfjm1 = var.Glob_DFjm1(num_elem)(0);
}
else
{
// matrices J_i DF_i^{-1} on quadrature points
int j = 4;
Matrix2_2 jacob_dfm1 = var.Glob_DFjm1(num_elem)(j);
}


IsNewFace

This attributed is used to know if the normale is outgoing to the element. If IsNewFace(i)(j) is true, the normales stored in Glob_normale are outgoing to the element i.

Example :

    EllipticProblem<LaplaceEquation<Dimension2> > var;

var.InitIndices(100);
var.SetTypeEquation("LAPLACE");

// then the mesh and finite element are constructed
var.PerformOtherInitializations(); // other initializations

// quadrature points and jacobian matrices are computed
var.ComputeMassMatrix();

int num_elem = 5;
int num_loc = 2; // local boundary of element
bool new_face = var.IsNewFace(num_elem)(num_loc);
int num_face = var.mesh.Element(num_elem).numBoundary(num_loc);
int num_pt = 1;
R2 normale = var.Glob_normale(num_face)(num_pt);
// changing sign of normale to enforce an outgoing normale to element num_elem
if (!new_face)
normale = -normale;