# # Module quanscient

This is the python binding for the 'PYSIMCORE' c++ library.

## # Functions

### #abs

def abs(
input: expression
) ‑> quanscient.expression

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

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

### #acos

def acos(
input: expression
) ‑> quanscient.expression

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

###### # Example
>>> expr = acos(0)
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

verbosity: int = 0
) ‑> bool

verbosity: int = 0
) ‑> bool

### #allcomputecohomology

def allcomputecohomology(
meshfile: str,
physregstocut: List[int],
subdomains: List[int] = []
) ‑> List[int]

### #allfinalize

def allfinalize() ‑> None

def allgather(
*args,
**kwargs
) ‑> None

### #allinitialize

def allinitialize(
verbosity: int = 1
) ‑> None

### #allintegrate

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

### #allmeasuredistance

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

### #allpartition

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

### #allsolve

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

### #andpositive

def andpositive(
exprs: List[expression]
) ‑> quanscient.expression

This returns an expression whose value is 1 for all evaluation points where the value of all input expression 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
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

###### # Example
>>> myarray = array1x1(1)

### #array1x2

def array1x2(
arg0: expression,
arg1: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

###### # Example
>>> myarray = array1x2(1,2)

### #array1x3

def array1x3(
arg0: expression,
arg1: expression,
arg2: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

arg2 : expression : term13

###### # Example
>>> myarray = array1x3(1,2,3)

### #array2x1

def array2x1(
arg0: expression,
arg1: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term21

###### # Example
>>> myarray = array2x1(1, 2)

### #array2x2

def array2x2(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

arg2 : expression : term21

arg3 : expression : term22

###### # Example
>>> myarray = array2x2(1,2, 3,4)

### #array2x3

def array2x3(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression,
arg4: expression,
arg5: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

arg2 : expression : term13

arg3 : expression : term21

arg4 : expression : term22

arg5 : expression : term23

###### # Example
>>> myarray = array2x3(1,2,3, 4,5,6)

### #array3x1

def array3x1(
arg0: expression,
arg1: expression,
arg2: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term21

arg2 : expression : term31

###### # Example
>>> myarray = array3x1(1, 2, 3)

### #array3x2

def array3x2(
arg0: expression,
arg1: expression,
arg2: expression,
arg3: expression,
arg4: expression,
arg5: expression
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

arg2 : expression : term21

arg3 : expression : term22

arg4 : expression : term31

arg5 : expression : term32

###### # 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
) ‑> quanscient.expression

This defines a vector or matrix operation of size . The array is populated in a row major way.

###### # Parameters

arg0 : expression : term11

arg1 : expression : term12

arg2 : expression : term13

arg3 : expression : term21

arg4 : expression : term22

arg5 : expression : term23

arg6 : expression : term31

arg7 : expression : term32

arg8 : expression : term33

###### # Example
>>> myarray = array3x3(1,2,3, 4,5,6, 7,8,9)

### #asin

def asin(
input: expression
) ‑> quanscient.expression

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

###### # Example
>>> expr = asin(sqrt(0.5))
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
) ‑> quanscient.expression

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

###### # Example
>>> expr = atan(0.57735)
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

### #barrier

def barrier() ‑> None

This is a collective MPI operation (must be called by all ranks). Processing will be waiting until all ranks have reached this call.

###### # Example
>>> import quanscient as qs
>>> qs.barrier()
###### # Notes

If this function is not called by all the ranks, then the other ranks that calls barrier() will wait indefinitely.

*args,
**kwargs
) ‑> None

### #cn

def cn(
n: float
) ‑> quanscient.expression

### #comp

def comp(
selectedcomp: int,
input: expression
) ‑> quanscient.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 is say 5 then an expression containing the entries of fifth row of 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)

### #compx

def compx(
input: expression
) ‑> quanscient.expression

This returns the first or x component of a column vector expression. For a matrix expression an expression containing the entries of first row of 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
) ‑> quanscient.expression

This returns the second or y component of a column vector expression. For a matrix expression an expression containing the entries of second row of 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
) ‑> quanscient.expression

This returns the third or z component of a column vector expression. For a matrix expression an expression containing the entries of third row of 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(
*args,
**kwargs
) ‑> List[quanscient.integration]

### #cos

def cos(
input: expression
) ‑> quanscient.expression

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

### #count

def count() ‑> int

This returns the number of processes in the universe.

###### # Example
>>> import quanscient as qs
>>> numranks = qs.count()

getrank()

### #crossproduct

def crossproduct(
a: expression,
b: expression
) ‑> quanscient.expression

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

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; top=3
>>> E = field("hcurl")
>>> n = normal(vol)
>>> cp = crossproduct(E, n)

### #curl

def curl(
input: expression
) ‑> quanscient.expression

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

###### # 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]
) ‑> quanscient.expression

### #dbtoneper

def dbtoneper(
toconvert: expression
) ‑> quanscient.expression

This converts a value to .

###### # Parameters

toconvert : expression : expression value in units.

###### # Example
>>> neperattenuation = dbtoneper(100.0)

### #determinant

def determinant(
input: expression
) ‑> quanscient.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

inverse()

### #detjac

def detjac() ‑> quanscient.expression

### #div

def div(
input: expression
) ‑> quanscient.expression

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

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1);
>>> div(u).write(vol, "divu.vtk", 1)

### #dof

def dof(
*args,
**kwargs
) ‑> quanscient.expression

### #doubledotproduct

def doubledotproduct(
a: expression,
b: expression
) ‑> quanscient.expression

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

###### # Example
>>> a = array2x2(1,2, 3,4)
>>> b = array2x2(11,12, 13, 14)
Expression size is 1x1
@ row 0, col 0 :
130

### #dt

def dt(
*args,
**kwargs
) ‑> quanscient.expression

### #dtdt

def dtdt(
*args,
**kwargs
) ‑> quanscient.expression

### #dtdtdt

def dtdtdt(
input: expression
) ‑> quanscient.expression

### #dtdtdtdt

def dtdtdtdt(
input: expression
) ‑> quanscient.expression

### #dx

def dx(
input: expression
) ‑> quanscient.expression

### #dy

def dy(
input: expression
) ‑> quanscient.expression

### #dz

def dz(
input: expression
) ‑> quanscient.expression

### #elementwiseproduct

def elementwiseproduct(
a: expression,
b: expression
) ‑> quanscient.expression

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

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

### #entry

def entry(
row: int,
col: int,
input: expression
) ‑> quanscient.expression

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

###### # Parameters

row : int : Row from which the entry is requested.

col : int : Column from which the entry is requested.

input : expression : Vector or matrix input 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

### #evaluate

def evaluate(
toevaluate: expression
) ‑> List[float]

def exchange(
*args,
**kwargs
) ‑> None

### #exp

def exp(
input: expression
) ‑> quanscient.expression

This returns an exponential function of base .

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

### #eye

def eye(
size: int
) ‑> quanscient.expression

This returns a size x size identity matrix

###### # Parameters

size : int : Size of the 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
) ‑> quanscient.expression

This returns an expression whose value is the interpolation order on each element for the provided field. The value is a constant on each element. When argument is set, the value returned is the lowest order required to include of the total shape function coefficient weight. An additional optional argument can be set to provide a minimum total weight below which the lowest possible field order is returned.

###### # Parameters

input : field : The field object whose interpolation order is requested.

alpha : double :

absthres : double :

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

### #finalize

def finalize() ‑> None

def gather(
*args,
**kwargs
) ‑> None

### #getcirculationsources

def getcirculationsources(
circulations: List[expression],
inittimederivatives: List[List[float]] = []
) ‑> List[quanscient.expression]

### #getepsilon0

def getepsilon0() ‑> float

Returns the value of vacuum permittivity .

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

### #getextrusiondata

def getextrusiondata() ‑> List[quanscient.expression]

This gives the relative depth in the extruded layer, the extrusion normal and tangents and . 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.write("diskpml.msh")
>>> pmldata = getextrusiondata()

mesh.extrude()

### #getharmonic

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

Returns the maximum number of threads allowed.

2

### #getmu0

def getmu0() ‑> float

Returns the value of vacuum permeability .

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

### #getpi

def getpi() ‑> float

Returns value of .

###### # Example
>>> pi = getpi()
3.141592653589793

### #getrandom

def getrandom() ‑> float

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

count()

### #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(
*args,
**kwargs
) ‑> List[float]

### #getversion

def getversion() ‑> int

### #getversionname

def getversionname() ‑> str

### #getx

def getx() ‑> quanscient.field

Returns the 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() ‑> quanscient.field

Returns the 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() ‑> quanscient.field

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

input: expression
) ‑> quanscient.expression

For a scalar input expression this is mathematically treated as the gradient of a scalar () and the output is a column vector with one entry per space derivative. For a vector input expression, this is mathematically treated as gradient of a vector () 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);

### #greenlagrangestrain

def greenlagrangestrain(
input: expression
) ‑> quanscient.expression

This defines the (nonlinear) Green-Lagrange strains in Voigt form . The input can either be the displacement field or its gradient.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> glstrain = greenlagrangestrain(u)
>>> glstrain.print()

### #grouptimesteps

def grouptimesteps(
*args,
**kwargs
) ‑> None

### #harm

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

### #ifpositive

def ifpositive(
condexpr: expression,
trueexpr: expression,
falseexpr: expression
) ‑> quanscient.expression

This returns a conditional expression. The expression value is trueexpr for all evaluation points where condexpr is larger or equal to zero. Otherwise, its value is falseexpr.

###### # Parameters

condexpr : expression : Expression the specifies the conditional argument.

trueexpr : expression : Expression value at points where condexpr evaluates to True.

falseexpr : expression : Expression value at points where condexpr evaluates to False.

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

### #initialize

def initialize() ‑> None

### #initializelogs

def initializelogs(
format: logformat
) ‑> None

INTERNAL -- do not call in script

Initializes the structured logging system

###### # Parameters

format : logformat : Format for the logs. One of logformat.plain or logformat.json

None

def insertrank(
name: str
) ‑> str

### #integral

def integral(
*args,
**kwargs
) ‑> quanscient.integration

### #inverse

def inverse(
input: expression
) ‑> quanscient.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)

determinant()

### #invjac

def invjac(
*args,
**kwargs
) ‑> quanscient.expression

### #isavailable

def isavailable() ‑> bool

### #isdefined

def isdefined(
physreg: int
) ‑> bool

This checks if a physical region is defined.

###### # Parameters

physreg : int : Physical region identifier.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; cirlce=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 is empty.

###### # Parameters

physreg : int : Physical region identifier.

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

### #isinside

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

This checks if a physical region is fully included in another region.

###### # Parameters

physregtocheck : int : Physical region to check.

physreg : int : Physical region with which physregtocheck is checked.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3; cirlce=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.

###### # Parameters

physregtocheck : int : Physical region to check.

physreg : int : Physical region with which physregtocheck is checked.

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

### #jac

def jac(
*args,
**kwargs
) ‑> quanscient.expression

### #linspace

def linspace(
a: float,
b: float,
num: int
) ‑> List[float]

meshfile: str
) ‑> List[List[quanscient.shape]]

filename: str,
delimiter: str = ',',
sizeincluded: bool = False
) ‑> List[float]

This loads from a file a list/vector. 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 list/vector length integer.

###### # Parameters

filename : str : The name of the file containing the entries of a list/vector.

delimiter : str, default=',' : Character separating each entry of the list/vector when writing to a file. E.g ',' or '\n'.

sizeincluded : bool, default=False : Must be set to True if the first number in the file is list/vector length integer.

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

writevector()

### #log

def log(
input: expression
) ‑> quanscient.expression

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

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

### #logspace

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

### #makeharmonic

def makeharmonic(
harms: List[int],
exprs: List[expression]
) ‑> quanscient.expression

### #max

def max(
*args,
**kwargs
) ‑> quanscient.expression

### #meshsize

def meshsize(
integrationorder: int
) ‑> quanscient.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.

###### # Parameters

integrationorder : int : Determines the accuracy of mesh size calculated. Higher the number, better the accuracy. Integration order can be negative.

###### # Raises

RuntimeError : If integration order < .

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

### #min

def min(
*args,
**kwargs
) ‑> quanscient.expression

### #mod

def mod(
input: expression,
modval: float
) ‑> quanscient.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

### #moveharmonic

def moveharmonic(
origharms: List[int],
destharms: List[int],
input: expression,
numfftharms: int = -1
) ‑> quanscient.expression

### #norm

def norm(
arg0: expression
) ‑> quanscient.expression

This gives the 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(
*args,
**kwargs
) ‑> quanscient.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(
*args,
**kwargs
) ‑> quanscient.expression

### #orpositive

def orpositive(
exprs: List[expression]
) ‑> quanscient.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
) ‑> List[quanscient.integration]

This returns

###### # Example
>>> ...
>>> u = field("h1xyz")
>>> ...
>>> elasticity = formulation()

This returns the formulation terms required to enforce on field a rotation or translation periodic condition between between boundary region and (meshes can be non-conforming).

### #pow

def pow(
base: expression,
exponent: expression
) ‑> quanscient.expression

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

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

dofp: expression,
tfp: expression,
soundspeed: expression,
neperattenuation: expression
) ‑> quanscient.expression

This function defines the equation for the Sommerfeld acoustic radiation condition

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

propogating in the direction perpendicular to the truncation boundary indeed satisfies the Sommerfeld radiation condition since

Zero artificial wave reflection at the truncation boundary happens only if it is perpendicular to the outgoing waves. In practical aplications 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 wavelength away).

An acoustic attenuation value can be provided (in ) in case of harmonic problems. For convenience use the function dbtoneper() to convert attenuation values to .

###### # 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
) ‑> quanscient.expression

This function defines the bi-directional coupling for acoustic-structure interaction at the medium interface. Field is the acoustic pressure and field is the mechanical displacement. Calling 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

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

To have a good matrix conditioning a scaling factor (e.g ) can be provided. In this case the pressure source is divided by and, to compensate, the pressure force is multiplied by . 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 ) in case of harmonic problems. For convenience use the function dbtoneper() to convert attenuation values to .

###### # 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(
*args,
**kwargs
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

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

An acoustic attenuation value can be provided (in ) in case of harmonic problems. For convenience use the function dbtoneper() to convert attenuation values to .

###### # Parameters

dofp : expression : dof of acoustic pressure field.

