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 $arccos$ or $cos^{-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

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

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

`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 $S_{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 ( $S_{11}, S_{12}, \ldots$ ), respectively.

For the definition of S-parameters, we assume that the input signal at port $i$ is a complex number $a_i$ times the reference mode of port $i$ . If port $i$ is driven, $a_i$ is the driving signal used for feeding; otherwise $a_i=0$ . The output signal at port $i$ is a complex number $b_i$ times the reference mode of port $i$ . The S-parameter $S_{ij}$ is defined such that the equation $S_{ij}=b_i/a_j$ holds when port $j$ 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)

`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 $S$ 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 $E$ and $H$ and the propagation constant $\gamma=\alpha+j\beta$ as a list of expressions. The returned modes are scaled such that the power into the waveguide,

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

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

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

`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

`allmeasuredistance`

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

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

\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

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

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

`array1x1`

def array1x1(arg0: expression) -> expression

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

\begin{bmatrix} arg0 \end{bmatrix}

`array1x2`

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

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

\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\times3$ . The array is populated in a row-major way.

\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\times1$ . The array is populated in a row-major way.

\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\times2$ . The array is populated in a row-major way.

\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\times3$ . The array is populated in a row-major way.

\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\times1$ . The array is populated in a row-major way.

\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\times2$ . The array is populated in a row-major way.

\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\times3$ . The array is populated in a row-major way.

\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 $arcsin$ or $sin^{-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

`atan`

def atan(input: expression) -> expression

This returns an expression that is the $arctan$ or $tan^{-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

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

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

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

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

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

`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 $u_1 = u_2$ between boundary region $\Gamma_1$ and $\Gamma_2$ (with $\Gamma_1 \subseteq \Gamma_2$ , meshes can be non-matching). In case $\Gamma_2$ is larger than $\Gamma_1$ ( $\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 $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 $\Gamma_1$ is the rotor-stator interface on the (non-moving) stator side while the boundary $\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 $x$ axis.

The condition is based on a Lagrange multiplier of the same type and the same harmonic content as the field $u_1$ and $u_2$ . The mortar finite element method is used to link the unknown $dof$ field on $\Gamma_1$ and $\Gamma_2$ so that there is no restriction on the mesh used for both regions.

`cos`

def cos(input: expression) -> expression

This returns an expression that is the $cos$ 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

`crossproduct`

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

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

\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.

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

`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 $dB$ units to $Nepers$ .

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

`detjac`

def detjac() -> expression

This returns the determinant of the Jacobian matrix.

Example

>>> detJ = detjac()

`div`

def div(input: expression) -> expression

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

s = \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)

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

`doubledotproduct`

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

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

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

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

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

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

`dx`

def dx(input: expression) -> expression

This returns the $x$ space derivative expression.

Example

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

`dy`

def dy(input: expression) -> expression

This returns the $y$ space derivative expression.

Example

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

`dz`

def dz(input: expression) -> expression

This returns the $z$ space derivative expression.

Example

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

`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.

C_{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 $e$ :

e^{input}

Example

>>> expr = exp(2)
>>> expr.print()
Expression size is 1x1
 @ row 0, col 0 :
7.38906

`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 $\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 $\boldsymbol{n}$ and tangents $\boldsymbol{t_1}$ and $\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()

`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 $\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()

`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

`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 $x$ and $z$ components are returned and the $y$ component includes a $2 \pi$ 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)

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

`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 $x$ 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)

`gety`

def gety() -> field

This returns the $y$ 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

`getz`

def getz() -> field

This returns the $z$ 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)

`grad`

def grad(input: expression) -> expression

For a scalar input expression, this is mathematically treated as the gradient of a scalar ( $\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 ( $\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)

`greenlagrangestrain`

def greenlagrangestrain(input: expression) -> expression

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

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

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

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

`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 \times 3$ Jacobian matrix.

>>> Ji = invjac()

`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

`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

`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

`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

`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 \times 3$ Jacobian matrix.

>>> J = jac()

`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:

\frac{b - a}{num -1}

Example

>>> vals = linspace(0.5, 2.5, 5)
>>> printvector(vals)
0.5 1 1.5 2 2.5

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

`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

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

`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 < $0$ 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)

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

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

\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 $u_1^{\prime}$ , $u_2^{\prime}$ and ${\lambda}^{\prime}$ . The example below illustrates the formulation terms needed to implement the Lagrange multiplier between interfaces $\Gamma_1$ and $\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 $\Gamma_1$ without a relative in $\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 $\Gamma_1$ is without a relative in $\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 $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)).

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

