Skip to content

Script API

This is the API definition of quanscient python module.

Functions

abs

def abs(input: expression) -> expression

This returns an expression that is the absolute value of the input expression.

Example
>>> expr = abs(-2.5)
>>> expr.print();
Expression size is 1x1
@ row 0, col 0 :
2.5

acos

def acos(input: expression) -> expression

This returns an expression that is the arccosarccos or cos1cos^{-1} of input. The output expression is in radians.

Example
>>> expr = acos(0)
>>> expr.print(); # in radians
Expression size is 1x1
@ row 0, col 0 :
1.5708
>>>
>>> (expr * 180/getpi()).print(); # in degrees
Expression size is 1x1
@ row 0, col 0 :
90
See Also

sin(), cos(), tan(), asin(), atan()

adapt

def adapt(verbosity: int = 0) -> bool

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.

Example
>>> all = 1
>>> q = shape("quadrangle", all, [0,0,0, 1,0,0, 1.2,1,0, 0,1,0], [5,5,5,5])
>>> mymesh = mesh([q])
>>> x = field("x"); y = field("y")
>>> criterion = 1 + sin(10*x)*sin(10*y)
>>>
>>> mymesh.setadaptivity(criterion, 0, 5)
>>>
>>> for i in range(5):
... criterion.write(all, f"criterion_{100+i}.vtk", 1)
... adapt(1)
See Also

mesh.setadaptivity(), field.setorder(), alladapt()

alladapt

def alladapt(verbosity: int = 0) -> bool

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().

Example
>>> ...
>>> alladapt()
See Also

adapt()

allcomputesparameters

def allcomputesparameters(
portsphysregs: list[int],
E: field,
Esources: list[expression],
Esols: list[vec],
integrationorder: int = 5
) -> list[list[float]]
def allcomputesparameters(
portsphysregs: list[int],
sources: list[expression],
sols: list[vec],
modeprojects: list[field]
) -> list[list[float]]
def allcomputesparameters(
portsphysregs: list[int],
sources: list[expression],
sols: list[vec],
modeprojectsorprimals: list[expression],
Rrefs: list[float]
) -> list[list[float]]

This function extracts the S-parameters SijS_{ij} 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,S_{11}, S_{12}, \ldots), respectively.

For the definition of S-parameters, we assume that the input signal at port ii is a complex number aia_i times the reference mode of port ii. If port ii is driven, aia_i is the driving signal used for feeding; otherwise ai=0a_i=0. The output signal at port ii is a complex number bib_i times the reference mode of port ii. The S-parameter SijS_{ij} is defined such that the equation Sij=bi/ajS_{ij}=b_i/a_j holds when port jj is driven and the others are not.

Example
>>> # Suppose a mesh has been loaded with physical regions (port1, port2, waveguide, pec_boundary) and E, mur, epsr_real, and eps_imag have been defined so that the following calls make sense:
>>> port1mode0 = qs.alleigenport(port1, 1, 0, 1000, E, mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, [0.0, 0.0, 1.0], 1e-06, 1000)[0]
>>> port2mode0 = qs.alleigenport(port2, 1, 0, 1000, E, mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, [0.0, 0.0, -1.0], 1e-06, 1000)[0]
>>>
>>> form = qs.formulation()
>>>
>>> # Prepare the lumpfield constraints:
>>> lumpfields = form.lump([port1, port2], [2, 3])
>>>
>>> # Electromagnetic waves
>>> form += qs.integral(waveguide, 3, qs.predefinedemwave(qs.dof(E), qs.tf(E), mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, "oo2"))
>>>
>>> # Drive and absorb the mode of port1 (feeding term in block number 1):
>>> form += qs.modedrive(port1, qs.sn(1), E, port1mode0[0], port1mode0[1], port1mode0[4], port1mode0[5], [0.0, 0.0, 1.0], lumpfields[0], 1)
>>>
>>> # Drive and absorb the mode of port2 (feeding term in block number 2):
>>> form += qs.modedrive(port2, qs.sn(1), E, port2mode0[0], port2mode0[1], port2mode0[4], port2mode0[5], [0.0, 0.0, -1.0], lumpfields[1], 2)
>>>
>>> sols = form.allsolve(relrestol=1e-06, maxnumit=1000, nltol=1e-05, maxnumnlit=-1, relaxvalue=-1, rhsblocks=[[0, 1], [0, 2]])
>>>
>>> sparams = qs.allcomputesparameters([port1, port2], [qs.sn(1), qs.sn(1)], sols, lumpfields)
See Also

alleigenport(), modedrive(), formulation.allsolve(), formulation.lump(), printsparameters()

alleigenport

def alleigenport(
portphysreg: int,
numeigenvalues: int,
targetbeta: float,
E: field,
mu: expression,
eps: expression,
portnormal: list[float],
ddmrelrestol: float,
ddmmaxnumit: int,
drivingfrequency: float = -1.0,
eigentol: float = 1e-06,
eigenmaxnumits: int = 1000,
integrationorder: int = 5,
verbosity: int = 1
) -> list[list[expression]]
def alleigenport(
portphysreg: int,
numeigenvalues: int,
targetreal: float,
targetimag: float,
E: field,
mu_real: expression,
mu_imag: expression,
eps_real: expression,
eps_imag: expression,
sigma_real: expression,
sigma_imag: expression,
portnormal: list[float],
ddmrelrestol: float,
ddmmaxnumit: int,
drivingfrequency: float = -1.0,
eigentol: float = 1e-06,
eigenmaxnumits: int = 1000,
integrationorder: int = 5,
verbosity: int = 1,
errorifzanisotropic: bool = True
) -> list[list[expression]]

This function computes the eigenmodes of a waveguide whose cross section SS is the given physical region portphysreg. The material parameters μ\mu, ϵ\epsilon, and σ\sigma 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 EE and HH and the propagation constant γ=α+jβ\gamma=\alpha+j\beta as a list of expressions. The returned modes are scaled such that the power into the waveguide,

P=12SRe(E×H)dS,P = \frac12\int_{S}\textup{Re}(E\times H^{*})\cdot \textup{d}S,

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.

Example
>>> # Load a mesh with physical regions:
>>> wg1_crosssection = 9; wg2_crossection = 10;
>>> portphysreg = 11; full_waveguide = 12;
>>> portboundary = 13; skin = 14;
>>> mymesh = qs.mesh()
>>> mymesh.selectskin(portboundary, portphysreg)
>>> mymesh.selectskin(skin)
>>> mymesh.partition()
>>> mymesh.load("gmsh:waveguide.msh", skin, 1, 1)
>>>
>>> # Set input parameter values:
>>> E = qs.field("hcurl")
>>> E.setorder(full_waveguide, 2)
>>> E.setconstraint(portboundary)
>>> epsr = qs.parameter(); mur = qs.parameter();
>>> epsr.setvalue(wg1_crosssection, 2.25)
>>> epsr.setvalue(wg2_crossection, 1.0)
>>> mur.setvalue(portphysreg, 1.0);
>>>
>>> # Attempt to find 4 modes, targetting the eigenvalue 3e6 * i:
>>> modes = qs.alleigenport(portphysreg, 4, 0, 3e6, E, mur * qs.getmu0(), 0.0, epsr * qs.getepsilon0(), 0.0, 0.0, 0.0, [0.0, 0.0, 1.0], 1e-6, 100, 1.4314e14, 1e-8)
>>>
>>> for i in range(len(modes)):
>>> # The fields are returned in the following order:
>>> modes[i][0].write(portphysreg, 'Etreal_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][1].write(portphysreg, 'Etimag_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][2].write(portphysreg, 'Elreal_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][3].write(portphysreg, 'Elimag_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][4].write(portphysreg, 'Htreal_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][5].write(portphysreg, 'Htimag_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][6].write(portphysreg, 'Hlreal_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> modes[i][7].write(portphysreg, 'Hlimag_rank' + str(qs.getrank()) + '_mode' + str(i) + '.vtu', 2)
>>> # Evaluate the eigenvalue expression to access the value:
>>> if (qs.getrank() == 0):
>>> alpha = modes[i][8].evaluate()
>>> beta = modes[i][9].evaluate()
>>> print('Mode ' + str(i) + ' has eigenvalue ' + str(alpha) + ' + ' + str(beta) + ' * i.')
See Also

modedrive(), rectangularport()

allintegrate

def allintegrate(
physreg: int,
expr: expression,
integrationorder: int
) -> float

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1"); x=field("x")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, 12*x)
>>>
>>> integratedvalue = allintegrate(vol, v, 4)
See Also

allinterpolate(), allprobe()

allinterpolate

def allinterpolate(
physreg: int,
expr: expression,
xyzcoord: list[float]
) -> float

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1"); x=field("x")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, 12*x)
>>> xyzcoord = [0.5,0.6,0.05]
>>>
>>> interpolated = allinterpolate(vol, v, xyzcoord)
>>> interpolated
[5.954696754335098]
See Also

allintegrate(), allprobe()

alllaplace

def alllaplace(
physreg: int,
regionv1: int,
regionv0: int,
coef: expression,
vorder: int
) -> field

allmax

def allmax(
physreg: int,
expr: expression,
refinement: int
) -> float

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1"); x = field("x")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, 12*x)
>>> maxdata = allmax(vol, v, 1)
>>> maxdata
12.0
See Also

allmin()

allmeasuredistance

def allmeasuredistance(
a: vec,
b: vec,
c: vec
) -> float

This returns a relative L2 norm according to the below formula:

abc\frac{|\boldsymbol{a} - \boldsymbol{b}|}{|\boldsymbol{c}|}
Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>>
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v), 0, 2)
>>> projection += integral(vol, dt(v)*tf(v) - 2*tf(v))
>>> projection += integral(vol, dtdt(dof(v))*tf(v))
>>>
>>> v_prev = vec(projection) # previous iteration solution
>>> v_curr = vec(projection) # current iteration solution
>>>
>>> relativeerror = 1.0
>>> while(relativeerror < 1e-5):
>>> projection.allsolve(1e-6, 500)
>>> v_curr.setdata()
>>> relativeerror = allmeasuredistance(v_curr, v_prev, v_curr)
>>> v_prev = v_curr.copy()