tfp : expression : test function of acoustic pressure field.

soundspeed: speed of sound in . neperattenuation: attenuation in . pmlterms : List[expression] : List of pml terms.

precondtype : str : Type of precondition. Default is an empty string "".

###### # Examples

Example 1: predefinedacousticwave(dofp:expression, tfp:expression, soundspeed:expression, neperattenuation:expression, precondtype:str="")

In the illustrative example below a highly-attenuated acustic 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 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))

doff: expression,
tff: expression,
v: expression,
alpha: expression,
beta: expression,
gamma: expression,
isdivvzero: bool = True
) ‑> quanscient.expression

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

where is the scalar quantity of interest and is the velocity that the quantity is moving with. With and set to unit, the classical advection-diffusion equation with diffusivity tensor is obtained. Set isdivvzero to True if 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)

predefineddiffusion()

### #predefinedaml

def predefinedaml(
k: expression,
shifted: bool = False
) ‑> List[quanscient.expression]

### #predefinedboxpml

def predefinedboxpml(
pmlreg: int,
innerreg: int,
k: expression,
shifted: bool = False
) ‑> List[quanscient.expression]

### #predefineddiffusion

def predefineddiffusion(
doff: expression,
tff: expression,
alpha: expression,
beta: expression
) ‑> quanscient.expression

This defines the weak formulation for the generalized diffusion equation:

where is the scalar quantity of interest. With set to unit, the classical diffusion equation with diffusivity tensor 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(
*args,
**kwargs
) ‑> quanscient.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:

where,

• is the displacement field vector
• is the cauchy stress tensor in Voight notation
• is the strain tensor in Voight notation
• is the order elasticity/stiffness tensor
• is the volumetric body force vector
• is the mass density

This is used when the material is isotropic. u is the mechanical displacement, Eyoung is the Young's modulus [Pa] and nupoisson is the Poisson's ratio. In 2D the option string must be either set to "planestrain" or "planestress" for a plane strain or plane stress assumption respectively.

>>> 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 . Set the prestress expressiom to 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 .

>>> 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, 150e9, 0.3, prestress))
>>> elasticity += integral(vol, -2300*dtdt(dof(u)) * tf(u)); # 2300 is mass density
>>>
>>> ... # nonlinear iteration loop

### #predefinedelasticwave

def predefinedelasticwave(
dofu: expression,
tfu: expression,
rho: expression,
H: expression,
neperattenuation: expression,
pmlterms: List[expression],
precondtype: str
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

### #predefinedelectrostaticforce

def predefinedelectrostaticforce(
input: expression,
E: expression,
epsilon: expression
) ‑> quanscient.expression

This function defines the weak formulation term for electrostatic forces. The first argument is 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 [] the electrostatic Maxwell stress tensor:

where is the electric permittivity, is the electric field and is the identity matrix. The electrostatic force density is N/m^3 so that the loading for a mechanical problem can be obtained by adding the following term:

where is the mechanical displacement. The term can be rewritten in the form that is provided by this function:

where is the infinitesimal strain tensor. This is identical to what is obtained using the virtual work principle. For details refer to 'Domain decomposition techniques for the nonlinear, steady state, finite element simulation of MEMS ultrasonic transducer arrays', page 40.

In this function a region should be provided to the test function argument to compute the force only for the degrees of freedom associated to that specific region (in the example below with tf(u, top) the force only acts on surface region 'top'. In any case a correct force calculation requires to include 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))

predefinedmagnetostaticforce()

### #predefinedelectrostatics

def predefinedelectrostatics(
dofv: expression,
tfv: expression,
epsilon: expression,
precondtype: str = ''
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

### #predefinedemwave

def predefinedemwave(
*args,
**kwargs
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

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

where is the electric field, is the magnetic permeability, is the electric permiitivity and is the electric conductivity. The real and imaginary parts of each material property can be provided.

###### # Parameters

dofE : expression : dof of electric field.

tfE : expression : test function of electric field.

mur : expression : real part of the magnetic permeability .

mui : expression : imaginary part of the magnetic permeability .

epsr : expression : real part of the electric permittivity .

epsi : expression : imaginary part of the electric permittivity .

sigr : expression : real part of the electric conducitivty .

sigi : expression : imaginary part of the electric conducitivty .

pmlterms : List[expression] : List of pml terms.

precondtype : str : Type of precondition. Default is an empty string "".

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

### #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
) ‑> quanscient.expression

2-Poisson wall distance equation ∇.∇ <φ<inf+ 1 = 0, φ ∈ Ω, φ = 0, φ ∈ 𝜕Ω --> DBC dφ/dn = 0, φ ∈ 𝜕Ω_N --> NBC where, φ is the approximate wall distance function

The system is linear and hence a linear solver is used. After calculating the field 'φ', a better approximation of distance field is obtained as follows: d(φ) = -|∇φ| ± sqrt(|∇φ|² + 2φ) --> Eq. (9)

The above formula is specific to Poisson parameter p=2. A more generic p-poisson equation and the corresponding improved approximation is given in http://samofar.eu/wp-content/uploads/2018/10/Bakker_Jelle_BSc-thesis_2018.pdfopen in new window and is implemented in wd::predefinednonlinearpoissonwalldistance

Reference: Computations of Wall Distances Based on Differential Equations Paul G. Tucker, Chris L. Rumsey, Philippe R. Spalart, Robert E. Bartels, and Robert T. Biedron AIAA Journal 2005 43:3, 539-549, https://doi.org/10.2514/1.8626open in new window

Excerpts from the above paper: "The derivation of Eq. (9) assumes extensive (infinite) coordinates in the non-normal wall directions. Hence, unlike the Eikonal, 'd' from Eq. (9) is only accurate close to walls. However, turbulence models only need 'd' accurate close to walls."

"The numerically benign nature the Poisson had comparable accuracy to the Eikonal. Hence this approach is recommended."

### #predefinedmagnetostaticforce

def predefinedmagnetostaticforce(
input: expression,
H: expression,
mu: expression
) ‑> quanscient.expression

This function defines the weak formulation term for magnetostatic forces. The first argument is mechanical displacement test function or its gradient, the second is the magnetic field expression and the third argument is the magnetic permeability ( must be a scalar).

Let us call [] the magnetostatic Maxwell stress tensor:

where is the magnetic permeability, is the magnetic field and is the identity matrix. The magnetostatic force density is N/m^3 so that the loading for a mechanical problem can be obtained by adding the following term:

where is the mechanical displacement. The term can be rewritten in the form that is provided by this function:

where is the infinitesimal strain tensor. This is identical to what is obtained using the virtual work principle. For details refer to 'Domain decomposition techniques for the nonlinear, steady state, finite element simulation of MEMS ultrasonic transducer arrays', page 40.

In this function a region should be provided to the test function argument to compute the force only for the degrees of freedom associated to that specific region (in the example below with tf(u, top) the force only acts on surface region 'top'. In any case a correct force calculation requires to include 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
>>> phi=field("h1"); u=field("h1xyz")
>>> phi.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(phi), 4*getpi()*1e-7))

predefinedelectrostaticforce()

### #predefinedmagnetostatics

def predefinedmagnetostatics(
*args,
**kwargs
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

### #predefinednavierstokes

def predefinednavierstokes(
dofv: expression,
tfv: expression,
v: expression,
dofp: expression,
tfp: expression,
mu: expression,
rho: expression,
dtrho: expression,
includetimederivs: bool = False,
isdensityconstant: bool = True,
isviscosityconstant: bool = True,
precondtype: str = ''
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

This defines the weak formulation for the general (nonlinear) flow of Newtonian fluids:

where,

• is the fluid density
• is the dynamic viscosity of the fluid
• is the pressure
• is the flow velocity

The formulation is provided in a form leading to a quadratic (Newton) convergence when solved iteratively in a loop. This formulation is only valid to simulate laminar as well as turbulent flows. Using it to simulate turbulent flows leads to a so- called DNS method (direct numerical simulation). DNS does not requires any turbulence model since it takes into account the whole range of spatial and temporal scales of the turbulence. It therefore requires a spatial and time refinement that for industrial applications typically exceeds the computing power of the most advanced supercomputers. As an alternative, RANS and LES method can be used for turbulent flow simulation.

The transition from a laminar to a turbulent flow is linked to a threshold value of the Reynolds number. For a flow in pipes typical Reynolds number below which the flow is laminar is about .

Arguments dtrho and gradrho are respectively the time derivative and the gradient of the density while includetimederivs gives the option to include or not the time-derivative terms in the formulation. In case the density constant argument is set to True, the fluid is supposed incompressible and the Navier-Stokes equations are further simplified since the divergence of the velocity is zero. If the viscosity in space (it does not have to constant in time) the constant viscosity argument can be set to True. By default the density and viscosity are supposed constant and the time-derivative terms are not included. Please note that to simulate the Stokes flow the LBB condtion has to be satisfied. This is achieved by using nodal (h1) type shape functions with an interpolation order of at least one higher for the velocity field than for the pressure field. Alternatively, an additional isotropic diffusive term or other stabilization techniques can be used to over the LBB limitation.

###### # Example
>>> mymesh = mesh("microvalve.msh")
>>> fluid = 2   # physical region
>>>
>>> v=field("h1xy"); p=field("h1")
>>> v.setorder(fluid, 2)
>>> p.setorder(fluid, 1)    # Satisfies the LBB condition
>>>
>>> laminar = formulation()
>>> laminar += integral(fluid, predefinednavierstokes(dof(v), tf(v), v, dof(p), tf(p), 8.9e-4, 1000, 0, 0))

predefinedstokes()

### #predefinednonlinearpoissonwalldistance

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

Generic p-Posison wall distance equation: ∇.(|∇φ|^(p-2) ∇φ) + 1 = 0, φ ∈ Ω, φ = 0, φ ∈ 𝜕Ω --> DBC dφ/dn = 0, φ ∈ 𝜕Ω_N --> NBC where, φ is the approximate distance field p is poisson parameter.

The p-Poisson parameter 𝑝 must be larger or equal to 2. Higher the parameter better the distance field approximation. The term |∇φ|^(p-2) represents an apparent diffusion coefficient, i.e D = |∇φ|^(p-2)

The above system is non-linear and hence an iterative Newton solver is used. After calculating the field 'φ', a better approximation of distance field is obtained as follows: d(φ) = -|∇φ|^(p-1) + (φ*p/(p-1) + |∇φ|^p)^((p-1)/p) --> Eq. (2.12)

Reference: Wall-Distance Calculation for Turbulence Modelling J. C. Bakker, Delft University of Technology http://samofar.eu/wp-content/uploads/2018/10/Bakker_Jelle_BSc-thesis_2018.pdfopen in new window

### #predefinedreciprocalwalldistance

def predefinedreciprocalwalldistance(
physreg: int,
wallreg: int,
reflength: float,
interpolationorder: int = 1,
relresddmtol: float = 1e-12,
maxnumddmit: int = 500,
relerrnltol: float = 1e-06,
maxnumnlit: int = 100,
verbosity: int = 1
) ‑> quanscient.expression

The proposed equation for the recirpocal wall distance (G=1/d) function reads

∇.(G.∇G) + (σ-1)G(∇.∇G) - (1+2σ)G⁴ = 0
DBC --> G = 1/xref, G ∈ 𝜕Ω
NBC --> dG/dn = 0,  G ∈ 𝜕Ω_N
IC  --> Ginit = [0 .. 1e-3]*(1/xref),   G ∈ Ω

Reference: Fares, E., and W. Schröder. "A differential equation for approximate wall distance." International journal for numerical methods in fluids 39.8 (2002): 743-762.

Excerpts from the above paper: "The desired smoothing is controlled by the value of the smoothing parameter 'sigma'. Larger the sigma value means a stronger smoothing at sharp edges, (but large deviation from exact distances). The value of wall boundary condition Gwall influences the smoothing too.

Reference length 'xref' is relevant in the definition of the initial and boundary conditions.
For one sided wall, xref does not play a role, since the solution is the exact distance for all
'σ' and 'Gwall'.

The formulation promises an enhancement of turbulence models at strong curved surfaces."

HINTS: For better approximations of distances --> smaller σ and larger Gwall If the above does not converge, lower the Gwall value.

### #predefinedstabilization

def predefinedstabilization(
stabtype: str,
delta: expression,
f: expression,
v: expression,
diffusivity: expression,
residual: expression
) ‑> quanscient.expression

### #predefinedstokes

def predefinedstokes(
dofv: expression,
tfv: expression,
dofp: expression,
tfp: expression,
mu: expression,
rho: expression,
dtrho: expression,
includetimederivs: bool = False,
isdensityconstant: bool = True,
isviscosityconstant: bool = True,
precondtype: str = ''
) ‑> Tuple[quanscient.expression, quanscient.preconditioner]

This defines the weak formulation for the Stokes (creeping) flow, a linear form of Navier-Stokes where the advective term is ignored as the inertial forces are smaller compared to the viscous forces:

where,

• is the fluid density
• is the dynamic viscosity of the fluid
• is the pressure
• is the flow velocity

This formulation is only valid to simulate the flow of Newtonian fluids (air, water, ...) with a very small Reynolds number ():

where is the characteristic length of the flow. Low flow velocities, high viscosities or small dimensions can lead to a valid Stokes flow approximation. Flows in microscale devices such as microvalves are also good candidates for Stokes flow simulations.

Arguments dtrho and gradrho are respectively the time derivative and the gradient of the density while includetimederivs gives the option to include or not the time-derivative terms in the formulation. In case the density constant argument is set to True, the fluid is supposed incompressible and the Navier-Stokes equations are further simplified since the divergence of the velocity is zero. If the viscosity in space (it does not have to constant in time) the constant viscosity argument can be set to True. By default the density and viscosity are supposed constant and the time-derivative terms are not included. Please note that to simulate the Stokes flow the LBB condtion has to be satisfied. This is achieved by using nodal (h1) type shape functions with an interpolation order of at least one higher for the velocity field than for the pressure field. Alternatively, an additional isotropic diffusive term or other stabilization techniques can be used to over the LBB limitation.

###### # Example
>>> mymesh = mesh("microvalve.msh")
>>> fluid = 2   # physical region
>>>
>>> v=field("h1xy"); p=field("h1")
>>> v.setorder(fluid, 2)
>>> p.setorder(fluid, 1)    # Satisfies the LBB condition
>>>
>>> stokesflow = formulation()
>>> stokesflow += integral(fluid, predefinedstokes(dof(v), tf(v), dof(p), tf(p), 8.9e-4, 1000, 0, 0))

predefinednavierstokes()

def printonrank(
rank: int,
toprint: str
) ‑> None

### #printtotalforce

def printtotalforce(
*args,
**kwargs
) ‑> List[float]

def printvector(
*args,
**kwargs
) ‑> None

### #printversion

def printversion() ‑> None

*args,
**kwargs
) ‑> None

def scatter(
*args,
**kwargs
) ‑> None

### #scatterwrite

def scatterwrite(
filename: str,
xcoords: List[float],
ycoords: List[float],
zcoords: List[float],
compxevals: List[float],
compyevals: List[float] = [],
compzevals: List[float] = []
) ‑> None

This writes to the output file a scalar or vector values at given coordinates. If atleast one of the compyevals or compzevals is not empty then the values saved are vectors and not scalars. For scalars, only compxevals must be provided.

###### # Raises

RuntimeError : if length of all the list arguments are not identical.

###### # Example
>>> Define coordinates of three points: (0.0,0.0,0.0), (1.0,1.0,0.0) and (2.0,2.0,0.0)
>>> coordx = [0.0, 1.0, 2.0]
>>> coordy = [0.0, 1.0, 2.0]
>>> coordz = [0.0, 0.0, 0.0]
>>>
>>> vals = [10, 20, 30]
>>> scatterwrite("scalarvalues.vtk", coordx, coordy, coordz, vals)
>>>
>>> xvals = [10, 20, 30]
>>> yvals = [40, 50, 60]
>>> scatterwrite("vectorvalues.vtk", coordx, coordy, coordz, xvals, yvals)

### #selectall

def selectall() ‑> int

Returns a new or an existing physical region that covers the entire domain.

###### # Reutrns

int Physical region that covers the entire domain.

###### # Example
>>> rega = 1; regb = 2
>>> qa = shape("quadrangle", rega, {0,0,0, 1,0,0, 1,1,0, 0,1,0}, {5,5,5,5})
>>> qb = shape("quadrangle", regb, {1,0,0, 2,0,0, 2,1,0, 1,1,0}, {5,5,5,5})
>>> mymesh = mesh([qa, qb])
>>> wholedomain = selectall()
>>> mymesh.write("mesh.msh")

### #selectintersection

def selectintersection(
physregs: List[int],
intersectdim: int
) ‑> int

Returns a new or an existing physical region that is the intersection of physical regions passed. The intersectdim argument determines the dimensional data from the intersection that would be utilized in subsequent operation such as setting constraint on a physical region or writing a field or expression to the physical region. For, intersectdim=3: from the intersection, only volumes are used in subsequent operations. intersectdim=2: from the intersection, only surfaces are used in subsequent operations. intersectdim=1: from the intersection, only lines are is used in subsequent operations. intersectdim=0: from the intersection, only points are is used in subsequent operations. This is useful in isolating only those dimensional data that might be necessary and ignoring the others arising from the intersection. Note that, the intersected region itself is not affected by intersectdim. This argument only determines which dimensional data is utilized when the intersected physical region is used in calculations or operations.

###### # Parameters

physregs : list : List of physical regions

intersectdim : int : 0, 1, 2 or 3. This determines which dimensional data from the intersected physical region is returned.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3              # physical regions
>>>
>>> # Influence of <code>intersectdim</code>
>>> intersectdim = 3
>>> intersectedreg = selectintersection([vol, vol], intersectdim)
>>> expression(1).write(intersectedreg, "out_3D.vtk", 1)    # uses only volume from the intersectedreg
>>>
>>> intersectdim = 2
>>> intersectedreg = selectintersection([vol, vol], intersectdim)
>>> expression(1).write(intersectedreg, "out_2D.vtk", 1)    # uses only surfaces from the intersectedreg
>>>
>>> intersectdim = 1
>>> intersectedreg = selectintersection([vol, vol], intersectdim)
>>> expression(1).write(intersectedreg, "out_1D.vtk", 1)    # uses only lines from the intersectedreg
>>>
>>> intersectdim = 0
>>> intersectedreg = selectintersection([vol, vol], intersectdim)
>>> expression(1).write(intersectedreg, "out_0D.vtk", 1)    # uses only points from the intersectedreg

### #selectnooverlap

def selectnooverlap() ‑> int

Returns a new or an existing physical region that covers no-overlap domain in case of overlap DDM and the entire domain otherwise.

### #selectunion

def selectunion(
physregs: List[int]
) ‑> int

Returns a new or an existing physical region that is the union of physical regions passed.

###### # Parameters

physregs : list : List of physical regions

###### # Example
>>> mymesh = mesh("disk.msh")
>>> sur=2; top=3                    # physical regions
>>> surandtop = selectunion([2,3])  # unioned physical region

def send(
*args,
**kwargs
) ‑> None

### #setaxisymmetry

def setaxisymmetry() ‑> None

This call should be placed at the very beginning of the code. After the call everything will be solved assuming axisymmetry (works for 2D meshes in the xy plane only). All equations should be written in their 3D form.

Please note that in order to correctly take into account the cylindrical coordinate change, the appropriate space derivative operators should be used. For example, the gradient of a vector operator required in the mechanical strain calculation to compute the gradient of mechanical displacement should not be defined manually using dx(), dy() and dz() space derivatives. The grad() operator should instead be called on the mechanical displacement vector.

###### # Example
>>> setaxisymmetry()

def setdata(
invec: vec
) ‑> None

### #setfundamentalfrequency

def setfundamentalfrequency(
f: float
) ‑> None

This defines the fundamental frequency (in ) required for multiharmonic problems.

###### # Example
>>> setfundamentalfrequency(50)

mnt: int
) ‑> None

