Data files in Montjoie

TypeEquation

Syntax :

TypeEquation = name_equation

This line specifies which equation is solved by the code. Actually this field is used only for some targets. As you noticed, for most of targets, the equation is already uniquely associated with the name of the target.

Example :

# to use Discontinuous Galerking for acoustics :
TypeEquation = ACOUSTIC_DG

Location :

Adimensionalization

Syntax :

Adimensionalization = YES
Adimensionalization = NO
Adimensionalization = ONE

This line is used to specify how equations are adimensionalized. The avaible choices are the following ones:

NO : no adimensionalization is performed, the user must provide physical values for all parameters (time step, frequency, etc).
YES:an adimensionalization is performed, but there should not be any impact on the user. The user must provide physical values, these values are changed inside the code, and outputs are modified in order to contain physical values.
ONE : the adimensionalization is performed by the user who provides adimensionalized values. With this mode, the velocity and wavelength at infinity are usually assumed to be equal to 1.

The default adimensionalization mode is ONE.

Example :

Adimensionalization = YES

Location :

Seed

Syntax :

Seed = n

This line should be the first line of the input file if you want to specify a random seed that will be common to all processors. This seed is used for random Gaussian sources.

Example :

Seed = 103

Location :

NumberPhysicalMedia

Syntax :

NumberPhysicalMedia = N

The number N is the maximal number for the volume references in the mesh. In Gmsh, it would correspond to the maximum number of Physical Surface (in 2-D) and Physical Volume (in 3-D). The physical indexes (such as epsilon, mu, sigma, etc) are stored in an array whose dimension is equal to this number. The default value is 20.

Example :

# if you have a physical volume of reference 100 (maximum value):
NumberPhysicalMedia = 100

Location :

Timer

Syntax :

Timer = Name

You can select the timer to use for the measurement of computation time. You can choose between Accurate, Real or Basic. The default timer is Real if you whave compiled with the flag MONTJOIE_WITH_REAL_TIMING, Basic otherwise. Accurate timer is based on the C function getrusage, Real timer is based on C function clock_gettime while Basic timer is based on C function clock.

Example :

# if you want to use real timing
Timer = Real

Location :

Frequency

Syntax :

Frequency = freq omega

You can specify either a frequency freq or a pulsation omega. The resulting pulsation is equal to freq*2*pi + omega. In case of adimensionalization, this frequency is multiplied by L/c where L is a characteristic wavelength and c the speed of light.

Example :

# For a frequency equal to 1 (then pulsation = 2 pi)
Frequency = 1.0  0

# For a pulsation equal to 10 (then frequency = 10/(2 pi) )
Frequency = 0  10.0

Location :

PhysicalFrequency

Syntax :

PhysicalFrequency = freq

For Helmholtz and Maxwell equations, we usually perform an adimensionalization so that a frequency equal to 1 corresponds in the physical space to a frequency approximatively equal to 300 Mhz. Therefore in this field, you can write the physical frequency (in Hz), and the conversion to adimensionalized frequency will be performed by dividing by the speed of light in free space.

Example :

# For a physical frequency equal to 2.5 Ghz
PhysicalFrequency = 2.5e9

Location :

Wavelength, WaveLength

Syntax :

  WaveLength = lambda
  Wavelength = lambda

You can specify the wavelength, the frequency is therefore equal to L/lambda where L is the characteristic wavelength specified in the parameter WavelengthAdim.

Example :

# for a wavelength of 200nm :
Wavelength = 200e-9

Location :

WavelengthAdim

Syntax :

  WavelengthAdim = L

You can specify the characteristic wavelength in this field. It represents the unit used in the mesh (a distance of 1 in the mesh represents L meters).

Example :

# if you have constructed your mesh with units of micrometers
Adimensionalization = YES
WavelengthAdim = 1e-6
  
# for a wavelength of 200nm :
Wavelength = 200e-9

Location :

IncidentAngle

Syntax :

IncidentAngle = teta phi

When the selected source is an incident plane wave, the incident angle is specified by these two angles in 3-D, and only by an angle in 2-D. The values are expressed in degrees. In 2-D the wave vector is equal to:

where k is the wave number. In 3-D the wave vector is given as:

Example :

# In 3-D (or axisymmetric), you enter theta and phi
IncidentAngle = 90 90

Location :

Polarization

Syntax :

Polarization = p0 p1 p2 ... pk

When the spatial dependence of the source is scalar, the polarization is the direction of the source, i.e. the constant vector appearing in the source. The number of components of the polarization depends on the solved equation. For scalar equations like Helmholtz equation, this value has no meaning, whereas it is useful for vectorial equations (like Maxwell equations). The default value is the constant vector (1,0,0...,0). When a total or diffracted field is computed, this polarization is the polarization of the incident wave.

Example :

# For maxwell equations, three components
# and you can specify polarization along ey
Polarization = 0 1.0 0

Location :

UseWarburtonTrick

Syntax :

UseWarburtonTrick = YES

This trick can be used to reduce the storage for discontinuous Galerkin. Basis functions are transformed as :

such that the mass matrix no longer depends on the geometry. As a result, the inverse of the mass matrix is not stored. This trick is interesting for curved tetrahedral elements.

Example :

UseWarburtonTrick = YES

Location :

Harmonic/VarAxisymProblem.cxx
Harmonic/VarGeometryProblem.cxx

UseSameDofsForPeriodicCondition

Syntax :

UseSameDofsForPeriodicCondition = type

For periodic condition, you can choose to use same dof numbers for periodic dofs. With this choice, the finite element matrix is naturally symmetric. But you cannot have a quasi-periodic condition (used for PERIODIC_X boundary condition). For this kind of condition, it is better to use a strong formulation or a weak formulation. The choices are the following:

YES : periodic dofs have the same numbers
NO : periodic dofs have different numbers, and the relation u_i = u_j is set strongly in the finite element matrix (instead of rows). The finite element matrix is non-symmetric.
WEAK : periodic dofs have different number and the periodic condition is enforced in the weak formulation. The finite element matrix is symmetric.

Example :

UseSameDofsForPeriodicCondition = WEAK

Location :

Mesh/NumberMesh.cxx

OriginePhase

Syntax :

OriginePhase = x0 y0 z0

For incident waves, the phase is equal to 1 at the origin O. You can set this origin at a different value from the default value (0,0,0).

Example :

# If you want the origin at (5, 0, 3)
OriginePhase = 5 0 3

Location :

ReferenceInfinity

Syntax :

ReferenceInfinity = ref

For incident waves, the infinite homogeneous media is often considered with physical indices (epsilon and mu) equal to 1. You can give the reference (in the mesh) for this infinite homogeneous media such that the wave number is correctly computed.

Example :

# if elements at infinity have reference 2
ReferenceInfinity = 2

Location :

OrderGeometry

Syntax :

OrderGeometry = order

By default, the same order is used both for geometry and finite element (isoparametric elements). However, you may want to specify different orders. Moreover, when you adopt a variable order strategy, you need to specify the order used to approximate geometry. Since the mesh uses classical nodal elements, a fixed order is used to describe curved elements.

Example :

# for example you use a variable order
OrderDiscretization = MEAN_EDGE AUTO 0.5

# you need to specify the order for the geometry
# otherwise a first-order approximation is used (linear elements)
OrderGeometry = 6

OrderDiscretization

Syntax :

OrderDiscretization = n

The order of approximation should be chosen depending on the mesh size and the wavelength. Usually, we will take 10 points per wavelength for a Q2 scheme and 7 points per wavelength for a Q7 scheme, i.e. the mesh size will be equal to the wavelength for an order of approximation equal to 7. Variable orders can also be specified (for Discontinuous Galerkin formulation or with hierarchical elements such as TRIANGLE_HIERARCHIC/TETRAHEDRON_HIERARCHIC). The attribution of an order to an element is purely geometric and based on this rule of thumb "10 points per wavelength". For each edge, an order r is attributed if the length of the edge is comprised between :

Here λ denotes the wavelength of the media (on this edge). When AUTO is provided, we take the following values for h_r

For example, it means that when there are more than 50 points per wavelength (the inverse of 0.02), a first-order is considered accurate enough and is chosen. Between 13.33 points and 50 points, a second-order is chosen, etc. You can replace AUTO with your own sequence of sizes h_r.

Example :

# the same order 3 will be used for all elements
OrderDiscretization = 3

# a variable order will be used
# for faces/volumes, we take the maximum among edge lengths to determine the order
# the last parameter if here a coefficient used to refine in order when it is decreased
# 1.0 is the nominal value for 10 points per wavelength
OrderDiscretization = MAX_EDGE AUTO 0.5

# you can define your own sequence of h_r
# MEAN_EDGE : we take the mean edge length
# the coefficient is here placed at the beginning for a manual sequence
OrderDiscretization = MEAN_EDGE 0.5   0.01 0.1 0.2 0.3 0.4 0.55 0.7 0.9 1.1 1.3   

# For high order, nodal triangles have the drawback that interpolation
# nodes are not present in the code
# You can either use a geometric order below 12 : (the geometry order is the second argument)
OrderDiscretization = 20 10 

# if you are sure that the mesh only contains quad/hexa, you can inform the code with the keyword QUAD
OrderDiscretization = 20 QUAD

Location :

Mesh/NumberMesh.cxx

FileMesh

Syntax :

FileMesh = nom_fichier
FileMesh = nom_fichier REFINED n

# 2-D regular mesh
FileMesh = REGULAR nx xmin xmax ymin ymax ref ref1 ref2 ref3 ref4
FileMesh = REGULAR_ANISO nx ny xmin xmax ymin ymax ref ref1 ref2 ref3 ref4

# 3-D regular meshes
FileMesh = REGULAR nx xmin xmax ymin ymax zmin zmax ref ref1 ref2 ref3 ref4 ref5 ref6
FileMesh = REGULAR_ANISO nx ny nz xmin xmax ymin ymax zmin zmax ref ref1 ref2 ref3 ref4 ref5 ref6
FileMesh = SPHERE rmin rmax nb_points 
FileMesh = SPHERE rmin rmax nb_points r0 r1 ... rk
FileMesh = SPHERE_CUBE rmin rmax nb_points
FileMesh = BALL rmin rmax nb_points ref REGULAR
FileMesh = BALL_CUBE rmin rmax nb_points ref REGULAR
FileMesh = EXTRUDE nb_layers z0 z1 ... zn nom_fichier AUTO \
           nb_ref_surf ref0 ref1 ref2 ... ref2n \
           nb_ref_vol ref0 ref1 ref2 ... refn

The file mesh can be either a file or a primitive (like a cube, sphere) for simple shapes. nom_fichier denotes the file name, the optional parameters "REFINED n" allows the user to split the initial mesh, each edge being divided into n small edges. This splitting procedure is implemented for hybrid meshes in 2-D and 3-D.

In 2-D, you can specify a rectangular domain [xmin, xmax] x [ymin, ymax] ( REGULAR ) with equally distributed points, nx is the number of points along x, the number of points along y is automatically computed. ref is the reference of the 2-D domain while ref1, ref2, ref3 and ref4 are the references of the four boundaries, respectively y = ymin, x = xmax, y = ymax and x = xmin. By default, the code computes automatically the number of points along y (by trying to have approximatively the same space step in x and y). If you want to enter manually nx and ny, you can use the keyword REGULAR_ANISO instead of REGULAR.

The two types of regular meshes in 2-D (REGULAR) with specifying TypeMesh = TRI or TypeMesh = QUAD

In 3-D, you can specify a parallelepiped [xmin, xmax] x [ymin, ymax] x [zmin, zmax], again you select only the number of points along x. ref is the reference of the 3-D domain, while ref1, ref2, ref3, ref4, ref5, ref6 are the references of the six boundaries, respectively x = xmin, y = ymin, z = zmin, z = zmax, y = ymax, x = xmax. The parameter SPHERE corresponds to spherical layers between an intern sphere of radius rmin (strictly positive) and extern sphere of radius rmax. nb_points is the number of points per unit, optional parameters r0, r1, ... rk are radius of intermediary spheres. The parameter SPHERE_CUBE provides same functionality as SPHERE except that the extern boundary is a cube, no intermediary layer is possible. You can also have the inside of the sphere meshed if you specify BALL (and BALL_CUBE if the extern boundary is a cube).

Several patterns for meshing an elementary cell (REGULAR).

Two left meshes have been generated by using SPHERE, while right mesh has been generated by using SPHERE_CUBE.

Meshes generated by using BALL, with TypeMesh = HEXA and TypeMesh = PYRAMID.

Finally, you can create a mesh by extruding a 2-D mesh along axis z (it will create hexahedra and wedges), intermediary layers can be inserted at z-coordinate z0, z1, z2 ..., zn. The initial 2-D mesh is read in the file whose name is nom_fichier. "AUTO" means here that the number of cells along z will be computed by evaluating the mesh step in 2-D, and trying to have the same mesh step for z. During the extrusion, you can change the references of lateral faces, by selecting appropriate ref0 ref1 ref2 ... refn, and the reference of top faces by chosing appropriate refn+1 ... ref2n. If a reference is set to 0, no referenced faces will be created at this place. For volume domains, you can also specify the reference for each domain, if a reference is set to 0, no volume is created, therefore you have created a hole.

Meshes generated by using EXTRUDE.

If you want to extrude in x or y (instead of z), EXTRUDE should be replaced by EXTRUDE_X or EXTRUDE_Y. Depending on the keyword TypeMesh, the initial mesh will be split. For pure hexahedral finite element, tetrahedral meshes will be split into hexahedra, for example.

Example :

# standard use -> a single file name
# the file format is detected thanks to the extension
FileMesh = cube.mesh
# you can also divide the initial mesh (each edge is subdivided into three edges here)
FileMesh = cube.mesh REFINED 3

# 2-D regular mesh of [-1,1]^2
FileMesh = REGULAR 11 -1 1 -1 1   1   1 2 3 4

# 2-D regular mesh of [-1,1]^2 with anisotropic space step
# here 11 points in x and 51 points in y
FileMesh = REGULAR_ANISO 11 51 -1 1 -1 1   1   1 2 3 4

# 3-D regular mesh of [-2,2]^3
FileMesh = REGULAR 21 -2 2 -2 2 -2 2   1   1 2 3 4 5 6

# 3-D regular mesh of [-2,2]^3 with anisotropic space step
# here 11 points in x, 21 in y and 5 in z
FileMesh = REGULAR_ANISO 11 21 5  -2 2 -2 2 -2 2   1   1 2 3 4 5 6

# regular mesh between sphere of radius 1 and sphere of radius 2
FileMesh = SPHERE 1.0 2.0 5.0

# regular mesh between sphere of radius 2 and cube [-3,3]^3
FileMesh = SPHERE_CUBE 2.0 3.0  2.5

# mesh of a ball of radius 2 with an internal boundary at r = 1
FileMesh = BALL 1.0 2.0 4.5 REGULAR

# mesh of a cylinder, base is a mesh of a disk, z belongs to [-1,1]
# we consider that the 2-D mesh contains only one reference
# (nb_ref_surf = nb_ref_volume = 1)
# therefore 2 will be the reference of lateral section and 3 of the
# top section
FileMesh = EXTRUDE 1 -1.0 1.0 disque.mesh AUTO \
           1 2 3 \
           1 1