allmin

def allmin(
physreg: int,
expr: expression,
refinement: int
) -> float

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1"); x = field("x")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, 12*x)
>>> mindata = allmin(vol, v, 1)
>>> mindata
-2.0
See Also

allmax()

allpartition

def allpartition(meshfile: str, skinphysreg: int = -1) -> str

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 function will be deprecated. Use mesh.partition() instead.

Example
>>> initialize()
>>> partfilename = allpartition("disk.msh")
>>> finalize()

allprobe

def allprobe(physreg: int, expr: expression) -> float

This functions returns the value of a scalar expression expr at the point region physreg.

Example
>>> vol=1; anynode=12
>>> mymesh = mesh()
>>> mymesh.selectanynode(anynode)
>>> mymesh.load("disk.msh")
>>>
>>> x = field("x")
>>> y = field("y")
>>> probedvalue = allprobe(anynode, 2*x+y)

allsolve

def allsolve(
relrestol: float,
maxnumit: int,
nltol: float,
maxnumnlit: int,
relaxvalue: float,
formuls: list[formulation],
verbosity: int = 1
) -> int

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.

Example
>>> ...
>>> allsolve(1e-8, 500, 1e-4, 1.0, [electrostatics])
See Also

solve(), formulation.solve(), formulation.allsolve()

andpositive

def andpositive(exprs: list[expression]) -> expression

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x=field("x"); y=field("y"); z=field("z")
>>> expr = andpositive([x,y]) # At points, where x>=0 and y>=0, value is 1, otherwise -1.
>>> expr.write(vol, "andpositive.vtk", 1)
See Also

ifpositive(), orpositive()

array1x1

def array1x1(arg0: expression) -> expression

This defines a vector or matrix operation of size 1×11\times1. The array is populated in a row-major way.

[arg0]\begin{bmatrix} arg0 \end{bmatrix}

array1x2

def array1x2(arg0: expression, arg1: expression) -> expression

This defines a vector or matrix operation of size 1×21\times2. The array is populated in a row-major way.

[arg0arg1]\begin{bmatrix} arg0 & arg1 \end{bmatrix}
Example
>>> myarray = array1x2(1,2)

array1x3

def array1x3(
arg0: expression,
arg1: expression,
arg2: expression
) -> expression

This defines a vector or matrix operation of size 1×31\times3. The array is populated in a row-major way.

[arg0arg1arg2]\begin{bmatrix} arg0 & arg1 & arg2 \end{bmatrix}
Example
>>> myarray = array1x3(1,2,3)

array2x1

def array2x1(arg0: expression, arg1: expression) -> expression

This defines a vector or matrix operation of size 2×12\times1. The array is populated in a row-major way.

[arg0arg1]\begin{bmatrix} arg0 \cr arg1 \end{bmatrix}
Example
>>> myarray = array2x1(1, 2)

array2x2

def array2x2(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression
) -> expression

This defines a vector or matrix operation of size 2×22\times2. The array is populated in a row-major way.

[arg0arg1arg2arg3]\begin{bmatrix} arg0 & arg1 \cr arg2 & arg3 \end{bmatrix}
Example
>>> myarray = array2x2(1,2, 3,4)

array2x3

def array2x3(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression,
arg4: expression,
arg5: expression
) -> expression

This defines a vector or matrix operation of size 2×32\times3. The array is populated in a row-major way.

[arg0arg1arg2arg3arg4arg5]\begin{bmatrix} arg0 & arg1 & arg2 \cr arg3 & arg4 & arg5 \end{bmatrix}
Example
>>> myarray = array2x3(1,2,3, 4,5,6)

array3x1

def array3x1(
arg0: expression,
arg1: expression,
arg2: expression
) -> expression

This defines a vector or matrix operation of size 3×13\times1. The array is populated in a row-major way.

[arg0arg1arg2]\begin{bmatrix} arg0 \cr arg1 \cr arg2 \end{bmatrix}
Example
>>> myarray = array3x1(1, 2, 3)

array3x2

def array3x2(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression,
arg4: expression,
arg5: expression
) -> expression

This defines a vector or matrix operation of size 3×23\times2. The array is populated in a row-major way.

[arg0arg1arg2arg3arg4arg5]\begin{bmatrix} arg0 & arg1 \cr arg2 & arg3 \cr arg4 & arg5 \end{bmatrix}
Example
>>> myarray = array3x2(1,2, 3,4, 5,6)

array3x3

def array3x3(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression,
arg4: expression,
arg5: expression,
arg6: expression,
arg7: expression,
arg8: expression
) -> expression

This defines a vector or matrix operation of size 3×33\times3. The array is populated in a row-major way.

[arg0arg1arg2arg3arg4arg5arg6arg7arg8]\begin{bmatrix} arg0 & arg1 & arg2 \cr arg3 & arg4 & arg5 \cr arg6 & arg7 & arg8 \end{bmatrix}
Example
>>> myarray = array3x3(1,2,3, 4,5,6, 7,8,9)

asin

def asin(input: expression) -> expression

This returns an expression that is the arcsinarcsin or sin1sin^{-1} of input. The output expression is in radians.

Example
>>> expr = asin(sqrt(0.5))
>>> expr.print(); # in radians
Expression size is 1x1
@ row 0, col 0 :
0.785398
>>>
>>> (expr * 180/getpi()).print(); # in degrees
Expression size is 1x1
@ row 0, col 0 :
45
See Also

sin(), cos(), tan(), acos(), atan()

atan

def atan(input: expression) -> expression

This returns an expression that is the arctanarctan or tan1tan^{-1} of input. The output expression is in radians.

Example
>>> expr = atan(0.57735)
>>> expr.print(); # in radians
Expression size is 1x1
@ row 0, col 0 :
0.523599
>>>
>>> (expr * 180/getpi()).print(); # in degrees
Expression size is 1x1
@ row 0, col 0 :
30
See Also

sin(), cos(), tan(), asin(), acos()

bode

def bode(reals: list[float], imags: list[float]) -> list[list[float]]

This function returns the magnitudes (unit: dB) and angles (unit: degrees) of the given complex numbers.

Example
>>> magnitudes = qs.bode([1, 1, 0], [0, 1, 1])[0]
>>> angles = qs.bode([1, 1, 0], [0, 1, 1])[1]
>>> angles[1]
45.0

cn

def cn(n: float) -> expression

This function takes as an argument the fundamental frequency multiplier. It is a shortform for cos(n2πft)cos(n * 2\pi f t).

Example
>>> f = 15.5 * 1e+9 # fundamnetal frequency
>>> setfundamentalfrequency(f)
>>> drivingsignal = cn(2) # same as cos(n* 2*getpi*f*t())
See Also

sn()

comp

def comp(selectedcomp: int, input: expression) -> expression

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> u.setvalue(vol, array3x1(10,20,30))
>>> vecexpr = 2*(u+u)
>>> comp(0, vecexpr).write(vol, "comp.vtk", 1)
See Also

compx(), compy(), compz()

complexdivision

def complexdivision(a: list[expression], b: list[expression]) -> list[expression]

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
>>> quot = qs.complexdivision([-54, 23], [2, 7])
>>> print(str(quot[0].evaluate()) + " + i * " + str(quot[1].evaluate()))
1.0 + i * 8.0

complexinverse

def complexinverse(a: list[expression]) -> list[expression]

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
>>> inv = qs.complexinverse([1, 2])
>>> print(str(inv[0].evaluate()) + " + i * " + str(inv[1].evaluate()))
0.2 + i * -0.4

complexproduct

def complexproduct(a: list[expression], b: list[expression]) -> list[expression]

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
>>> prod = qs.complexproduct([1, 8], [2, 7])
>>> print(str(prod[0].evaluate()) + " + i * " + str(prod[1].evaluate()))
-54.0 + i * 23.0

computeradiationpattern

def computeradiationpattern(
skinregion: int,
region: int,
E: field,
mu: float,
epsilon: float,
numelems: int
) -> iodata
def computeradiationpattern(
skinregion: int,
region: int,
p: field,
c: float,
numelems: int
) -> iodata

compx

def compx(input: expression) -> expression

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().

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> u.setvalue(vol, array3x1(10,20,30))
>>> vecexpr = 2*(u+u)
>>> compx(vecexpr).write(vol, "compx.vtk", 1)
See Also

comp(), compy(), compz()

compy

def compy(input: expression) -> expression

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().

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> u.setvalue(vol, array3x1(10,20,30))
>>> vecexpr = 2*(u+u)
>>> compx(vecexpr).write(vol, "compx.vtk", 1)
See Also

comp(), compx(), compz()

compz

def compz(input: expression) -> expression

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().

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> u.setvalue(vol, array3x1(10,20,30))
>>> vecexpr = 2*(u+u)
>>> compx(vecexpr).write(vol, "compx.vtk", 1)
See Also

comp(), compx(), compy()

continuitycondition