Sets the maximum number of threads allowed.

###### # Parameters

mnt : int : input value specifying the maximum number of threads allowed.

2

### #setphysicalregionshift

def setphysicalregionshift(
shiftamount: int
) ‑> None

### #settime

def settime(
t: float
) ‑> None

This sets the time variable t.

###### # Example
>>> settime(1e-3)

### #settimederivative

def settimederivative(
*args,
**kwargs
) ‑> None

### #sin

def sin(
input: expression
) ‑> quanscient.expression

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

###### # Example
>>> expr = sin(getpi()/2)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
1
>>>
>>> expr = sin(2)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
0.909297

### #sn

def sn(
n: float
) ‑> quanscient.expression

### #solve

def solve(
*args,
**kwargs
) ‑> List[quanscient.vec]

### #sqrt

def sqrt(
input: expression
) ‑> quanscient.expression

This returns an expression that is square root of input expression.

###### # Example
>>> expr = sqrt(2)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
1.41421

### #strain

def strain(
input: expression
) ‑> quanscient.expression

This defines the (linear) engineering strains in Voigt form . The input can either be the displacement field or its gradient.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> engstrain = strain(u)
>>> engstrain.print()

def sum(
*args,
**kwargs
) ‑> None

### #t

def t() ‑> quanscient.expression

This gives the time variable in the form of an expression. The evaluation gives a value equal to gettime().

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setconstraint(vol, sin(2*t()))

### #tan

def tan(
input: expression
) ‑> quanscient.expression

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

###### # Example
>>> expr = tan(getpi()/4)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
1
>>>
>>> expr = tan(1)
>>> expr.print()
Expression size is 1x1
@ row 0, col 0 :
1.55741

### #tangent

def tangent() ‑> quanscient.expression

This defines a tangent vector with unit norm.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> top=3
>>> tangent().write(top, "tangent.vtk", 1)

### #tf

def tf(
*args,
**kwargs
) ‑> quanscient.expression

### #trace

def trace(
a: expression
) ‑> quanscient.expression

This computes the trace of a square matrix expression. The returned expression is a scalar.

###### # Example
>>> a = array2x2(1,2, 3,4)
>>> tracea = trace(a)
>>> tracea.print()
Expression size is 1x1
@ row 0, col 0 :
5

### #transpose

def transpose(
input: expression
) ‑> quanscient.expression

This returns an expression that is the transpose of a vector or matrix expression.

###### # Example
>>> colvec = array3x1(1,2,3)
>>> rowvec = transpose(colvec)
>>> rowvec.print()
Expression size is 1x3
@ row 0, col 0 :
1
@ row 0, col 1 :
2
@ row 0, col 2 :
3
>>> matexpr = expression(3,3, [1,2,3, 4,5,6, 7,8,9])
>>> transposed = transpose(matexpr)

### #vonmises

def vonmises(
stress: expression
) ‑> quanscient.expression

This returns the von Mises stress expression coreesponding to the 3D stress tensor provided as argument. The stress tensor should be provided in Voigt form .

##### #setrelaxationfactor
def setrelaxationfactor(
self,
relaxfact: float
) ‑> None

This sets the relaxation factor for the fixed-point nonlinear iteration performed at every timestep for nonlinear problems. If the relaxation factor is not set, the default value of is set. If is the solution obtained at a current iteration, is solution at previous teration, then the new solution at the current iteration is updated as

where is the relaxation factor.

##### #settimederivative
def settimederivative(
self,
sol: List[vec]
) ‑> None

This sets the current solution for the time derivatives and to sol[0] and sol[1] respectively.

genalpha.gettimederivative()

##### #settimestep
def settimestep(
self,
timestep: float
) ‑> None

This sets the current timestep.

##### #settolerance
def settolerance(
self,
nltol: float
) ‑> None

This sets the tolerance for the fixed-point nonlinear iteration performed at every timestep for nonlinear problems. If the tolerance is not set, the default value of is considered.

##### #setverbosity
def setverbosity(
self,
verbosity: int
) ‑> None

This sets the verbosity level. For debugging, higher verbosity is recommended. If the verbosity is not set, the default value of is considered.

### #impliciteuler

class impliciteuler(
formul: formulation,
dtxinit: vec,
verbosity: int = 3,
isrhskcconstant: List[bool] = [False, False, False]
)

This defines the impliciteuler object to solve in time the formulation formul with the fields' state as initial solution for and dtxinit for . The isrhskcconstant list can be used to specify whether the RHS vector, K matrix and C matrix are constant in time (default is False). If isrhskcconstant[i]=True, the corresponding vector/matrix are generated once and then reused. In rhskcconstant[i]:

• corresponds to the rhs vector
• corresponds to the K matrix
• corresponds to the C matrix

If the K and C matrix are constant in time, the factorization of the algebraic problem is also reused.

The impliciteuler object allows to perform an implicit (backward) Euler time resolution for a problem of the form , be it linear or nonlinear. The solutions for as well as are made available. For nonlinear problems a fixed-point iteration is performed at every timestep until the relative error (norm of relative solution vector change) is less than the prescribed tolerance.

#### # Notes

Even if the rhs vector can be reused the Dirichlet constraints will neverthless be recomputed at each timestep.

#### # Methods

##### #allnext
def allnext(
*args,
**kwargs
) ‑> None

This is a collective MPI operation and hence must be called by all ranks. It is similar to the impliciteuler.next() function but the resolution is performed on all ranks using DDM. This method runs the implicit Euler algorithm for one timestep. After the call, the time and field values are updated on all the DDM ranks. This method can be used for both linear and nonlinear problems.

###### # Examples

Example 1: impliciteuler.allnext(relrestol: double, maxnumit: int, timestep: double) This is used for linear problems. The relrestol and maxnumit arguments correspond to the stopping criteria for DDM solver. The DDM iterations stop if either relative residual tolerance is less than relrestol or if the number of DDM iteration reaches maxnumit. Use timestep=-1 for automatic time-adaptivity.

Example 2: impliciteuler.next(relrestol: double, maxnumit: int, timestep: double, maxnumnlit:int) This is used for nonlinear problems. The relrestol and maxnumit arguments correspond to the stopping criteria for DDM solver. The DDM iteration stops if either relative residual tolerance is less than relrestol or if the number of DDM iteration reaches maxnumit. A nonlinear fixed-point iteration is performed for at most maxnumnlit or until the relative error (norm of relative solution vector change) is smaller than the tolerance prescribed in impliciteuler.settolerance(). Set timestep=-1 for automatic time-adaptivity and maxnumnlit=-1 for unlimited nonlinear iterations. This method returns the number of nonlinear iterations performed.

impliciteuler.next()

##### #count
def count(
self
) ‑> int

This counts the total number of steps computed.

##### #gettimederivative
def gettimederivative(
self
) ‑> quanscient.vec

This returns the current solution for the first time derivative .

impliciteuler.settimederivative()

##### #gettimes
def gettimes(
self
) ‑> List[float]

This returns all the time values stepped through.

##### #gettimestep
def gettimestep(
self
) ‑> float

This returns the current timestep.

##### #next
def next(
*args,
**kwargs
) ‑> None

This runs the implicit Euler algorithm for one timestep. After the call, the time and field values are updated. This method can be used for both linear and nonlinear problems depending on the number of arguments passed.

###### # Examples

Example 1: implicit.next(timestep: double) This is used for linear problems. \boldsymbol{Use for automatic time-adaptivity}.

Example 2: implicit.next(timestep: double, maxnumnlit:int) This is used for nonlinear problems. A nonlinear fixed-point iteration is performed for at most maxnumnlit or until the relative error (norm of relative solution vector change) is smaller than the tolerance prescribed in impliciteuler.settolerance(). Set timestep=-1 for automatic time-adaptivity and maxnumnlit=-1 for unlimited nonlinear iterations. This method returns the number of nonlinear iterations performed.

##### #postsolve
def postsolve(
self,
formuls: List[formulation]
) ‑> None

This defines the set of formulations that must be solved every resolution of the formulation provided to the impliciteuler constructor. The formulations provided here must lead to a system of the form (no damping or mass matrix allowed).

impliciteuler.presolve()

##### #presolve
def presolve(
self,
formuls: List[formulation]
) ‑> None

This defines the set of formulations that must be solved every resolution of the formulation provided to the impliciteuler constructor. The formulations provided here must lead to a system of the form (no damping or mass matrix allowed).

impliciteuler.postsolve()

self,
tol: float,
mints: float,
maxts: float,
reffact: float = 0.5,
coarfact: float = 2.0,
coarthres: float = 0.5
) ‑> None

This sets the configuration for automatic time-adaptivity. The timestep will be adjusted between the minimum timestep () and maximum timesteps () to reach the requested relative error tolerance (). The relative error is defined as