Location :

AdditionalMesh

Syntax :

AdditionalMesh = nom_fichier

Actually, the syntax is the same as FileMesh. The aim of this keyword is to load additional meshes after the initial mesh. For almost all the simulations, these additional meshes won't be read... The initial aim of this entry was to deal with non-conformal meshes and load each part of the mesh separately, but non-conformal meshes are not currently handled by Montjoie .

Example :

# First, you specify the initial mesh
FileMesh = first.mesh

# then you can add as many meshes as you want
AdditionalMesh = second.mesh
AdditionalMesh = third.mesh

Location :

Harmonic/DistributedProblem.cxx

NbProcessorsPerMode

Syntax :

NbProcessorsPerMode = n

This entry is still experimental and should not be used.

Example :

NbProcessorsPerMode = 8

Location :

SaveSplitDomain

Syntax :

SaveSplitDomain = YES file_name

If you add this keyword, the partioning of the mesh is saved through the array Epart. The Epart array is an array that gives for each element the processor number. This array is saved in binary in the specified file.

Example :

# partitioning saved in the file Epart.dat
SaveSplitDomain = YES Epart.dat

Location :

Harmonic/DistributedProblem.cxx

SplitDomain

Syntax :

SplitDomain = type_partitioning

This entry specifies the partioning tool used to split the mesh into several sub-domains for execution on parallel architecture. The default partitioning is Metis, you have the following choices:

Metis : Metis partioning
Scotch : Scotch partioning (available if Montjoie has been compiled with Scotch)
Concentric : the mesh is split in concentric layers
User : the user provides directly the Epart array
Layered : the user provides a mesh with different physical areas

The Epart array is an array that gives for each element the processor number. If User partioning is selected, the user must provide a file containing this array of integers (in binary format). For concentric layers, the user provides a set of radii r_n. All the elements whose center satisfy the condition :

belong to the processor n+1. The processor 0 takes elements such that :

The last processor takes elements such that:

This partioning creates circular layers in 2-D, and spherical layers in 3-D.

Example :

# partitioning with Scotch
SplitDomain = Scotch

# user-partitioning 
# the user provides the array epart in a binary file
SplitDomain = User file_with_epart.dat

# concentric partioning
# the user provides radii, and elements between r_n and r_{n+1}
# belong to processor n
# for two processors, you have to specify only one radius
SplitDomain = Concentric 1.5

# Layered partitioning
# each reference (Physical Volume in 3-D, Physical Surface in 2-D)
# belongs to a different processor
SplitDomain = Layered AUTO

# you can manually set the references for each processor
# for each processor, you write nb_ref ref0 ref1 ... ref_n
# In this example, processor 0 is associated with references 1 and 5
# processor 1 is associated with references 2, 4 and 6
SplitDomain = Layered 2  1 5   3   2 4 6

Location :

Harmonic/DistributedProblem.cxx

TypeMesh

Syntax :

TypeMesh = type

You can specify the type of elements to use in the case of regular meshes generated by the code as described in the description of the keyword FileMesh. This type is also used to determine if the original mesh is split into small elements to satisfy the requirements. In 2-D, you have the following types:

TRI : the mesh contains only triangles. If the original mesh contains quadrilaterals, each quadrilateral is split into two triangles.
QUAD : the mesh contains only quadrilaterals. If the original mesh contains triangles, each triangle is split into three quadrilaterals and each quadrangle is split into four quadrilaterals

In 3-D, you have the following types:

TETRA : the mesh contains only tetrahedra. If the original contains other elements, they are split into tetrahedra.
HEXA : the mesh contains only hexahedra. If the mesh is purely tetrahedral, each tetrahedra is split into four hexahedra.
PYRAMID : the regular mesh contains only pyramids. No splitting is performed if the original mesh contains other elements.
WEDGE : the regular mesh contains only wedges (triangular prisms). No splitting is performed if the original mesh contains other elements.
HYBRID : the regular mesh contains pyramids and tetrahedra.

Below, we show some illustrations of regular meshes created by Montjoie

The two types of regular meshes in 2-D (REGULAR) with specifying TypeMesh = TRI or TypeMesh = QUAD

Several patterns for meshing an elementary cell (REGULAR).

Example :

# In 2-D, you can use triangular or quadrilateral mesh
TypeMesh = TRI
TypeMesh = QUAD

# In 3-D, more choices are available
# HYBRID corresponds to the case where each cube is made of 2 tets and
# 2 pyramids
TypeMesh = TETRA
TypeMesh = PYRAMID
TypeMesh = HYBRID
TypeMesh = WEDGE
TypeMesh = HEXA

Location :

Mesh/Mesh2D.cxx
Mesh/Mesh3D.cxx

MeshPath

Syntax :

MeshPath = directory

Often the meshes are always in the same directory. As a result, you can specify this directory, and the meshes will be searched inside this directory. The default directory is the current directory.

Example :

# You can use current directory and enter the relative path of the mesh in FileMesh
MeshPath = ./
FileMesh = MyDirectory/MyMeshes/cube.mesh

# Or you can separate the directory from the file
MeshPath = MyDirectory/MyMeshes/
FileMesh = cube.mesh

Location :

IrregularMesh

Syntax :

IrregularMesh = YES
IrregularMesh = NO

If you want bad quadrilateral/hexahedral meshes, you can insert this entry, and it will generate tetrahedral meshes split into hexahedra. In 2-D, the regular mesh will be made of triangles split into quadrilaterals. The default value is NO.

Example :

IrregularMesh = YES

Location :

SlightModificationOnRegularMesh

Syntax :

SlightModificationOnRegularMesh = YES 
SlightModificationOnRegularMesh = RANDOM

This functionality is available in 3-D only and modifies the regular mesh (if FileMesh = REGULAR is present) such that there are non-affine elements. In order to generate non-affine elements, one vertex out of 2 (in each direction) is translated:

If YES is selected each odd vertex is translated with this same quantity. If RANDOM is selected, each odd vertex is vertex by this quantity multiplied by random numbers (between -1 and 1) for each component.

Example :

# non-affine elements with a fixed translation vector
SlightModificationOnRegularMesh = YES

# non-affine elements with a random translation
SlightModificationOnRegularMesh = RANDOM

Location :

Mesh/Mesh3D.cxx

RefinementVertex

Syntax :

RefinementVertex = x y level ratio

For 2-D meshes, you can ask a local refinement around the vertex of coordinates (x,y). level is the number of refinements performed, and ratio the quotient of coarse edge to fine edge. There is no local refinement applied by default.

Local refinement around vertex (0, 0) with two levels and ratio equal to 3.0.

Example :

# around point (1,0), you locally refine 5 times with a ratio of 2.
RefinementVertex = 1 0 5 2.0

Location :

AddVertex

Syntax :

AddVertex = x y

For 2-D meshes, you can ask to create a new vertex. The mesh will be modified in order to include this vertex while being conformal. There can be triangles generated by the process.

Example :

# you want to create a new point (-2.5, 0.3)
AddVertex = -2.5 0.3

# then you can ask a local refinement around this vertex
RefinementVertex = -2.5 0.3  5  2.0

Location :

ThresholdMesh

Syntax :

ThresholdMesh = eps

You can specify a given threshold for 2-D (or 3-D points). The default value is 1e-10. If the distance between two points is below this threshold, they are assumed to be the same.

Example :

# depending on the characteristic length, you may modify the threshold for points of the mesh :
ThresholdMesh = 1e-12

Location :

Harmonic/VarProblemBase.cxx
Elliptic/PhysicalProperty.cxx

TypeCurve

Syntax :

# curves in 2-D
# a circle of center (x0, y0) and radius r0 
TypeCurve = ref1 ... refn CIRCLE x0 y0 r0
# an ellipse (of axis x and z)
TypeCurve = ref1 ... refn ELLIPSE x0 y0 a b
# a C1 spline
TypeCurve = ref1 ... refn SPLINE
# local spline (i.e. on each edge, we use a third order approximation
by using neighboring points, but this technique does not lead to a C1
curve)
TypeCurve = ref1 ... refn LOCAL_SPLINE


# curves in 3-D
# a sphere of center (x0, y0, z0, r0)
TypeCurve = ref1 ... refn SPHERE x0 y0 z0 r0
# a cylinder of radius r0 and axis (u0, v0, w0) and point (x0, y0, z0)
# belongs to the axi 
TypeCurve = ref1 ... refn CYLINDER x0 y0 z0 u0 v0 w0  r0
# a cone of center (x0, y0, z0), of axis (u0, v0, w0) and angle teta
TypeCurve = ref1 ... refn CONE x0 y0 z0 u0 v0 w0 teta

You can specify curved surfaces, but a limited choice of curves is proposed. A more general solution is to directly provide a curved mesh (Gmsh produces curved meshes for instance). In 2-D, the use of splines is quite general, but provides only a third order approximation of the geometry. By default, no curves are present (elements are straight). In 2-D the following curves are proposed:

CIRCLE x0 y0 radius: the curve is an arc of circle. The parameters can be given, if not the code will find the correct parameters automatically.
ELLIPSE x0 y0 a b: the curve is an arc of ellipse.
PEANUT x0 y0 a b c: the curve is a part of a peanut.
SPLINE : the curve is approximated by cubic splines
LOCAL_SPLINE : the curve is approximated by cubic polynomials (computed with only neighbor points).

In 3-D, the following curves are proposed:

SPHERE x0 y0 z0 radius: the curved surface is a part of a sphere. The parameters can be given, if not the code will find the correct parameters automatically.
CYLINDER x0 y0 z0 u0 v0 w0 radius: the curved surface is a part of a cylinder whose axis is directed along (u0, v0, w0) and (x0, y0, z0) is a point on the axis. The parameters can be given, if not the code will find the correct parameters automatically.
CONE x0 y0 z0 u0 v0 w0 alpha: the curved surface is a part of a cone whose axis is directed along (u0, v0, w0) and (x0, y0, z0) is a point on the axis, alpha is the angle of the cone. The parameters can be given, if not the code will find the correct parameters automatically.

Example :

# If you don't specify the parameters of the circle, the
# code finds them automatically
TypeCurve = 1 CIRCLE 

# For ellispes, you need to fill parameters ...
TypeCurve = 2 ELLIPSE 0.0 0.0  3.0 1.0

# For SPLINE, no bifurcation is allowed (a curve with two branches).
# It is allowed for LOCAL_SPLINE, but the result should be ugly ...
TypeCurve = 4 SPLINE


# For spheres, cylinders and cones, you don't need to specify parameters, they are automatically found by an optimization procedure
TypeCurve = 3 SPHERE
TypeCurve = 5 CYLINDER
TypeCurve = 6 CONE

Location :

Mesh/MeshBoundaries.cxx

MateriauDielec

Syntax :

MateriauDielec = ref ISOTROPE rho mu sigma

This entry is useful to specify physical properties of a media, such as dielectric permittivity epsilon and magnetic permeability mu for Maxwell's equations. Since it depends on the considered equation, we refer the reader to the sections devoted to each equation : Helmholtz, Maxwell, Aeroacoustics, Elastic. The default values are 1 (and 0 for damping indexes). Each index can be constant (e.g. 2.4, 0.5, etc) or variable. In order to specify a variable coefficient, we expect one of the following value :

RANDOM : the index is given through its values on a regular grid. Fourth-order interpolation is performed to compute the values on quadrature points
SINUS : the index is a product of sine functions
MESH : the index is given on a mesh (different from the computational mesh), currently this functionality has been disabled
SAME_MESH : the index is given on the same mesh as the computational mesh
PERIODIC : it is similar to RANDOM, the index is periodized outside the bounding box
QUASI_PERIODIC : it is similar to RANDOM, the index is periodized outside the bounding box (except on the center cell where the index is constant)
RADIAL : the index is assumed to be radial (it depends only on r). In a file the user gives the radii and associated values of the index. Cubic spline interpolation is used to interpolate on quadrature points.
USER : the index is given analytically by the user in the file UserSource.cxx (function ComputeVariableUserIndex). In that case, no interpolation is performed.
VELOCITY : the index is equal to 1/c² where c is the given constant
SQUARE : the index is equal to c² where c is the given constant
ANGLE : the index is equal to teta π/180 where teta is the given constant
any other value : the index is constant and equal to the given value

For each keyword, additional parameters are expected. As a result, each coefficient rho, mu, sigma indicated by the line :

MateriauDielec = ref ISOTROPE rho mu sigma

can be replaced by one of the following choice :

constant_value
ANGLE teta
SQUARE c
VELOCITY c
USER offset amplitude coef_infinity
RADIAL file_index
RANDOM xmin xmax ymin ymax (zmin zmax) coefx coefy (coefz) offset amplitude file_index PRECISION coef_infinity
PERIODIC xmin xmax ymin ymax (zmin zmax) coefx coefy (coefz) offset amplitude file_index PRECISION coef_infinity
QUASI_PERIODIC xmin xmax ymin ymax (zmin zmax) coefx coefy (coefz) offset amplitude file_index PRECISION coef_infinity
MESH mesh_file offset amplitude data_file order coef_infinity
SAME_MESH mesh_file offset amplitude data_file order coef_infinity
SINUS xmin xmax ymin ymax (zmin zmax) coefx coefy (coefz) offset amplitude kx ky (kz) coef_infinity

The fields inside parenthesis are needed in 3-D only. coef_infinity is needed when the inhomogeneous media encounters an absorbing boundary condition or PML layers, because they need the knowledge of the value of physical indices at infinity. All these parameters can lead to a lengthy line, usually we put the character '\' to separate each index. Another alternative is to use the keyword PhysicalMedia such that you specify each index on a separated line. Below we describe the different choices of varying media.

SINUS is used when you want to have a sinusoidal variation of coefficient on a bounding box [xmin, xmax] x [ymin, ymax]. Outside the bounding box, the coefficient will be equal to offset. The formula used for SINUS is :

MESH is used when you know the coefficient on a mesh. There will be problems if the mesh provided does not contain the mesh where the simulation is performed, that often occurs when curved boundaries are present. As a result, this functionality has been disabled. If you use the same mesh for the simulation and the definition of the coefficient, you have to write SAME_MESH since the interpolation is obvious and is not time consuming. In order to produce the data file, you have to start from an initial data file such as:

TypeElement = QUADRANGLE_LOBATTO
TypeEquation = TIME_REISSNER_MINDLIN

# path where the meshes are stored
MeshPath = ./

# frequency of the temporal source
# Frequency = a b means that the pulsation omega = 2*pi*a + b
Frequency = 5.0 0.0

# file where the mesh is stored
FileMesh = plaque.mesh

# In practice all the triangles are split into quadrilaterals
# so you have to write this line
TypeMesh = QUAD

# order of space approximation
OrderDiscretization = 4

# other fields are not very important for the next step

The only requirement in this data file is that all lines related to the mesh are present (such as RefinementVertex, TypeMesh, etc) and the order of discretization. Nodal points on the resulting mesh are generated by write_nodal. If write_nodal has not been compiled, you can compile it by writing :

make write_nodal

Nodal points (that are written in the file nodal_points.dat) and the resulting mesh (written in the file test.mesh) are generated by executing write_nodal:

./write_nodal.x time_carre_plate.ini

The second argument is the name of the data file. A third optional argument is the order of approximation. If this argument is provided, eg.:

./write_nodal.x time_carre_plate.ini 10

It will ignore the value given in the data file and use this value (here 10). It can be useful if you want to use different orders of approximation, one for the finite element computation and one for the definition of the index. It is advised to copy the resulting mesh in another file :

cp test.mesh test_plaque.mesh

Then, you can implement the computation of the index from coordinates. You can start from the file src/Program/Test/write_index.cc, copy it in another file and write the implementation of your index. write_index.cc evaluates the index as a continuous index. You can provide piecewise continuous indexes, with continuous index for each reference of the mesh. If you want to specify a piecewise continuous index, you can start from the file src/Program/Test/write_index_plaque.cc and modify the function ComputeIndex:

// x, y : coordinates of the 2-D point where the index must be evaluated
// ref : reference of the element where the point is
// n : index number
// delta : value of the index n at this point
void ComputeIndex(const Real_wp& x, const Real_wp& y, int ref, int n, Complexe& delta)
{
  switch (ref)
    {
    case 1: delta = 0.002*exp(-(x*x+y*y)); break;
    case 2 : delta = 0.005; break;
    case 5 : delta = 0.01*exp(-0.3*x*x); break;
    case 6 : delta = 0.007*exp(-0.5*y*y); break;
    }
}

The file containing the values of the index is then generated by compiling and executing the file you have modified :

./write.x nodal_points.dat test_plaque.mesh index.don

You can specify several file names after index.don if there are several indexes to evaluate. In that case, n will refer to the index number. Once the index file has been generated, you can fill your line in the data file

# SAME_MESH maillage offset amplitude file_index order coef_infinity
MateriauDielec = 1 ISOTROPE SAME_MESH test_plaque.mesh 0.0 1.0 index.don 4 1.0
MateriauDielec = 2 ISOTROPE SAME_MESH test_plaque.mesh 0.0 1.0 index.don 4 1.0

The order provided is the order used to generate the file nodal_points.dat. It can be different from the order used in the simulation. Here you have to repeat the same line for each reference. The final value of rho will be equal to :

offset + amplitude * delta

where delta corresponds to values contained in the file index.don.

RANDOM is used when you want to have a variation of coefficient specified by a file on a box [xmin, xmax] x [ymin, ymax] (x [zmin, zmax] in 3-D). Outside the box, the coefficient will be equal to offset. The formula used for RANDOM is :

The function h(x,y) is computed with fourth-order interpolation of values of h on a regular grid (these values are directly given in a file). You can write these values in binary format from a dense matrix nu with the Matlab functions write_full_?.m (contained in the folder MATLAB, with the Blas convention c, d, s, z for complex float, double, float, and complex double). You can also write the indices with Python, for example if you execute the following lines:

X = linspace(-5.0, 5.0, 400);
Y = linspace(-3.0, 3.0, 240);
[XI, YI] = meshgrid(X, Y);
rho = cos(0.5*pi*XI)*sin(pi*YI);
write_full("../rho_var.don", transpose(rho));

The transpose is here due to the meshgrid function and is inherent to 2-D fields. In Matlab write_full has to be replaced by write_full_d in this case (double precision). In Python, write_full handles both real and complex data, a third argument optional is used to switch to single precision :

write_full("../rho_var.don", transpose(rho), False);

In this last case, rho_var.don is written in single precision. In the data file, you have to specify the type of data (FLOAT, DOUBLE, COMPLEX or COMPLEX_DBLE). For example, you can write

Example :

# MateriauDielec = 1 ISOTROPE RANDOM xmin xmax ymin ymax coefx coefy offset amplitude file_index PRECISION coef_infinity
# the two last values here are mu and sigma
MateriauDielec = 1 ISOTROPE RANDOM -2 2 -2 2 0 0 0.5 0.1 rho_var.don DOUBLE 1.0 \
                                   2.5 0.5

PERIODIC uses the same formula as for RANDOM. The difference is that the point x is translated into the specified bounding box before applying the formula. By construction, the index will be periodic with a period given by the bounding box.

QUASI_PERIODIC is almost equivalent to PERIODIC. The only difference is that if the point is in the central cell (not translated), the index will be set to the given offset. As a result the physical index is uniform in the central cell, and periodic in all other cells.

USER will use directly the function ComputeVariableUserIndex to compute the variable index. In this function, the index number refers to a given index (0 for rho, 1 for mu, 2 for sigma, etc) and the component number to a component of the index (for vectors or tensors). In that case, the user must modify the file Source/UserSource.cxx to enter the desired expression, recompile the code and launch it. The data file will have a line of the type

Example :

# if you want a variable rho and mu
# USER offset amplitude coef_infinity for each index
MateriauDielec = 1 ISOTROPE USER 0.5 1.0 1.0 USER 0.0 2.0 1.0

RADIAL specify an index that depends on the norm of x (i.e. the radius). The values of the index are given in a text file with two columns:

r0 index0
r1 index1
r2 index2

The first column contains the different radii, the second column the associated values. A cubic spline is contructed from this data to interpolate the index on any point.

If you have a print level greater or equal to 7, and if you are running the code in sequential, varying indexes will be written in files :

rho_ref1_G0_U0.dat
rho_ref1_G1_U0.dat
mu_ref1_G0_U0.dat
rho_ref2_G0_U0.dat
...

The file names begin with the name of the physical index (rho, sigma, mu, epsilon, etc) then the reference number follows, then the grid number and finally the component of the index.

Location :

PhysicalMedia

Syntax :

PhysicalMedia = name_media ref ISOTROPE parameters

This entry is useful to specify physical properties of a media. The media can be uniform, in that case parameters will contain the numerical values of the media. It can also be non-uniform, and each component of the physical media has to be specified as detailed in MateriauDielec. The keyword PhysicalMedia is more convenient, because you specify a different kind of anisotropy for each media whereas the type of anisotropy is common in MateriauDielec. The names of physical media are detailed in the section devoted to the desired equation.

Example :

# For Maxwell's equations, there are three physical indexes : epsilon, mu, sigma
# here we specify different anisotropy for the three indexes (which are matrices)
PhysicalMedia = epsilon 1 ISOTROPE 2.0
PhysicalMedia = mu 1 ORTHOTROPE 1.0 1.5 1.8
PhysicalMedia = sigma 1 ANISOTROPE 0.1 0.02 0.01  0.2 0.005  0.3

ExactIntegration

Syntax :

ExactIntegration = YES
ExactIntegration = NO
ExactIntegration = YES order
ExactIntegration = NO order

For discontinuous Galerkin formulation, if you say YES, you will use an unsplit formulation, which will be probably stable if integrals are exactly computed (i.e. you are using Gauss points), if you say NO, a split formulation will be used, slower but stable for approximate integration. Details are avaible in the paper of N.Castel, G. Cohen and M. Duruflé. This value is important only for aeroacoustics and default value is YES. The optional parameter order is the over-integration that can be used.

Example :

ExactIntegration = NO

Location :

Harmonic/VarProblemBase.cxx
Elliptic/PhysicalProperty.cxx

MixedFormulation

Syntax :

MixedFormulation = YES_OR_NO

If the user specifies a mixed formulation, the second-order equation such as Helmholtz equation :

is transformed in a first-order system, e.g :

This kind of transformation is implemented for Helmholtz equation, Maxwell's equations and elastodynamics. When continuous (or edge elements for Maxwell's equations) finite elements are used to solve one of these equations, the additional unknown v is chosen discontinuous. The advantage of this first-order system is that it is well adapted to time-domain simulations (with the adjunction of PML layers). This mixed formulation is automatically chosen if the selected time scheme is a scheme adapted to first-order formulation (such as Runge-Kutta schemes). For time-harmonic simulations, there is no interest in considering this mixed formulation. As a result, the mixed formulation is not chosen by default, you can force to use it by inserting this line in the data file. This formulation is interesting for the computation of eigenmodes, since the eigenvalue problem is linear with respect to omega even with PML or first-order absorbing boundary condition.

Example :

MixedFormulation = YES

Location :

PenalizationDG

Syntax :

PenalizationDG = alpha delta

For Discontinuous Galerkin formulation and linear problems, additional stabilizing terms need to be incorporated into the discontinuous Galerkin formulation (we can call them penalty terms as well). This is somehow equivalent to choose upwind fluxes instead of centered fluxes, which correspond to take alpha and delta equal to 0. Default choice is upwind fluxes, negative values provide more stability (positive values are instable).

Example :

# Classical choice of parameters -1, -1
PenalizationDG = -1.0 -1.0

# upwind fluxes :
PenalizationDG = Upwind

# centered fluxes :
PenalizationDG = 0.0 0.0

Location :

CoefficientPenalization

Syntax :

CoefficientPenalization = AUTO

For Discontinuous Galerkin formulation and linear problems, additional stabilizing terms need to be incorporated into the discontinuous Galerkin formulation (we can call them penalty terms as well). For SIPG (Symmetric Interior Penalty Galerkin), the coefficient depends on the order and the mesh size as:

By default, this coefficient is chosen and multiplied by the first coefficient given in PenalizationDG. If you do not want to start from this value, you can put this line with MANUAL (instead of AUTO) and the coefficient is then directly given in PenalizationDG.

Example :

# if you put auto, the penalization term is equal to
#  beta \int_{\partial K} [u] [v] (for Helmholtz equation)
# beta will be equal to 2.0 * alpha where alpha = r(r+1) / h
# usually this value works pretty well for any kind of mesh, you can
# modify the 2.0 in Penalization DG (the negative sign is mandatory)
CoefficientPenalization = AUTO
PenalizationDG = -2.0 0.0

# if you want to set beta directly :
CoefficientPenalization = MANUAL
PenalizationDG = -100.0 0.0

Location :

PointDirichlet

Syntax :

  PointDirichlet = x0 y0 z0

You can specify dirichlet conditions on vertices for continuous finite elements by entering the coordinates of the vertices. When it is used to recover uniqueness of the solution for a Poisson equation with only Neumann boundary conditions, it doesn't work as nicely as expected since the point where Dirichlet condition is enforced is obvious to detect when you look at the solution. There is no default value.

Example :

# in 3-D :
PointDirichlet = 1.0 -2.0 3.0

# you can add as many lines you want
# each line will appends a new Dirichlet node
PointDirichlet = 0.0 5.0 2.4

Location :

OrderAbsorbingBoundaryCondition

Syntax :

  OrderAbsorbingBoundaryCondition = order

When you write ConditionReference = 1 ABSORBING, an absorbing boundary condition is set on edges of reference 1. By default, this condition is a first order absorbing boundary condition. For Helmholtz equation only, some high order absorbing boundary conditions are available.

Example :

OrderAbsorbingBoundaryCondition = 1

# For Helmholtz equation, 2 is the classical second-order ABC with Laplace-Beltrami operator (Joly-Vacus)
OrderAbsorbingBoundaryCondition = 2

# other boundary conditions has been developed by Chabassier, Bergot, Barucq
OrderAbsorbingBoundaryCondition = 2 GRAZING

# PARAMETERS gamma theta zeta
OrderAbsorbingBoundaryCondition = 2 PARAMETERS 0.25 0.5 0.5

# PARAMETERS gamma theta zeta GRAZING
OrderAbsorbingBoundaryCondition = 2 PARAMETERS 0.25 0.5 0.5 GRAZING

Location :

ModifiedImpedance

Syntax :

ModifiedImpedance = CURVE
ModifiedImpedance = coef

For first-order absorbing boundary condition and Helmholtz equation, you can multiply the default impedance value (-i omega sqrt(rho mu) ) by a coefficient. This can be useful for circular boundaries, where the modification of this impedance will reduce reflection. For curved boundary, you can write the keyword CURVE, such that the usual factor h is added (h is the mean curvature).

Example :

# first-order absorbing boundary condition
OrderAbsorbingBoundaryCondition = 1

# usually if the boundary is curved, we have du/dn = (ik + h) u
# where h is the mean curvature (k1+k2)/2
# if you want to take into account this term, you put CURVE
ModifiedImpedance = CURVE

# another possibility is to set directly a coef such that
# du/dn = ik coef u
# you write the coefficient coef below:
ModifiedImpedance = (1.002,0.002)

Location :

OrderHighConductivityBoundaryCondition

Syntax :

OrderHighConductivityBoundaryCondition = order

Objects with a high conductivity can be replaced by an equivalent boundary condition. Order 0 is the classical perfectly conducting boundary condition. Order 1, 2 and 3 are implemented and provide more accurate results. The implementation is achieved only for Helmholtz equation and axisymmetric Maxwell equations. Default value is set to 0.

Example :

  OrderHighConductivityBoundaryCondition = 3

Location :

Exit_IfNo_BoundaryCondition

Syntax :

Exit_IfNo_BoundaryCondition = YES
Exit_IfNo_BoundaryCondition = NO

When an isolated boundary (belonging to only one element) is found, and there is no associated boundary condition, the code will stop the simulation. If you want to continue the simulation, you have to insert this entry.

Example :

Exit_IfNo_BoundaryCondition = NO

Location :

ConditionReference

Syntax :

ConditionReference = n1 n2 ... nk DIRICHLET
ConditionReference = n1 n2 ... nk NEUMANN
ConditionReference = n1 n2 ... nk ABSORBING
ConditionReference = n1 n2 ... nk HIGH_CONDUCTIVITY epsilon
ConditionReference = n1 n2 ... nk IMPEDANCE lambda
ConditionReference = n1 n2 PERIODICITY

Each entry of this type has to be inserted for each considered boundary condition. Dirichlet, Neumann, absorbing, impedance or periodic boundary conditions are classical and implemented for all the equations solved by Montjoie. High conductivity boundary condition is only available for Helmholtz and time-harmonic Maxwell's equations. For each condition, you can specify an arbitrary number of references (the integers n1, n2, ..., nk). These reference numbers are the physical line numbers defined in the mesh (for Gmsh). For periodic conditions, two and only two references are required, the reference of the two parallel hyperplanes.

Example :

# accepted types : DIRICHLET, SUPPORTED, NEUMANN, ABSORBING, HIGH_CONDUCTIVITY, IMPEDANCE
# SUPPORTED is used for Reissner-Mindlin or elastodynamics
# such that you specify Dirichlet condition for some components
# and Neumann condition for other components
ConditionReference = 1 2 3 DIRICHLET
ConditionReference = 4 NEUMANN
ConditionReference = 5 6 ABSORBING
ConditionReference = 5 6 HIGH_CONDUCTIVITY 0.1
ConditionReference = 7 8 9 IMPEDANCE 2.5

# For periodic conditions, you can choose among
# PERIODICITY, PERIODIC_X, PERIODIC_Y, PERIODIC_Z, CYCLIC
# PERIODIC_X, PERIODIC_Y, PERIODIC_Z, CYCLIC will induce that 
# Fourier series of the source are computed such that
# the user gives only a generating part of the geometry
ConditionReference = 10 13 PERIODICITY

Location :

ForceDirichletSymmetry

Syntax :

ForceDirichletSymmetry = YES