def continuitycondition(
gamma1: int,
gamma2: int,
u1: field,
u2: field,
errorifnotfound: bool = True,
lagmultorder: int = 0
) -> list[integration]
def continuitycondition(
gamma1: int,
gamma2: int,
u1: field,
u2: field,
rotcent: list[float],
rotangz: float,
angzmod: float,
factor: float,
lagmultorder: int = 0
) -> list[integration]

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)

>>> ...
>>> u1 = field("h1xyz")
>>> u2 = field("h1xyz")
>>>
>>> elasticity = formulation()
>>> ...
>>> elasticity += continuitycondition(gamma1, gamma2, u1, u2)

This returns the formulation terms required to enforce u1=u2u_1 = u_2 between boundary region Γ1\Gamma_1 and Γ2\Gamma_2 (with Γ1Γ2\Gamma_1 \subseteq \Gamma_2, meshes can be non-matching). In case Γ2\Gamma_2 is larger than Γ1\Gamma_1 (Γ1Γ2\Gamma_1 \subset \Gamma_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)

>>> ...
>>> # Rotor-stator interface
>>> rotorside=11; statorside=12
>>> ...
>>> # Rotor rotation around z axis
>>> alpha = 30.0
>>> ...
>>> az = field("h1")
>>> ...
>>> magentostatics = formulation()
>>> ...
>>> magnetostatics += continuitycondition(statorside, rotorside, az, az, [0,0,0], alpha, 45.0, -1.0)

This returns the formulation terms required to enforce field continuity across an angzmodangzmod degrees slice of a rotor-stator interface where the rotor geometry is rotated by rotangzrotangz around the zz axis with rotation center at rotcentrotcent. 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 -11 for antiperiodicity. Boundary Γ1\Gamma_1 is the rotor-stator interface on the (non-moving) stator side while the boundary Γ2\Gamma_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 xx axis.


The condition is based on a Lagrange multiplier of the same type and the same harmonic content as the field u1u_1 and u2u_2. The mortar finite element method is used to link the unknown dofdof field on Γ1\Gamma_1 and Γ2\Gamma_2 so that there is no restriction on the mesh used for both regions.

See Also

periodicitycondition(), symmetrycondition()

cos

def cos(input: expression) -> expression

This returns an expression that is the coscos of input. The input expression is in radians.

Example
>>> expr = cos(getpi())
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
-1
>>>
>>> expr = cos(1)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
0.540302
See Also

sin(), tan(), asin(), acos(), atan()

crossproduct

def crossproduct(a: expression, b: expression) -> expression

This computes the cross-product of two vector expressions. The returned expression is a vector.

c=a×b\boldsymbol{c} = \boldsymbol{a} \times \boldsymbol{b}
Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; top=3
>>> E = field("hcurl")
>>> n = normal(vol)
>>> cp = crossproduct(E, n)

curl

def curl(input: expression) -> expression

This computes the curl of a vector expression. The returned expression is a vector.

w=×u\boldsymbol{w} = \nabla \times \boldsymbol{u}
Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1);
>>> curl(u).write(vol, "curlu.vtk", 1)
See Also

grad(), div()

d

def d(
expr: expression,
var: expression,
flds: list[field],
deltas: list[expression]
) -> expression

This functions approximates the derivative of an expression when its variable is moved by a delta.

Example

dbtoneper

def dbtoneper(toconvert: expression) -> expression

This converts the expression toconvert from a dBdB units to NepersNepers.

Example
>>> neperattenuation = dbtoneper(100.0)

determinant

def determinant(input: expression) -> expression

This returns the determinant of a square matrix.

Example
>>> matexpr = expression(3,3, [1,2,3, 6,5,4, 8,9,7])
>>> detmat = determinant(matexpr)
>>> detmat.print()
Expression size is 1x1
@ row 0, col 0 :
21
See Also

inverse()

detjac

def detjac() -> expression

This returns the determinant of the Jacobian matrix.

Example
>>> detJ = detjac()
See Also

jac(), invjac()

div

def div(input: expression) -> expression

This computes the divergence of a vector expression. The returned expression is a scalar.

s=us = \nabla \cdot \boldsymbol{u}
Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1);
>>> div(u).write(vol, "divu.vtk", 1)
See Also

grad(), curl()

dof

def dof(input: expression) -> expression
def dof(input: expression, physreg: int) -> expression

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.

Examples

Example 1: dof(input:expression)

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))

Example 2: dof(input:expression, physreg:int)

>>> ...
>>> projection += integral(vol, dof(v, vol)*tf(v) - 2*tf(v))
See Also

tf()

doubledotproduct

def doubledotproduct(a: expression, b: expression) -> expression

This computes the double-dot product of two matrix expressions. The returned expression is a scalar.

A:B=i,jAijBij\boldsymbol{A:B} = \sum_{i,j} A_{ij} B_{ij}
Example
>>> a = array2x2(1,2, 3,4)
>>> b = array2x2(11,12, 13, 14)
>>> addotb = doubledotproduct(a, b)
>>> addotb.print()
Expression size is 1x1
@ row 0, col 0 :
130

dt

def dt(input: expression) -> expression
def dt(
input: expression,
initdt: float,
initdtdt: float
) -> expression

This returns the first-order time derivative expression.

Examples

Example 1: dt(input:expression)

>>> ...
>>> setfundamentalfrequency(50)
>>> vmh = field("h1", [2,3])
>>> vmh.setorder(vol, 1)
>>> dt(vmh).write(vol, "dtv.vtk", 1)

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.

>>> ...
>>> dtapprox = dt(t()*t(), 0, 2)
See Also

dtdt(), dtdtdt(), dtdtdtdt()

dtdt

def dtdt(input: expression) -> expression
def dtdt(
input: expression,
initdt: float,
initdtdt: float
) -> expression

This returns the second-order time derivative expression.

Examples

Example 1: dtdt(input:expression)

>>> ...
>>> setfundamentalfrequency(50)
>>> vmh = field("h1", [2,3])
>>> vmh.setorder(vol, 1)
>>> dtdt(vmh).write(vol, "dtdtv.vtk", 1)

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.

>>> ...
>>> dtapprox = dtdt(t()*t(), 0, 2)
See Also

dt(), dtdtdt(), dtdtdtdt()

dtdtdt

def dtdtdt(input: expression) -> expression

This returns the third-order time derivative expression.

Example
>>> ...
>>> setfundamentalfrequency(50)
>>> vmh = field("h1", [2,3])
>>> vmh.setorder(vol, 1)
>>> dtdtdt(vmh).write(vol, "dtdtdtv.vtk", 1)
See Also

dt(), dtdt(), dtdtdtdt()

dtdtdtdt

def dtdtdtdt(input: expression) -> expression

This returns the fourth-order time derivative expression.

Example
>>> ...
>>> setfundamentalfrequency(50)
>>> vmh = field("h1", [2,3])
>>> vmh.setorder(vol, 1)
>>> dtdtdtdt(vmh).write(vol, "dtdtdtdtv.vtk", 1)
See Also

dt(), dtdt(), dtdtdt()

dx

def dx(input: expression) -> expression

This returns the xx space derivative expression.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> dx(v).write(vol, "dxv.vtk", 1)
See Also

dy(), dz()

dy

def dy(input: expression) -> expression

This returns the yy space derivative expression.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> dy(v).write(vol, "dyv.vtk", 1)
See Also

dx(), dz()

dz

def dz(input: expression) -> expression

This returns the zz space derivative expression.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> dz(v).write(vol, "dzv.vtk", 1)
See Also

dx(), dy()

elasticwavespeed

def elasticwavespeed(
rho: expression,
E: expression,
nu: expression
) -> list[expression]
def elasticwavespeed(rho: expression, H: expression) -> list[expression]

elementwiseproduct

def elementwiseproduct(a: expression, b: expression) -> expression

This computes the element-wise product of two matrix expressions a and b. The returned expression has the same size as the two input expressions.

Cij=AijBijC_{ij} = A_{ij} B_{ij}
Example
>>> a = array2x2(1,2, 3,4)
>>> b = array2x2(11,12, 13, 14)
>>> aelwb = elementwiseproduct(a, b)
>>> aelwb.print()
Expression size is 2x2
@ row 0, col 0 :
11
@ row 0, col 1 :
24
@ row 1, col 0 :
39
@ row 1, col 1 :
56

emwavespeed

def emwavespeed(mur: expression, epsr: expression) -> expression

entry

def entry(
row: int,
col: int,
input: expression
) -> expression

This gets the (row, col) entry in the input vector or matrix expression.

Example
>>> u = array3x2(1,2, 3,4, 5,6)
>>> arrayentry_row2col0 = entry(2, 0, u) # entry from third row (index=2), first column (index=0)
>>> arrayentry_row2col0.print()
Expression size is 1x1
@ row 0, col 0 :
5

exp

def exp(input: expression) -> expression

This returns an exponential function of base ee:

einpute^{input}
Example
>>> expr = exp(2)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
7.38906
See Also

pow()

eye

def eye(size: int) -> expression

This returns a size x size identity matrix.

Example
>>> II = eye(2)
>>> II.print()
Expression size is 2x2
@ row 0, col 0 :
1
@ row 0, col 1 :
0
@ row 1, col 0 :
0
@ row 1, col 1 :
1

fieldorder

def fieldorder(
input: field,
alpha: float = -1.0,
absthres: float = 0.0
) -> expression

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
>>> mymesh = mesh("disk.msh")
>>> vol = 1; sur = 2; top = 3
>>> v = field("h1")
>>> v.setorder(vol, 2)
>>> fo = fieldorder(v)
>>> fo.write(vol, "fieldorder.vtk", 1)

getc0

def getc0() -> float