to measure the relative deviation from a constant time derivative.

Arguments and gives the factor to use when the time step is refined or coarsened respectively. The timestep is refined when the relative error is above or when the maximum number of nonlinear iterations is reached. The timestep is coarsened when the relative error is below the product and the nonlinear loop has converged in less than the maximum number of iterations.

##### #setrelaxationfactor
def setrelaxationfactor(
self,
relaxfact: float
) ‑> None

This sets the relaxation factor for the fixed-point nonlinear iteration performed at every timestep for nonlinear problems. If the relaxation factor is not set, the default value of is set. If is the solution obtained at a current iteration, is solution at previous teration, then the new solution at the current iteration is updated as

where is the relaxation factor.

##### #settimederivative
def settimederivative(
self,
sol: vec
) ‑> None

This sets the current solution for the first time derivatives to sol[0].

impliciteuler.gettimederivative()

##### #settimestep
def settimestep(
self,
timestep: float
) ‑> None

This sets the current timestep.

##### #settolerance
def settolerance(
self,
nltol: float
) ‑> None

This sets the tolerance for the fixed-point nonlinear iteration performed at every timestep for nonlinear problems. If the tolerance is not set, the default value of is considered.

##### #setverbosity
def setverbosity(
self,
verbosity: int
) ‑> None

This sets the verbosity level. For debugging, higher verbosity is recommended. If the verbosity is not set, the default value of is considered.

### #indexmat

class indexmat(
**kwargs
)

The indexmat object stores a row-major array of integers that corresponds to a dense matrix. For storing array of doubles, see densemat object.

#### # Examples

There are number of ways of instatiating an indexmat object. There are listed below:

Example 1: indexmat(numberofrows:int, numberofcolumns:int) The following creates a matrix with 2 rows and 3 columns. The entries may be undefined.

>>> B = indexmat(2,3)

Example 2: indexmat(numberofrows:int, numberofcolumns:int, initvalue:int) This creates a matrix with 2 rows and 3 columns. All entries are assigned the value initvalue.

>>> B = indexmat(2,3, 12)
>>> B.print()
Matrix size is 2x3
12 12 12
12 12 12

Example 3: indexmat(numberofrows:int, numberofcolumns:int, valvec:List[int]) This creates a matrix with 2 rows and 3 columns. The entries are assigned the values of valvec. The length of valvec is expected to be equal to the total count of entries in the matrix. So for creating a matrix of size , length of valvec must be 6.

>>> B = indexmat(2,3, [1,2,3,4,5,6])
>>> B.print()
Matrix size is 2x3
1 2 3
4 5 6

Example 4: indexmat(numberofrows:int, numberofcolumns:int, init:int, step:int) This creates a matrix with 2 rows and 3 columns. The first entry is assigned the value init and the consecutive entries are assigned values that increase by steps of step.

>>> B = indexmat(2,3, 0, 1)
>>> B.print()
Matrix size is 2x3
0 1 2
3 4 5

Example 5: indexmat(input:List[indexmat]) This creates a matrix that is the vertical concatenation of input matrices. Since, the concatenation occurs vertically, the number of columns in all the input matrices must match.

>>> A = indexmat(2,3, 0)
>>> B = indexmat(1,3, 2)
>>> AB = indexmat([A,B])
>>> AB.print()
Matrix size is 3x3
0 0 0
0 0 0
2 2 2

#### # Methods

##### #count
def count(
self
) ‑> int

This counts and returns the total number of entries in the dense matrix.

###### # Example
>>> B = indexmat(2,3)
>>> B.count()
6
##### #countcolumns
def countcolumns(
self
) ‑> int

This counts and returns the number of columns in the dense matrix.

###### # Example
>>> B = indexmat(2,3)
>>> B.countcolumns()
3
##### #countrows
def countrows(
self
) ‑> int

This counts and returns the number of rows in the dense matrix.

###### # Example
>>> B = indexmat(2,3)
>>> B.countrows()
2
##### #print
def print(
self
) ‑> None

This prints the entries of the dense matrix.

###### # Example
>>> B = indexmat(2,3, 0,1)
>>> B.print()
Matrix size is 2x3
0 1 2
3 4 5
##### #printsize
def printsize(
self
) ‑> None

This prints the size of the dense matrix.

###### # Example
>>> B = indexmat(2,3)
>>> B.printsize()
Matrix size is 2x3

### #integration

class integration

### #logformat

class logformat(
value: int
)

Holds possible log formats

Members -----=

plain

json

#### # Class variables

##### #json

Type: quanscient.logformat

##### #plain

Type: quanscient.logformat

### #mat

class mat(
**kwargs
)

The mat object holds a sparse algebraic square matrix.

#### # Raises

RuntimeError : If a mesh object is not available prior to creating a mat object.

#### # Notes

Before creating a mat object, ensure that a mesh object is available. If a mesh object is not already available, create an empty mesh object.

#### # Examples

There are number of ways of instatiating an indexmat object. There are listed below:

Example 1: mat(matsize:int, rowaddresses:indexmat, coladdresess:indexmat, vals:densemat) This creates a sparse matrix object of size matsizematsize. The rowaddresses and coladdresess provide the location (row,col) of non-zero values in the sparse matrix. The non-zero values are provided in the dense matrix vals. Note that a mesh object must already be available before instatiating mat object.

>>> rows = indexmat(7,1, [0,0,1,1,1,2,2])
>>> cols = indexmat(7,1, [0,1,0,1,2,1,2])
>>> vals = densemat(7,1, [11,12,13,14,15,16,17])
>>>
>>> mymesh = mesh()
>>> A = mat(3, rows, cols, vals)
>>> A.print()
A block 3x3:

Mat Object: 1 MPI processes type: seqaij row 0: (0, 11.) (1, 12.) row 1: (0, 13.) (1, 14.) (2, 15.) row 2: (1, 16.) (2, 17.)

D block 3x0:

Mat Object: 1 MPI processes type: seqaij row 0: row 1: row 2:

Example 2: mat(myformulation:formulation, rowaddresses:indexmat, coladdresses:indexmat, vals:densemat) This creates a sparse matrix object whose dof() structure is the one in the formulation projection. The rowaddresses and coladdresess provide the location (row,col) of non-zero values in the sparse matrix. The non-zero values are provided in the dense matrix vals.

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>>
>>> addresses = indexmat(numberofrows=projection.countdofs(), numberofcolumns=1, init=0, step=1)
>>> vals = densemat(numberofrows=projection.countdofs(), numberofcolumns=1, init=12)
>>>

#### # Methods

##### #copy
def copy(
self
) ‑> quanscient.mat

This creates a full copy of the matrix. Only the values are copied. (E.g: the mat.reusefactorization() is set back to the default no resue.)

###### # Example
>>> A =
>>> copiedmat = A.copy()
##### #countcolumns
def countcolumns(
self
) ‑> int

This counts and returns the number of columns in the matrix.

###### # Example
>>> numcols = A.countcolumns()
##### #countnnz
def countnnz(
self
) ‑> int

This counts and returns the number of non-zero entries in the matrix which is the sub-matrix of with eliminated Dirichlet constraints. Refer mat.getainds().

If the requested information is not available, then is returned.

###### # Example
>>> numnnz = A.countnnz()
##### #countrows
def countrows(
self
) ‑> int

This counts and returns the number of rows in the matrix.

###### # Example
>>> numrows = A.countrows()
##### #getainds
def getainds(
self
) ‑> quanscient.indexmat

Let us call dinds the set of unknowns that have a Dirichlet constraint and ainds the remaining unknowns. The mat object holds sub-matrices and such that

where is a square matrix equal to with eliminated Dirichlet constraints. is an all zero matrix and is the square identity matrix of all Dirichlet constraints. Matrices and are stored with their local indexing. The methods mat.getainds() and mat.getdinds() gives the global indexing (i.e index in ) of each local index in and .

###### # Example
>>> ainds = A.getainds()

mat.getdinds()

##### #getdinds
def getdinds(
self
) ‑> quanscient.indexmat

This outputs dinds.

###### # Example
>>> dinds = A.getdinds()

mat.getainds()

##### #print
def print(
self
) ‑> None

This prints the matrix size and values.

>>> A.print()
##### #reusefactorization
def reusefactorization(
self
) ‑> None

The matrix factorization will be reused in allsolve().

### #mesh

class mesh(
**kwargs
)

The mesh object holds the finite element mesh of the geometry.

#### # Examples

A mesh object based on a mesh file can be created through native reader or via the gmsh api. To get more information on the physical regions of the mesh, the verbosity argument can be set ot .

>>> # Creating a mesh object with the native reader:
>>> mymesh = mesh("disk.msh")
>>>
>>> # Creating a mesh object with GMSH API:
>>> mymesh = mesh("gmsh:disk.msh")

In the domain decomposition framework, creating a mesh object requires two additional arguments: globalgeometryskin and numoverlaplayers. Furthermore, the mesh is treated as a part of a global mesh. Each MPI rank owns only a part of the global mesh and all ranks must perform the call collectively. The argument globalgeometryskin is the part of the global mesh skin that belongs to the current rank. It can only hold elements of dimension one lower than geometry dimension (use -1 if empty). The global mesh skin cannot intersect itself. The mesh parts are overlapped by the number of overlap layers requested. More than one overlap layer cannot be guaranteed everywhere as the overlapping is limited to the direct neighbouring domains. mesh(filename:str, globalgeometryskin:int, numoverlaplayers:int, verbosity:int=1)

In the above examples, the mesh object were created based on a mesh file. Similarly, mesh objects can be created based on shape objects.

>>> # define physical regions
>>> faceregionnumber=1; lineregionnumber=2
>>>
>>> # define a quadrangle shape object
>>>
>>> # get the leftline from the contour of the quad shape object
>>> contourlines = quadface.getsons()  # returns a list
>>> leftline = contourlines[3]
>>> leftline.setphysicalregion(lineregionnumber)
>>>
>>> # create mesh object based on the quadrangle shape and its left side line

Creating a mesh object from the shape object can be carried out also in the domain decomposition framework using the following syntax: mesh(inputshapes:List[shape], globalgeometryskin:int, numoverlaplayers:int, verbosity:int=1)

It is also possible to combine multiple meshes. Elements shared by the input meshes can either be merged or not by setting the bool value for argument mergeduplicates. For every input mesh a new physical region containing all elements is created. Set verbosity equal to 2 to get information on physical regions in the mesh.

>>> mymesh = mesh("disk.msh")
>>> mymesh.shift(2,1,0)     # all entities are shifted by x,y,z amount
>>> mymesh.write("shifted.msh")
>>>
>>> mergedmesh = mesh(True, ["disk.msh", "shifted.msh"])
>>> mergedmesh.write("merged.msh")
>>>
>>> mergedmesh = mesh("merged.msh", 2)

#### # Methods

##### #extrude
def extrude(
self,
newphysreg: int,
newbndphysreg: int,
bnd: int,
extrudelens: List[float]
) ‑> None

This extrudes the boundary region bnd. After the extrusion process, newphysreg will contain the extruded region and newbndphysreg will contain the extrusion end boundary.

###### # Parameters

newphysreg : int : The physical region that will contain the extruded region.

newbndphysreg : int : The physical region that will contain the end boundary of extrusion.

bnd : int : The boundary region that is extruded.

extrudelens : List[float] : List specifiying the size of each layer in the extrusion. Length of list determines the number of mesh layers in the extrusion. If is given as the extrusion length for each layer, an optimal value is automatically calculated.

###### # Example

We use the disk.msh for the example here.

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

getextrusiondata()

##### #getdimension
def getdimension(
self
) ‑> int

This returns the dimension of the highest dimension element in the mesh 0D, 1D, 2D or 3D.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> dim = mymesh.getdimension()
>>> dim
3
##### #getdimensions
def getdimensions(
self
) ‑> List[float]

This returns the x, y and z mesh dimensions in meters.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> dims = mymesh.getdimensions()
>>> dims
[2.0, 2.0, 0.1]
##### #getphysicalregionnumbers
def getphysicalregionnumbers(
self,
dim: int = -1
) ‑> List[int]

This returns all physical region numbers of a given dimension. Use or no argument to get the regions of all dimensions.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> allphysregs = mymesh.getphysicalregionnumbers()
[4, 2, 3, 1]
*args,
**kwargs
) ‑> None

This method allows an empty mesh object to be populated with mesh data. It takes in the same corresponding arguments as required in instantiating a mesh object directly. The only difference with direct instantiation is that this methid requires that an empty mesh object is already created. If this method is called by a non-empty mesh object any existing mesh data are lost.

###### # Examples
>>> # Create an empty mesh object
>>> mymesh = mesh()
>>>
>>>
>>> # Load a mesh file with GMSH API:
>>> mymesh = mesh("gmsh:disk.msh", 2)

>>> # define physical regions
>>> faceregionnumber=1; lineregionnumber=2
>>>
>>> # define a quadrangle shape object
>>>
>>> # get the leftline from the contour of the quad shape object
>>> contourlines = quadface.getsons()  # returns a list
>>> leftline = contourlines[3]
>>> leftline.setphysicalregion(lineregionnumber)
>>>
>>> # Load a mesh from the quadrangle shape and its left side line
>>> mymesh = mesh()     # creates an empty mesh object

Combing multiple meshes with load method

>>> mymesh = mesh("disk.msh")
>>> mymesh.shift(2,1,0)     # all entities are shifted by x,y,z amount
>>> mymesh.write("shifted.msh")
>>>
>>> mergedmesh = mesh()
>>> mergedmesh.write("merged.msh")
>>>
>>> mergedmesh = mesh("merged.msh", 2)

Notes: cohomologycut, split, partition, skin, selectskin, selectbox, selectsphere, selectlayer, selectexclusion, seleetanynode

##### #move
def move(
*args,
**kwargs
) ‑> None

This moves the whole or part of the mesh object by the x, y and z component of expression u in the x, y and z direction.

###### # Examples
>>> vol = 1
>>> mymesh = mesh("disk.msh")
>>> x=field("x"); y=field("y")
>>>
>>> # Move the whole mesh object
>>> mymesh.move(array3x1(0,0, sin(x*y)))
>>> mymesh.write("moved.msh")
>>>
>>> # Move only the mesh part on physical region 'vol'
>>> mymesh.move(vol, array3x1(0,0, sin(x*y)))
>>> mymesh.write("moved.msh")
def partition(
*args,
**kwargs
) ‑> None
##### #printdimensions
def printdimensions(
self
) ‑> List[float]