By writing YES, the user asks to treat Dirichlet condition in the same way for symmetric or unsymmetric matrices. By default, if the matrix is unsymmetric, rows associated with Dirichlet dofs are replaced by zeros with one on the diagonal. This operation "unsymmetrizes" the matrix. If the matrix is symmetric, rows and columns associated with Dirichlet dofs are erased except one on the diagonal. Columns are stored such that the right-hand side is modified with a hetereogeneous Dirichlet condition. By writing YES, this treatment (usually performed for symmetric matrices) is performed whatever is the type of the finite element matrix. The solution should the same, it may improve the condition number of the matrix.

Example :

ForceDirichletSymmetry = YES

Location :

DirichletCoefMatrix

Syntax :

DirichletCoefMatrix = alpha

By default, if the matrix is unsymmetric, rows associated with Dirichlet dofs are replaced by zeros with one on the diagonal. Instead of one, the user can provide its own coefficient.

Example :

DirichletCoefMatrix = 1e5

Location :

StorageModes

Syntax :

StorageModes = YES

When a cyclic condition is imposed (or PERIODIC_X, PERIODIC_Y, PERIODIC_Z), the solution is computed through its Fourier series. Each mode can be stored, and the final solution is obtained by applying a FFT. However, the user can ask that the modes are not stored (to reduce the memory usage), in that case the solution is incremented (FFT can no longer be used).

Example :

# you can write YES or NO
StorageModes = YES

Location :

NbModesPeriodic

Syntax :

NbModesPeriodic = nx ny nz

When periodic conditions are imposed (PERIODIC_X, PERIODIC_Y and/or PERIODIC_Z), the solution is computed through its Fourier series. You need to give the number of Fourier modes along each direction. For CYCLIC, it is not needed to specifiy this number, since the code will detect automatically the angle between the two cyclic boundaries and deduce the number of modes.

Example :

# if you have periodicity in y and z for example :
ConditionReference = 1 3 PERIODIC_Y
ConditionReference = 2 4 PERIODIC_Z

# Number of modes for each direction (here one mode for x)
NbModesPeriodic = 1 8 16

Location :

TransverseDampingPML

Syntax :

TransverseDampingPML = ratio

You can add a damping in the transverse direction (in y if the PML is in x). By doing that, the PML is no longer PML, but it can eliminate instabilities (e.g. in elastodynamics).

Example :

# this ratio should be small compared to 1
TransverseDampingPML = 0.01

Location :

DampingPML

Syntax :

DampingPML = alpha
DampingPML = fct parameters

You can increase damping coefficient inside PML in order to reduce the reflection on the outside boundary. Default value is set to 1.0, for which a reflection of 0.1 % should be observed for a thickness of one wavelength. This keyword can also be used to specify a different damping function (different from the usual parabolic function). For a parabolic damping we consider the following expression

where a is the thickness of the PML, v_max is the maximal velocity in the PMLs and x₀ the starting position of PML. The coefficient alpha corresponds to the value given in DampingPML. For a constant function, we consider :

SHIFTED_PARABOLE corresponds to the following expression

Example :

# 2.0 should provide an error of 1e-6
DampingPML = 2.0

# if you want to define other damping functions
# for example a constant damping : 
# available choices are CONSTANT, PARABOLE and SHIFTED_PARABOLE
DampingPML = CONSTANT 1.0

# you can indicate a linear profile 
DampingPML = LINEAR 1.0

# for shifted parabole, you specify a and b (in the expression a x^2 + b)
DampingPML = SHIFTED_PARABOLE 2.0 0.5

Location :

AddPML

Syntax :

AddPML = YES axis thickness nb_layers
AddPML = ref type_pml thickness nb_layers

PML layers are not generated during the meshing phase, they are directly added by Montjoie thanks to that entry. You ask the code to add PML layers around the provided mesh. You have to specify the thickness of PML layers (often we take one wavelength), and the number of layers (e.g. one layer is usual with a Q7 approximation). For this last parameter, you can enter AUTO so that the number of layers is automatically computing depending on the thickness and the mesh step. the parameter "axis" is a combination of X, Y, Z, X+, X-, Y+, Y-, Z-, Z+ such that you can add PML layers in any direction you want. By default, no PML layers are added.

Example :

# If you want PML in X and Y, of thickness 1, number of layers is automatically chosen
AddPML = YES XY 1.0 AUTO

# If you want only PML in the right side, and bottom side :
AddPML = YES X+Y- 1.0 2

# If you want PML only for z > 0
AddPML = YES Z+ 2.0 3

# for radial PML : PML_R
AddPML = YES PML_R 1.0 AUTO

# for a mix of radial and cartesian PML, we 
# extrude PML from a reference 
AddPML = 2 PML_R 0.4 2

Location :

Mesh_Boundaries.cxx

EigenvalueTolerance

Syntax :

EigenvalueTolerance = tolerance

You can specify here the stopping criterion used for the determination of eigenvalues. The default value is 1e-6.

Example :

EigenvalueTolerance = 1e-6

Location :

EigenvalueMaxNumberIterations

Syntax :

EigenvalueMaxNumberIterations = nb_iter

You can specify here the maximum number of iterations for the eigenvalue solver. The default value is 1000.

Example :

EigenvalueMaxNumberIterations = 1000

Location :

UseCholeskyForEigenvalue

Syntax :

UseCholeskyForEigenvalue = YES

The Cholesky factorization of the mass matrix is computed such that the generalized eigenvalue problem

is transformed into a standard eigenvalue problem:

where L is the Cholesky factor of the mass matrix M.

Example :

UseCholeskyForEigenvalue = YES

Location :

EigenvalueSolver

Syntax :

EigenvalueSolver = name

You can choose the solver to use (between Arpack, Feast or Slepc). The default solver is Arpack.

Example :

// if you want to choose Slepc
EigenvalueSolver = Slepc

Location :

FileEigenvalue

Syntax :

FileEigenvalue = file_name

You can specify here the file where the eigenpulsations are written. By default, the eigenpulsations are written in the file Omega.dat. They are written in binary format (loadable with load1D in Python).

Example :

FileEigenvalue = Omega.dat

Location :

Eigenvalue

Syntax :

Eigenvalue = YES mode shift nb_eigenvalues

You can ask the computation of eigenvalues of the following generalized eigenvalue problem:

where K is the stiffness matrix, M the mass matrix. For Laplace equation (and continuous finite elements), these matrices are given as:

This eigenvalue problem is solved by calling an eigenvalue solver interfaced in Seldon (i.e. Arpack, Anasazi or Feast). The different modes of eigenvalues computation are the following:

REGULAR : regular mode, the largest eigenvalues of M^-1 K are searched. This modes does not involve a factorization of a linear system.
SHIFTED : shifted mode, the largest eigenvalues of are computed in order to obtain eigenvalues close to the shift sigma.
COMPLEX_SHIFT : complex-shifted mode, the largest eigenvalues of where the shift sigma is complex, only the imaginary part of the factorized matrix is used here (whereas the real part is used for SHIFTED).
BUCKLING : Buckling mode, the largest eigenvalues of are computed in order to obtain eigenvalues close to the shift sigma.
CAYLEY : Cayley mode, the largest eigenvalues of are computed in order to obtain eigenvalues close to the shift sigma.
INVERT : this mode consists of transforming a generalized eigenvalue problem into a standard eigenvalue problem by solving mass matrix.

REGULAR mode is usually to avoid since it requires a lot of iterations in order to have all the needed eigenvalues. Therefore, for mass lumped elements, we advise to use SHIFTED, and for other elements, CAYLEY or INVERT mode can be used. shift is the frequency around which you want to know a given number of eigenvalues. shift can be equal to the following values:

LARGE : largest eigenvalues are searched
SMALL : smallest eigenvalues are searched
COMPLEX : clustered eigenvalues are searched close to a complex shift
shift : clustered eigenvalues are searched close to the given shift

Example :

# if you want 20 eigenvalues near frequency 0.5 for mass lumped elements
Eigenvalue = YES SHIFTED 0.5 20

# for 10 largest eigenvalues:
Eigenvalue = YES REGULAR LARGE 10

# for a real matrix and complex shift
Eigenvalue = YES COMPLEX_SHIFT COMPLEX 0.0 2.0 30

Location :

Solver/Preconditioner.cxx

Smoother

Syntax :

Smoother = type_smoother param_smoother

In this field, you specify the smoother used in multigrid preconditioning. Classical smoothers are available (JACOBI, BLOCK_JACOBI, SSOR), some equations can require a specific smoother like Maxwell's equations. The default smoother is Jacobi.

Example :

# Damped Jacobi with omega = 0.5
Smoother = JACOBI 0.5
# you can ask 2 iterations of Jacobi
Smoother = JACOBI 0.5 2
 
# similar for Block-Jacobi or Sor
Smoother = BLOCK_JACOBI 0.5
Smoother = SSOR 1.0

Location :

DampingParameters

Syntax :

DampingParameters = alpha delta

In this field, you specify the coefficient alpha used as damping for time-harmonic equations. delta is used in multigrid preconditioning if you want to decrease frequency on coarser grids. Default value of delta is 1. Default value of alpha is 1 for real problems and (1,0.5) for complex problems (advised value).

Example :

# example where the frequency will be multiplied by 0.7 for each coarser level
DampingParameters = (1.0,0.5) 0.7

Location :

Solver/Preconditioner.cxx

DirectSolver

Syntax :

DirectSolver = MUMPS
DirectSolver = ILUT threshold type_ilu level_fill droptol permtol alpha_ilu

In this field, you specify the direct solver used by the coarsest level in multigrid preconditioning. Usually, you have to choose between a complete factorization (using Mumps for instance) or an incomplete factorization. Incomplete factorization can lead to a substantial gain in storage. The default value is MUMPS. threshold is the value above which all the values of the initial matrix are removed (in absolute value), droptol is the threshold used by incomplete factorization (in relative value).

Example :

# for a complete factorization
DirectSolver = MUMPS

# for an incomplete factorization with epsilon = 0.01 (common threshold)
DirectSolver = ILUT 1e-5 0 100000 1e-2 0.1 1

Location :

Solver/Preconditioner.cxx

NumberMaxIterations

Syntax :

NumberMaxIterations = nb_max_iter restart

You can specify here the maximum number of iterations allowed for the iterative solution of finite element matrix, and the restart parameter (used by some Krylov solvers like GMRES and GCR). Default values are 1000 iterations allowed and a restart parameter equal to 10.

Example :

NumberMaxIterations = 1000 10

Location :

Solver/SolveSystem.cxx

Tolerance

Syntax :

Tolerance = tolerance

You can specify here the stopping criterion used for the iterative solution of finite element matrix. The default value is 1e-6.

Example :

Tolerance = 1e-6

Location :

Solver/SolveSystem.cxx

TypeResolution

Syntax :

TypeResolution = type_solver type_precond param_precond

You can specify here the iterative solver to use, preconditioning technique and optional parameters related to the preconditioner. In the table below, we present an exhaustive list of iterative solvers. We specify if the solver is able to solve symmetric complex, hermitian or unsymmetric matrices (all the solvers are of course able to solve symmetric real linear systems). Besides, some algorithms need the matrix-vector product with the transpose matrix. In the data file the name of the solver is written in capital letters.

Name of the solver	Symmetric complex	Hermitian	Unsymmetric	Transpose matrix
BiCg	CoCg	Cg	Yes	Yes
BiCgcr	Yes	No	No	No
BiCgStab	Yes	Yes	Yes	No
BiCgStabl	Yes	Yes	Yes	No
Cg	No	Yes	No	No
Cgne	Yes	Yes	Yes	Yes
Cgs	Yes	No	Yes	No
CoCg	Yes	No	No	No
Gcr	Yes	Yes	Yes	No
Gmres	Yes	Yes	Yes	No
Lsqr	Yes	No	No	No
MinRes	Yes	No	No	No
QCgs	Yes	No	Yes	No
Qmr	Qmr_Sym	Yes	Yes	Yes
Qmr_Sym	Yes	No	No	No
Symmlq	Yes	No	No	No
TfQmr	Yes	Yes	Yes	No

In this table, gray cells indicate which solver should be preferred for a symmetric or hermitian matrix. For example, if the matrix is symmetric, it is better to use CoCg instead of BiCg in order to save memory usage. Several preconditioning techniques are available (parameters are given in italic) :

IDENTITY : the preconditioning is the identity matrix. If selected, the iterative algorithm will be performed without preconditioning.
JACOBI omega nb_iter : the precondioning is damped Jacobi solver, omega is the relaxation parameter, nb_iter (optional parameter equal to 1) is the number of iterations
BLOCK_JACOBI size_block omega nb_iter : the precondioning is damped block-Jacobi solver. size_block is the size of blocks (dofs that are close are gathered), omega is the relaxation parameter, nb_iter (optional parameter equal to 1) is the number of iterations
SOR omega nb_iter : the preconditioning is S.O.R (Successive Over-Relaxations). omega is the relaxation parameter, nb_iter (optional parameter equal to 1) is the number of iterations
SSOR omega nb_iter : the preconditioning is S.S.O.R (Symmetric Successive Over-Relaxations). omega is the relaxation parameter, nb_iter (optional parameter equal to 1) the number of iterations.
ILUT type_facto lvl droptol pivot diag_coef : the preconditioning is Incomplete Factorization with threshold (ILUT). type_facto is equal to ILU_0, ILU_D, ILU_K or ILUT_K. lvl is the level k used by ILU(k) method. droptol is the dropping threshold, pivot the pivot threshold, diag_coef is used in ILU_D.
DIRECT : the preconditioning is a direct solver.
LOW_ORDER : the preconditioning is a direct solver for a finite element matrix constructed with first order and a subdivided mesh
MULTIGRID nb_cycles rmin nb_smooth_iter: the preconditioning is an iteration of a multigrid solution. nb_cycles is the number of cycles (1 for a V-cycle, 2 for a W-cycle). rmin is the minimal order for which a direct solver will be used. If the coarse order is less or equal to this order, a direct solver is used. nb_smooth_iter is the number of smoothing iterations.
SUBDOMAIN : the preconditioning is a basic subdomain preconditioning, each subdomain is associated with a processor.

Other preconditionings may be implemented for each specific solved equation. We have listed only preconditionings common for all equations.

Example :

# By default one iteration of Jacobi = diagonal preconditioning
TypeResolution = COCG JACOBI 1.0
# You can select more iterations (JACOBI omega nb_iter)
TypeResolution = COCG JACOBI 0.5 4
# or use Block-Jacobi algorithm (BLOCK_JACOBI size_block omega nb_iter)
TypeResolution = QMR_SYM BLOCK_JACOBI 10 0.5 1

# By default one iteration of SOR (SOR omega)
TypeResolution = BICGCR SOR 1.0
# You can select more iterations (SOR omega nb_iter)
TypeResolution = BICGCR SOR 1.5 4
# You can select more iterations or Symmetric SOR
TypeResolution = BICGCR SSOR 0.5 2

# Incomplete Factorization with threshold
# (ILUT threshold 0 max_elt droptol permtol alpha)
# max_elt is the maximum number of entries on a single row
# threshold is the value below which entries of original matrix are
# removed. droptol is the actual threshold used by incomplete factorization.
TypeResolution = GMRES ILUT 1e-5 0 100000 1e-2 0.1 1 

