This function is used to perform a h/p/hp adaptation according to the defined h/p/hp adaptivity settings. To define the
h-adaptivity use function mesh.setadaptivity(). To define a p-adaptivity for a field use function field.setorder().
The function returns True if the mesh or any field order was changed by the adaption and returns False if no changes were
made.
This is a collective MPI operation and hence must be called by all ranks. It replaces the adapt() function in the DDM
framework. This function is used to perform a h/p/hp adaptation according to the defined h/p/hp adaptivity settings. To define the
h-adaptivity use function mesh.setadaptivity(). To define a p-adaptivity for a field use function field.setorder().
This function extracts the S-parameters Sij corresponding to the physical regions portphysregs from the solution list sols.
The parameters sources and modeprojects should contain what was used as drivesignal and lumpfield with modedrive function.
The function returns two lists that contain the real and imaginary parts of the S-parameters in row-major order (S11,S12,…), respectively.
For the definition of S-parameters, we assume that the input signal at port i is a complex number ai times the reference mode of port i.
If port i is driven, ai is the driving signal used for feeding; otherwise ai=0.
The output signal at port i is a complex number bi times the reference mode of port i.
The S-parameter Sij is defined such that the equation Sij=bi/aj holds when port j is driven and the others are not.
This function computes the eigenmodes of a waveguide whose cross section S is the given physical region portphysreg.
The material parameters μ, ϵ, and σ are required as input arguments; these can be anisotropic, and both real and imaginary parts should be given.
The surface normal portnormal is assumed to point in the direction of propagation.
The boundary condition (PEC) and the approximation order is extracted from the input field E.
For each found mode, the function returns the real and imaginary parts of the transverse and longitudal components of the fields E and H and the propagation constant γ=α+jβ as a list of expressions.
The returned modes are scaled such that the power into the waveguide,
P=21∫SRe(E×H∗)⋅dS,
equals one watt; integrationorder determines the order of the integration rule used to compute this integral.
By default, the function throws an error if material parameters are anisotropic in the direction of propagation.
This is because the computed modes can not be used with modedrive function in that case.
To compute the modes anyway, the error check can be disabled by setting errorifzanisotropic = false.
This integrates the expression expr over the physical region physreg. The integration is exact up to the order of
polynomials specified in the argument integrationorder.
This interpolates the expression expr value at a single point whose [x,y,z] coordinate is provided as an argument.
The flattened interpolated field values are returned if the point was found in the elements of the
physical region physreg.
This returns the maximum value of the expression expr obtained over the geometric region physreg by splitting all
elements refinement times in each direction. Increasing the refinement will thus lead to a more accurate maximum value,
but at an increased computational cost. The maximum value is exact when the refinement nodes added to the elements correspond
to the position of maximum. For a first-order nodal shape function interpolation, on a mesh that is not curved, the maximum
is always exact to machine precision.
This returns a relative L2 norm according to the below formula:
∣c∣∣a−b∣
Example
allmin
This returns the minimum value of the expression expr obtained over the geometric region physreg by splitting all
elements refinement times in each direction. Increasing the refinement will thus lead to a more accurate minimum value,
but at an increased computational cost. The minimum value is exact when the refinement nodes added to the elements correspond
to the position of minimum. For a first-order nodal shape function interpolation, on a mesh that is not curved, the minimum
is always exact to machine precision.
This is a collective MPI operation and hence must be called by all ranks. This function partitions the requested mesh into a
number of parts equal to the number of ranks. All parts are saved to disk and the part file name for each rank is returned.
The GMSH API must be available to partition into more than one part. A physical region containing the global skin can be
created with the last argument.
This functions returns the value of a scalar expression expr at the point region physreg.
Example
allsolve
This is a collective MPI operation and hence must be called by all ranks. It solves across all the ranks a nonlinear problem
with a fixed-point iteration.
This returns an expression whose value is 1 for all evaluation points where the value of all the input expressions is larger
or equal to zero. Otherwise, its value is -1.
This returns the selected component of a column vector expression. For a column vector expression, selectedcomp is 0 for
the first, 1 for the second component and 2 for the third component respectively. For a matrix expression, the whole corresponding
row is returned in the form of an expression. Thus, if selectedcomp, for example is 5, then an expression containing the entries of the
fifth row of the matrix is returned.
Complexdivision computes the quotient a / b of two complex-valued scalar expressions a and b, which are presented as lists of length two (containing the real and imaginary parts).
Example
complexinverse
Complexinverse computes the inverse 1 / a of a complex-valued scalar expression a, which is presented as a list of length two (containing the real and imaginary parts).
Example
complexproduct
Complexproduct computes the product a * b of two complex-valued expressions a and b, which are presented as lists of length two (containing the real and imaginary parts).
Example
computeradiationpattern
compx
This returns the first or x component of a column vector expression. For a matrix expression, an expression containing the entries
of the first row of the matrix is returned. This is equivalent to setting selectedcomp=0 in comp().
This returns the second or y component of a column vector expression. For a matrix expression, an expression containing the entries
of the second row of the matrix is returned. This is equivalent to setting selectedcomp=1 in comp().
This returns the third or z component of a column vector expression. For a matrix expression, an expression containing the entries
of the third row of the matrix is returned. This is equivalent to setting selectedcomp=2 in comp().
This returns the formulation terms required to enforce field continuity.
Examples
Example 1: continuitycondition(gamma1:int, gamma2:int. u1:field, u2:field, errorifnotfound:bool=True, lagmultorder:int=0)
This returns the formulation terms required to enforce u1=u2 between boundary region Γ1 and Γ2 (with
Γ1⊆Γ2, meshes can be non-matching). In case Γ2 is larger than Γ1 (Γ1⊂Γ2) the boolean flag errorifnotfound must be set to false.
Example 2:continuitycondition(gamma:int, gamma2:int, u1:field, u2:field, rotcent:List[double], rotangz:double, angzmod:double, factor:double, lagmultorder:int=0)
This returns the formulation terms required to enforce field continuity across an angzmod degrees slice of a rotor-stator
interface where the rotor geometry is rotated by rotangz around the z axis with rotation center at rotcent. This
situation arises for example in electric motor simulations when (anti)periodicity can be considered and thus only a slice
of the entire 360 degrees needs to be simulated. Use a factor of -1 for antiperiodicity. Boundary Γ1 is the
rotor-stator interface on the (non-moving) stator side while the boundary Γ2 is the interface on the rotor side. In
the unrotated position the bottom boundary of the stator and rotor slice must be aligned with the x axis.
The condition is based on a Lagrange multiplier of the same type and the same harmonic content as the field u1 and u2. The
mortar finite element method is used to link the unknown dof field on Γ1 and Γ2 so that there is no
restriction on the mesh used for both regions.
This declares an unknown field (dof\ for degree of freedom). The dofs are defined only on the region physreg which when
not provided is set to the element integration region.
This computes the double-dot product of two matrix expressions. The returned expression is a scalar.
A:B=i,j∑AijBij
Example
dt
This returns the first-order time derivative expression.
Examples
Example 1: dt(input:expression)
Example 2: dt(input:expression, initdt:double, initdtdt)
This gives the transient approximation of the first-order time derivative of a space-independent expression. The
initial values must be provided when using generalized alpha (genalpha) and are ignored otherwise.
This returns the second-order time derivative expression.
Examples
Example 1: dtdt(input:expression)
Example 2: dtdt(input:expression, initdt:double, initdtdt:double)
This gives the transient approximation of the second-order time derivative of a space-independent expression. The
initial values must be provided when using generalized alpha (genalpha) and are ignored otherwise.
This returns an expression whose value is the interpolation order on each element for the provided input field. The value is a
constant on each element. When the argument alpha is set, the value returned is the lowest order required to include
alpha percentage of the total shape function coefficient weight. An additional optional argument absthres can be set to provide a
minimum total weight below which the lowest possible field order is returned.
Example
getc0
getcirculationports
getcirculationsources
getdimension
This returns the x, y and z mesh dimensions.
Example
getepsilon0
This returns the value of vacuum permittivity ϵ0=8.854187812813e−12Fm−1.
Example
getextrusiondata
This gives the relative depth δ in the extruded layer, the extrusion normal n and tangents t1
and t2. This is useful when creating Perfectly Matched Layers (PMLs).
This returns a single harmonic from a multi-harmonic expression. Set a positive last argument to use an FFT to compute the
harmonic. The returned expression is on harmonic 1.
Example
getmu0
This returns the value of vacuum permeability μ0=1.2566370621219e−06NA−2.
Example
getpi
This returns value of π.
Example
getrandom
This returns a random value uniformly distributed between 0.0 and 1.0.
This returns the components of the total magnetostatic/electrostatic force acting on a given region. In the axisymmetric
case zero x and z components are returned and the y component includes a 2π factor to provide the force acting
on the corresponding 3D shape. Units are ”N per unit depth” in 2D and N in 3D and 2D axisymmetry.
Examples
Example 1: gettotalforce(physreg:int, EorH:expression, epsilonormu:expression, extraintegrationorder:int=0)
Example 2: gettotalforce(physreg:int, meshdeform:expression, EorH:expression, epsilonormu:expression, extraintegrationorder:int=0)
This is similar to the above function but the total force is computed on the mesh deformed by the field u.
For a scalar input expression, this is mathematically treated as the gradient of a scalar (∇v) and the output is a
column vector with one entry per space derivative. For a vector input expression, this is mathematically treated as the gradient
of a vector (∇u) and the output has one row per component of the input and one column per space derivative.
This defines the (nonlinear) Green-Lagrange strains in Voigt form (ϵxx,ϵyy,ϵzz,γyz,γxz,γxy). The input can either be the displacement field or its gradient.
This writes a .pvd ParaView file to group a set of .vtu files that are time solutions at the time values provided in
timevals.
Examples
Example 1: grouptimesteps(filename::str, filestogroup:List[str], timevals:List[double])
Example 2: grouptimesteps(filename::str, fielprefix:str, firstint:int, timevals:List[double])
This is similar to the previous function except that the full list of file names to group does not have to be provided. The
file names are constructed from the file prefix with an appended integer starting from ‘firstint’ by steps of 1. The filenames
are ended with .vtu.
harm
ifpositive
This returns a conditional expression. The argument condexpr specifies the conditional argument. The expression value is
trueexpr for all evaluation points where condexpr is larger or equal to zero. Otherwise, its value is falseexpr.
This function loads a mesh file to shapes. The output holds a shape for every physical region of dimension d (0D, 1D, 2D, 3D)
defined in the mesh file. The loaded shapes can be edited (extruded, deformed, …) and grouped with other shapes to create a
new mesh. Note that the usage of loaded shapes might be more limited than other shapes.
Example
loadvector
This loads a list from the file filename. The delimiter specfies the character separating each entry in the file.
The sizeincluded must be set to True if the first number in the file is the size of the list.
This returns a multi-harmonic expression whose harmonic numbers and expressions are provided as arguments. The argument
expressions must be on harmonic 1.
Example
max
This returns an expression whose value is the maximum of the two input arguments a and b.
This returns an expression whose value is the length/area/volume for each 1D/2D/3D mesh element respectively. The value is
constant on each mesh element. The integrationorder determines the accuracy of mesh size calculated. Higher the number,
the better the accuracy. Integration order cannot be negative and if the integration order < 0 a RuntimeError is raised.
Example
min
This returns an expression whose value is the minimum of the two input arguments a and b.
This is a modulo function. This returns an expression equal to the remainder resulting from the division
of input by modval.
Example
modedrive
Modedrive defines a boundary condition in the physical region portphysreg by feeding and absorbing the mode found by alleigenport (parameters Etreal, Etimag, Htreal, and Htimag).
The input field E is constrained to be a complex multiple of the mode; lumpfield defines this multiplier and must have been prepared into the formulation beforehand.
If modedrive is applied with different modes in the same region, E is constrained to be their linear combination.
The mode is fed according to drivesignal, which should be expressed as sin(ωt+shift) or cos(ωt+shift), where ω is the driving frequency and shift is an optional phase shift.
The parameter blocktag is the block number of the feeding term (for the absorbing term this is always 0).
The feeding term only affects the right hand side, so this information can be used to solve the same formulation with multiple right hand sides, e.g. when computing S-parameters.
This returns an expression equal to the input expression with a selected and moved harmonic content. Set a positive last
argument to use an FFT to compute the harmonics of the input expression.
Example
norm
This gives the L2 norm of an expression input.
Example
normal
This defines a normal vector with unit norm. If a physical region is provided as an argument then the normal points out of it.
if no physical region is provided then the normal can be flipped depending on the element orientation in the mesh.
Examples
on
This function allows to use fields, unknown dof\ fields or general expressions across physical regions with possibly non-
matching meshes by evaluating the expression argument using a (x, y, z) coordinate interpolation. It makes it straightforward
to setup the mortar finite element method to enforce general relations, such as field equality u1=u2, at
the interface Γ of non-matching meshes. This can for example be achieved with a Lagrange multiplier λ such that
∫Γ(λ(u1′−u2′)+λ′(u1−u2))dΓ=0
holds for any appropriate field u1′, u2′ and λ′. The example below illustrates the
formulation terms needed to implement the Lagrange multiplier between interfaces Γ1 and Γ2.
Examples
Example 1: on(physreg:int, expr:expression, errorifnotfound:bool=True)
physreg is the physical region across which the expression expr is evaluated. The expr argument can be fields or dof\
fields or any general expressions.
When setting the flag errorifnotfound=False, any point in Γ1 without a relative in Γ2 (and vice versa)
does not contribute to the assembled matrix. The default value is True for which an error is raised if a point in Γ1 is without
a relative in Γ2 (and vice versa).
The case where there is no unknown dof\ term in the expression argument is described below:
With the on function, the (x, y, z) coordinates corresponding to each Gauss point of the integral are first calculated then the
dz(z)\ expression is evaluated through interpolation at these (x, y, z) coordinates on region vol. With on here dz(z)\ is correctly evaluated as 1 because the z-derivative calculation is performed on the volume region vol. Without the
on operator the z-derivative would be wrongly calculated on the top face of the disk (a plane perpendicular to the z-axis).
If a requested interpolation point cannot be found (because it is outside of physreg or because the interpolation algorithm
fails to converge, as can happen on curved 3D elements) then an error occurs unless errorifnotfound is set to False. In the
latter case, the value returned at any non-found coordinate is zero, without raising an error.
Example 2: on(physreg:int, expr:coordshift, expr:expression, errorifnotfound:bool=True)
This is similar to the previous example but here the (x, y, z) coordinates at which to interpolate the expression are shifted
by (x+compx(coordshift), y+compy(coordshift), z+compz(coordshift)).
orpositive
This returns an expression whose value is 1 for all evaluation points where at least one input expression has a value
larger or equal to zero. Otherwise, its value is -1.
This returns the formulation terms required to enforce on field u a rotation or translation periodic condition between
boundary region Γ1 and Γ2 (meshes can be non-conforming). A factor different than 1 can be
provided to scale the field on Γ2 (use -1 for antiperiodicity).
Example
The condition is based on a Lagrange multiplier of the same type and the same harmonic content as the field u1 and u2. The
mortar finite element method is used to link the unknown dof field on Γ1 and Γ2 so that there is no
restriction on the mesh used for both regions.
More advanced periodic conditions can be implemented easily using on() function.
Zero artificial wave reflection at the truncation boundary happens only if it is perpendicular to the outgoing waves. In
practical applications however the truncation boundary is not at an infinite distance from the acoustic source and the wave
amplitude is not constant and thus, some level of artificial reflection cannot be avoided. To minimize this effect the
truncation boundary should be placed as far as possible from the acoustic source (at least a few wavelengths away).
An acoustic attenuation value can be provided (in Neper/m) in case of harmonic problems. For convenience use the function
dbtoneper() to convert dB/m attenuation values to Np/m.
Example
predefinedacousticstructureinteraction
This function defines the bi-directional coupling for acoustic-structure interaction at the medium interface. Field p is
the acoustic pressure and field u is the mechanical displacement. Calling n the normal to
the interface pointing out of the solid region, the bi-directional coupling is obtained by adding the fluid pressure loading
to the structure
fpressure=−pn
as well as linking the structure acceleration to the fluid pressure normal derivative using Newton’s law:
∂np=−ρfluid∂t2∂2u⋅n
To have a good matrix conditioning a scaling factor s (e.g s=1e10) can be provided. In this case, the pressure source
is divided by s and, to compensate, the pressure force is multiplied by s. This leads to the correct membrane deflection
but the pressure field is divided by the scaling factor.
An acoustic attenuation value can be provided (in Neper/m) in case of harmonic problems. For convenience use the function
dbtoneper() to convert dB/m attenuation values to Np/m.
Example
predefinedacousticwave
This function defines the equation for (linear) acoustic wave propagation:
∇2p−c21∂t2∂2p=0
An acoustic attenuation value can be provided (in Neper/m) in case of harmonic problems. For convenience use the function
dbtoneper() to convert attenuation values from dB/m to Np/m.
The arguments have the following meaning:
dofpis the dof of the acoustic pressure field.
tfp is the test function of the acoustic pressure field.
soundspeed is the speed of sound in m/s.
neperattenuation is the attenuation in Np/m.
pmlterms is the list of pml terms.
precondtype is the type of precondition.
Examples
Example 1: predefinedacousticwave(dofp:expression, tfp:expression, soundspeed:expression, neperattenuation:expression, precondtype:str="")
In the illustrative example below, a highly-attenuated acoustic wave propogation in a rectangular 2D box is simulated.
Example 2: predefinedacousticwave(dofp:expression, tfp:expression, soundspeed:expression, neperattenuation:expression, pmlterms:List[expression], precondtype:str="")
This is the same as the previous example but with PML boundary conditions.
predefinedadmittancebcincompressible
This function returns the weak formulation of admittance boundary condition for incompressible flows:
where p is the pressure, V the velocity of the flow, n the normal vector to the boundary (pointing outward),
and Y(ω) the frequency-dependent admittance in harmonic form. It is assumed that the tabulated values of Y(ω) are provided
as function of frequency by the user.
This defines the weak formulation for the generalized advection-diffusion equation:
β∂t∂c−∇⋅(α∇c)+γ∇⋅(cv)=0
where c is the scalar quantity of interest and v is the velocity that the quantity is moving with. With
β and γ set to unit, the classical advection-diffusion equation with diffusivity
tensor α is obtained. Set isdivvzero to True if ∇⋅v is zero (for
incompressible flows).
This is a collective MPI operation and hence must be called by all the ranks. This function returns
[ℜ(∣D∣)ℑ(∣D∣)ℜ(D)ℑ(D)ℜ(D−1)ℑ(D−1)]
where D is the PML transformation matrix for a square box in a square box. A hyperbolic or shifted hyperbolic
PML can be selected with the last argument.
Example
predefineddiffusion
This defines the weak formulation for the generalized diffusion equation:
β∂t∂c−∇⋅(α∇c)=0
where c is the scalar quantity of interest. With β set to unit, the classical diffusion equation with diffusivity
tensor α is obtained.
This defines a classical linear elasticity formulation.
Examples
Example 1: predefinedelasticity(dofu:expression, tfu:expression, Eyoung:expression, nupoisson:Expression, myoption:str="")
This defines a classical linear isotropic elasticity formulation whose strong form is:
⎩⎨⎧∇⋅σ−ρu¨+F=0σ=H:ϵϵ=21[∇u+(∇u)T]
where,
u is the displacement field vector
σ is the Cauchy stress tensor in Voigt notation
ϵ is the strain tensor in Voigt notation
H is the 4th order elasticity/stiffness tensor
F is the volumetric body force vector
ρ is the mass density
This is used when the material is isotropic. u is the mechanical displacement, Eyoung is the Young’s modulus [Pa] and
nupoisson is the Poisson’s ratio. In 2D the option string must be either set to “planestrain” or “planestress” for a
plane strain or plane stress assumption respectively.
Example 2: predefinedelasticity(dofu:expression, tfu:expression, elasticitymatrix:expression, myoption:str="")
This extends the previous definition (Example 1) to general anisotropic materials. The elasticity matrix [Pa] must be provided such
that it relates the stress and strain in Voigt notation.
Example 3: predefinedelasticity(dofu:expression, tfu:expression, u:field, Eyoung:expression, nupoisson:expression, prestress:expression, myoption:str="")
This defines an isotropic linear elasticity formulation with geometric nonlinearity taken into account (full-Lagrangian
formulation using the Green-Lagrange strain tensor). Problems with large displacements and rotations can be simulated with
this equation but strains must always remain small. Buckling, snap-through and the likes or eigenvalues of prestressed
structures can be simulated with the above equation in combination with a nonlinear iteration loop.
The prestress vector [Pa] must be provided in Voigt notation (σxx,σyy,σzz,σyz,σxz,σxy).
Set the prestress expression to 0.0 for no prestress.
Example 4: predefinedelasticity(dofu:expression, tfu:expression, u:field, elasticitymatrix:expression, prestress:expression, myoption:str="")
This extends the previous definition (Example 3) to general anisotropic materials. The elasticity matrix [Pa] must be provided such
that it relates the stress and strain in Voigt notation. Similarly, the prestress vector [Pa] must be also provided in Voigt
notation (σxx,σyy,σzz,σyz,σxz,σxy).
predefinedelasticradiation
predefinedelectrostaticforce
This function defines the weak formulation term for electrostatic forces. The first argument is the mechanical displacement test
function or its gradient, the second is the electric field expression and the third argument is the electric permittivity (
must be a scalar).
Let us call T [N/m2] the electrostatic Maxwell stress tensor:
T=ϵE⊗E−21ϵ(E⋅E)I
where ϵ is the electric permittivity, E is the electric field and I is the identity matrix.
The electrostatic force density is ∇⋅T[N/m^3] so that the loading for a mechanical problem can
be obtained by adding the following term:
∫Ω(∇⋅T)⋅u′dΩ
where u is the mechanical displacement. The term can be rewritten in the form that is provided by this function:
−∫ΩTϵ′dΩ
where ϵ is the infinitesimal strain tensor. This is identical to what is obtained using the virtual
work principle. For details refer to ‘Domain decomposition techniques for the nonlinear, steady state, finite element
simulation of MEMS ultrasonic transducer arrays’, page 40.
In this function, a region should be provided to the test function argument to compute the force only for the degrees of
freedom associated to that specific region (in the example below with tf(u, top) the force only acts on the surface region ‘top’.
In any case, a correct force calculation requires including in the integration domain all elements in the region where
the force acts and in the element layer around it (in the example below ‘vol’ includes all volume elements touching
surface ‘top’).
This defines the equation for (linear) electromagnetic wave propagation:
∇×(μ1∇×E)+σ∂t∂E+ϵ∂t2∂2E=0
where E is the electric field, μ is the magnetic permeability, ϵ is the electric permiitivity
and σ is the electric conductivity. The real and imaginary parts of each material property can be provided.
The argument have the following meaning:
dofE is the dof of the electric field.
tfE is the test function of the electric field.
mur and mui is the real and imaginary part of the magnetic permeability μ.
epsr and epsi is the real and imaginary part of the electric permittivity ϵ.
sigr and sigi is the real and imaginary part of the electric conductivity σ.
where v is the velocity, p the pressure, u the displacement of the structure, Ω the fluid domain, and Γ the whole boundary of the fluid domain. The dashed variables are the test functions in the weak formulation.
The no-slip boundary condition at the fluid-structure interface is applied by using a Lagrange multiplier (λ) method as
∫Γfsi((v−u˙)⋅λ′+λ⋅v′)dΓfsi=0
The Lagrange multiplier λ naturally contains the viscous forces acting on the fluid, which can be applied directly to the structure with an inverted sign along with the normal pressure force as
∫Γfsi(pn−λ)⋅u′dΓfsi=0
where Γfsi the boundary interface between fluid and structure domains.
This function returns the weak formulation of impedance boundary condition for incompressible flows:
p(ω)=Z(ω)(V⋅n)
where p is the pressure, V the velocity of the flow, n the normal vector to the boundary (pointing outward),
and Z(ω) the frequency-dependent impedance in harmonic form. It is assumed that the tabulated values of Z(ω) are provided
as function of frequency by the user.