This prints and returns the x, y and z mesh dimensions.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> mymesh.printdimensions()
Mesh dimensions:
x: 2 m
y: 2 m
z: 0.1 m
[2.0, 2.0, 0.1]
##### #rotate
def rotate(
*args,
**kwargs
) ‑> None

This rotates the whole or part of the mesh object first by ax degrees around x axis followed by ay degrees around y axis and then by az degrees around z axis.

###### # Examples
>>> vol = 1
>>> mymesh = mesh("disk.msh")
>>>
>>> # rotate the whole mesh object
>>> mymesh.rotate(20,60,90)
>>> mymesh.write("rotated.msh")
>>>
>>> # rotate only the mesh part on physical region 'vol'
>>> mymesh.rotate(vol, 20,60,90)
>>> mymesh.write("rotated.msh")
##### #scale
def scale(
*args,
**kwargs
) ‑> None

This scales the whole or part of the mesh object first by a factor x, y and z respectively in the x, y and z direction.

###### # Examples
>>> vol = 1
>>> mymesh = mesh("disk.msh")
>>>
>>> # scale the whole mesh object
>>> mymesh.scale(0.1,0.2,1.0)
>>> mymesh.write("scaled.msh")
>>>
>>> # scale only the mesh part on physical region 'vol'
>>> mymesh.scale(vol, 0.1,0.2,1.0)
>>> mymesh.write("scaled.msh")
##### #selectanynode
def selectanynode(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains a single node arbitrarily chosen in the region physregtoselectfrom. If no region is selected (i.e. if physregtoexcludefrom is empty) or if the argument physregtoexcludefrom is not provided, then the arbitraray node is chosen considering the whole domain. The new region newphysreg is created when the mesh.load() method is called on the mesh object.

###### # Examples
>>> vol=1; anynode=12
>>> mymesh = mesh()
>>> mymesh.selectanynode(excluded, vol, [box])   # a node is chosen from 'vol' region
>>> # (OR)
>>> # mymesh.selectanynode(excluded, [box])      # a node is chosen from the whole domain
>>> mymesh.write("out.msh")
##### #selectbox
def selectbox(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains elements of the region physregtoselectfrom that are in the box delimited by [,, ,, ,] given in boxlimit. If no region is selected (i.e. if physregtoselectfrom is empty) or if the argument physregtoselectfrom is not provided, then the box region is created considering the whole domain.

The new region newphysreg is created when the mesh.load() method is called on the mesh object. The elements populated in the new region newphysreg are of dimension selecteddim.

###### # Examples
>>> vol=1; boxregion=12
>>> mymesh = mesh()
>>> mymesh.selectbox(boxregion, vol, 3, [0,1,0, 1,0,0.1])   # select box region from the 'vol' region
>>> # (OR)
>>> # mymesh.selectbox(boxregion, 3, [0,1,0, 1,0,0.1])      # select box region from the whole domain
>>>
>>> v=field("h1"); x=field("x"); y=field("y"); z=field("z")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, x*y*z)
>>> v.write(vol, "v.vtk", 1)
>>> v.write(boxregion, "vboxregion.vtk", 1)
##### #selectexclusion
def selectexclusion(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains the elements of the region physregtoexcludefrom that are not in physregtoexclude. If no region is selected (i.e. if physregtoexcludefrom is empty) or if the argument physregtoexcludefrom is not provided, then the new region is created considering the whole domain. The new region newphysreg is created when the mesh.load() method is called on the mesh object.

###### # Examples
>>> vol=1; sur=2; top=3; box=11; excluded=12
>>> mymesh = mesh()
>>> mymesh.selectbox(box, vol, 3, [0,2, -2,2, -2,2])
>>> mymesh.selectexclusion(excluded, vol, [box])   # physregtoexcludefrom = 'vol'
>>> # (OR)
>>> # mymesh.selectexclusion(excluded, [box])      # physregtoexcludefrom = whole domain
>>> mymesh.write("out.msh")
##### #selectlayer
def selectlayer(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains the layer of elements of the region physregtoselectfrom that touches the region physregtostartgrowth. If no region is selected (i.e. if physregtoselectfrom is empty) or if the argument physregtoselectfrom is not provided, then the layer region is created considering the whole domain. When multiple layers are requested through the argument numlayers, they are grown on top of each other. The new region newphysreg is created when the mesh.load() method is called on the mesh object.

###### # Examples
>>> vol=1; sur=2; top=3; layerregion=12
>>> mymesh = mesh()
>>> mymesh.selectlayer(layerregion, vol, sur, 1)   # select layer region from the 'vol' region
>>> # (OR)
>>> # mymesh.selectlayer(layerregion, sur, 1)      # select layer region from the whole domain
>>> mymesh.write("out.msh")
##### #selectskin
def selectskin(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains elements that forms the skin of the selected physical regions. If no region is selected (i.e. if physregtoselectfrom is empty) or if the argument physregtoselectfrom is not provided, then the skin region is created considering the whole domain.

The skin region newphysreg is created when the mesh.load() method is called on the mesh object. The dimension of the skin region is always one dimension less than that of the physical regions selected. Note that space derivatives or 'hcurl' field evaluations on a surface do not usually lead to the same values as a volume evaluation.

###### # Examples
>>> vol=1; skin=12
>>> mymesh = mesh()
>>> mymesh.selectskin(skin, vol)    # select skin region from the 'vol' region
>>> # (OR)
>>> # mymesh.selectskin(skin)         # select skin region from the whole domain
>>>
>>> v=field("h1"); x=field("x"); y=field("y"); z=field("z")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, x*y*z)
>>> v.write(vol, "v.vtk", 1)
>>> v.write(skin, "vskin.vtk", 1)
##### #selectsphere
def selectsphere(
*args,
**kwargs
) ‑> None

This tells the mesh object to create a new physical region newphysreg that contains elements of the region physregtoselectfrom that are in the sphere of prescribed radius and of center [, ,] as given in centercoords. If no region is selected (i.e. if physregtoselectfrom is empty) or if the argument physregtoselectfrom is not provided, then the sphere region is created considering the whole domain.

The new region newphysreg is created when the mesh.load() method is called on the mesh object. The elements populated in the new region newphysreg are of dimension selecteddim.

###### # Examples
>>> vol=1; sphereregion=12
>>> mymesh = mesh()
>>> mymesh.selectsphere(sphereregion, vol, 3, [1,0,0], 1)   # select sphere region from the 'vol' region
>>> # (OR)
>>> # mymesh.selectsphere(sphereregion, 3, [1,0,0], 1)      # select sphere region from the whole domain
>>>
>>> v=field("h1"); x=field("x"); y=field("y"); z=field("z")
>>> v.setorder(vol, 1)
>>> v.setvalue(vol, x*y*z)
>>> v.write(vol, "v.vtk", 1)
>>> v.write(sphereregion, "vsphereregion.vtk", 1)
self,
criterion: expression,
lownumsplits: int,
highnumsplits: int
) ‑> None

Each element in the mesh will be adapted (refined/coarsened) based on the value of a positive criterion (h-adaptivity). The max range of the criterion is split into a number of intervals equal to the number of refinement levels in range lownumsplits and highnumsplits. All intervals have the same size. The barycenter value of the criterion on each element is considered to select the interval, and therefore the corresponding refinement of each mesh element. As an example, for a criterion with highest value 900 over the entire domain and a low/high refinement level requested of 1/3 the refinement on mesh elements with criterion value in range 0 to 300, 300 to 600, 600 to 900 will be 1, 2, 3 levels respectively.

###### # 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)
>>>
>>>
>>> for i in range(5):
...     criterion.write(all, f"criterion_{100+i}.vtk", 1)
##### #setcohomologycuts
def setcohomologycuts(
*args,
**kwargs
) ‑> None

This makes the mesh object aware of the cohomology cut regions.

###### # Example
>>> mymesh = mesh()
>>> mymesh.setcohomologycuts([chreg1, chreg2])
##### #setphysicalregions
def setphysicalregions(
self,
dims: List[int],
nums: List[int],
geometryentities: List[List[int]]
) ‑> None
##### #shift
def shift(
*args,
**kwargs
) ‑> None

This translates the whole or part of the mesh object by x, y and z amount in the x, y and z direction.

###### # Examples
>>> vol = 1
>>> mymesh = mesh("disk.msh")
>>> x=field("x"); y=field("y")
>>>
>>> # shift/translate the whole mesh object
>>> mymesh.shift(1.0, 2.0, 3.0)
>>> mymesh.write("shifted.msh")
>>>
>>> # shift/translate only the mesh part on physical region 'vol'
>>> mymesh.shift(vol, 1.0, 2.0, 3.0)
>>> mymesh.write("shifted.msh")
##### #split
def split(
self,
n: int = 1
) ‑> None

This splits each element in the mesh n times. Element quality is maximized and element curvature is taken into account. Each element is split recursively n times as follows:

• point \leftarrow 1 point
• line \leftarrow 2 lines
• triangle \leftarrow 4 triangles
• tetrahedron \leftarrow 8 tetrahedra
• hexahedron \leftarrow 8 hexahedra
• prism \leftarrow 8 prisms
• pyramid \leftarrow 6 pyramids + 4 tetrahedra
###### # Example
>>> mymesh = mesh()
>>> mymesh.split()
>>> mymesh.write("splitdisk.msh")
##### #use
def use(
self
) ‑> None

This allows to select which mesh to use in case multiple meshes are available. This call invalidates all objects that are based on the previously selected mesh for as long as the latter is not selected again.

###### # Example
>>> finemesh = mesh()
>>> finemesh.split(2)
>>> coarsemesh = mesh("disk.msh")
>>> finemesh.use()
##### #write
def write(
*args,
**kwargs
) ‑> None

This writes the mesh object to a given input filename.

###### # Examples
>>> # mesh data
>>> mymesh = mesh("disk.msh")
>>> vol=1; top=2; sur=3

If a physical region is passed in the first argument, then only part of the mesh object included in that physical region is written:

>>> mymesh.write(vol, "out.msh")

If only the file name is provided as an argument, then all the physical regions of the mesh are written.

>>> mymesh.write("out.msh)
>>> # or equivalently:
>>> mymesh.write("out1.msh", physregs=[-1], option=1)

The argument physregs is the list of physical regions that will be written if the argument option=1. If option=-1 then all the physical regions except the ones in the list physregs will be written. The default value for physregs=-1 which is equivalent to considering all the physical regions (the option argument is ignored when physregs=-1).

>>> mymesh.write("out2.msh", [-1], 1)       # all physical regions are written
>>> mymesh.write("out3.msh", [-1], -1)      # all physical region will be written, ignores 'option' argument
>>> mymesh.write("out4.msh", [1,2], 1)      # physical regions 1 and 2 will be written
>>> mymesh.write("out5.msh", [1,2], -1)     # all physical regions except 1 and 2 will be written

### #parameter

class parameter(
**kwargs
)

The parameter object can hold different expression objects on different geometric regions.

#### # Examples

A parameter object can be a scalar. The following creates an empty object.

>>> mymesh = mesh("disk.msh")
>>> E = parameter()

A parameter object can also be a 2D array. The following creates an empty object.

>>> E = parameter(3,3)  # a 3x3 parameter matrix

#### # Methods

self,
physreg: int,
input: expression
) ‑> None
##### #allintegrate
def allintegrate(
*args,
**kwargs
) ‑> float

This is a collective MPI operation and hence must be called by all ranks. This integrates a parameter across all the DDM ranks.

###### # Examples

...

>>> myparameter = parameter()
>>> myparameter.setvalue(vol, 12.0)
>>> integralvalue = myparameter.allintegrate(vol, 4)
>>>
>>> # allintegrate on a mesh deformed by field 'u'
>>> integralvalueondeformedmesh = myparameter.allintegrate(vol, u, 4)
##### #allinterpolate
def allinterpolate(
*args,
**kwargs
) ‑> List[float]

This is a collective MPI operation and hence must be called by all ranks. Its functionality is as described in pararmeter.interpolate but considers the physical region partitioned across the DDM ranks. The xyz coordinate argument must be the same for all ranks.

###### # Examples
>>> ...
>>> p = parameter(3,1)
>>> p.setvalue(vol, array3x1(x,y,z))
>>> interpolated = p.allinterpolate (vol, {0.5,0.6,0.05});
>>>
>>> # allinterpolate on a mesh deformed by field 'u'
>>> interpolated = p.allinterpolate (vol, u, {0.5,0.6,0.05});
##### #allmax
def allmax(
*args,
**kwargs
) ‑> List[float]

This is a collective MPI operation and hence must be called by all ranks. It computes the max value across all the DDM ranks. The evaluation of the max value can be restricted to a box or evaluated on a deformed mesh similar to as described in parameter.max().

##### #allmin
def allmin(
*args,
**kwargs
) ‑> List[float]

This is a collective MPI operation and hence must be called by all ranks. It computes the min value across all the DDM ranks. The evaluation of the min value can be restricted to a box or evaluated on a deformed mesh similar to as described in parameter.min().

##### #atbarycenter
def atbarycenter(
self,
physreg: int,
onefield: field
) ‑> quanscient.vec

This outputs a vec object whose structure is based on the field argument onefield and which contains the parameter evaluated at the barycenter of each reference element of physical region physreg. The barycenter of the reference element might not be identical to the barycenter of the actual element in the mesh (for curved elements, for general quadrangles, hexahedra and prisms). The evaluation at barycenter is constant on each mesh element.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x = field("x"); f = field("one")
>>>
>>> # Evaluating the parameter
>>> p = parameter()
>>> p.setvalue(vol, 12*x)
>>> p.write(vol, "parameter.vtk", 1)
>>>
>>>> # Evaluating the same parameter at barycenter
>>> myvec = p.atbarycenter(vol, f)
>>> f.setdata(vol, myvec)
>>> f.write(vol, "barycentervalues.vtk", 1)
##### #countcolumns
def countcolumns(
self
) ‑> int

This returns the number of columns in the parameter.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> E = parameter(2,3)
>>> E.countcolumns()
3
##### #countrows
def countrows(
self
) ‑> int

This returns the number of rows in the parameter.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> E = parameter(2,3)
>>> E.countrows()
2
##### #integrate
def integrate(
*args,
**kwargs
) ‑> float

This integrates a parameter over the physical region physreg. The integration is exact upto the order of polynomials specified in the argument integrationorder.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> myparameter = parameter()
>>> myparameter.setvalue(vol, 12.0)
>>> integralvalue = myparameter.integrate(vol, 4)
>>>
>>> # integrate on a mesh deformed by field 'u'
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> integralvalueondeformedmesh = myparameter.integrate(vol, u, 4)
##### #interpolate
def interpolate(
*args,
**kwargs
) ‑> List[float]

This inteprolates the parameter at a single point whose [x,y,z] coordinate is provided as an argument. The flattened interpolated parameter values are returned if the point was found in the elements of the physical region physreg. If not found an empty list is returned. A parameter can also be interpolated on a deformed mesh by passing its corresponding field.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x=field("x"); y=field("y"); z=field("z")
>>> xyzcoord = [0.5,0.6,0.05]
>>>
>>> p = parameter(3,1)
>>> p.setvalue(vol, array3x1(x,y,z))
>>> interpolated = p.interpolate(vol, xyzcoord)
>>> interpolated
[0.5, 0.6, 0.05]
>>>
>>> # interpolation on the mesh deformed by field 'u'
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> interpolated = p.interpolate(vol, u, xyzcoord)
##### #max
def max(
*args,
**kwargs
) ‑> List[float]

This gives the max value of the parameter over the geometric region physreg. The max value is obtained by splitting all elements refinement times in each direction. Increasing the refinement will thus lead to a more accurate max value, but at an increasing computational cost. The max value is exact when the refinement nodes added to the elements correspond to the position of max. For a first order nodal shape function interpolation, on a mesh that is not curved, the max is always exact to machine precision.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x = field("x")
>>> p = 2*x
>>> maxdata = p.max(vol, 1)
>>> maxdata[0]
2.0

The search of the max value can be restricted to a box delimited by the last argument whose form is [xboxmin,xboxmax, yboxmin, yboxmin, zboxmax, zboxmin]. The output returned is a list of form [maxvalue, xcoordmax, ycoordmax, zcoordmax] or an empty list if the physical region argument is empty or is not in the box provided. If the argument defining the box is not provided, then the whole geometric region is considered for evaluating the max value.

>>> maxdatainbox = p.max(vol, 5, [-2,0, -2,2, -2,2])

The max value can also be evaluated on the geometry deformed by a field (possibly a curved mesh). The max location and the delimiting box are on the undeformed mesh.

>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> maxdataondeformedmesh = p.max(vol, u, 1)

parameter.min()

##### #min
def min(
*args,
**kwargs
) ‑> List[float]

This gives the min value of the parameter over the geometric region physreg. The min value is obtained by splitting all elements refinement times in each direction. Increasing the refinement will thus lead to a more accurate min value, but at an increasing computational cost. The min value is exact when the refinement nodes added to the elements correspond to the position of min. For a first order nodal shape function interpolation, on a mesh that is not curved, the min is always exact to machine precision.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> x = field("x")
>>> p = 2*x
>>> mindata = p.min(vol, 1)
>>> mindata[0]
-2.0

The search of the min value can be restricted to a box delimited by the last argument whose form is [xboxmin,xboxmax, yboxmin, yboxmin, zboxmax, zboxmin]. The output returned is a list of form [maxvalue, xcoordmax, ycoordmax, zcoordmax] or an empty list if the physical region argument is empty or is not in the box provided. If the argument defining the box is not provided, then the whole geometric region is considered for evaluating the min value.

>>> mindatainbox = p.min(vol, 5, [-2,0, -2,2, -2,2])

The min value can also be evaluated on the geometry deformed by a field (possibly a curved mesh). The min location and the delimiting box are on the undeformed mesh.

>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> u.setorder(vol, 1)
>>> mindataondeformedmesh = p.min(vol, u, 1)

parameter.max()

##### #print
def print(
self
) ‑> None

This prints the information on the parameter to the console.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> E = parameter()
>>> E.setvalue(vol, 150e9)
>>> E.print()
##### #setvalue
def setvalue(
self,
physreg: int,
input: expression,
op: str = 'set'
) ‑> None

This sets the parameter ìnput on the physical region physreg.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> E = parameter()
>>> E.setvalue(vol, 150e9)
##### #write
def write(
*args,
**kwargs
) ‑> None

This evaluates a parameter in the physical region physreg and writes it to the file filename. The lagrangeorder is the order of interpolation for evaluation of the parameter.

###### # Examples
>>> # setup
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> u = field("h1xyz")
>>> v = field("h1", [1,2,3])
>>> u.setorder(vol, 1)
>>> v.setorder(vol, 1)
>>>
>>> # interpolation order for writing a parameter
>>> p = parameter()
>>> p.setvalue(1e8*u)
>>> p.write(vol, "uorder1.vtk", 1)    # interpolation order is 1
>>> p.write(vol, "uorder3.vtk", 3)    # interpolation order is 3

### #port

class port(
**kwargs
)

The port object represents a scalar lumped quantity.

#### # Examples

A port object with initial zero value is created as:

>>> V = port()

A mulitharmonic port object wiht initial zero value can be created by passing a list of harmonic numbers. Refer the multiharmonic field constructor for the meaning of the harmonic numbers.

>>> V = port([2,3])

#### # Methods

##### #cos
def cos(
self,
freqindex: int
) ‑> quanscient.port

This gets a port that is the harmonic at freqindex times the fundamental frequency in port V.

###### # Example
>>> V = port([1,2,3,4,5])
>>> Vs = V.cos(0)
>>> Vs.getharmonics()
1

port.sin()

##### #getharmonics
def getharmonics(
self
) ‑> List[int]

This returns the list of harmonics of the port object.

###### # Example
>>> V = port([1,2,3])
>>> harms = V.getharmonics()
>>> harms
[1, 2, 3]
##### #getname
def getname(
self
) ‑> str

This gets the name of the port object.

###### # Example
>>> V = port()
>>> V.setname("LumpedMass")
>>> V.getname()
'LumpedMass'

port.setname()

##### #getvalue
def getvalue(
self
) ‑> float

This returns the value of the port object.

###### # Example
>>> V = port()
>>> V.setvalue(-1.2)
>>> val = V.getvalue()
>>> val
-1.2

port.setvalue()

##### #harmonic
def harmonic(
*args,
**kwargs
) ‑> quanscient.port

This returns a port that is the harmonic/list of harmonics of the port object.

###### # Example
>>> V = port([1,2,3])
>>> V23 = V.harmonic([2,3])
>>> V23.getharmonics()
[2, 3]
##### #print
def print(
self
) ‑> None

This prints the information of the port object.

###### # Example
>>> V = port([2,3])                 # create a multiharmonic port object
>>> V.harmonic(2).setvalue(1)       # set the value of 2nd harmonic
>>> V.harmonic(3).setvalue(0.5)     # set the value of 3rd harmonic
>>> V.print()
Port harmonic 2 has value 1
Port harmonic 3 has value 0.5
##### #setname
def setname(
self,
name: str
) ‑> None

This sets the name for the port object.

###### # Example
>>> V = port()
>>> V.setname("LumpedMass")
>>> V.print()
Port LumpedMass has value 0

port.getname()

##### #setvalue
def setvalue(
self,
portval: float
) ‑> None

This sets the value of the port object.

###### # Examples
>>> V = port()
>>> V.setvalue(10.0)
>>> V.print()
Port has value 10

To set the value of mulitharmonic port:

>>> V = port([2,3])                 # create a multiharmonic port object
>>> V.harmonic(2).setvalue(1)       # set the value of 2nd harmonic
>>> V.harmonic(3).setvalue(0.5)     # set the value of 3rd harmonic
>>> V.print()
Port harmonic 2 has value 1
Port harmonic 3 has value 0.5
##### #sin
def sin(
self,
freqindex: int
) ‑> quanscient.port

This gets a port that is the harmonic at freqindex times the fundamental frequency in port V.

###### # Example
>>> V = port([1,2,3,4,5])
>>> Vs = V.sin(2)
>>> Vs.getharmonics()
4

port.cos()

### #preconditioner

class preconditioner

### #shape

class shape(
**kwargs
)

The shape objects are meshed geometric entities. The mesh created based on shapes can be written in .msh format at any time for visualization in GMSH. It might needed to change the color and visibility options in the menu Tools > Options > Mesh of GMSH.

#### # Examples

Depending on the number and type of arguments, different shape objects can be created for different purposes.

• Creates a shape with the coordinates of all nodes provided as input:
• Example 2: for points and lines: myshape = shape(shapename:str, physreg:int, coords:List[double])
• Creates a shape based on the coordinates of the corner nodes in the shape
• Example 3: for lines and arcs: myshape = shape(shapename:str, physreg:int, coords:List[double], nummeshpts:int)
• Example 4: for triangles and quadrangles: myshape = shape(shapename:str, physreg:int, coords:List[double], nummeshpts:List[int])
• Creates a shape based on sub-shapes provided.
• Example 5: for lines and arcs: myshape = shape(shapename:str, physreg:int, subshapes:List[shape], nummeshpts:int)
• Example 6: for straight-edged triangles and quadrangles: myshape = shape(shapename:str, physreg:int, subshapes:List[shape], nummeshpts:List[int])
• Example 7: for curved triangles and quandrangles. Also, union of several shapes: myshape = shape(shapename:str, physreg:int, subshapes:List[shape])
• Creates a disk shape.
• Example 8: myshape = shape(shapename:str, physreg:int, centercoords:List[double], radius:double, nummeshpts:int)
• Example 9: myshape = shape(shapename:str, physreg:int, centerpoint:shape, radius:double, nummeshpts:int)

Example 1: Creating an empty shape object

>>> myshape = shape()

Example 2: myshape = shape(shapename:str, physreg:int, coords:List[double]). This can be used to create a line going through a list of nodes whose x,y,z coordinates are provided. A physical region number is also provided to have access to the geometric regions of interest in the finite element simulation.

>>> linephysicalregion  = 1
>>> myline = shape("line", linephysicalregion, [0,0,0, 0.5,0.5,0, 1,1,0, 1.5,1,0, 2,1,0])
>>> mymesh = mesh([myline])
>>> mymesh.write("meshed.msh")

If the nodes in the mesh need to be accessed, a point shape object can be created with corresponding nodal coordinates. The nodes can then be accessed through the physical region provided.

>>> pointphysicalregion = 2
>>> point1_coords = [0,0,0]
>>> point2_coords = [2,1,0]
>>> mypoint1 = shape("point", pointphysicalregion, point1_coords)
>>> mypoint2 = shape("point", pointphysicalregion, point2_coords)
>>> # The point 1 and 2 are now available in the pointphysicalregion=2.
>>>
>>> p = field("h1")
>>> p.setorder(linephysicalregion, 1)
>>> p.setconstraint(pointphysicalregion, 2) # Dirichlet boundary constraint will be applied on point 1 and 2

Example 3: myshape = shape(shapename:str, physreg:int, coords:List[double], nummeshpts:int) This can be used to create: a straight line between the first (x1,y1,z1) and last point (x2,y2,z2) provided. a circualr arc between the first (x1,y1,z1) and second point (x2,y2,z2) whose center is the third point (x3,y3,z3). The nummeshpts argument corresponds to the number of nodes in the meshed shape. At least two nodes are expected.

>>> linephysicalregion=1; arcphysicalregion=1
>>> myline = shape("line", linephysicalregion, [0,0,0, 1,-1,1], 10)     # creates a line mesh with 10 nodes
>>> myarc = shape("arc", arcphysicalregion, [1,0,0, 0,1,0, 0,0,0], 8)   # creates an arc mesh with 8 nodes
>>> mymesh = mesh([myline, myarc])
>>> mymesh.write("meshed.msh")

Example 4: myshape = shape(shapename:str, physreg:int, coords:List[double], nummeshpts:List[int]) This can be used to create: a straight-edge quadrangle with a full quandrangle structured mesh. a straight-edge triangle with s structured mesh made of triangles along the edge linking the second and third node and quadrangles everywhere else.

>>> mytriangle = shape("triangle", trianglephysicalregion, [1,0,0, 2,0,0, 1,1,0], [10,10,10])
>>> mymesh.write("meshed.msh")
The <code>coords</code> arguments provides the <code>x,y,z</code>coordinates of the corner nodes. E.g. (0,0,0), (1,0,0), (1,1,0) and (-0.5,1,0) for the quadrangle.
The <code>nummeshpts</code> argument specifies the number of nodes to mesh each of the contour lines. At least two nodes are expected for each contour line.
All contour lines must have the same number of nodes for the triangle shape while for the quadrangle shape the contour lines facing each other
must have the same number of nodes.

Example 5: myshape = shape(shapename:str, physreg:int, subshapes:List[shape], nummeshpts:int) This can be used to create the following shapes from the list of subshapes provided: a straight line between the first (x1,y1,z1) and last point (x2,y2,z2) provided. a cricle arc between the first (x1,y1,z1) and second point (x2,y2,z2) whose center is the third point. The nummeshpts argument corresponds to the number of nodes in the meshed shape.

>>> # Point subshapes
>>> point1 = shape("point", -1, [0,0,0])
>>> point2 = shape("point", -1, [1,0,0])
>>> point3 = shape("point", -1, [0,1,0])
>>> point4 = shape("point", -1, [1,-1,1])
>>>>
>>> # Creating line and arc shapes from point subshapes
>>> linephysicalregion=1; arcphysicalregion=2
>>> myline = shape("line", linephysicalregion, [point1, point4], 10)
>>> myarc = shape("arc", arcphysicalregion, [point2, point3, point1], 8)
>>> mymesh = mesh([myline, myarc])
>>> mymesh.write("meshed.msh")

Example 6: myshape = shape(shapename:str, physreg:int, subshapes:List[shape], nummeshpts:List[int]) This can be used to create the following shapes from the list of subshpaes provided: a straight-edge quadrangle with a full quadrangle structured mesh a straight-edge triangle with a structured mesh made of triangles along the edge linking the second and third node and quadrangles everywhere. The subshapes argument provides the list of corner point shapes. The nummeshptsargument gives the number of nodes to mesh each of the contour lines. At least two nodes are expected for each contour line. All contour lines must have the same number of nodes for the triangle shape while for quadrangle shape the contour lines facing each other must have the same number of nodes.

>>> # Point subshapes
>>> point1 = shape("point", -1, [0,0,0])
>>> point2 = shape("point", -1, [1,0,0])
>>> point3 = shape("point", -1, [1,1,0])
>>> point4 = shape("point", -1, [0,1,0])
>>> point5 = shape("point", -1, [2,0,0])
>>>
>>> # Creating traingle and quadrangle shape from subshapes of corner points
>>> mytraingle = shape("triangle", trianglephysicalregion, [point2, point5, point3], [8,8,8])
>>> mymesh.write("meshed.msh")

Example 7: myshape = shape(shapename:str, physreg:int, subshapes:List[shape]) This can be used to create: a curved quardrangle with full quadrangle structured mesh. a curved triangle with structured mesh made of triangles along the edge linking the second and third node and quandrangles everywhere. a shape that is the union of several shapes of the same dimension. The subshapes argument provides the contour shapes (clockwise or anti-clockwise). All contour lines must have the same number of nodes for triangle shape while for quadrangle shape the contour lines facing each other must have the same number of nodes.

>>> # Creating subshapes
>>> line1 = shape("line", -1, [-1,-1,0, 1,-1,0], 10)
>>> arc2 = shape("arc", -1, [1,-1,0, 1,1,0, 0,0,0], 12)
>>> line3 = shape("line", -1, [1,1,0, -1,1,0], 10)
>>> line4 = shape("line", -1, [-1,1,0, -1,-1,0], 12)
>>> line5 = shape("line", -1, [1,-1,0, 3,-1,0], 12)
>>> arc6 = shape("arc", -1, [3,-1,0, 1,1,0, 1.6,-0.4,0], 12)
>>>
>>> mytriangle = shape("triangle", trianglephysicalregion,[line5, arc6, arc2])
>>> myunion = shape("union", unionphysicalregion, [line1, arc2, line3, line4])
>>>
>>> mymesh = mesh([myquadrangle, mytraingle, myunion])
>>> mymesh.write("meshed.msh")

Example 8: myshape = shape(shapename:str, physreg:int, centercoords:List[double], radius:double, nummeshpts:int) This is used to create a 2D disk with structured mesh centered around centercoords. The nummeshptsargument corresponds to the number of nodes in the contour circle of the disk. Since the disk has a structured mesh, the number of mesh nodes must be a multiple of 4. The radiusarguments provides the radius of the disk.

>>> diskphysicalregion=1
>>> mydisk = shape("disk", diskphysicalregion, [1,0,0], 2, 40)
>>> mymesh = mesh([mydisk])
>>> mymesh.write("meshed.msh")

Example 9: myshape = shape(shapename:str, physreg:int, centerpoint:shape, radius:double, nummeshpts:int) This is used to create a 2D disk with structured mesh centered around point shape centerpoint. The nummeshpts argument corresponds to the number of nodes in the contour circle of the disk. Since the disk has a structured mesh, the number of mesh nodes must be a multiple of 4. The radiusarguments provides the radius of the disk.

>>> diskphysicalregion=1
>>> centerpoint = shape("point", -1, [1,0,0])
>>> mydisk = shape("disk", diskphysicalregion, centerpoint, 2, 40)
>>> mymesh = mesh([mydisk])
>>> mymesh.write("meshed.msh")

mesh

#### # Methods

##### #duplicate
def duplicate(
self
) ‑> quanscient.shape

This outputs a shape that is duplicate of the initial shape. All subshapes are duplicated recursively as well but the object equality relations between subshapes are identical between a shape and its duplicate.

##### #extrude
def extrude(
*args,
**kwargs
) ‑> quanscient.shape

A given shape is extruded in the requested direction (z by defalt) to form a higher dimensional shape. The extrude function works for 0D, 1D and 2D shapes.

###### # Parameters

physreg : int : Physical region to which the extruded shape is set.

height : double : Height of extrusion in the requested direction.

numlayers : int : Value specifying the number of node layers the extruded mesh should contain.

extrudedirection : List[int] : Unit vector specifying the direction of extrusion.

###### # Examples

Example 1: myshape = shape.extrude(physreg:int, height:double, numlayers:int, extrudedirection:List[double])

>>> volumephysicalregion = 100
>>> myvolume = myquadrangle.extrude(volumephysicalregion, 1.4, 6, [0,0,1])
>>> mymesh = mesh([myvolume])
>>> mymesh.write("meshed.msh")

Example 2: myshape = shape.extrude(physreg:List[int], height:List[double], numlayers:[int], extrudedirection:List[double]). This extends the extrude function to multiblock extrusion.

>>> mytriangle = shape("triangle", 1, [0,0,0, 1,0,0, 0,1,0], [6,6,6])
>>>
>>> '''
>>> Creating multiblock extrusion:
>>> ---------|----------|----------|-----------
>>> block    | physreg  |  height  |  numlayers
>>> ---------|----------|----------|-----------
>>> Block 1: |    11    |   0.5    |    3
>>> Block 2: |    12    |   0.3    |    5
>>> ---------|----------|----------|-----------
>>> In block 1, the initial shape is extruded to a height of 0.5 and contains 3 node layers and the extruded shape is set to physical region 11.
>>> In block 2, the initial shape is extruded to a height of 0.3 (starting from height 0.5) and contains 5 node layers and the extruded shape is set to physical region 12.
>>> '''
>>> myvolumes = mytriangle.extrude([11,12], [0.5,0.3], [3,5], [0,0,1])    # creates two extruded shapes
>>> mymesh = mesh(myvolumes)
>>> mymesh.write("meshed.msh")
##### #getcoords
def getcoords(
self
) ‑> List[float]

This returns the coordinates of all nodes in the shape mesh.

###### # Examples
[0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0]
##### #getcurvatureorder
def getcurvatureorder(
self
) ‑> int

This returns the curvature order of a given shape.

###### # Example
>>> q = shape("quadrangle", 1, [0,0,0, 1,0,0, 1,1,0, 0,1,0], [2,2,2,2])
>>> q.getcurvatureorder()
1
##### #getdimension
def getdimension(
self
) ‑> int

This gives the shape dimension (0D, 1D, 2D or 3D).

###### # Examples
>>> mypoint = shape("point", 111, [0,0,0])
>>> mypoint.getdimension()
0
>>>
>>> myline = shape("line", 222, [0,0,0, 0.5,0.5,0, 1,1,0, 1.5,1,0, 2,1,0])
>>> myline.getdimension()
1
>>>
>>> myarc = shape("arc", 333, [1,0,0, 0,1,0, 0,0,0], 8)
>>> myarc.getdimension()
1
>>>
>>> mytriangle = shape("triangle", 444, [1,0,0, 2,0,0, 1,1,0], [10,10,10])
>>> mytriangle.getdimension()
2
>>>
2
>>>
>>> mylines[0].getdimension()
1
##### #getname
def getname(
self
) ‑> str

This returns the name of the shape.

###### # Examples
>>> mypoint = shape("point", 111, [0,0,0])
>>> mypoint.getname()
'point'
>>>
>>> myline = shape("line", 222, [0,0,0, 0.5,0.5,0, 1,1,0, 1.5,1,0, 2,1,0])
>>> myline.getname()
'line'
>>>
>>> myarc = shape("arc", 333, [1,0,0, 0,1,0, 0,0,0], 8)
>>> myarc.getname()
'arc'
>>>
>>> mytriangle = shape("triangle", 444, [1,0,0, 2,0,0, 1,1,0], [10,10,10])
>>> mytriangle.getname()
'triangle'
>>>
>>>
>>> mylines[0].getname()
'line'
##### #getphysicalregion
def getphysicalregion(
self
) ‑> int

This gives the physical region number for a given shape. The physical region is used in the finite element simulation to identify a region.

###### # Returns

int: : Returns -1 if the physical region was not set, else a corresponding positive integer.

###### # Examples
>>>
>>> myline1 = mylines[0]
>>> myline1.setphysicalregion(2)
>>> myline1.getphysicalregion()
2
>>>
>>> myline2 = mylines[1]
>>> myline1.getphysicalregion()
-1
##### #getsons
def getsons(
self
) ‑> List[quanscient.shape]

This returns a list containing the direct subshapes of the shape object. For a quadrangle, its 4 contour lines are returned. For a triangle, its 3 contour lines are returned.

###### # Example
>>> for myline in mylines:
...     myline.setphysicalregion(2)
...
>>> mymesh = mesh(mylines)
>>> mymesh.write("meshed.msh")
##### #move
def move(
self,
u: expression
) ‑> None

This moves the shape (and all its subshapes recursively) in the x,y and z direction by a value provided in the 3x1 expression array. When moving multiple shapes that share common subshapes, ensure that subshapes are not moved multiple times.

###### # Parameter

u: expression 3x1 array expression that specifies that values by which shape is moved in x,y,z direction

###### # Example
>>> x=field("x"); y=field("y"); z=field("z")
>>> mymesh.write("meshed.msh")
##### #rotate
def rotate(
self,
alphax: float,
alphay: float,
alphaz: float
) ‑> None

This rotates the shape (and all its subshapes recursively) first by alphax degress around the x-axis, then alphay degress around the y-axis and finally alphaz degress around the z-axis. When rotating multiple shapes that share common subshapes make sure that subshapes are not rotated multiple times.

###### # Parameters

alphax : double : Degrees by which the shape is rotated around x-axis.

alphay : double : Degrees by which the shape is rotated around y-axis.

alphaz : double : Degrees by which the shape is rotated around z-axis.

###### # Example
>>> mymesh.write("meshed.msh")
##### #scale
def scale(
self,
scalex: float,
scaley: float,
scalez: float
) ‑> None

This scales the shape (and all its subspaces recursively) in the x,y and z direction by a give factor. Factor of 1 keeps the shape unchanged. When scaling multiple shapes that share common subshapes make sure the subshapes are not scaled multiple times.

###### # Parameters

scalex : double : Value by which the shape is scaled in x-direction.

scaley : double : Value by which the shape is scaled in y-direction.

scalez : double : Value by which the shape is scaled in z-direction.

###### # Example
>>> mymesh.write("meshed.msh")
##### #setphysicalregion
def setphysicalregion(
self,
physreg: int
) ‑> None

This sets the physical region number for a given shape. Subshapes are not affected. The physical region is used in the finite element simulation to identify a region.

###### # Parameters

physreg : int : An integer identifying the physical region.

###### # Example
>>> myline = myquadrangle.getsons()[0]  # the shape 'myline' is not associated with any physical region yet.
>>> myline.setphysicalregion(linephysicalregion)
>>> mymesh.write("meshed.msh")
##### #shift
def shift(
self,
shiftx: float,
shifty: float,
shiftz: float
) ‑> None

This shifts the shape (and all its subshapes recursively) in the x,y, and z direction by a double value. When shifting multiple shapes that share common subshapes make sure the subspaces are not shifted multiple times.

###### # Parameters

shiftx : double : Value by which the shape is moved in x-direction.

shifty : double : Value by which the shape is moved in y-direction.

shiftz : double : Value by which the shape is moved in x-direction.

###### # Example
>>> mymesh.write("meshed.msh")

### #spanningtree

class spanningtree(
physregs: List[int]
)

The spanningtree object holds a spanning tree whose edges go through all nodes in the mesh without forming a loop. The tree growth starts on the physical regions provided as argument.

#### # Parameters

physregs : List[int] : List of physical regions where the spanning tree is first fully grown before being extended everywhere.

#### # Example

A spanning tree object is created by passing the physical regions 'sur' and 'top'. Hence, here the tree is first fully grown on face regions 'sur' and 'top' before extending everywhere.

>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3;    # physical regions
>>> spantree = spannintree([sur, top])

#### # Methods

##### #countedgesintree
def countedgesintree(
self
) ‑> int

This returns the number of edges in the spanning tree.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3;    # physical regions
>>> spantree = spannintree([sur, top])
>>> spantree.countedgesintree()
1859
##### #write
def write(
self,
filename: str
) ‑> None

This writes the tree to a file for visualization.

###### # Parameters

filename : str : Name of the file to which the data samples are written.

###### # Example
>>> mymesh = mesh("disk.msh")
>>> vol=1; sur=2; top=3;    # physical regions
>>> spantree = spannintree([sur, top])
>>> spantree.write("spantree.vtk")

### #spline

class spline(
**kwargs
)

The spline object allows interpolation in a discrete data set using cubic (natural) splines.

#### # Raises

RuntimeError : If a mesh object is not available prior to creating a spline object.

#### # Notes

Before creating a spline object, ensure that a mesh object is available. If a mesh object is not already available, create an empty mesh object.

#### # Examples

Say the data samples are in a text file, with each data separated by "," as shown below:

>>> 273,5e9,300,4e9,320,2.5e9,340,1e9

A spline object can then be created by reading the the x-y data contained in the text file.

>>> mymesh = mesh()     # if a mesh object is not already available.
>>> spl1 = spline(filename="measured.txt", delimiter="\n")

The x-y data samples can also be provided in two separate lists or tuples. In that case, a spline object is created as follows:

>>> mymesh = mesh()     # if a mesh object is not already available.
>>> temperature = (273, 300, 320, 340)
>>> youngsmodulus = (5e9, 4e9, 2.5e9, 1e9)
>>> spl2 = spline(xin=temperature, yin=youngsmodulus)

#### # Notes

The ordering of the samples provided does not matter. Internally they are always sorted in the ascending order of data.

#### # Methods

##### #evalat
def evalat(
*args,
**kwargs
) ‑> quanscient.densemat

Interpolates at given point(s) that falls within the original data range provided.

###### # Parameters

input : double or List[double] or Tuple(double) : An point or a list/tuple of points for interpolation.

###### # Returns

double : a interpolated scalar value if the input was a scalar.

List[double] : a list of interpolated values if the input was a list/tuple.

###### # Raises

RuntimeError : If at least one of the input for interpolation is outside the range of the original data samples provided.

###### # Examples
>>> mymesh = mesh()     # if a mesh object is not already available.
>>>
>>> temperature = (320, 273, 340, 300)          # original x input
>>> youngsmodulus = (2.5e9, 5e9, 1e9, 4e9)      # original y input
>>> spl = spline(temperature, youngsmodulus)

Example 1:

>>> spl.evalat(303)
4199390995.630462

Example 2:

>>> spl.evalat([298,304,275])
[3912277862.548358, 4278701622.971286, 4835693423.344276]

Example 3:

>>> spl.evalat(250)
RuntimeError: Error in 'spline' object: data requested in interval (250,250) is out of the provided data range (273,340)

Example 4:

>>> spl.evalat([290, 310, 400])
RuntimeError: Error in 'spline' object: data requested in interval (290,400) is out of the provided data range (273,340)
##### #getderivative
def getderivative(
self
) ‑> quanscient.spline

This returns the derivative of the spline.

The spline polynomial and its derivative :

###### # Returns

spline : spline object that holds the data samples corresponding to the derivative of the spline.

###### # Example
>>> mymesh = mesh()     # if a mesh object is not already available.
>>>
>>> temperature = (273, 300, 320, 340)
>>> youngsmodulus = (5e9, 4e9, 2.5e9, 1e9)
>>>
>>> spl = spline(temperature, youngsmodulus)
>>> spl.write("spline_data.txt", 0, ",")
273,5000000000,300,4000000000,320,5000000000,340,1000000000
>>>
>>> dspl = spl.getderivative()
>>> dspl.write("splinederivative_data.txt", 0, ",")
273,-82402205.576362908,300,53693300.041614674,320,-58198085.726175644,340,-270900957.13691229
##### #getxmax
def getxmax(
self
) ‑> float

This returns the maximum value that input takes in the original data provided.

###### # Example
>>> mymesh = mesh()     # if a mesh object is not already available.
>>>
>>> temperature = (320, 273, 340, 300)          # original x input
>>> youngsmodulus = (2.5e9, 5e9, 1e9, 4e9)      # original y input
>>> spl = spline(temperature, youngsmodulus)
>>> spl.getxmax()
340

spline.getxmin()

##### #getxmin
def getxmin(
self
) ‑> float

This returns the minimum value that input takes in the original data provided.

###### # Example
>>> mymesh = mesh()     # if a mesh object is not already available.
>>>
>>> temperature = (320, 273, 340, 300)          # original x input
>>> youngsmodulus = (2.5e9, 5e9, 1e9, 4e9)      # original y input
>>> spl = spline(temperature, youngsmodulus)
>>> spl.getxmin()
273

spline.getxmax()

##### #set
def set(
self,
xin: List[float],
yin: List[float]
) ‑> None
##### #write
def write(
self,
filename: str,
numsplits: int,
delimiter: str = '\n'
) ‑> None

This writes to file a refined version of the original data samples with data sorted in the ascending order. It can be used to visualize the interpolation obtained with cubic splines.

###### # Parameters

filename : str : Name of the file to which the data samples are written.

numsplits : int : The number of additional points between two successive input considered for evaluation and subsequent writing. Minimum value required is .

delimiter : str, optional : String specifying the delimiter. The default value is "\n".

###### # Raises

RuntimeError : If numsplits < 0. Minimum value required is .

###### # Examples

Create a spline object:

>>> mymesh = mesh()     # if a mesh object is not already available.
>>>
>>> temperature = (320, 273, 340, 300)          # original x input
>>> youngsmodulus = (2.5e9, 5e9, 1e9, 4e9)      # original y input
>>> spl = spline(temperature, youngsmodulus)

Example 1: If numsplits = 0, no additional points are considered and the original data samples are written.

>>> numsplits = 0
>>> spl.write("spline_data.txt", numsplits, ",")
273,5000000000,300,4000000000,320,5000000000,340,1000000000

Example 2: If numsplits = 1, between two successive inputs, one additional point is considered for evaluation and subsequent writing.

>>> numsplits = 1
>>> spl.write("spline_data.txt", numsplits, "\n")
273
5000000000
286.5
4040677668.5393257
300
4000000000
310
4779728464.4194756
320
5000000000
330
3531757178.5268416
340
1000000000

Similarly, if numsplits = 2, between two successive inputs,two additional points are considered for evaluation and subsequent writing.

class universe

Type: list

mnt: int
) ‑> None

### #vec

class vec(
**kwargs
)

The vec object holds a vector, be it the solution vector of an algebraic problem or its right handside.

#### # Examples

Example 1: vec(formul:formulation) This creates an all zero vector whose structure and size is the one of formulation projection.

>>> 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))
>>>
>>> b = vec(projection)