`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 $u$ a rotation or translation periodic condition between boundary region $\Gamma_1$ and $\Gamma_2$ (meshes can be non-conforming). A factor different than $1$ can be provided to scale the field on $\Gamma_2$ (use - $1$ 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)

More advanced periodic conditions can be implemented easily using on() function.

`pow`

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

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

{base}^{exponent}

Example

>>> expr = pow(2, 5)
>>> expr.print()
Expression size is 1x1
 @ row 0, col 0 :
32

`predefinedacousticradiation`

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

This function defines the equation for the Sommerfeld acoustic radiation condition

\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(\omega t - kr)

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

\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/m$ ) in case of harmonic problems. For convenience use the function dbtoneper() to convert $dB/m$ attenuation values to $Np/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 $p$ is the acoustic pressure and field $\boldsymbol{u}$ is the mechanical displacement. Calling $\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

\boldsymbol{f}_{pressure} = -p \boldsymbol{n}

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

\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 $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

>>> ...
>>> 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:

\nabla^2 p - \frac{1}{c²} \frac{\partial^2 p}{\partial t^2} = 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.

>>> 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 $p$ is the pressure, $\bf{V}$ the velocity of the flow, $\mathbf{n}$ the normal vector to the boundary (pointing outward), and $Y(\omega)$ the frequency-dependent admittance in harmonic form. It is assumed that the tabulated values of $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))

`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:

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

where $c$ is the scalar quantity of interest and $\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 $\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))

`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

\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 $\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:

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

where $c$ 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))

`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:

\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,

$\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
$\boldsymbol{H}$ is the $4^{th}$ order elasticity/stiffness tensor
$\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 $(\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy})$ . Set the prestress expression to $0.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 $(\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 $\boldsymbol{T}$ [ $N/m^2$ ] the electrostatic Maxwell stress tensor:

\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, $\boldsymbol{E}$ is the electric field and $\boldsymbol{I}$ is the identity matrix. The electrostatic force density is $\nabla \cdot \boldsymbol{T} [$ N/m^3 $]$ so that the loading for a mechanical problem can be obtained by adding the following term:

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

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

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

`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:

\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 $\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

\int_{\Omega} \frac{\partial \rho}{\partial t} \cdot p^{\prime} d\Omega + \int_{\Omega} (\rho \nabla \cdot \mathbf{v}) {p^{\prime}} d\Omega = 0

\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 $\mathbf{v}$ is the velocity, $p$ the pressure, $\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

\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

\int_{\Gamma_{\text{fsi}}} (p \mathbf{n} -\boldsymbol{\lambda})\cdot \mathbf{u}^{\prime} d\Gamma_{\text{fsi}} = 0

where $\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)))

`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(\omega) = Z(\omega) (\mathbf{V} \cdot \mathbf{n})

where $p$ is the pressure, $\mathbf{V}$ the velocity of the flow, $\mathbf{n}$ the normal vector to the boundary (pointing outward), and $Z(\omega)$ the frequency-dependent impedance in harmonic form. It is assumed that the tabulated values of $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))

`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:

\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*}

Script API

Functions

abs

Example

acos

Example

See Also

adapt

Example

See Also

alladapt

Example

See Also

allcomputesparameters

Example

See Also

alleigenport

Example

See Also

allintegrate

Example

See Also

allinterpolate

Example

See Also

alllaplace

allmax

Example

See Also

allmeasuredistance

Example

allmin

Example

See Also

allpartition

Example

allprobe

Example

allsolve

Example

See Also

andpositive

Example

See Also

array1x1

array1x2

Example

array1x3

Example

array2x1

Example

array2x2

Example

array2x3

Example

array3x1

Example

array3x2

Example

array3x3

Example

asin

Example

See Also

atan

Example

See Also

bode

Example

cn

Example

See Also

comp

Example

See Also

complexdivision

Example

complexinverse

Example

complexproduct

`abs`

`acos`

`adapt`

`alladapt`

`allcomputesparameters`

`alleigenport`

`allintegrate`

`allinterpolate`

`alllaplace`

`allmax`

`allmeasuredistance`

`allmin`

`allpartition`

`allprobe`

`allsolve`

`andpositive`

`array1x1`

`array1x2`

`array1x3`

`array2x1`

`array2x2`

`array2x3`

`array3x1`

`array3x2`

`array3x3`

`asin`

`atan`

`bode`

`cn`

`comp`

`complexdivision`

`complexinverse`

`complexproduct`

`computeradiationpattern`

`compx`

`compy`

`compz`

`continuitycondition`

`cos`

`crossproduct`

`curl`

`d`

`dbtoneper`

`determinant`

`detjac`

`div`

`dof`

`doubledotproduct`

`dt`

`dtdt`

`dtdtdt`

`dtdtdtdt`

`dx`

`dy`

`dz`

`elasticwavespeed`

`elementwiseproduct`

`emwavespeed`

`entry`

`exp`