getcirculationports

def getcirculationports(harmonicnumbers: list[int] = []) -> tuple[list[port], list[expression]]

getcirculationsources

def getcirculationsources(circulations: list[expression], inittimederivatives: list[list[float]] = []) -> list[expression]

getdimension

def getdimension(physreg: int) -> int

This returns the x, y and z mesh dimensions.

Example
>>> mymesh = mesh("disk.msh")
>>> dims = mymesh.getdimensions()

getepsilon0

def getepsilon0() -> float

This returns the value of vacuum permittivity ϵ0=8.854187812813e12Fm1\epsilon_0 = 8.854187812813e-12 Fm^{-1}.

Example
>>> eps0 = getepsilon0()
8.854187812813e-12

getextrusiondata

def getextrusiondata() -> list[expression]

This gives the relative depth δ\delta in the extruded layer, the extrusion normal n\boldsymbol{n} and tangents t1\boldsymbol{t_1} and t2\boldsymbol{t_2}. This is useful when creating Perfectly Matched Layers (PMLs).

Example

We use the disk.msh for the example here.

>>> vol=1; sur=2; top=3; circle=4 # physical regions defined in disk.msh
>>> mymesh = mesh()
>>>
>>> # predefine extrusion
>>> volextruded=5; bndextruded=6; # physical regions that will be utilized in extrusion
>>> mymesh.extrude(volextruded, bndextruded, sur, [0.1,0.05])
>>> mymesh.load("disk.msh") # extrusion is performed when the mesh is loaded.
>>> mymesh.write("diskpml.msh")
>>> pmldata = getextrusiondata()
See Also

mesh.extrude()

getharmonic

def getharmonic(
harmnum: int,
input: expression,
numfftharms: int = -1
) -> expression

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
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>> v = field("h1", [2])
>>> v.setorder(vol, 1)
>>> v.harmonic(2).setvalue(vol, 1)
>>> constcomp = getharmonic(1, abs(v), 10)
>>> constcomp.write(vol, "constcomp.pos", 1)

getmu0

def getmu0() -> float

This returns the value of vacuum permeability μ0=1.2566370621219e06NA2\mu_0 = 1.2566370621219e-06 NA^{-2}.

Example
>>> mu0 = getmu0()
1.2566370621219e-06

getpi

def getpi() -> float

This returns value of π\pi.

Example
>>> pi = getpi()
3.141592653589793

getrandom

def getrandom() -> float

This returns a random value uniformly distributed between 0.0 and 1.0.

Example
>>> rnd = getrandom()

getrank

def getrank() -> int

This returns the rank of the current process.

Example
>>> import quanscient as qs
>>> rank = qs.getrank()
See Also

count

getsstkomegamodelconstants

def getsstkomegamodelconstants() -> list[float]

getsubversion

def getsubversion() -> int

gettime

def gettime() -> float

This gets the value of the time variable t.

Example
>>> settime(1e-3)
>>> gettime()
0.001
See Also

settime(), t()

gettotalforce

def gettotalforce(
physreg: int,
EorH: expression,
epsilonormu: expression,
extraintegrationorder: int = 0
) -> list[float]
def gettotalforce(
physreg: int,
meshdeform: expression,
EorH: expression,
epsilonormu: expression,
extraintegrationorder: int = 0
) -> list[float]

This returns the components of the total magnetostatic/electrostatic force acting on a given region. In the axisymmetric case zero xx and zz components are returned and the yy component includes a 2π2 \pi factor to provide the force acting on the corresponding 3D shape. Units are ”NN per unit depth” in 2D and NN in 3D and 2D axisymmetry.

Examples

Example 1: gettotalforce(physreg:int, EorH:expression, epsilonormu:expression, extraintegrationorder:int=0)

>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> phi = field("h1")
>>> phi.setorder(vol, 2)
>>>
>>> mu0 = 4 * getpi() * 1e-7
>>> mu = parameter()
>>> mu.setvalue(vol, mu0)
>>> totalforce = gettotalforce(vol, -grad(phi), mu)

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.

>>> ...
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> totalforce = gettotalforce(vol, u, -grad(phi), mu)
See Also

printtotalforce()

getturbulentpropertiessstkomegamodel

def getturbulentpropertiessstkomegamodel(
v: expression,
rho: expression,
viscosity: expression,
walldistance: expression,
kp: expression,
omegap: expression,
logomega: expression
) -> list[expression]

getversion

def getversion() -> int

getversionname

def getversionname() -> str

getx

def getx() -> field

This returns the xx coordinate.

Example

An expression for distance can be calculated as follows:

>>> x = getx()
>>> y = gety()
>>> z = getz()
>>> d = sqrt(x*x+y*y+z*z)
See Also

gety(), getz()

gety

def gety() -> field

This returns the yy coordinate.

Example

In CFD applications, a parabolic inlet velocity profile can be prescribed as follows:

>>> y = gety() # y coordinate
>>> U = 0.3 # maximum velocity
>>> h = 0.41 # height of the domain
>>> u = 4*U*y(h-y)/(h*h) # parabolic inlet velocity profile expression
See Also

getx(), getz()

getz

def getz() -> field

This returns the zz coordinate.

Example

An expression for distance can be calculated as follows:

>>> x = getx()
>>> y = gety()
>>> z = getz()
>>> d = sqrt(x*x+y*y+z*z)
See Also

getx(), gety()

grad

def grad(input: expression) -> expression

For a scalar input expression, this is mathematically treated as the gradient of a scalar (v\nabla{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\nabla{\boldsymbol{u}}) and the output has one row per component of the input and one column per space derivative.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1);
>>> grad(v).write(vol, "gradv.vtk", 1)
See Also

div(), curl()

greenlagrangestrain

def greenlagrangestrain(input: expression) -> expression

This defines the (nonlinear) Green-Lagrange strains in Voigt form (ϵxx,ϵyy,ϵzz,γyz,γxz,γxy)(\epsilon_{xx},\epsilon_{yy},\epsilon_{zz},\gamma_{yz}, \gamma_{xz},\gamma_{xy}). The input can either be the displacement field or its gradient.

ϵij=12[uiXj+ujXi+uiXiuiXj+ujXiujXj+ukXiukXj] \epsilon_{ij} = \frac{1}{2} \left[ \frac{\partial u_i}{\partial X_j} + \frac{\partial u_j}{\partial X_i} + \frac{\partial u_i}{\partial X_i} \frac{\partial u_i}{\partial X_j} + \frac{\partial u_j}{\partial X_i} \frac{\partial u_j}{\partial X_j} + \frac{\partial u_k}{\partial X_i} \frac{\partial u_k}{\partial X_j} \right]
Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> glstrain = greenlagrangestrain(u)
>>> glstrain.print()
See Also

strain()

grouptimesteps

def grouptimesteps(
filename: str,
filestogroup: list[str],
timevals: list[float]
) -> None
def grouptimesteps(
filename: str,
fileprefix: str,
firstint: int,
timevals: list[float]
) -> None

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])

>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> v.write(vol, "v1.vtu", 1)
>>> v.write(vol, "v2.vtu", 1)
>>> grouptimesteps("v.pvd", ["v1.vtu", "v2.vtu"], [0.0, 10.0])

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.

>>> ...
>>> grouptimesteps("v.pvd", "v", 1, [0.0, 10.0])

harm

def harm(
harmnum: int,
input: expression,
numfftharms: int = -1
) -> expression

ifpositive

def ifpositive(
condexpr: expression,
trueexpr: expression,
falseexpr: expression
) -> expression

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1
>>> x=field("x"); y=field("y"); z=field("z")
>>> expr = ifpositive(x+y, 1, -1) # At points, where x+y>=0, value is 1, otherwise -1.
>>> expr.write(vol, "ifpositive.vtk", 1)
See Also

andpositive(), orpositive()

insertrank

def insertrank(name: str) -> str

integral