Example 2: vec(vecsize:int, addresses:indexmat, vals:densemat) This creates a vector with given values at given addresses.

>>> allinitialize()
>>> vals = densemat(3,1, [5,10,20])
>>>
>>> b = vec(3, addresses, vals)
>>> b.print()

Vec Object: 1 MPI processes type: seq 5. 10. 20.

>>> allfinalize()

#### # Methods

##### #copy
def copy(
self
) ‑> quanscient.vec

This creates a full copy of the vector object.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> copiedvec = sol.copy()
##### #getallvalues
def getallvalues(
self
) ‑> quanscient.densemat

This gets the values of all the entries of the vector in the sequential order.

###### # Returns

densemat : A column matrix with number of rows equal to the vector object size.

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                 # creates a zero vector
>>>
>>> vals = densemat(myvec.size(),1, 12)
>>> myvec.setallvalues(vals)                # all the entries now contain a value of 12
>>>
>>> vecvals = myvec.getallvalues()
##### #getvalue
def getvalue(
*args,
**kwargs
) ‑> float

This gets the value of the vector object at the given address. The address provides the index at which the entry is requested.

###### # Parameters

address : int : Index at which the entry of a vector is requested.

###### # Returns

densemat : A matrix of size .

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                 # creates a zero vector
>>>
>>> vals = densemat(myvec.size(),1, 12)
>>> myvec.setallvalues(vals)                # all the entries now contain a value of 12
>>>
>>> vecvals = myvec.getvalue(2)
##### #getvalues
def getvalues(
self,
) ‑> quanscient.densemat