# preconditioning by low-order elements (ok for interpolatory basis)
TypeResolution = BICGSTAB LOW_ORDER

# p-multigrid preconditioning (MULTIGRID gamma pmin nb_iter)
# gamma = 1 -> V-cycle, 2 -> W-cycle
# pmin is the minimal order used when decreasing order for coarse meshes
TypeResolution = QMR MULTIGRID 2 2 1

Location :

Solver/SolveSystem.cxx
Solver/Preconditioner.cxx

TypeSolver

Syntax :

TypeSolver = AUTO
TypeSolver = DIRECT
TypeSolver = ITERATIVE
TypeSolver = MULTIGRID
TypeSolver = ILUT

This field overrides the field TypeResolution, i.e. this last field won't be taken into account if you have defined TypeSolver in another part of the datafile. The aim of this field is to simplify how the user selects the solver he wants to use, because TypeResolution is more complex to fill. AUTO will select the "best" linear solver known for this problem, well it is not yet correctly implemented. DIRECT will select a direct solver, Mumps if available. ITERATIVE will select an iterative solver compatible with the type of matrix but without preconditioning, while MULTIGRID will use a multigrid preconditioning with ad hoc parameters. ILUT proposes an incomplete factorization as preconditioning.

Example :

TypeSolver = AUTO

Location :

PivotThreshold

Syntax :

PivotThreshold = eps

This field is used if you want to change the threshold to use for partial pivoting in the direct solver.

Example :

PivotThreshold = 1e-3

Location :

NbThreadsPerNode

Syntax :

NbThreadsPerNode = n

This field is used if you want to change the number of threads that will be launched (usually during numerical factorization of the finite element matrix) for each MPI process.

Example :

NbThreadsPerNode = 8

Location :

EstimationConditionNumber

Syntax :

EstimationConditionNumber = YES nb_iter

If you are using a direct solver, you can estimate the condition number of the finite element matrix. We use the formula lambda_max/lambda_min even if the matrix is not symmetric. The estimate of extrema eigenvalues is achieved by using a power iteration method on the matrix and its factorization.

Example :

EstimationConditionNumber = YES 20

Location :

StaticCondensation

Syntax :

StaticCondensation = YES

If you are using a direct solver, you can decrease the required storage by performing a static condensation (internal degrees of freedom are eliminated). This is particularly efficient for higher method and in 2-D.

Example :

StaticCondensation = YES

Location :

ScalingMatrix

Syntax :

ScalingMatrix = YES

Before the factorization (or the iterative solution), the finite element matrix can be scaled if this field is set to YES. The default value is NO, no scaling is performed if not asked. It may improve the convergence and accuracy (due to round-off errors).

Example :

ScalingMatrix = YES

Location :

MumpsMemoryCoefficient

Syntax :

MumpsMemoryCoefficient = coef_init coef_max inc_coef

You can change the parameters used in Mumps to factorize the finite element matrix. coef_init is the initial ratio used to store the factorization (compared to the estimate). A value of 1.0 means that you will use the memory as estimated. Usually because of pivoting, the factors may take much more memory than initially estimated, that's why a larger initial ratio may save computational time. coef_max is the maximal ratio allowed (the computation will be stopped if the memory requirements are larger). If the memory is too small, Mumps will be called again while the ratio is lower than this value, the ratio will be multiplied by inc_coef at each iteration.

Example :

MumpsMemoryCoefficient = 2.0 100.0 2.0

Location :

Source/DefineSourceElliptic.cxx

TypeSource

Syntax :

TypeSource = type_source type_function param1 param2 ... paramn
TypeSource = num_source type_source type_function param1 param2 ... paramn

The first parameter is the type of the source to be chosen between the following types:

SRC_VOLUME param_f: the source will include a volume term of the type
where the function f is specified in the remaining parameters.
SRC_SURFACE ref param_f: for a Dirichlet condition, it will induce a hetereogeneous Dirichlet condition :
For other boundary conditions, the following surface term will be added to the source:
the surface Γ is specified through its reference (parameter ref after SRC_SURFACE) and the function f is described by remaining parameters.
SRC_DIRAC x0 y0 z0 Polarization p0 p1 ... pN: the source will include a Dirac term of the type (it is equivalent to the term in SRC_VOLUME where f is a Dirac function)
where p is the polarization of the Dirac. The origin of the Dirac x_s is set to the point given in OriginePhase if no parameters follow SRC_DIRAC. If there are parameters after SRC_DIRAC, they are assumed to store the coordinates of the source. If the parameter Polarization is present, it is followed by the polarization of the Dirac. If not present, p is initialized by the polarization given in the keyword Polarization
SRC_USER param_f : the source is provided in the file UserSource.cxx (class UserDefinedSource). If you modify this file, you will have to recompile Montjoie. The parameters are stored in the attribute source_space_param present in the class.
SRC_MODE ref param_mode: the source on a surface Gamma is a mode (e.g. computed by solving an eigenvalue problem). Since this type of source is very dependent on the solved equation, it is detailed for each solved equation.
SRC_TOTAL_FIELD param_incident: the total field is computed. There will be source terms in PML layers or in absorbing boundary conditions.
SRC_DIFFRACTED_FIELD param_incident: the diffracted field is computed. There will be source terms in hetereogeneities or on surfaces where boundary conditions are set (such as Dirichlet, Neumann or Impedance boundary conditions). Only absorbing boundary conditions do not include a source term.

SRC_TOTAL_FIELD or SRC_DIFFRACTED_FIELD are followed by parameters describing the incident field. The following incident fields are implemented:

PLANE_WAVE : incident plane wave
HANKEL_WAVE : incident Hankel wave
where r is the distance to the origin and h₀¹ is the Hankel function in 2-D and the spherical Hankel function in 3-D. In 3-D, we have :
GAUSSIAN_BEAM: incident gaussian beam
where k is the wave number and
where w is the waist. Here the wave vector is oriented along Oz. You can specify any incident angle (the expressions are just rotated).
LAYERED_PLANE_WAVE nb_layers d0 b0 a0 d1 b1 a1 ... dN : it is an incident plane wave adapted to a layered medium (with different wave numbers). It is implemented in 2-D only:
On each layer, the incident wave satisfies a Helmholtz equation
The coefficients A_j, B_j are obtained by solving a linear system (obtained with transmission conditions on each interface, and imposing a purely transmitted wave). The coefficients d_j, b_j and a_j are given in the parameters.

Some parameters can be provided after the description of the incident field, they are starting with a keyword and values:

Theta Theta : the incident angle θ, by default the value is the value given in the keyword IncidentAngle.
Phi Phi : the incident angle φ (3-D only), by default the value is the value given in the keyword IncidentAngle.
Origin x0 y0 z0 : the origin of the plane, by default the value is the value given in the keyword OriginePhase. In 2-D, only two values x0 and y0 must be provided.
Pulsation omega : the wave number, by default the value is the value given in the keyword Frequency.
Waist w : the waist of the gaussian, by default the value is equal to 0. As a result, this keyword is needed for gaussian beams (it is not used for other incident fields).

Example :

# angles theta and phi
IncidentAngle = 60.0 40.0

# Plane wave incident field
TypeSource = SRC_DIFFRACTED_FIELD PLANE_WAVE

# origin of the incident field
OriginePhase = 0.0 2.0 0.0

# polarization of the incident field
Polarization = 0.0 1.0 0.0

# incident field with Hankel function
TypeSource = SRC_DIFFRACTED_FIELD HANKEL_WAVE Origin -5.0 0.0 0.0

# incident field with a gaussian beam
TypeSource = SRC_TOTAL_FIELD GAUSSIAN_BEAM Waist 6.4

# Layered incident wave 
# LAYERED_PLANE_WAVE nb_layers d0 b0 a0 d1 b1 a1 ... dN
# here two layers at with d = [-0.5, 0.2, 2.5]
TypeSource = SRC_TOTAL_FIELD LAYERED_PLANE_WAVE 2 -0.5  2.0 0.8  0.2 3.0 1.2  2.5

If you choose SRC_TOTAL_FIELD, SRC_DIFFRACTED_FIELD or SRC_USER, only a line TypeSource is necessary to describe the source. For other sources (ie SRC_VOLUME, SRC_SURFACE, SRC_DIRAC and SRC_MODE) , several lines TypeSource can describe the same source. We consider that terms are added

For SRC_VOLUME, the function f can be chosen among the following list:

GAUSSIAN x0 y0 (z0) r1 r2 : a gaussian centered at point (x0, y0, z0) (z0 is not present in 2-D) and with a distribution radius equal to r1. The gaussian is equal to:
where r is the distance to the center. The parameter alpha is set such that
r2 is a cut-off radius. If not provided, it will be set to r1. The parameter beta is given as
UNIFORM : an uniform source (f=1). The polarization can be set by writing Polarization followed by values. The source f is constant only on the specified reference, such that you can specify a piecewise constant source if the mesh contains different physical references.

VARIABLE_BINARY : a variable source. The source is provided directly on quadrature points of the mesh with a binary file. Quadrature points must be written by specifying WriteQuadraturePoints=YES. Then, you can write the source in a binary as follows (in Python, by including visuND package):

# first you read the points contained in file quadrature_points.dat
PT = loadtxt('../quadrature_points.dat')
x = PT[:, 0]
y = PT[:, 1]

# then you evaluate your source on these points
f = exp(-5.0*(x*x + y*y))

# and write the source in a binary file
write_full("../binary_source.dat", f)

# if the source has different components, they are written in the same file
g = zeros(nb_components*len(x))
f1 = sin(x)*sin(y)
f2 = cos(x)*cos(y)
...
g[0:len(x)] = f1
g[len(x):2*len(x)] = f2
...
write_full("../binary_source.dat", g)