def integral(
physreg: int,
tointegrate: expression,
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> integration
def integral(
physreg: int,
meshdeform: expression,
tointegrate: expression,
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> integration
def integral(
physreg: int,
numcoefharms: int,
tointegrate: expression,
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> integration
def integral(
physreg: int,
numcoefharms: int,
meshdeform: expression,
tointegrate: expression,
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> integration
def integral(
physreg: int,
tointegrate: tuple[expression, preconditioner],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
meshdeform: expression,
tointegrate: tuple[expression, preconditioner],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
numcoefharms: int,
tointegrate: tuple[expression, preconditioner],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
numcoefharms: int,
meshdeform: expression,
tointegrate: tuple[expression, preconditioner],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
tointegrate: list[tuple[expression, int, preconditioner]],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
meshdeform: expression,
tointegrate: list[tuple[expression, int, preconditioner]],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
numcoefharms: int,
tointegrate: list[tuple[expression, int, preconditioner]],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]
def integral(
physreg: int,
numcoefharms: int,
meshdeform: expression,
tointegrate: list[tuple[expression, int, preconditioner]],
integrationorderdelta: int = 0,
blocknumber: int = 0
) -> list[tuple[integration, preconditioner]]

inverse

def inverse(input: expression) -> expression

This returns the inverse of a square matrix.

Example
>>> matexpr = array2x2(1,2, 3,4)
>>> invmat = inverse(matexpr)
>>> invmat.print()
Expression size is 2x2
@ row 0, col 0 :
-2
@ row 0, col 1 :
1
@ row 1, col 0 :
1.5
@ row 1, col 1 :
-0.5
>>>
>>> matexpr = expression(3,3, [1,2,3, 6,5,4, 8,9,7])
>>> invmat = inverse(matexpr)
See Also

determinant()

invjac

def invjac(row: int, col: int) -> expression
def invjac() -> expression

Example 1: invjac(row:int, col:int)

This returns the inverse of the Jacobian matrix at the entry (row, col).

>>> Ji = invjac(1,2)

Example 2: invjac()

This returns the whole 3×33 \times 3 Jacobian matrix.

>>> Ji = invjac()
See Also

jac(), detjac()

isdefined

def isdefined(physreg: int) -> bool

This checks if a physical region physreg is defined.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; circle=4;
>>> isdefined(sur) # returns True as the physical region sur=2 is defined
True
>>>
>>> isdefined(5) # returns False as the mesh loaded has no physical region 5
>>> False
See Also

isempty(), isinside(), istouching()

isempty

def isempty(physreg: int) -> bool

This checks if a physical region physreg is empty.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; circle=4;
>>> isempty(sur)
False
See Also

isdefined(), isinside(), istouching()

isinside

def isinside(physregtocheck: int, physreg: int) -> bool

This checks if a physical region is fully included in another region. The physreg is the physical region with which physregtocheck is checked.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; circle=4;
>>> isinside(sur, vol)
True
>>>
>>> isinside(vol, sur)
False
See Also

isdefined(), isempty(), istouching()

istouching

def istouching(physregtocheck: int, physreg: int) -> bool

This checks if a region is touching another region. The physreg is the physical region with which physregtocheck is checked.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; circle=4;
>>> istouching(sur, top)
True
See Also

isdefined(), isempty(), isinside()

jac

def jac(row: int, col: int) -> expression
def jac() -> expression

Example 1: jac(row:int, col:int)

This returns the inverse of the Jacobian matrix at the entry (row, col).

>>> J = jac(1,2)

Example 2: jac()

This returns the whole 3×33 \times 3 Jacobian matrix.

>>> J = jac()
See Also

invjac(), detjac()

linspace

def linspace(
a: float,
b: float,
num: int
) -> list[float]

This gives a vector of num equally spaced values from a to b. The space between each values is calculated as:

banum1\frac{b - a}{num -1}
Example
>>> vals = linspace(0.5, 2.5, 5)
>>> printvector(vals)
0.5 1 1.5 2 2.5
See Also

logspace()

loadshape

def loadshape(meshfile: str) -> list[list[shape]]

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
>>> diskshapes = loadshape("disk.msh")
>>>
>>> # Add a thin slice on top of the disk (diskshapes[2][1] is the top face of the disk)
>>> thinslice = diskshapes[2][1].extrude(5, 0.02, 2)
>>>
>>> mymesh = mesh([diskshapes[2][0], diskshapes[2][1], diskshapes[3][0], thinslice])
>>> mymesh.write("editeddisk.msh")

loadvector

def loadvector(
filename: str,
delimiter: str = ',',
sizeincluded: bool = False
) -> list[float]

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.

Example
>>> v = [2.4, 3.14, -0.1]
>>> writevector("vecvals.txt", v, '\n', True)
>>> vloaded = loadvector("vecvals.txt", \n', True)
>>> vloaded
[2.4, 3.14, -0.1]
See Also

printvector(), writevector()

log

def log(input: expression) -> expression

This returns an expression that is the natural logarithm of the input expression.

Example
>>> expr = log(2);
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
0.693147

log10

def log10(input: expression) -> expression

logspace

def logspace(
a: float,
b: float,
num: int,
basis: float = 10.0
) -> list[float]

This returns the basis to the power of each of the num values in the linspace.

Example
>>> vals = logspace(1, 3, 3) # resulting linspace: 1, 2, 3
>>> printvector(vals) # 10^1, 10^2, 10^3
10 100 1000
See Also

linspace()

makeharmonic

def makeharmonic(harms: list[int], exprs: list[expression]) -> expression
def makeharmonic(
harms: list[int],
expr: expression,
numfftharms: int
) -> expression

This returns a multi-harmonic expression whose harmonic numbers and expressions are provided as arguments. The argument expressions must be on harmonic 1.

Example
>>> ...
>>> harmexpr = makeharmonic([1,2,4], [11, v.harmonic(2), 14])
>>> harmexpr.write(vol, "harmexpr.pos", 1)

max

def max(a: expression, b: expression) -> expression
def max(a: field, b: field) -> expression
def max(a: parameter, b: parameter) -> expression

This returns an expression whose value is the maximum of the two input arguments a and b.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x=field("x"); y=field("y"); z=field("z")
>>> max(x,y).write(vol, "max.pos", 1)
See Also

min()

meshsize

def meshsize(integrationorder: int) -> expression

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 < 00 a RuntimeError is raised.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1; sur = 2; top = 3
>>> h = meshsize(0)
>>> h.write(top, "meshsize.vtk", 1)

min

def min(a: expression, b: expression) -> expression
def min(a: field, b: field) -> expression
def min(a: parameter, b: parameter) -> expression

This returns an expression whose value is the minimum of the two input arguments a and b.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x=field("x"); y=field("y"); z=field("z")
>>> min(x,y).write(vol, "min.pos", 1)
See Also

max()

mod

def mod(input: expression, modval: float) -> expression

This is a modulo function. This returns an expression equal to the remainder resulting from the division of input by modval.

Example
>>> expr = mod(10, 9)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
1
>>> expr = mod(99, 100)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
99
>>> expr = mod(2.55, 0.6)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
0.15

modedrive

def modedrive(
portphysreg: int,
drivesignal: expression,
E: field,
Etreal: expression,
Etimag: expression,
Htreal: expression,
Htimag: expression,
portnormal: list[float],
lumpfield: field,
blocktag: int,
extraintegrationorder: int = 0
) -> list[integration]

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)\sin(\omega t + \text{shift}) or cos(ωt+shift)\cos(\omega t + \text{shift}), where ω\omega is the driving frequency and shift\text{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.

Example
>>> # Suppose a mesh has been loaded with physical regions (port1, port2, waveguide, pec_boundary) and E, mur, epsr_real, and eps_imag have been defined so that the following calls make sense:
>>> port1mode0 = qs.alleigenport(port1, 1, 0, 1000, E, mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, [0.0, 0.0, 1.0], 1e-06, 1000)[0]
>>> port2mode0 = qs.alleigenport(port2, 1, 0, 1000, E, mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, [0.0, 0.0, -1.0], 1e-06, 1000)[0]
>>>
>>> form = qs.formulation()
>>>
>>> # Prepare the lumpfield constraints:
>>> lumpfields = form.lump([port1, port2], [2, 3])
>>>
>>> # Electromagnetic waves
>>> form += qs.integral(waveguide, 3, qs.predefinedemwave(qs.dof(E), qs.tf(E), mur * qs.getmu0(), 0, epsr_real * qs.getepsilon0(), epsr_imag * qs.getepsilon0(), 0, 0, "oo2"))
>>>
>>> # Drive and absorb the mode of port1 (feeding term in block number 1):
>>> form += qs.modedrive(port1, qs.sn(1), E, port1mode0[0], port1mode0[1], port1mode0[4], port1mode0[5], [0.0, 0.0, 1.0], lumpfields[0], 1)
>>>
>>> # Drive and absorb the mode of port2 (feeding term in block number 2):
>>> form += qs.modedrive(port2, qs.sn(1), E, port2mode0[0], port2mode0[1], port2mode0[4], port2mode0[5], [0.0, 0.0, -1.0], lumpfields[1], 2)
>>>
>>> sols = form.allsolve(relrestol=1e-06, maxnumit=1000, nltol=1e-05, maxnumnlit=-1, relaxvalue=-1, rhsblocks=[[0, 1], [0, 2]])
See Also

alleigenport(), formulation.allsolve(), formulation.lump(), allcomputesparameters()

moveharmonic

def moveharmonic(
origharms: list[int],
destharms: list[int],
input: expression,
numfftharms: int = -1
) -> expression

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
>>> ...
>>> movedharm = moveharmonic([1,2], [5,3], 11+v)
>>> movedharm.write("vol", "movedharm.pos", 1)

norm

def norm(arg0: expression) -> expression

This gives the L2L2 norm of an expression input.

Example
>>> myvector = array3x1(1,2,3)
>>> normL2 = norm(myvector)
>>> normL2.print()
Expression size is 1x1
@ row 0, col 0 :
3.74166

normal

def normal() -> expression
def normal(pointoutofphysreg: int) -> expression

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
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3 # physical regions
>>> normal(vol).write(sur, "normal.vtk", 1)

on

def on(
physreg: int,
expression: expression,
errorifnotfound: bool = True
) -> expression
def on(
physreg: int,
coordshift: expression,
expression: expression,
errorifnotfound: bool = True
) -> expression

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=u2u_1 = u_2, at the interface Γ\Gamma of non-matching meshes. This can for example be achieved with a Lagrange multiplier λ\lambda such that

Γ( λ(u1u2)+λ(u1u2) )dΓ=0\int_{\Gamma} \left( \ \lambda (u_1^{\prime} - u_2^{\prime}) + {\lambda}^{\prime} (u_1 - u_2) \ \right) d \Gamma = 0

holds for any appropriate field u1u_1^{\prime}, u2u_2^{\prime} and λ{\lambda}^{\prime}. The example below illustrates the formulation terms needed to implement the Lagrange multiplier between interfaces Γ1\Gamma_1 and Γ2\Gamma_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.

>>> ...
>>> u1=field("h1xyz"); u2=field("h1xyz"); lambda=field("h1xyz")
>>> ...
>>> elasticity = formulation()
>>> ...
>>> elasticity += integral(gamma1, dof(lambda)*tf(u1) )
>>> elasticity += integral(gamma2, -on(gamma1, dof(lambda)) * tf(u2) )
>>> elasticity += integral(gamma1, (dof(u1) - on(gamma2, dof(u2)))*tf(lambda) )

When setting the flag errorifnotfound=False, any point in Γ1\Gamma_1 without a relative in Γ2\Gamma_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\Gamma_1 is without a relative in Γ2\Gamma_2 (and vice versa).

The case where there is no unknown dof \ term in the expression argument is described below:

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>> x=field("x"); y=field("y"); z=field("z"); v=field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(top, dof(v)*tf(v) - on(vol, dz(z))*tf(v))
>>> projection.solve()
>>> v.write(top, "dzz.pos", 1)

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 11 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)).

>>> ...
>>> projection += integral(top, dof(v)*tf(v) - on(vol, array3x1(2*x,2*y,2*z), dz(z))*tf(v))

orpositive

def orpositive(exprs: list[expression]) -> expression

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.

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x=field("x"); y=field("y"); z=field("z")
>>> expr = orpositive([x,y]) # At points, where x>=0 or y>=0, value is 1, otherwise -1.
>>> expr.write(vol, "orpositive.vtk", 1)
See Also

ifpositive(), andpositive()

periodicitycondition

def periodicitycondition(
gamma1: int,
gamma2: int,
u: field,
dat1: list[float],
dat2: list[float],
factor: float,
lagmultorder: int = 0
) -> list[integration]

This returns the formulation terms required to enforce on field uu a rotation or translation periodic condition between boundary region Γ1\Gamma_1 and Γ2\Gamma_2 (meshes can be non-conforming). A factor different than 11 can be provided to scale the field on Γ2\Gamma_2 (use -11 for antiperiodicity).

Example
>>> ...
>>> u = field("h1xyz")
>>> ...
>>> elasticity = formulation()
>>> ...
>>> # In case gamma2 is gamma1 rotated by (ax, ay, az) degrees around first the x, then y and then z-axis.
>>> # Rotation center is (cx, cy, cz)
>>> cx=0; cy=0; cz=0
>>> ax=0; ay=0; az=60
>>>
>>> elasticity += periodicitycondition(gamma1, gamma2, u, [cx,cy,cz], [ax,ay,az], 1.0)
>>>
>>> # In case gamma2 is gamma1 translated by a distance d in direction (nx, ny, nz)
>>> d=0.8
>>> nx=1; ny=0; nz=0
>>> elasticity += peridocitycondition(gamma1, gamma2, u, [nx,ny,nz], [d], 1.0)

The condition is based on a Lagrange multiplier of the same type and the same harmonic content as the field u1u_1 and u2u_2. The mortar finite element method is used to link the unknown dofdof field on Γ1\Gamma_1 and Γ2\Gamma_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.

See Also

continuitycondition(), symmetrycondition()

pow

def pow(base: expression, exponent: expression) -> expression

This is a power function. This returns an expression equal to base to the power exponent:

baseexponent{base}^{exponent}
Example
>>> expr = pow(2, 5)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
32
See Also

exp()

predefinedacousticradiation

def predefinedacousticradiation(
dofp: expression,
tfp: expression,
soundspeed: expression,
neperattenuation: expression
) -> expression

This function defines the equation for the Sommerfeld acoustic radiation condition

np+1cpt=0\partial_{\boldsymbol{n}} p + \frac{1}{c} \frac{\partial p}{\partial t} = 0

which forces the outgoing pressure waves at infinity: a pressure field of the form

p(r,t)=P cos(ωtkr)p(r, t) = P \ cos(\omega t - kr)

propagating in the direction er\boldsymbol{e_r} perpendicular to the truncation boundary indeed satisfies the Sommerfeld radiation condition since

np=erp=k P sin(ωtkr)=wc P sin(ωtkr)=1c(P cos(ωtkr))t \partial_{\boldsymbol{n}} p = \partial_{\boldsymbol{e}_r} p = k \ P \ sin(\omega t - kr) = \frac{w}{c} \ P \ sin(\omega t - kr) = - \frac{1}{c} \frac{\partial (P \ cos(\omega t -kr))}{\partial t}

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/mNeper/m) in case of harmonic problems. For convenience use the function dbtoneper() to convert dB/mdB/m attenuation values to Np/mNp/m.

Example
>>> ...
>>> acoustics += integral(sur, predefinedacousticradiation(dof(p), tf(p), 340, dbtoner(500)))

predefinedacousticstructureinteraction

def predefinedacousticstructureinteraction(
dofp: expression,
tfp: expression,
dofu: expression,
tfu: expression,
soundspeed: expression,
fluiddensity: expression,
normal: expression,
neperattenuation: expression,
scaling: float = 1.0
) -> expression

This function defines the bi-directional coupling for acoustic-structure interaction at the medium interface. Field pp is the acoustic pressure and field u\boldsymbol{u} is the mechanical displacement. Calling n\boldsymbol{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\boldsymbol{f}_{pressure} = -p \boldsymbol{n}

as well as linking the structure acceleration to the fluid pressure normal derivative using Newton’s law:

np=ρfluid2ut2n\partial_{\boldsymbol{n}} p = - \rho_{fluid} \frac{\partial^2 \boldsymbol{u}}{\partial t^2} \cdot \boldsymbol{n}

To have a good matrix conditioning a scaling factor ss (e.g s=1e10s = 1e10) can be provided. In this case, the pressure source is divided by ss and, to compensate, the pressure force is multiplied by ss. 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/mNeper/m) in case of harmonic problems. For convenience use the function dbtoneper() to convert dB/mdB/m attenuation values to Np/mNp/m.

Example
>>> ...
>>> u = field("h1xy", [2,3])
>>> u.setorder(sur, 2)
>>> acoustics += integral(left, predefinedacousticstructureinteraction(dof(p), tf(p), dof(u), tf(u), 340, 1.2, array2x1(1,0), dbtoneper(500), 1e10)

predefinedacousticwave

def predefinedacousticwave(
dofp: expression,
tfp: expression,
soundspeed: expression,
neperattenuation: expression,
precondtype: str = ''
) -> list[tuple[expression, int, preconditioner]]
def predefinedacousticwave(
dofp: expression,
tfp: expression,
soundspeed: expression,
neperattenuation: expression,
pmlterms: list[expression],
precondtype: str = ''
) -> list[tuple[expression, int, preconditioner]]

This function defines the equation for (linear) acoustic wave propagation:

2p1c22pt2=0\nabla^2 p - \frac{1}{c²} \frac{\partial^2 p}{\partial t^2} = 0

An acoustic attenuation value can be provided (in Neper/mNeper/m) in case of harmonic problems. For convenience use the function dbtoneper() to convert attenuation values from dB/mdB/m to Np/mNp/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/sm/s.
  • neperattenuation is the attenuation in Np/mNp/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.

>>> sur=1; left=2; wall=3
>>> h=10e-3; l=50e-3
>>> q = shape("quadrangle", sur, [0,0,0, l,0,0, l,h,0, 0,h,0], [250,50,250,50])
>>> ll = q.getsons()[3]
>>> ll.setphysicalregion(left)
>>> lwall = shape("union", wall, [q.getsons()[0], q.getsons()[1], q.getsons()[2]])
>>>
>>> mymesh = mesh([q, ll, lwall])
>>>
>>> setfundamentalfrequency(40e3)
>>>
>>> # Wave propogation requires both the in-phase (2) and quadrature (3) harmonics:
>>> p=field("h1", [2,3]); y=field("y")
>>> p.setorder(sur, 2)
>>> # In-phase only pressure source
>>> p.harmonic(2).setconstraint(left, y*(h-y)/(h*h/4))
>>> p.harmonic(3).setconstraint(left, 0)
>>> p.setconstraint(wall)
>>>
>>> acoustics = formulation()
>>> acoustics += integral(sur, predefinedacousticwave(dof(p), tf(p), 340, dbtoneper(500)))
>>> acoustics.solve()
>>> p.write(sur, "p.vtu", 2)
>>>
>>> # Write 50 timesteps for a time visualization
>>> # p.write(sur, "p.vtu", 2, 50)

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.

>>> ...
>>> pmlterms = [detDr, detDi, Dr, Di, invDr, invDi]
>>> acoustics += integral(sur, predefinedacousticwave(dof(p), tf(p), 340, 0, pmlterms))

predefinedadmittancebcincompressible

def predefinedadmittancebcincompressible(
physreg: int,
dofv: expression,
tfv: expression,
v: field,
dofp: expression,
yharm: expression
) -> expression

This function returns the weak formulation of admittance boundary condition for incompressible flows:

p(\omega) = \frac{1.0}{Y(\omega)} (\mathbf{V} \cdot {\mathbf{n})

where pp is the pressure, V\bf{V} the velocity of the flow, n\mathbf{n} the normal vector to the boundary (pointing outward), and Y(ω)Y(\omega) the frequency-dependent admittance in harmonic form. It is assumed that the tabulated values of Y(ω)Y(\omega) are provided as function of frequency by the user.

Example
>>> # Fundamental frequency
>>> f = 100.0
>>> setfundamentalfrequency(f)
>>>
>>> # physical regions
>>> fluid=1, inlet=2, admittanceoutlet=3
>>>
>>> # Pressure field
>>> p = field("h1", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> p.setorder(all, 1)
>>>
>>> # Velocity field
>>> V = field("h1xyz", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> V.setorder(all, 2)
>>>
>>> p.harmonic(1).setconstraint(admittanceoutlet)
>>>
>>> form = formulation()
>>>
>>> omega = 2*getpi()*f
>>>
>>> yharm = makeharmonic([1,2,3,4,5,6,7,8,9,10,11], [yreal(0), yreal(f), yimag(f), yreal(2*f), yimag(2*f), yreal(3*f), yimag(3*f), yreal(4*f), yimag(4*f), yreal(5*f), yimag(5*f)])
>>>
>>> form += integral(admittanceoutlet, 5, predefinedimpedancebcincompressible(fluid, V, dof(V), tf(V), dof(p), yharm))
See Also

predefinedimpedancebcincompressible()

predefinedadvectiondiffusion

def predefinedadvectiondiffusion(
doff: expression,
tff: expression,
v: expression,
alpha: expression,
beta: expression,
gamma: expression,
isdivvzero: bool = True
) -> expression

This defines the weak formulation for the generalized advection-diffusion equation:

βct(αc)+γ(cv)=0\beta \frac{\partial c}{\partial t} - \nabla \cdot \left( \boldsymbol{\alpha} \nabla c\right) + \gamma \nabla \cdot \left( c \boldsymbol{v} \right)= 0

where cc is the scalar quantity of interest and v\boldsymbol{v} is the velocity that the quantity is moving with. With β\beta and γ\gamma set to unit, the classical advection-diffusion equation with diffusivity tensor α\boldsymbol{\alpha} is obtained. Set isdivvzero to True if v\nabla \cdot \boldsymbol{v} is zero (for incompressible flows).

Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> T = field("h1")
>>> v = field("h1xyz")
>>> T.setorder(vol, 1)
>>> v.setorder(vol, 1)
>>> advdiff = formulation()
>>> advdiff += integral(vol, predefinedadvectiondiffusion(dof(T), tf(T), v, 1e-4, 1.0, 1.0, True))
See Also

predefineddiffusion()

predefinedaml

def predefinedaml(c: expression, shifted: bool = False) -> list[expression]

predefinedboxpml

def predefinedboxpml(
pmlreg: int,
innerreg: int,
c: expression,
shifted: bool = False
) -> list[expression]

This is a collective MPI operation and hence must be called by all the ranks. This function returns

[(D)(D)(D)(D)(D1)(D1)] \Big[ \quad \Re (|\boldsymbol{D}|) \quad \Im (|\boldsymbol{D}|) \quad \Re (\boldsymbol{D}) \quad \Im (\boldsymbol{D}) \quad \Re (\boldsymbol{D}^{-1}) \quad \Im (\boldsymbol{D}^{-1}) \quad \Big]

where D\boldsymbol{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
>>> ...
>>> k = 2*getpi()*freq/c # wave number
>>> Dterms = predefinedboxpml(pmlreg, innerreg, k)

predefineddiffusion

def predefineddiffusion(
doff: expression,
tff: expression,
alpha: expression,
beta: expression
) -> expression

This defines the weak formulation for the generalized diffusion equation:

βct(αc)=0\beta \frac{\partial c}{\partial t} - \nabla \cdot \left( \boldsymbol{\alpha} \nabla c\right) = 0

where cc is the scalar quantity of interest. With β\beta set to unit, the classical diffusion equation with diffusivity tensor α\boldsymbol{\alpha} is obtained.

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; top=3
>>>
>>> # Temperature field in [K]:
>>> T = field("h1")
>>> T.setorder(vol, 1)
>>> T.setconstraint(top, 298)
>>>
>>> # Material properties:
>>> k = 237.0 # thermal conductivity of aluminium [W/mK]
>>> cp = 897.0 # heat capacity of aluminium [J/kgK]
>>> rho = 2700 # density of aluminium [Kg/m^3]
>>>
>>> heatequation = formulation()
>>> heatequation += integral(vol, predefineddiffusion(dof(T), tf(T), k, rho*cp))
See Also

predefinedadvectiondiffusion()

predefinedelasticity

def predefinedelasticity(
dofu: expression,
tfu: expression,
Eyoung: expression,
nupoisson: expression,
myoption: str = ''
) -> expression
def predefinedelasticity(
dofu: expression,
tfu: expression,
elasticitymatrix: expression,
myoption: str = ''
) -> expression
def predefinedelasticity(
dofu: expression,
tfu: expression,
u: field,
Eyoung: expression,
nupoisson: expression,
prestress: expression,
myoption: str = ''
) -> expression
def predefinedelasticity(
dofu: expression,
tfu: expression,
u: field,
elasticitymatrix: expression,
prestress: expression,
myoption: str = ''
) -> expression

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 ⁣: ⁣ϵϵ=12[u+(u)T]\begin{cases} \nabla \cdot \boldsymbol{\sigma} - \rho \boldsymbol{\ddot{u}} + \boldsymbol{F} = 0 \\[5pt] \boldsymbol{\sigma} = \boldsymbol{H} \! : \! \boldsymbol{\epsilon} \\[5pt] \boldsymbol{\epsilon} = \frac{1}{2} \Bigl[\nabla \boldsymbol{u} + \bigl(\nabla \boldsymbol{u}\bigr)^T \Bigr] \end{cases}

where,

  • u\boldsymbol{u} is the displacement field vector
  • σ\boldsymbol{\sigma} is the Cauchy stress tensor in Voigt notation
  • ϵ\boldsymbol{\epsilon} is the strain tensor in Voigt notation
  • H\boldsymbol{H} is the 4th4^{th} order elasticity/stiffness tensor
  • F\boldsymbol{F} is the volumetric body force vector
  • ρ\rho 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.

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>>
>>> u = field("h1xyz")
>>> u.setorder(vol, 2)
>>> u.setconstraint(sur)
>>>
>>> elasticity = formulation()
>>> elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), 150e9, 0.3))
>>>
>>> # elasticity += integral(vol, -2300*dtdt(dof(u)) * tf(u)); # 2300 is mass density
>>>
>>> # Atmospheric pressure load (volumetric force) on top face deformed by field u (might require a nonlinear iteration)
>>> elasticity += integral(top, u, -normal(vol)*1e5 * tf(u))
>>>
>>> elasticity.solve()
>>> u.write(top, "u.vtk", 2)

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.

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>>
>>> u = field("h1xyz")
>>> u.setorder(vol, 2)
>>> u.setconstraint(sur)
>>>
>>> # It is enough to only provide the lower triangular part of the elasticity matrix as it is symmetric
>>> H = expression(6,6, [195e9, 36e9,195e9, 64e9,64e9,166e9, 0,0,0,80e9, 0,0,0,0,80e9, 0,0,0,0,0,51e9])
>>>
>>> elasticity = formulation()
>>> elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), H))
>>>
>>> # elasticity += integral(vol, -2300*dtdt(dof(u)) * tf(u)); # 2300 is mass density
>>>
>>> # Atmospheric pressure load (volumetric force) on top face deformed by field u (might require a nonlinear iteration)
>>> elasticity += integral(top, u, -normal(vol)*1e5 * tf(u))
>>>
>>> elasticity.solve()
>>> u.write(top, "u.vtk", 2)

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)(\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy}). Set the prestress expression to 0.00.0 for no prestress.

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>>
>>> u = field("h1xyz")
>>> u.setorder(vol, 2)
>>>
>>> ... # boundary conditions
>>>
>>> prestress = expression(6,1, [10e6,0,0,0,0,0])
>>>
>>> elasticity = formulation()
>>> elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), u, 150e9, 0.3, prestress))
>>> elasticity += integral(vol, -2300*dtdt(dof(u)) * tf(u)); # 2300 is mass density
>>>
>>> ... # nonlinear iteration loop

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)(\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy}).

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3
>>>
>>> u = field("h1xyz")
>>> u.setorder(vol, 2)
>>>
>>> ... # boundary conditions
>>>
>>> # It is enough to only provide the lower triangular part of the elasticity matrix as it is symmetric
>>> H = expression(6,6, [195e9, 36e9,195e9, 64e9,64e9,166e9, 0,0,0,80e9, 0,0,0,0,80e9, 0,0,0,0,0,51e9])
>>>
>>> prestress = expression(6,1, [10e6,0,0,0,0,0])
>>>
>>> elasticity = formulation()
>>> elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), u, H, prestress))
>>> elasticity += integral(vol, -2300*dtdt(dof(u)) * tf(u)); # 2300 is mass density
>>>
>>> ... # nonlinear iteration loop