This gets the values in the vector that are at the indices given in addresses.

###### # Parameters

addresses : indexmat : A column matrix storing the indices at which the entries of a vector is requested.

###### # Returns

densemat : A column matrix with number of rows equal to the length of the addresses.

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                 # creates a zero vector
>>>
>>> vals = densemat(myvec.size(),1, 12)
>>> myvec.setallvalues(vals)                # all the entries now contain a value of 12
>>>
self,
filename: str
) ‑> None

This loads the data of a vector object from a file (either .bin or .txt ASCII format). This only works correctly if the dof structure of the calling vector object is exactly same as that of the vector object from which the file was written to the disk. In other words, the same set of formulation contributions must be defined in the same order.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>> sol.write("vecdata.bin")        # writes the vector data to a file
>>>

vec.write()

##### #noautomaticupdate
def noautomaticupdate(
self
) ‑> None

After this call the vector object will not have its value automatically updated after hp-adaptivity. If the automatic update is not needed then this call is recommended to avoid a possibly costly update to the vector values.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> sol.noautomaticupdate()
##### #norm
def norm(
self,
type: str = '2'
) ‑> float

This returns the , or norm of the vector. Default is the norm.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> vals = densemat(sol.size(),1, 12)
>>> sol.setallvalues(vals)                # all the entries now contain a value of 12
>>>
>>> sol.norm()
108.6646