GRADIENT_BINARY : a variable source to be integrated against gradient of basis functions (or curl for Maxwell's equations), i.e:
This source has to be evaluated on quadrature points of the mesh in binary (as explained for the previous item VARIABLE_BINARY.

For SRC_SURFACE, the function f can be chosen among the following list:

GAUSSIAN x0 y0 (z0) r1 r2 : a gaussian centered at point (x0, y0, z0) (z0 is not present in 2-D) and with a distribution radius equal to r1. The gaussian is the same as described above.
UNIFORM : an uniform source (f=1). The polarization can be set by writing Polarization followed by values.
INCIDENT_WAVE : the function f is an incident wave, the parameters are the same as those used to described an incident field (see above).
VARIABLE file_points file_value : the function f is variable and provided through a list of points and a list of values. The points are written in a text file while values are written in binary. Below, we list Python commands to write these two files (by including visuND package) :
```
# exemple of a source on boundary y = -2
PT = zeros([N, 2])
PT[:, 0] = linspace(-2.0, 2.0, N)
PT[:, 1] = -2.0*ones(N)
savetxt("Points.txt", PT)

A = zeros([N, 1]) + 1j*zeros([N, 1])
A[:,0] = exp(1j*PT[:,0])
write_full("Values.dat", A)
```
the matrix A contains values. Each column corresponds to a different component of f. For scalar equations (such as Helmholtz equation), there is only one component to provide.

Example :

# Volume source : GAUSSIAN x0 y0 radius (in 2-D)
TypeSource = SRC_VOLUME GAUSSIAN 0 0 1.0

# Volume source : f = 1 on reference 1 and 0 otherwise
TypeSource = SRC_VOLUME 1 UNIFORM Polarization 1.0

# Volume source : f given directly on quadrature points
# syntax SRC_VOLUME ref VARIABLE_BINARY values.dat
# non null only on reference 2 here
TypeSource = SRC_VOLUME 1 VARIABLE_BINARY binary_source.dat

# case where we compute integral against gradient of basis functions
TypeSource = SRC_VOLUME 1 GRADIENT_BINARY binary_source.dat

# Surface source on reference 2 :
# Polarization can be used to set the amplitude for scalar equations
# TypeSource = SRC_SURFACE 2 GAUSSIAN x0 y0 r1 r2 Polarization p0 ... pN
TypeSource = SRC_SURFACE 2 GAUSSIAN -2.0 0.4  1.0 1.5 Polarization 1.0

# Other surface source on reference 3
TypeSource = SRC_SURFACE 3 UNIFORM Polarization 2.0

# other source for reference 4
TypeSource = SRC_SURFACE 4 INCIDENT_WAVE PLANE_WAVE Theta 60.0 Polarization 0.4

# variable source for reference 5
TypeSource = SRC_SURFACE 5 VARIABLE Points.txt Values.dat

# the right hand side computed is obtained by summing all contributions

Location :

ThresholdRhs

Syntax :

ThresholdRhs = epsilon

The threshold for right hand sides is used to avoid computing the solution of linear systems when the right hand side is considered as negligible. For such linear systems, the solution will be set to 0. The default value is the machine precision (1e-16 for double precision).

Example :

ThresholdRhs = 1e-100

Location :

Source/DefineSourceElliptic.cxx

InitialCondition

Syntax :

InitialCondition = USER param1 param2 ... paramn

No default initial condition is implemented in Montjoie, except null initial conditions. However, you can enter the initial condition by modifying the file UserSource.cxx, but you will need to recompile Montjoie for each modification of the initial condition. Nevertheless, you can provide parameters after USER, so that you can reuse those parameters in the file UserSource.cxx.

Example :

# after USER you can give any number of parameters
InitialCondition = USER 2.0

Location :

TemporalSource

Syntax :

TemporalSource = type_source parameters

Most often in Montjoie, we consider source with separated variables, i.e. obtained by multiplying a function F(t) depending on time with a function g(x) depending on space. For example, for the wave equation, we consider the first-order system:

The space source g(x) is defined by the keyword TypeSource while the time source F(t) is defined by the keyword TemporalSource. F(t) is usually the primitive of a function f(t), since the equivalent second-order formulation of the wave equation reads:

As a result, the obtained solution will be the same with schemes adapted to the first-order system or with schemes adapted to the second-order equation. Some predefined functions are available to specify F(t), they use the frequency defined by the keyword Frequency. You can also provide your own time source, either with an ASCII file containing a list of values, or by modifying the file UserSource.cxx. There is no default value, null source will be computed if you don't insert this field. The following time sources are available (parameters are written in italic) :

RICKER : a Ricker (i.e. second-derivative of a gaussian) :
The frequency f₀ is the central frequency of the Ricker
MODIFIED_RICKER t0 : a Ricker with a center time t₀
The center time t₀ is provided as the first parameter.
GAUSSIAN w0 : a gaussian with a center plateau of width w₀ (initial time t₀ is equal to 0)
HARMONIC T0 : a sine function with smooth cut-off, T₀ is the final time of the source.
SINUS_GAUSSIAN alpha : a gaussian (with parameter alpha) modulated by a sine function:
The time t₀ is initialized such that
FILE file_name t0 tf : the function F(t) is given with values contained in a file. The first parameter is the file name, the two other parameters are the initial time and final time. The values of the source are assumed to be known at equidistributed times :
where
N is the number of values stored in the file. The file contains an unique column with only values (in text file) :
F(t) is then computed for any time t by a cubic interpolation.
USER parameters : F(t) is given by the user in the file UserSource.cxx (method Evaluate of the class TimeUserSource). The parameters can be used in the constructor of the class TimeUserSource.

For the wave equation (targets acous2D and acous3D), an incident plane wave can be specified (TypeSource = SRC_DIFFRACTED_FIELD or SRC_TOTAL_FIELD). In that case, the source is not written in separated variables. The incident plane wave is given as

where F is the time source (as specified by the keyword TemporalSource), is the wave vector (given through the frequency and the incident angles), c₀ is the velocity of the media at infinity. The parameter alpha is equal to:

if AUTO has been selected in TypeSource:

TypeSource = SRC_TOTAL_FIELD AUTO

otherwise it is equal to the parameter given by the user (instead of AUTO). The source is deduced from the definition of the incident wave (for a total field, it implies a source in PML or ABC).

Example :

# a gaussian with a plateau of width=0.2
TemporalSource = GAUSSIAN 0.2

# exp(-alpha (t-t0)^2) * sin (2 pi f t)
TemporalSource = SINUS_GAUSSIAN 0.2

# values given on a file
TemporalSource = FILE source.dat 0 10

Location :

TimeStep

Syntax :

TimeStep = dt

The time step is provided in this field, otherwise the default value is equal to 0.01. If dt is set to AUTO, the CFL will be computed by evaluating the maximal eigenvalue of M^-1 K (where M is the mass matrix and K the stiffness matrix). This operation is expensive, but sometimes it can be interesting to have a good estimate of the CFL. Once the CFL is evaluated, the time step is then chosen close to the CFL limit (0.999 CFL).

Example :

TimeStep = 0.1

# automatic choice of TimeStep by evaluating CFL
TimeStep = AUTO

# you can also choose the time duration and the number of iterations
# the time step will be equal to   Time Duration / Number iterations
# TimeStep = DIVISION Tfinal nb_iterations
TimeStep = DIVISION 10.0  200

Location :

TimeInterval

Syntax :

TimeInterval = t0 tf

You can specify the final time tf and initial time t0 of the experience. Default values are set to 0.

Example :

TimeInterval = 0 100

Location :

RandomInitialCondition

Syntax :

RandomInitialCondition = Small

You can specify that the initial condition is set with small random values (usually of 1e-270). It happens that the simulation can be slow at the beginning because of a null initial condition. Indeed, denormal values appear and slow the computation. To overcome this problem, you can start from a non-null initial condition with these small random values. Another possibility is to compile with some options in order to prohibit denormal values:

// in your program, you add these lines
#include <xmmintrin.h>

#define CSR_FLUSH_TO_ZERO         (1 << 15)

// in the main program, you have to add these lines
int main()
{
   unsigned csr = __builtin_ia32_stmxcsr();
   csr |= CSR_FLUSH_TO_ZERO;
   __builtin_ia32_ldmxcsr(csr);

}

Example :

RandomInitialCondition = Small

Location :

DisplayRate

Syntax :

   DisplayRate = N
   DisplayRate = N YES

The norm of the solution is display each N iterations. You can set N with this field. By default, N depends on the print level chosen. If the second optional parameter is present, the computational time is also displayed.

Example :

# each ten iterations, the norm of the solution is displayed
DisplayRate = 10  

# you can also display computational time
DisplayRate = 10 YES

Location :

LoadReprise

Syntax :

LoadReprise = YES

If asked, the code will try to read the data saved in a previous simulation and resume the simulation from the last iterate saved.

Example :

LoadReprise = YES

Location :

SaveReprise

Syntax :

SaveReprise = nb_iter

If present, the code will save data needed to resume the simulation later. The data will be saved each nb_iter iterations. By default, the code does not save anything, and the simulation can not be resumed.

Example :

# data is written on the disk each 100 iterations
SaveReprise = 100

Location :

PathReprise

Syntax :

PathReprise = directory

The user can provide the directory where datas needed to resume a simulation should be written. It can be interesting to specify a local directory in order to have a fast writing.

Example :

# temporary data is written on /scratch
PathReprise = /scratch

Location :

NormeMaxSolution

Syntax :

NormeMaxSolution = value

The user can provide a maximal norm allowed for the solution. This value is used to determine if the time scheme is stable or not. If the norm of the solution is greather than this value, the simulation is stopped (in order to avoid slow computations with NaNs). By default, the maximal norm allowed is set to 1e30.

Example :

# If you want to specifiy a different maximal norm
NormeMaxSolution = 100.0

Location :

Instationary/TimeSchemes.cxx

OrderTimeScheme

Syntax :

OrderTimeScheme = order type_scheme

Here you provide the time scheme and order. Available time schemes are :

Scheme	Available orders	First-order Formulation	ABC/Damping	Stable
LEAP_FROG_TRUE	2	Yes	No	CFL
LEAP_FROG	2 and 4	No	Yes	CFL
OPTIMAL_LEAP_FROG	2 r	No	Yes	CFL
LEAP_FROG_PML	2 and 4	Yes	Yes	CFL
OPTIMAL_LEAP_FROG_PML	2 r	Yes	No	CFL
NYSTROM	2 r	Yes	No	No
RUNGE_KUTTA_NYSTROM	from 2 to 12	No	No	CFL
SYMMETRIC_MULTISTEP	2 r	No	No	CFL
IMPLICIT_SYMMETRIC_MULTISTEP	2 r	No	No	CFL
ADAMS_BASHFORTH	r	Yes	Yes	CFL
ADAMS_BASHFORTH_MOULTON	r	Yes	No	No
TALEZER	r	Yes	Yes	CFL
RUNGE_KUTTA	from 1 to 8	Yes	Yes	CFL
RUNGE_KUTTA_LOW	4	Yes	Yes	CFL
GAUSS_RUNGE_KUTTA	2 r	Yes	Yes	Always
THETA_SCHEME	2	No	Yes	Always
SDIRK	From 1 to 5	Yes	Yes	Always
DIRK	4	Yes	Yes	Always
CRANK_NICOLSON	2	Yes	Yes	Always
ADAMS_IMPLICIT	r	Yes	Yes	No
MILNE_SIMPSON	r	Yes	Yes	No
BDF	r	Yes	Yes	No
PADE	2 r	Yes	Yes	Always

Some schemes can be unavailable for a given equation. The default value is second order leap frog scheme. Some time schemes need additional parameters, time schemes are detailed in this section.

Example :

TimeScheme = 2 LEAP_FROG_PML
TimeScheme = 4 LEAP_FROG
TimeScheme = 4 RUNGE_KUTTA_LOW

# you need to specify teta here
TimeScheme = 2 THETA_SCHEME 0.25

Location :

ParametersRCS

Syntax :

ParametersRCS = YES ref AUTO MONOSTATIC
ParametersRCS = YES ref order BISTATIC

The computation of radar cross section can be done for Helmholtz and time-harmonic Maxwell equations. Monostatic or bistatic radar cross sections are available. The surface where the computation of the integral is performed is defined by the number ref. This number is not a reference number, but a body number, that can include several reference numbers. order is the order of integration used to evaluate integrals. If AUTO is selected, the order will be taken equal to the order of approximation. The default value is NO (no radar cross section is computed).

Example :

# Body is defined with TypeBody
# TypeBody = ref num_body
# here body 1 is formed with references 2 and 3
TypeBody = 2 1
TypeBody = 3 1

# ParameterRCS = YES num_body AUTO type_rcs
ParametersRCS = YES 1 AUTO BISTATIC

# AngleRCS = teta0 teta1 nb_angles
AngleRCS = 0.0 180.0 181

# file where rcs will be stored
FileRCS = RcsDisque.dat

Location :

AngleRCS

Syntax :

AngleRCS = teta0 teta1 nb_points

Usually radar cross section is computed between two angles teta0 and teta1 and with a given number of angles (nb_points). The angles should be expressed in degrees. The default values are 0, 180 degrees and 181 points (so that the RCS is computed at every degree).

Example :

AngleRCS = 0.0 180.0 361

Location :

FileRCS

Syntax :

FileRCS = file_name
FileRCS = file_name file_far_field

The radar cross section will be saved in this file. The default value is "rcs.dat" . If a second argument is specified, the far field is also stored (the far field is the complete vector for Maxwell's equations).

Example :

# by default only radar cross section will be stored
FileRCS = Rcs.dat

# you can give a file name for the far field
FileRCS = Rcs.dat FarField.dat

Location :

TransparencyCondition

Syntax :

TransparencyCondition = YES ref AUTO
TransparencyCondition = YES ref order
TransparencyCondition = YES ref order ref_domain

For Hemholtz and time-harmonic Maxwell equations, it is possible to use a transparent condition (instead of PML or ABC), which is exact but relies on the representation formula to iterate on the solution. This transparent condition is detailed in M. Durufle's thesis (and in a paper of C. Hazard and M. Lenoir). For the representation formula, a direct integration is used, no multipole strategy is exploited. ref is a body number, you can specify an order of integration, AUTO correspond to selecting the order of approximation. ref_domain is an optional parameter that gives the reference of the volume contained between the artificial boundary and the extern boundary.

Example :

# Body is defined with TypeBody
# TypeBody = ref num_body
# here body 1 is formed with references 2 and 3
TypeBody = 2 1
TypeBody = 3 1

# TransparencyCondition = YES num_body order
TransparencyCondition = YES 1 AUTO

Location :

ParamResolutionTransparency

Syntax :

ParamResolutionTransparency = type_resolution nb_iterations tolerance restart

The use of transparent condition requires an iterative solver. type_resolution informs the type of iterative solver, nb_iterations the maximum number of iterations, tolerance the stopping criterion and restart the restart parameter (for GMRES for example). Be careful to select an iterative solver which works for unsymmetric linear systems. The default solver is GMRES(10) with a maximum of 100 iterations and a stopping criterion equal to 1e-6. The number of iterations is usually small, especially when the distance between the extern boundary and the intern boundary (on which representation formula is computed) is large.

Example :

# ParamResolutionTransparency = algo nb_iter_max  epsilon  restart_parameter
ParamResolutionTransparency = GMRES 100 1e-6 10

Location :

TypeBody

Syntax :

TypeBody = ref body1 body2 ... bodyn

You can group surfaces of different references into a single body (with a unique number). This body number is then used for the computation of radar cross section or for the transparent condition. By default, no body is created.

Example :

# Body 2 includes surfaces of reference 2 and 3
# Body 1 includes surfaces of reference 1 and 3
TypeBody = 1  1
TypeBody = 2  2
TypeBody = 3  1 2

Location :

Mesh/MeshBoundaries.cxx

ExplicitMatrixFEM

Syntax :

ExplicitMatrixFEM = YES_NO_or_AUTO

If the value is YES, the finite element matrix will be effectively computed even though it is not needed. If the value is NO, and if an iterative solver (for time-harmonic simulations) is selected or for explicit time-stepping (for time-domain simulations), the finite element matrix is not computed. The matrix-vector product is made with a matrix-free implementation. Usually only the jacobian matrices DF_i are stored and used for this matrix-vector product. If the value is AUTO, the code will choose if the finite element matrix is stored or not.

Example :

# if you want that the finite element matrix is computed
ExplicitMatrixFEM = YES

Location :

StorageMatrix

Syntax :

StorageMatrix = YES nom_fichier

During the computation, the finite element matrix can be stored in a file, so that you can check that the computation is correctly done. It is better to use a direct solver if you want to store the matrix in a file, because if you choose an iterative solver, this will probably perform a matrix vector product will all canonical vectors (1,0,...0), (0,1,0, ...,0), ..., (0,0, ...,1). And this can be very long if the matrix is large. But this is also a way to check that the matrix-free implementation is correctly implemented for the finite element and equation you consider. If the file name ends with mtx, the matrix will be stored in the MatrixMarket format. The default value is set to NO.

Example :

#if you want to store the matrix in the file mat.dat
StorageMatrix = YES mat.dat

# to write the FEM matrix in MatrixMarket format, use the extension mtx
StorageMatrix = YES mat.mtx

Location :

ThresholdMatrix

Syntax :

ThresholdMatrix = tolerance

If you hope to obtain a sparser matrix, you can set a threshold, and all the entries of the matrix smaller than this threshold won't be added to the matrix during assembly phase.

Example :

ThresholdMatrix = 1e-6

Location :

PrintLevel

Syntax :

PrintLevel = level

This entry can be modified if you want to increase or decrease the level of displayed messages. A value of -1 will display nothing (silent mode). The default value is 0. The maximum value is 10 but displays a lot of informations and writes stuff on the disk. An intermediate value is 6 and is better in a normal use.

Example :

PrintLevel = 2

Location :

WriteQuadraturePoints

Syntax :

WriteQuadraturePoints = YES/NO
WriteQuadraturePoints = YES NO_PML_POINTS

If you specify YES, quadrature points will be written in the file quadrature_points.dat. If the simulation is launched with several processors, the quadrature points of each processor will be written in a separated file : quadrature_points_P0.dat, quadrature_points_P1.dat, etc. The second argument NO_PML_POINTS is optional. If present, quadrature points of elements in PML elements will not be written.

Example :

WriteQuadraturePoints = YES

Location :

Elliptic/AeroAcoustic/AeroAcoustic.cxx

EnergyConservingAeroacousticModel

Syntax :

  EnergyConservingAeroacousticModel = type_model

This entry sets which model is used for aeroacoustics (target aero2D, harmonic_aero2D). The following models are available:

YES : conservative model (an energy is conserved)
BogeyBaillyJuve : Bogey-Bailly-Juve model
Galbrun : simplified Galbrun's equation
NO : simplified Linearized Euler's equation

Example :

EnergyConservingAeroacousticModel = Galbrun

Location :

DropUnstableTerms

Syntax :

  DropUnstableTerms = Convective
  DropUnstableTerms = Convective alpha
  DropUnstableTerms = NonUniform

This entry can be used to drop some terms for target galbrun_axi

Example :

DropUnstableTerms = Convective 0.2

Location :

Elliptic/AeroAcoustic/AxiSymGalbrun.cxx

ApplyConvectiveCorrectionSource

Syntax :

ApplyConvectiveCorrectionSource  = YES

This entry can be used to apply convective derivative to the volume source:

Example :

ApplyConvectiveCorrectionSource  = YES

Location :

Elliptic/AeroAcoustic/AxiSymGalbrun.cxx

ForceFlowNeumann

Syntax :

  ForceFlowNeumann = YES/NO

If specified to YES, the Neumann condition will be enforced. Usually, Neumann condition will be weakly satisfied by the steady flow (computed by a Discontinous Galerkin method for example). As a result, Neumann condition is approximately verified, and this can produce numerical instabilites, that are removed by enforcing a strongly Neumann condition. If no entry is provided, this operation is not performed, default value is NO.

Example :

ForceFlowNeumann = YES

Location :

Elliptic/AeroAcoustic/HarmonicGalbrun.cxx

NumberModes

Syntax :

  NumberModes = AUTO tolerance
  NumberModes = AUTO tolerance nb_min_modes
  NumberModes = SINGLE num_mode
  NumberModes = nb_modes

For the finite element solution in an axisymmetric configuration, each Fourier mode is solved separately. The selection of Fourier modes to solve can be automatically performed by specifying tolerance. When the norm of the right hand side is lower than tolerance (in relative value), the computation is stopped. You can also specify a maximum number of modes nb_modes, and all the modes between -nb_modes and nb_modes will contribute to the solution. Default value is automatic selection of number of mondes with tolerance equal to 1e-6.

Example :

# computation will be stopped when the norm of x_m will be smaller than 1e-12 max || x_i ||
NumberModes = AUTO 1e-12

# you can enforce the computation of first modes (in order to compute correctly the maximum)
NumberModes = AUTO 1e-12 8

# you can specify the number of modes to compute
NumberModes = 20

# you can ask the solution for a single mode
NumberModes = SINGLE 1

Location :

Harmonic/VarAxisymProblem.cxx

ForceDiagonalMass

Syntax :

ForceDiagonalMass = YES

For axisymmetric Helmholtz equation, the mass matrix can be lumped by using Gauss-Radau points. This functionality can be triggered with this keyword.

Example :

# diagonal mass matrix
ForceDiagonalMass = YES

Location :

Harmonic/VarAxisymProblem.cxx

AddFlowTerm

Syntax :

AddFlowTerm = YES

For Helmholtz equation, a convection term can be added if this keyword is present. If the value is set to YES, Helmholtz equation is given as:

The default value is NO (the flow M is assumed to be zero).

Example :

# first, you need to say that a convection term is present
AddFlowTerm = YES

# then you can specify a flow for the different domains
PhysicalMedia = M 1 ISOTROPE 0.5 0.0

Location :

Elliptic/Helmholtz/VarHelmholtz.cxx

AddSlot

Syntax :

  AddSlot = x0 y0 x1 y1 epsilon N

For 2-D Helmholtz equation, if MONTJOIE_WITH_THIN_SLOT_MODEL has been defined in the compiled file (or with faster compilation), you can specify a thin slot in the computational domain. The thin slot extends between points [x0, y0] and [x1,y1], with a thickness epsilon.

First-order models are used, a Dirichlet-to-Neumann model or a 1-D model (which avoids the problem of resonant frequencies). For more details about the used models, you can read S. Tordeux's thesis. N denotes the number of points to discretized the 1-D solution along the slot (used in the case of 1-D model). There is no default value since no slot is included at the beginning.

Example :

# specify a slot of extremities [0,0] [2,0] and of thickness 0.1
AddSlot = 0 0 2 0 0.1 200 
# you can add as many slots as you want with different parameters
AddSlot = -1 -2 1 -2 0.2 100

Location :

Elliptic/Helmholtz/VarHelmholtz.cxx

ModelSlot

Syntax :

  ModelSlot = DTN/MESH_1D

For 2-D Helmholtz equation, if MONTJOIE_WITH_THIN_SLOT_MODEL has been defined in the compiled file, you can specify a thin slot in the computational domain. First-order models are used, a Dirichlet-to-Neumann model or a 1-D model (which avoids the problem of resonant frequencies). For more details about the used models, you can read S. Tordeux's thesis. DTN corresponds to the Dirichlet-to-Neumann model while MESH_1D corresponds to the 1-D model. Default value is MESH_1D.

Example :

ModelSlot = DTN

Location :

Elliptic/Helmholtz/VarHelmholtz.cxx

CalculEnveloppe

Syntax :

  CalculEnveloppe = YES

For axisymmetric Helmholtz equation (target helm_axi), you can ask to compute the enveloppe, i.e. the field v such that :

It can be interesting if the wave propagates mainly in z-direction

Example :

CalculEnveloppe = YES

Location :

Elliptic/Helmholtz/AxiSymHelmholtz.cxx

FormulationAxisymmetric

Syntax :

  FormulationAxisymmetric = R3

For axisymmetric Helmholtz equation, the usual variational formulation (without flow term) is equal to :

There is no singularity because an homogeneous Dirichlet condition is set on the axis for m different from 0. If the formulation R3 is selected, we replace basis functions as:

It implies the following variational formulation

This formulation is denoted R3 because of the presence of r³ weight for the mass matrix. With this variational formulation, Dirichlet condition is no longer enforced, and we have only polynomials in r.

Example :

# if you want to use the r^3 formulation
FormulationAxisymmetric = R3

Location :

Elliptic/Helmholtz/AxiSymHelmholtz.cxx

TimeReversal

Syntax :

TimeReversal = DIRECT t0 tf dt fichier_maillage fichier_du_dn ref
TimeReversal = INVERSE t0 tf dt fichier_maillage fichier_du_dn ref

This entry gives the possibility to simulate time reversal experiences. It is currently only available for 2-D acoustic equation, but the extension would be straightforward for other equations. If you specify DIRECT, the values of du/dn are computed on a boundary of reference ref and stored in the file fichier_du_dn each dt. Then, you can use these datas (modify them if you want), to simulate a time reversal (INVERSE). The same parameters should must be kept (t0, tf, ...). fichier_maillage contains a mesh of the boundary. The reference can change, since the mesh can be different in direct and inverse experiences. There is no default value, time reversal is optional.

Example :

TimeReversal = DIRECT 10.0 20.0 0.1 boundary.mesh dudn.dat 2
TimeReversal = INVERSE 10.0 20.0 0.1 boundary.mesh dudn.dat 4

Location :

Hyperbolic/Acoustic/TimeReversal.cxx

FileCoefficientsQ

Syntax :

FileCoefficientsQ = file_name

This entry allows the user to compute coefficients Q_i:

for each element i. These coefficients are written in an output file. By default, such coefficients are not computed.

Example :

# if you want to compute q_i and store them:
FileCoefficientsQ = file_Qi.dat

Location :

Elliptic/Maxwell/3D/HarmonicMaxwell3D.cxx

ModifiedFormulation

Syntax :

ModifiedFormulation = YES

This entry enables the user to use a modified formulation for axisymmetric Maxwell's equations. Usually, this formulation is more accurate. By default, this formulation is not selected.

Example :

ModifiedFormulation = YES

Location :

Elliptic/Maxwell/Axi/AxiSymHcurlMaxwell.cxx

OutputHy

Syntax :

OutputHy = ref0 .. refN PARAM file_name n

By inserting this line, the user asks to write H_y on a set of referenced curves. n is the number of subdivisions for each interval.

Example :

# on edges of reference 2 and 3, Hy is computed with 10 subdivisions for each edge
OutputHy = 2 3 PARAM Hy.dat 10

Location :

Elliptic/Maxwell/Axi/MaxwellAxiSymHarmonic.cxx

DisplayStress

Syntax :

DisplayStress = NO

For elastodynamics equations, it may be interesting to display the stress σ :

when the gradient is required. If NO is selected, the gradient of u is displayed, otherwise the stress is displayed. The default value is YES, if this line is not present the stress is displayed.

Example :

# if you want to know the gradient of u, and not the stress :
DisplayStress = NO

Location :

Elliptic/Elastic/VarElastic.cxx

ThicknessPlate

Syntax :

ThicknessPlate = delta
ThicknessPlate = VARIABLE ref param_media

The thickness of the plate is specified, otherwise the default value is 0.1. It is only used for Reissner-Mindlin equations. For a variable thickness the parameters are detailed in the description of MateriauDielec.

Example :

ThicknessPlate = 0.1

# for a variable thickness
ThicknessPlate = VARIABLE 1 SAME_MESH test_plaque.mesh 0.0 1.0 delta.don 8 0.02

Location :

Hyperbolic/Elastic/ReissnerMindlin.cxx

SismoPlane

Syntax :

# for 2-D and axisymmetric targets
SismoPlane = AUTO nx ny
SismoPlane = xmin xmax ymin ymax nx ny

# 3-D targets
SismoPlane = x0 y0 z0 xA yA zA xB yB zB nx ny

Solution can be computed on a plane, on the rectangle [xmin, xmax] x [ymin, ymax] in 2-D, and on the parallelogram defined by the origin (x0, y0, z0) and the two extremities (xA, yA, zA), (xB, yB, zB) in 3-D. nx and ny are the number of points along direction x and direction y. By default, no output is done on planes. In 2-D, AUTO means that the code will use the bounding box of the physical domain (PML layers are excluded).

Example :

# for regular grids in 2-D
SismoPlane = AUTO 101 101
SismoPlane = -5 5 -5 5  101 101

# for 3-D plane {0} x [-1, 1] x [0 2]
SismoPlane = 0 -1  0   0 1 0   0 -1  2  101 101

Location :

SismoPlaneAxi

Syntax :

# axisymmetric targets
SismoPlaneAxi = x0 y0 z0 xA yA zA xB yB zB nx ny

Solution can be computed on a 3-D parallelogram defined by the origin (x0, y0, z0) and the two extremities (xA, yA, zA), (xB, yB, zB) in 3-D. nx and ny are the number of points along direction x and direction y.

Example :

# for output on a 3-D plane {0} x [-1, 1] x [0 2] for axisymmetric case
SismoPlaneAxi = 0 -1  0   0 1 0   0 -1  2  101 101

Location :

SismoGrille

Syntax :

# output on three planes intersecting at point (x0, y0, z0)
SismoGrille = x0 y0 z0 AUTO nx ny nz
SismoGrille = x0 y0 z0 xmin xmax ymin ymax zmin zmax nx ny nz

3-D solutions can be computed on three planes : {x0} x [ymin, ymax] x [zmin,zmax], and [xmin, xmax] x {y0} x [zmin, xmax] and [xmin, xmax] x [ymin,ymax] x {z0}. The number of points on each direction are nx, ny and nz. If AUTO is specified, the code uses the bounding box of the physical domain (PML layers are excluded).

Example :

# AUTO : bounding box is the bounding box of the mesh
SismoGrille = 0 0 0 AUTO 101 101 101

Location :

SismoLine

Syntax :

# output on the line joining A to B
SismoLine = xA xB yA yB nb_points

# For 3-D
SismoLine = xA xB yA yB zA zB nb_points

Output on lines of extremities (xA, yA, zA) and (xB, yB, zB) can be performed, nb_points is the number of points on the line.

Example :

# Line x = 2, y in [0, 10]
SismoLine = 2 2 0 10 201

Location :

SismoLineAxi

Syntax :

# output on a 3-D line
For 3-D
SismoLine = xA xB yA yB zA zB nb_points

Output on 3-D lines of extremities (xA, yA, zA) and (xB, yB, zB) can be performed, nb_points is the number of points on the line. This field is meaningful only for axisymmetric targets.

Example :

# Line x = 2, y in [0, 10], z in [-4, 4]
SismoLineAxi = 2 2 0 10 -4.0 4.0 201

Location :

SismoCircle

Syntax :

# 2-D circles
SismoCircle = x0 y0 radius nb_points

# 3-D circle
SismoCircle = x0 y0 z0 nx ny nz ra rb nb_points

In 2-D, you specify the center and radius of the circle. In 3-D, you specify the center (x0, y0, z0) of the circle, the normale (nx, ny, nz) of the plane containing the circle, and ra, rb the two radii (for an ellipse). nb_points is the number of points in the circle.

Example :

# output on a circle of center (0, 0) and radius 3 with 201 points
SismoCircle = 0 0  3.0  201

Location :

SismoCircleAxi

Syntax :

# 3-D circle
SismoCircleAxi = x0 y0 z0 nx ny nz ra rb nb_points

The solution can be computed on points on a 3-D ellipse. This field is meaningful only for axisymmetric targets. nx, ny, nz denotes the normale to the plane where the circle lies. ra and rb are two radii.

Example :

SismoCircleAxi = 0 0 0  1 0 0  3.0 5.0  201

Location :

SismoPoint

Syntax :

# 2-D point
SismoPoint = x0 y0

# 3-D point
SismoPoint = x0 y0 z0

You can compute the solution on a single point ... This is actually useful for time-domain simulations in order to obtain a seismogramm. In this case, the output is written in the SAME file and in ASCII format so that you can directly read the file in Matlab by using command load (or in gnuplot). If you insert several SismoPoint, the solution will be computed on all the points given, as explained in FileOutputPoint.

Example :

SismoPoint = 2 0

Location :

SismoPointAxi

Syntax :

# 3-D point
SismoPointAxi = x0 y0 z0

You can compute the solution on a single 3-D point ... This is actually useful for time-domain simulations in order to obtain a seismogramm. In this case, the output is written in the SAME file and in ASCII format so that you can directly read the file in Matlab by using command load (or in gnuplot). If you insert several SismoPoint, the solution will be computed on all the points given, as explained in FileOutputPoint. This field is meaningful only for axisymmetric targets.

Example :

SismoPointAxi = 2 0 0.8

Location :

SismoPointsFile

Syntax :

SismoPointsFile = nom_fichier

The solution can be computed on points defined by the user. In that case, you have to provide a file in ASCII format, like that :

x0 y0 z0
x1 y1 z1 
...

In 2-D, of course you have only two coordinates to give.

Example :

SismoPointsFile = points.dat

Location :

SismoPointsFileAxi

Syntax :

SismoPointsFileAxi = file_name

The solution can be computed on 3-D points defined by the user. In that case, you have to provide a file in ASCII format, like that :

x0 y0 z0
x1 y1 z1 
...

This field is meaningful for axisymmetric targets only.

Example :

SismoPointsFileAxi = points.dat

Location :

SismoGrille3D

Syntax :

# output on 3-D regular grid [xmin, xmax] x [ymin, ymax] x [zmin, zmax]
SismoGrille3D = AUTO nx ny nz
SismoGrille3D = xmin xmax ymin ymax zmin zmax nx ny nz

3-D solutions can be computed on the 3-D grid [xmin, xmax] x [ymin,ymax] x [zmin, zmax]. The number of points on each direction are nx, ny and nz. In our opinion, you have to be careful, because such output files are very expensive in memory. If you ask 200 points in each direction with complex double precision, a single file will take 122 Mo ... That's why it is often better to compute solution only on three planes by using the keyword SismoGrille. If AUTO is specified, the codes uses the bounding box of the physical domain (PML layers are excluded).

Example :

SismoGrille3D = AUTO 101 101 101

Location :

OutputFormat

Syntax :

OutputFormat = BINARY/ASCII FLOAT/DOUBLE

Actually only binary files are usable since in Matlab, there are no script files to read ascii files produced by Montjoie . Furthermore, seismogramms on points and radar cross sections are written in ASCII and directly readable with Matlab by using command load. The binary files can be written in simple precision (FLOAT) or double precision (DOUBLE), both formats are correctly interpreted in Matlab (by script loadND). You can also ask outputs in vtk format (readable by Paraview), BINARY has to be replaced by VTK and ASCII by ASCII_VTK. The default format is BINARY DOUBLE.

Example :

# outputs in "Matlab" format and medit for mesh outputs
OutputFormat = BINARY DOUBLE
OutputFormat = BINARY FLOAT

# For outputs in vtk format only :
OutputFormat = VTK DOUBLE

Location :

DirectoryOutput

Syntax :

DirectoryOutput = AUTO
DirectoryOutput = path

The user can specify the directory where output files will be written. If equal to AUTO, the output files will be written in [STIFFOUT]/num where num is a number that will be incremented for each next simulation. The default value is ./ (output files are written in the current directory).

Example :

# Directory where output files are written
DirectoryOutput = /home/user/Results/

Location :

ElectricOrMagnetic

Syntax :

ElectricOrMagnetic = component

This field is useful to select the component of the solution (and the gradient/curl) you want to know. If you wish to see all the components and gradients, you will put -1, whereas a positive integer will select a single component. If you want only components of the solution (and not gradients), put -2. The default component is 0 (so only the first component of solution). If you want to compute all the components of the solution, gradients and also hessians, you have to put -5. The hessians are computed through an interpolation on nodal points.

Example :

# all the components of the solution
ElectricOrMagnetic = -1

# all the components of the solution but not the gradients
ElectricOrMagnetic = -2

# only Ez for example
ElectricOrMagnetic = 2

# if you want to have all components, gradients and hessians
ElectricOrMagnetic = -5

Location :

FineMeshLobatto

Syntax :

FineMeshLobatto = YES

You can require to compute the values of the solution on the nodal values of the mesh. Instead of subdividing the mesh in regular subdivisions, you subdivide it on Gauss-Lobatto points, so that the vertices of the fine are coincident with the nodal points. The default value is NO.

Example :

FineMeshLobatto = YES

Location :

WriteSolutionQuadrature

Syntax :

WriteSolutionQuadrature = YES file_name_solution file_name_weights

You can require to compute the values of the solution on the quadrature points of the mesh. The second parameter is the file name where the solution will be written. The third parameter is the file name where the quadrature weights will be written. They can be read in Python by typing

  u = load1D('sol.dat');
  w = load1D_real('poids.dat');

Example :

WriteSolutionQuadrature = YES sol.dat poids.dat

Location :

MovePointsSurface

Syntax :

MovePointsSurface = YES component coef

You can require that outputs on surface meshes are produced with translated points. The second argument is the component number used to translate the point. Z-coordinate of the solution is modified:

where m is the component number. In 2-D computation, a 3-D mesh is generated with a z-coordinate.

Example :

MovePointsSurface = YES 1 0.2

Location :

FiniteElement/PointsReference.cxx

NonLinearSolver

Syntax :

NonLinearSolver = type_solver threshold nb_max_iteration

With this keyword, you can specify the non-linear solver used to invert transformation Fi. This transformation transforms any reference element into the real element. As a result, it is necessary to invert this transformation when points of output grids are localized in the mesh. The default solver is Newton solver, usually it works pretty well.

Example :

# if you want to specify a different threshold or different maximum number of iterations 
NonLinearSolver = NEWTON 1e-14 20

# other solvers are LVM (=Levenberg-Marquardt) and MINPACK
NonLinearSolver = LVM 1e-15 10
NonLinearSolver = MINPACK 1e-15 50

Location :

SismoMeshVolumetric

Syntax :

SismoMeshVolumetric = AUTO
SismoMeshVolumetric = nb_subdiv

You can require to compute the values of the solution on a subdivided mesh (with regular subdivisions). The number of subdivisions can be explicitely provided or equal to the order of approximation (AUTO). There is no default value. If FineMeshLobatto is enabled, the subdivisions are based on Gauss-Lobatto points instead of equidistributed points.

Example :

# Output on subdivided mesh with nb_subdiv = 3 for each edge
SismoMeshVolumetric = 3

Location :

SismoMeshSurfacic

Syntax :

SismoMeshSurfacic = AUTO BODY ref1 ref2 ... refn
SismoMeshSurfacic = AUTO CONDITION ref1 ref2 ... refn
SismoMeshSurfacic = AUTO REFERENCE ref1 ref2 ... refn
SismoMeshSurfacic = AUTO REFERENCE ALL

You can require to compute the values of the solution on the surfaces of the mesh (therefore in 3-D only). As for SismoMeshVolumetric, the solution is computed on a regular subdivision of the surfacic mesh. The number of subdivisions can be explicitely provided or equal to the order of approximation (AUTO). You can select all the surfaces of the mesh (REFERENCE ALL), or provide the references of the surface. There is no default value.

Example :

# surfacic mesh of all the surfaces
SismoMeshSurfacic = 3 REFERENCE ALL

# surfacic mesh of surface 1 and 5
SismoMeshSurfacic = AUTO REFERENCE 1 5

# you can ask output on bodies (defined with TypeBody)
SismoMeshSurfacic = AUTO BODY 2

# or on same boundary conditions (but you need to know numbers affected to each boundary condition)
SismoMeshSurfacic = AUTO CONDITION 1

Location :

SismoOutsidePoints

Syntax :

SismoOutsidePoints = file_points file_diffrac file_total dt
SismoOutsidePoints = CIRCLE x0 y0 z0 radius nb_points file_diffrac file_total dt
SismoOutsidePoints = CIRCLE_STEREO x0 y0 z0 radius nb_points dist file_diffrac file_total dt
SismoOutsidePoints = LINE x0 y0 z0 x1 y1 z1 nb_points file_diffrac file_total dt
SismoOutsidePoints = PLANE x0 y0 z0 x1 y1 z1 x2 y2 z2 Nx Ny file_diffrac file_total dt
SismoOutsidePoints = POINT x0 y0 z0 file_diffrac file_total dt

You can require to compute the values of the solution on points located outside the computational mesh (through an integral representation). This functionality is available only for wave equation and Helmholtz equation. You can provide a file containing a list of points or specify points computed on a circle, line or plane. CIRCLE_STEREO corresponds to pairs of points on a circle, the distance between the two points of a pair is given by the user. The last parameter is the time step dt (values are computed each dt), this parameter is not needed for time-harmonic simulations (e.g. Helmholtz equation). The output files are readable in Python/Matlab with the function ReadFarField in time-domain simulation and with loadND for time-harmonic simulation. Of course, in 2-D simulations, the z-coordinate should be omitted in the list of parameters. The surface where the integrals are computed (throught integral representation) is described in ParametersRCS. It is assumed that the media is homogeneous outside this surface (and no source is present). The values of the media at infinity are initialized if you have provided a line with ReferenceInfinity. When a scattered field is computed (diffracted field or total field), only the diffracted field will be computed.

Example :

# To set rho0 and mu0
# you give the reference of the homogeneous media (infinity)
ReferenceInfinity = 1

# the solution is computed through an integral over surface specified with TypeBody/ParametersRCS
TypeBody = 1 1
ParametersRCS = YES 1 AUTO

# you can write all the desired points in a file Points3D.txt
# SismoOutsidePoints = file_points file_diffrac file_total dt
SismoOutsidePoints = Points3D.txt diffracExt.dat totalExt.dat 0.02

# output on points located regularly on a circle
# SismoOutsidePoints = CIRCLE x0 y0 z0 radius nb_points file_diffrac file_total dt
SismoOutsidePoints = CIRCLE 0.0 0.0 0.0   4.0 101 diffracExt.dat totalExt.dat 0.02

# output on pairs of points located regularly on a circle
# SismoOutsidePoints = CIRCLE_STEREO x0 y0 z0 radius nb_points dist file_diffrac file_total dt
SismoOutsidePoints = CIRCLE_STEREO 0.0 0.0 0.0   4.0 101 0.02 diffracExt.dat totalExt.dat 0.02

# output on a single point
# SismoOutsidePoints = POINT x0 y0 z0 file_diffrac file_total dt
SismoOutsidePoints = POINT 0.0 0.0 0.0 diffracExt.dat totalExt.dat 0.02

# output on a 2-D line for a 2-D simulation (you omit the z-coordinates)
# time-harmonic simulation : dt is not needed
# you specify the two extremities of the line (x0, y0) and (x1, y1)
# SismoOutsidePoints = LINE x0 y0 x1 y1 N diffracFile totalFile
SismoOutsidePoints = LINE -3.0 -3.0   3.0 -3.0  501 diffracExt.dat totalExt.dat

# output on a 2-D plane (2-D simulation)
# (x0, y0) is the base of the parallelogramm
# (x1, y1) and (x2, y2) the two points adjacent to (x0, y0)
# Nx and Ny the number of points on each side of the parallelogramm
# SismoOutsidePoints = PLANE x0 y0  x1 y1  x2 y2  Nx Ny diffracFile totalFile
SismoOutsidePoints = PLANE -3.0 -3.0   3.0 -3.0   -3.0 3.0  201 201 diffracExt.dat totalExt.dat

Location :

FileOutputGrille, FileOutputGrille3D, FileOutputPlane, FileOutputLine, etc

Syntax :

FileOutputPlane = nom_fichier_diffrac nom_fichier_total

In this field, you specify the files where the diffracted field and total will be written. All the outputs on planes will be stored in the same file, all the outputs on lines will be in another same file, etc. The exceptions are FileOutputMeshVolumetric, FileOutputMeshSurfacic, which are required AFTER each SismoOutputMeshVolumetric (or SismoOutputMeshSurfacic). In this case, different files are used for the creation of volumetric or surfacic meshes. There is no default value, you need to specify file names.

Example :

# For Matlab files, extension .dat
FileOutputPlane = diffracTest.dat totalTest.dat

# For Medit files, extension .bb
FileOutputMeshVolumetric = diffracTest.bb totalTest.bb

Location :

FileOutputPoint

Syntax :

FileOutputPoint = nom_fichier_diffrac nom_fichier_total

In this field, you specify the files where the diffracted field and total will be written. In the case where there is no diffracted/total field distincition (a simple gaussian source for example), the output is always written in the total field, that is to say the second parameter. For time-domain computations, seismogramms are stored in this file. So the file name is not evolving (0001.dat, 0002.dat) and all the values are stored in the same file, contrary to all other outputs (FileOutputPlane, FileOutputGrille, etc). This seismogramm is usually in ascii format, and can be read by using load function in Matlab (loadtxt in Python). It will look like

t0  u(x0, t0) u(x1, t0) u(x2, t0) ... u(xn, t0)
t1  u(x0, t1) u(x1, t1) u(x2, t1) ... u(xn, t1)
...
tn  u(x0, tn) u(x1, tn) u(x2, tn) ... u(xn, tn)

If ElectricOrMagnetic has been set to a positive value, u(x0, t0) will be equal to the asked component of the solution on the point x0, u(x1, t0) on the point x1, etc. The points x0, x1, ..., xn are given by inserting as many lines "SismoPoint = xi yi" as you want. The order will be the same as the order specified in the data file. If ElectricOrMagnetic has been set to -1, the columns will be the following :

t0  u(x0, t0) u(x1, t0) u(x2, t0) ... u(xn, t0) du(x0, t0) du(x1, t0) du(x2, t0) ... du(xn, t0)
t1  u(x0, t1) u(x1, t1) u(x2, t1) ... u(xn, t1) du(x0, t1) du(x1, t1) du(x2, t1) ... du(xn, t1)
...
tn  u(x0, tn) u(x1, tn) u(x2, tn) ... u(xn, tn) du(x0, tn) du(x1, tn) du(x2, tn) ... du(xn, tn)

where du(x0, t0) will be the gradient of u at point x0. u may have several components and each component is written in a separate column, and du have more components, each component being written in a separate column. If ElectricOrMagnetic has been set to -2, only u is written with all its components.

Example :

# for time-domain computations, first name is not used
FileOutputPoint = diffrac.dat Sismo.dat

Location :

ParametersOutputGrille, ParametersOutputGrille3D, ParametersOutputPlane, ParametersOutputLine, etc

Syntax :

ParametersOutputPlane = t0 tf dt

This field is devoted to transient equations for which snapshots are written at time t in [t0, tf] at every dt. The exceptions are ParametersOutputMeshVolumetric, ParametersOutputMeshSurfacic, which are required AFTER each SismoOutputMeshVolumetric (or SismoOutputMeshSurfacic). In this case, different time intervals can be specified for each output. The default values are set to 0.

Example :

# output for t in [0, 10] for each dt = 0.1
ParametersOutputPlane  = 0.0  10.0  0.1

# for seismograms, you ask a binary output :
ParametersOutputPoint = 0.0 10.0 0.01 BINARY

Location :

DataFileExperiment

Syntax :

DataFileExperiment = input_file ref0 ref1 ... ref_n

This field is devoted for RTM and inverse problems. For these problems, experimental data is needed. This data is generated synthetically by running the case described in the file input_file. The references ref0, ref1 ... ref_n are the references where the solution is measured (on quadrature points).

Example :

# synthetical data is generated by "running" the input file INI/experiment.ini
# and measuring the solution on references 1, 2 and 3
DataFileExperiment = INI/experiment.ini 1 2 3

Location :

DataFileSimulation

Syntax :

DataFileSimulation = input_file ref0 ref1 ... ref_n

This field is devoted for RTM and inverse problems. For these problems, the simulated data and back-propagated data are generated by starting from the case detailed in the file input_file. The references ref0, ref1 ... ref_n are the references where the solution is measured and back-propagated (on quadrature points).

Example :

# simulated and back-propagated data are generated with INI/simu.ini
DataFileSimulation = INI/simu.ini 1 2 3

Location :

TypeSourceLine

Syntax :

TypeSourceLine = nb_source REF ref Gaussian xA yA  xB yB radius rmax
TypeSourceLine = nb_source REF ref PlaneWave tetaA tetaB
TypeSourceLine = nb_source IncidentWave tetaA tetaB

This field is devoted for RTM and inverse problems. It defines a sequence of sources (usually along a line). nb_source is the number of sources, ref is the reference (for surface sources) associated with the surface. Gaussian specifies surface Gaussian sources, and the center of gaussians describes the line between two points. PlaneWave speficies a surface plane wave source, the incident angle will describe the given interval [tetaA, tetaB]. Finally IncidentWave is used to specify a total field source with different incident angles.

Example :

# example with 2-D simulation and gaussian
# 20 sources between points (-2, 2) and (4, 2)
# radius of the Gaussian : 0.1
# 0.2 is a cut-off radius (the gaussian is truncated for r ≥ 0.2)
TypeSourceLine = 20 REF 1 Gaussian -2.0 2.0  4.0 2.0   0.1 0.2

# instead of a Gaussian, you can put a plane wave exp(i k x)
# where k = omega/c (cos teta, sin teta)
# 20 sources with teta between 30 and 120 degrees :
TypeSourceLine = 20 REF 1 PlaneWave 30.0 120.0

# incident field (the total field will be computed)
# the incident wave will be a plane wave
TypeSourceLine = 20 IncidentWave 30.0 120.0

Location :

TypeSourceCircle

Syntax :

TypeSourceCircle = nb_source REF ref Gaussian xA yA  radiusA tetaA tetaB radius rmax

This field is devoted for RTM and inverse problems. It defines a sequence of sources (usually along a circle). nb_source is the number of sources, ref is the reference (for surface sources) associated with the surface. Gaussian specifies surface Gaussian sources, the center of the gaussians describes a circle of center (xA, yA) and radius radiusA, with the two angles tetaA and tetaB. radius and rmax are the radii of the Gaussian.

Example :

# example with 2-D simulation and gaussian
# 20 sources in a circle of center (0,0) and radius 2
# half-circle => angles between 0 and 180
# radius of the Gaussian : 0.1
# 0.2 is a cut-off radius (the gaussian is truncated for r ≥ 0.2)
TypeSourceCircle = 20 REF 1 Gaussian 0.0 0.0 2.0  0.0  180.0   0.1 0.2

Location :