predefinedelasticradiation

def predefinedelasticradiation(
dofu: expression,
tfu: expression,
rho: expression,
H: expression
) -> expression

predefinedelectrostaticforce

def predefinedelectrostaticforce(
input: expression,
E: expression,
epsilon: expression
) -> expression

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\boldsymbol{T} [N/m2N/m^2] the electrostatic Maxwell stress tensor:

T=ϵ EE12ϵ(EE) I\boldsymbol{T} = \epsilon \ \boldsymbol{E} \otimes \boldsymbol{E} - \frac{1}{2} \epsilon \left( \boldsymbol{E} \cdot \boldsymbol{E} \right) \ \boldsymbol{I}

where ϵ\epsilon is the electric permittivity, E\boldsymbol{E} is the electric field and I\boldsymbol{I} is the identity matrix. The electrostatic force density is T[\nabla \cdot \boldsymbol{T} [N/m^3]] so that the loading for a mechanical problem can be obtained by adding the following term:

Ω(T)udΩ\int_{\Omega} \left( \nabla \cdot \boldsymbol{T} \right) \cdot \boldsymbol{u}^{\prime} d\Omega

where u\boldsymbol{u} is the mechanical displacement. The term can be rewritten in the form that is provided by this function:

ΩT ϵdΩ-\int_{\Omega} \boldsymbol{T} \ \boldsymbol{\epsilon}^{\prime} d\Omega