vec.sum()

##### #permute
def permute(
self,
rowpermute: indexmat,
invertit: bool = False
) ‑> None

This rearranges the vector in the order of indices prescribed in rowpermute. The inverse permutation is performed if the boolean flag invertit is set to True. The rowpermute decribes the mapping or inverse mapping function.

###### # Parameters

rowpermute : indexmat : A column matrix storing the order of vector indices for permutation.

invertit : bool, default: False : If set to True, inverse permutation is performed.

###### # Example
>>> rows = indexmat(6,1, [0,1,2,3,4,5])
>>> vals = densemat(6,1, [00,10,20,30,40,50])
>>>
>>> v = vec(6, rows, vals)
>>> permuterows = indexmat(6,1, [3,1,4,5,0,2])
>>> v.permute(permuterows, invertit=False)
>>> v.print()
Vec Object: 1 MPI processes
type: seq
30.
10.
40.
50.
0.
20.
>>>
>>> # Inverting the permutation on the above will give back the original order of the vector.
>>> v = vec(6, rows, vals)
>>> permuterows = indexmat(6,1, [3,1,4,5,0,2])
>>> v.permute(permuterows, invertit=True)
>>> v.print()
Vec Object: 1 MPI processes
type: seq
0.
10.
20.
30.
40.
50.
##### #print
def print(
self
) ‑> None

This prints the values of the vector object.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> sol.print()
##### #setallvalues
def setallvalues(
self,
valsmat: densemat,
op: str = 'set'
) ‑> None

This replaces all the entries in the vector object by the values in valsmat. The addresses of valsmat are assumed to be in sequential order. If op='set', the values are replaced and if op='add' the values are instead added to existing ones. This method works on all the enrties.

###### # Parameters

valsmat : densemat : A column matrix storing the values that are replaced in or added to the vector object.

op : str, default: 'set' : Equal to 'set' if the values must replaced. For adding values to exisitng ones use 'add' instead.

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                     # creates a zero vector
>>>
>>> vals = densemat(myvec.size(),1, 12)
>>> myvec.setallvalues(vals)
##### #setdata
def setdata(
*args,
**kwargs
) ‑> None

This sets to the vector the data from the fields and ports defined in the formulation.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)   # creates a zero vectorgiven op: str)`
This transfers the data of a field <code>myfield</code> to the addresses in vector object corresponding to the physical region <code>physreg</code>.
>>> sol.setdata(vol, v);    # populates the vector with the data from the field v

Example 2: vec.setdata() This transfers to the vector the data from all fields and ports defined in the associated formulation.

>>> ...
>>> sol.setdata();
##### #setvalue
def setvalue(
*args,
**kwargs
) ‑> None

This replaces the value in the vector object at the given address by the values in value. The 'address' provides the index at which the entry is replaced by value. If op='set', the value is replaced and if op='add' the value is added to existing one. This method works only on a given single entry.

###### # Parameters

address : int : Index at which the entry of a vector is replaced/added.

value : double : Value that is set/added in the vector object.

op : str, default: 'set' : Equal to 'set' if the value must replaced. For adding the value to exisitng one use 'add' instead.

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                     # creates a zero vector
>>>
>>> myvec.setvalue(2, 2.32)
##### #setvalues
def setvalues(
self,
valsmat: densemat,
op: str = 'set'
) ‑> None

This replaces the values in the vector object at the given addresses by the values in valsmat. If op='set', the values are replaced and if op='add' the values are instead added to existing ones. This method works only on entries given in the addresses.

###### # Parameters

addresses : indexmat : A column matrix storing the indices at which the entries of a vector are replaced/added.

valsmat : densemat : A column matrix storing the values that are replaced in or added to the vector object.

op : str, default: 'set' : Equal to 'set' if the values must replaced. For adding values to exisitng ones use 'add' instead.

###### # Example
>>> 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))
>>>
>>> myvec = vec(projection)                     # creates a zero vector
>>>
>>> vals = densemat(myvec.size(),1, 12)
>>>
>>> myvec.setvalues(addresses, vals)            # op is the default 'set'. All entries are replaced by value 12.
>>> myvec.setvalues(addresses, vals, 'set')     # All entries are replaced by value 12.
##### #size
def size(
self
) ‑> int

This returns the size of the vector object. If the vector was instatiated from a formulation, then the vector size is equal to the number of dofs in that formulation.

###### # Example
>>> 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))
>>>
>>> b = vec(projection)
>>> b.size()
82
##### #sum
def sum(
self
) ‑> float

This returns the sum of all the values in the vector.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> vals = densemat(sol.size(),1, 12)
>>> sol.setallvalues(vals)                # all the entries now contain a value of 12
>>>
>>> sol.sum()
984.0

vec.norm()

##### #updateconstraints
def updateconstraints(
self
) ‑> None

This updates the values of all Dirichlet constraint entries in the vector.

###### # Example
>>> 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))
>>>
>>> b = vec(projection)
>>>
>>> v.setconstraint(sur, 1)
>>> b.updateconstraints()
##### #write
def write(
self,
filename: str
) ‑> None

This writes all the data in the vector object to disk in a lossless and compact form. The file can be written in binary .bin format (extremely compact but less portable) or in ASCII .txt format (portable).

###### # Parameters

filename : str : The name of the file to which the data from the vector object is written.

###### # Examples
>>> mymesh = mesh("disk.msh")
>>> vol = 1
>>> v = field("h1")
>>> v.setorder(vol, 1)
>>> projection = formulation()
>>> projection += integral(vol, dof(v)*tf(v) - 2*tf(v))
>>> sol = vec(projection)
>>>
>>> sol.write("vecdata.txt")    # writes the vector data to a file

### #vectorfieldselect

class vectorfieldselect

#### # Methods

##### #setdata
def setdata(
self,
physreg: int,
myfield: field,
op: str = 'set'
) ‑> None

### #wallclock

class wallclock

This initializes the wall clock object.

#### # Methods

##### #pause
def pause(
self
) ‑> None

This pauses the clock. The wallclock.pause() and wallclock.resume() functions allow to time selected operations in loop.

###### # Example
>>> myclock = wallclock()
>>> myclock.pause()
>>> # Do something
>>> myclock.resume()
>>> myclock.print()

wallclock.resume()

##### #print
def print(
self,
toprint: str = ''
) ‑> None

This prints the time elapsed in the most appropriate format (, , or ). It also prints the message passed in the argument toprint (if any).

###### # Parameters

toprint : str, optional : message to print. The default value is an empty string "".

###### # Example
>>> myclock = wallclock()
>>> myclock.print("Time elapsed")   # or myclock.print()
##### #resume
def resume(
self
) ‑> None

This resumes the clock. The wallclock.pause() and wallclock.resume() functions allow to time selected operations in loop.

###### # Example
>>> myclock = wallclock()
>>> myclock.pause()
>>>
>>> for i in range (0, 10):
>>>     myclock.resume()
>>>     # Do something and time it
>>>     myclock.pause()
>>>     # Do something else
>>> myclock.print()

wallclock.pause()

##### #tic
def tic(
self
) ‑> None

This resets the clock.

###### # Example
>>> myclock = wallclock()
>>> myclock.tic()

wallclock.toc()

##### #toc
def toc(
self
) ‑> float

This returns the time elapsed (in ).

###### # Example
>>> myclock = wallclock()
>>> timeelapsed = myclock.toc()