where ϵ\boldsymbol{\epsilon} 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’).

Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; top=3
>>> v=field("h1"); u=field("h1xyz")
>>> v.setorder(vol,1)
>>> u.setorder(vol,2)
>>>
>>> elasticity = formulation()
>>> elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), 150e9, 0.3))
>>> elasticity += integral(vol, predefinedelectrostaticforce(tf(u,top), -grad(v), 8.854e-12))
See Also

predefinedmagnetostaticforce()

predefinedemwave

def predefinedemwave(
dofE: expression,
tfE: expression,
mur: expression,
mui: expression,
epsr: expression,
epsi: expression,
sigr: expression,
sigi: expression,
precondtype: str = ''
) -> tuple[expression, preconditioner]
def predefinedemwave(
dofE: expression,
tfE: expression,
mur: expression,
mui: expression,
epsr: expression,
epsi: expression,
sigr: expression,
sigi: expression,
pmlterms: list[expression],
precondtype: str = ''
) -> tuple[expression, preconditioner]

This defines the equation for (linear) electromagnetic wave propagation:

×(1μ×E)+σEt+ϵ2Et2=0 \nabla \times \left( \frac{1}{\mu} \nabla \times \boldsymbol{E} \right) + \sigma \frac{\partial \boldsymbol{E}}{\partial t} + \epsilon \frac{\partial^2 \boldsymbol{E}}{\partial t^2} = 0

where E\boldsymbol{E} is the electric field, μ\mu is the magnetic permeability, ϵ\epsilon is the electric permiitivity and σ\sigma 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 μ\mu.
  • epsr and epsi is the real and imaginary part of the electric permittivity ϵ\epsilon.
  • sigr and sigi is the real and imaginary part of the electric conductivity σ\sigma.
  • pmlterms is the list of pml terms.
  • precondtype is the type of precondition.
Examples

Example 1: predefinedemwave(dofE:expression, tfE:expression, mur:expression, mui:expression, epsr:expression, epsi:expression, sigr:expression, sigi:expression, precondtype:str="")

>>> ...
>>> E = field("hcurl", [2,3])
>>> ...
>>> maxwell += integral(sur, predefinedemwave(dof(E), tf(E), mu0,0, epsilon0,0, 0,0))

Example 2: predefinedemwave(dofE:expression, tfE:expression, mur:expression, mui:expression, epsr:expression, epsi:expression, sigr:expression, sigi:expression, pmlterms:List[expression], precondtype:str="")

This is the same as the previous example but with PML boundary conditions.

>>> ...
>>> pmlterms = [detDr, detDi, Dr, Di, invDr, invDi]
>>> maxwell += integral(sur, predefinedemwave(dof(E), tf(E), mu0,0, epsilon0,0, 0,0, pmlterms))

predefinedfluidstructureinteraction

def predefinedfluidstructureinteraction(
fsireg: int,
fluidreg: int,
dofv: expression,
tfv: expression,
v: field,
dofp: expression,
dofu: expression,
tfu: expression
) -> expression
def predefinedfluidstructureinteraction(
fsireg: int,
dofv: expression,
tfv: expression,
v: field,
dofu: expression,
tfu: expression
) -> expression

This function returns the formulation of fluid-structure interaction for incompressible flows:

The Navier-Stokes equation for incompressible flows in weak form is given by

ΩρtpdΩ+Ω(ρv)pdΩ=0 \int_{\Omega} \frac{\partial \rho}{\partial t} \cdot p^{\prime} d\Omega + \int_{\Omega} (\rho \nabla \cdot \mathbf{v}) {p^{\prime}} d\Omega = 0 ΩρvtvdΩ+Ωρ(vv)vdΩ=ΩρpvdΩΩμ(v+vT)vdΩ+Γμ(v+vT)nvdΓ\begin{align*} \int_{\Omega} \rho \frac{\partial \mathbf{v}}{\partial t} \cdot \mathbf{v^{\prime}} d\Omega + \int_{\Omega} \rho (\mathbf{v} \cdot \nabla \mathbf{v}) \cdot \mathbf{v^{\prime}} d\Omega = & - \int_{\Omega} \rho \nabla p \cdot \mathbf{v^{\prime}} d\Omega \\ & - \int_{\Omega} \mu (\nabla \mathbf{v} + \nabla \mathbf{v}^T)\cdot \nabla \mathbf{v^{\prime}} d\Omega \\ & + \int_{\Gamma} \mu (\nabla \mathbf{v} + \nabla \mathbf{v}^T)\cdot \mathbf{n} \cdot \mathbf{v^{\prime}} d\Gamma \\ \end{align*}

where v\mathbf{v} is the velocity, pp the pressure, u\mathbf{u} the displacement of the structure, Ω\Omega the fluid domain, and Γ\Gamma 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 (λ\boldsymbol{\lambda}) method as

Γfsi((vu˙)λ+λv)dΓfsi=0 \int_{\Gamma_{\text{fsi}}} \Bigl( \left(\mathbf{v} - \dot{\mathbf{u}} \right) \cdot \boldsymbol{\lambda^{\prime}} + \boldsymbol{\lambda} \cdot \mathbf{v}^{\prime} \Bigr) d\Gamma_{\text{fsi}} = 0

The Lagrange multiplier λ\boldsymbol{\lambda} 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λ)udΓfsi=0 \int_{\Gamma_{\text{fsi}}} (p \mathbf{n} -\boldsymbol{\lambda})\cdot \mathbf{u}^{\prime} d\Gamma_{\text{fsi}} = 0

where Γfsi\Gamma_{\text{fsi}} the boundary interface between fluid and structure domains.

Example
>>> mymesh = mesh("micropillar.msh")
>>> solid = 1 # physical region
>>> fluid = 2 # physical region
>>> fsinterface = 3 # Fluid-structure interface region
>>>
>>> v=field("h1xyz"); p=field("h1"); u=field("h1xyz")
>>> v.setorder(fluid, 2)
>>> p.setorder(fluid, 1) # Satisfies the LBB condition
>>>
>>> fsi = formulation()
>>> fsi += integral(solid, predefinedelasticity(dof(u), tf(u), H))
>>> fsi += integral(fluid, umesh, predefinednavierstokes(dof(v), tf(v), v, dof(p), tf(p), 8.9e-4, 1000, 0, 0))
>>> fsi += integral(fsinterface, umesh, predefinedfluidstructureinteraction(fsinterface, fluid, dof(v), tf(v), v, dof(p), dof(u), tf(u)))
See Also

predefinednavierstokes()

predefinedimpedancebcincompressible

def predefinedimpedancebcincompressible(
physreg: int,
dofv: expression,
tfv: expression,
v: field,
dofp: expression,
zharm: expression
) -> expression

This function returns the weak formulation of impedance boundary condition for incompressible flows:

p(ω)=Z(ω)(Vn) p(\omega) = Z(\omega) (\mathbf{V} \cdot \mathbf{n})

where pp is the pressure, V\mathbf{V} the velocity of the flow, n\mathbf{n} the normal vector to the boundary (pointing outward), and Z(ω)Z(\omega) the frequency-dependent impedance in harmonic form. It is assumed that the tabulated values of Z(ω)Z(\omega) are provided as function of frequency by the user.

Example
>>> # Fundamental frequency
>>> f = 100.0
>>> setfundamentalfrequency(f)
>>>
>>> # physical regions
>>> fluid=1, inlet=2, impedanceoutlet=3
>>>
>>> # Pressure field
>>> p = field("h1", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> p.setorder(all, 1)
>>>
>>> # Velocity field
>>> V = field("h1xyz", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> V.setorder(all, 2)
>>>
>>> p.harmonic(1).setconstraint(impedanceoutlet)
>>>
>>> form = formulation()
>>>
>>> omega = 2*getpi()*f
>>>
>>> zharm = makeharmonic([1,2,3,4,5,6,7,8,9,10,11], [zreal(0), zreal(f), zimag(f), zreal(2*f), zimag(2*f), zreal(3*f), zimag(3*f), zreal(4*f), zimag(4*f), zreal(5*f), zimag(5*f)])
>>>
>>> form += integral(impedanceoutlet, 5, predefinedimpedancebcincompressible(fluid, V, dof(V), tf(V), dof(p), zharm))
See Also

predefinedadmittancebcincompressible()

predefinedlinearpoissonwalldistance

def predefinedlinearpoissonwalldistance(
physreg: int,
wallreg: int,
interpolationorder: int = 1,
relresddmtol: float = 1e-12,
maxnumddmit: int = 500,
relerrnltol: float = 1e-06,
maxnumnlit: int = 100,
project: bool = True,
verbosity: int = 1
) -> expression
def predefinedlinearpoissonwalldistance(
physreg: int,
meshdeform: expression,
wallreg: int,
interpolationorder: int = 1,
relresddmtol: float = 1e-12,
maxnumddmit: int = 500,
relerrnltol: float = 1e-06,
maxnumnlit: int = 100,
project: bool = True,
verbosity: int = 1
) -> expression

This function calculates and returns distance values from the wall based on the linear Poisson wall distance equation:

(ϕ)+1=0on  Ωϕ=0on  Γudϕdn=0on  ΓN\begin{align*} \nabla \cdot (\nabla \phi) + 1 & = 0 \qquad &on \ \ &\Omega \\[5pt] \phi & = 0 \qquad &on \ \ &\Gamma_u \\[5pt] \frac{d \phi}{dn} & = 0 \qquad &on \ \ &\Gamma_{N} \end{align*}