Polar splines

Basic modules

class struphy.polar.basic.PolarDerhamSpace(derham, space_id)

Bases: VectorSpace

Derham space with polar basis in eta1-eta2.

Parameters:

• derham (struphy.feec.psydac_derham.Derham) – Discrete Derham complex.

• space_id (str) – Space identifier for the field (H1, Hcurl, Hdiv, L2 or H1vec).

property dtype

TODO

property comm

TODO

property space_id

TODO

property n_comps

TODO

property n

TODO

property d

TODO

property parent_space

The parent space (StencilVectorSpace or BlockVectorSpace) of which the PolarDerhamSpace is a sub-space.

property starts

TODO

property ends

TODO

property dimension

TODO

property n_polar

Number of polar basis functions in each component (tuple for vector-valued).

property n_rings

Number of rings to be set to zero in tensor-product basis (tuple for vector-valued).

property n2

Tuple holding total (global) number of basis function in eta2, for each component.

property n3

Tuple holding total (global) number of basis function in eta3, for each component.

property type_of_basis_3

Tuple holding type of spline basis (B-splines or M-splines), for each component.

property parallel

zeros()

Creates an element of the vector space filled with zeros.

axpy(a, x, y)

Increment the vector y with the a-scaled vector x, i.e. y = a * x + y, provided that x and y belong to the same vector space V (self). The scalar value a may be real or complex, depending on the field of V.

Parameters:

• a (scalar) – The scaling coefficient needed for the operation.

• x (Vector) – The vector which is not modified by this function.

• y (Vector) – The vector modified by this function (incremented by a * x).

inner(x, y)

Evaluate the inner vector product between two vectors of this space V.

If the field of V is real, compute the classical scalar product. If the field of V is complex, compute the classical sesquilinear product with linearity on the second vector.

TODO [YG 01.05.2025]: Currently, the first vector is conjugated. We want to reverse this behavior in order to align with the convention of FEniCS.

Parameters:

• x (Vector) – The first vector in the scalar product. In the case of a complex field, the inner product is antilinear w.r.t. this vector (hence this vector is conjugated).

• y (Vector) – The second vector in the scalar product. The inner product is linear w.r.t. this vector.

Returns:

The scalar product of the two vectors. Note that inner(x, x) is a non-negative real number which is zero if and only if x = 0.

Return type:

float | complex

class struphy.polar.basic.PolarVector(V)

Bases: Vector

Element of a PolarDerhamSpace.

An instance of a PolarVector consists of two parts:

1. a list of xp.arrays of the polar coeffs (not distributed)

2. a tensor product StencilVector/BlockVector of the parent space with inner rings set to zero (distributed).

Parameters:

V (PolarDerhamSpace) – Vector space which the polar vector to be created belongs to.

property space

TODO

property dtype

TODO

property pol

Polar coefficients as xp.array.

property tp

Tensor product Stencil-/BlockVector with inner rings set to zero.

property ghost_regions_in_sync

Whether ghost regions of tensor product part are up-to-date.

dot(v)

Scalar product with another instance of PolarVector.

set_vector(v)

In-place setter for polar + tensor product coeffiecients.

set_tp_coeffs_to_zero()

Sets inner tensor-product rings that make up the polar splines to zero.

toarray(allreduce=False)

Converts the polar vector to a 1d numpy array.

toarray_tp()

Converts the Stencil-/BlockVector to a 1d numpy array but NOT the polar part.

copy(out=None)

TODO

update_ghost_regions()

Update ghost regions before performing non-local access to vector elements (e.g. in matrix-vector product).

Parameters:

direction (int) – Single direction along which to operate (if not specified, all of them).

exchange_assembly_data()

Exchange assembly data before performing non-local access to vector elements (e.g. in matrix-vector product).

Parameters:

direction (int) – Single direction along which to operate (if not specified, all of them).

conjugate()

No need for complex conjugate

struphy.polar.basic.set_tp_rings_to_zero(v, n_rings)

Sets a certain number of rings of a Stencil-/BlockVector in eta_1 direction to zero.

Parameters:

• v (StencilVector | BlockVector) – The vector whose inner rings shall be set to zero.

• n_rings (tuple) – The number of rings that shall be set to zero (has length 1 for StencilVector and 3 for BlockVector).

Extraction operators

class struphy.polar.extraction_operators.PolarExtractionBlocksC1(domain, derham)

Bases: object

Class for operator blocks defining C1 basis extraction operators and DOF extraction operators.

/ / / / /

k = 1 / / / / /

k = 0 / / / / /

——————————- —————–

j = 0 | | | | | j = 1 | | | | |

| | | |

i = 0 | … | i = n_rings | i > n_rings |

| | | | /

| | | | /

j = n2-1 | | | | | /

polar rings | first tp ring | outer tp zone |

Fig.Indexing of 3d spline tensor-product (tp) basis functions/DOFs.

A linear combination of n_rings rings is used to construct n_polar polar basis functions/DOFs and “first tp ring” basis functions/DOFs.

Parameters:

• cx (array-like) – Control points defining the x-component of a 2D B-spline mapping.

• cy (array-like) – Control points defining the y-component of a 2D B-spline mapping.

property cx

property cy

property pole

property n0

property n1

property n2

property d0

property d1

property d2

property n_rings

property n_polar

property tau

property xi_0

property xi_1

property e_ten_to_pol

property p_ten_to_pol

property p_ten_to_ten

property grad_pol_to_pol

property grad_pol_to_ten

property grad_e3

property curl_pol_to_pol

property curl_pol_to_ten

property curl_e3

property div_pol_to_pol

property div_pol_to_ten

property div_e3

class struphy.polar.extraction_operators.PolarSplines_C0_2D(n0, n1)

Bases: object

class struphy.polar.extraction_operators.PolarSplines_C1_2D(cx, cy)

Bases: object

class struphy.polar.extraction_operators.PolarSplines(tensor_space, cx, cy)

Bases: object

Linear operators

class struphy.polar.linear_operators.PolarExtractionOperator(V, W, blocks_ten_to_pol=None, blocks_ten_to_ten=None, transposed=False)

Bases: LinOpWithTransp

Linear operator mapping from Stencil-/BlockVectorSpace (V) to PolarDerhamSpace (W).

For fixed third index k, the dot product maps tensor-product (tp) basis functions/DOFs a_ijk

1. with index i < n_rings (“polar rings”) to n_polar polar basis functions/DOFs,

2. with index i <= n_rings to n2 tp basis functions/DOFs (“first tp ring”),

and leaves the outer tp zone unchanged (identity map). For notation, see Fig. below.

/ / / / /

k = 1 / / / / /

k = 0 / / / / /

——————————- —————–

j = 0 | | | | | j = 1 | | | | |

| | | |

i = 0 | … | i = n_rings | i > n_rings |

| | | | /

| | | | /

j = n2-1 | | | | | /

polar rings | first tp ring | outer tp zone |

Fig. : Indexing of 3d spline tensor-product (tp) basis functions/DOFs/coefficients.

Parameters:

• V (StencilVectorSpace | BlockVectorSpace) – Domain of the operator (always corresponding to the case transpose=False).

• W (PolarDerhamSpace) – Codomain of the operator (always corresponding to the case transpose=False).

• blocks_ten_to_pol (list) –

  2D nested list with matrices that map inner most n_rings tp rings to n_polar polar basis functions/DOFs.

  • shape[m][n] = (n_polar[m], n_rings[n]*n2) if transposed=False.

  • shape[m][n] = (n_rings[m]*n2, n_polar[n]) if transposed=True.

• blocks_ten_to_ten (list) –

  2D nested list with matrices that map inner most n_rings + 1 tp rings to n2 tp basis functions/DOF on “first tp ring”.

  • shape[m][n] = (n2, (n_rings[n] + 1)*n2) if transposed=False.

  • shape[m][n] = ((n_rings[m] + 1)*n2, n2) if transposed=True.

• transposed (bool) – Whether the transposed extraction operator shall be constructed.

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

The data type of the coefficients of the linear operator, upon convertion to matrix.

property tosparse

Convert to a sparse matrix in any of the formats supported by scipy.sparse.

property toarray

Convert to Numpy 2D array.

property transposed

property blocks_ten_to_pol_shapes

property blocks_ten_to_pol

property blocks_ten_to_ten_shapes

property blocks_ten_to_ten

property blocks_e3

dot(v, out=None)

Dot product mapping from Stencil-/BlockVector to PolarVector (or vice versa in case of transposed)

Parameters:

• v (StencilVector | BlockVector | (PolarVector, if transposed)) – Input (domain) vector.

• out (PolarVector | (StencilVector | BlockVector, if transposed), optional) – Optional output vector the result will be written to in-place.

Returns:

out – Output (codomain) vector.

Return type:

PolarVector | (StencilVector | BlockVector, if transposed)

transpose(conjugate=False)

Returns the transposed operator.

class struphy.polar.linear_operators.PolarLinearOperator(V, W, tp_operator=None, blocks_pol_to_ten=None, blocks_pol_to_pol=None, blocks_e3=None, transposed=False)

Bases: LinOpWithTransp

Linear operator mapping from PolarDerhamSpace (V) to PolarDerhamSpace (W).

The dot product maps a PolarVector to a PolarVector with

1. outer tp zone to outer tp zone (stencil),

2. polar coeffs to polar coeffs (dense),

3. polar coeffs to “first tp ring” (dense).

“Polar rings” are alway zero for PolarVectors. For notation, see Fig. below.

/ / / / /

k = 1 / / / / /

k = 0 / / / / /

——————————- —————–

j = 0 | | | | | j = 1 | | | | |

| | | |

i = 0 | … | i = n_rings | i > n_rings |

| | | | /

| | | | /

j = n2-1 | | | | | /

polar rings | first tp ring | outer tp zone |

Fig. : Indexing of 3d spline tensor-product (tp) basis functions/DOFs/coefficients.

Parameters:

• V (PolarDerhamSpace) – Domain of the operator (always corresponding to the case transposed=False).

• W (PolarDerhamSpace) – Codomain of the operator (always corresponding to the case transposed=False).

• tp_operator (LinOpWithTransp) – Standard (stencil) linear operator on the outer tp zone.

• blocks_pol_to_ten (list) –

  2D nested list with matrices that map polar coeffs to “first tp ring”.

  • shape[m][n] = ((n_rings[m] + 1)*n2, n_polar[n]) if transposed=False.

  • shape[m][n] = (n_polar[m], (n_rings[n] + 1)*n2) if transposed=True.

• blocks_pol_to_pol (list) –

  2D nested list with matrices that map polar coeffs to polar coeffs.

  • shape[m][n] = (n_polar[m], n_polar[n]).

• blocks_e3 (list) – 2D nested list with matrices that map in eta_3 direction.

• transposed (bool) – Whether to create the transposed extraction operator.

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

The data type of the coefficients of the linear operator, upon convertion to matrix.

property tosparse

Convert to a sparse matrix in any of the formats supported by scipy.sparse.

property toarray

Convert to Numpy 2D array.

property transposed

property tp_operator

property blocks_pol_to_ten_shapes

property blocks_pol_to_ten

property tp_blocks_shapes

property tp_blocks

property blocks_pol_to_pol_shapes

property blocks_pol_to_pol

property blocks_e3_shapes

property blocks_e3

dot(v, out=None)

Dot product mapping from PolarVector to PolarVector

Parameters:

• v (struphy.polar.basic.PolarVector) – Input (domain) vector.

• out (struphy.polar.basic.PolarVector, optional) – If given, the result will be written into this vector in-place.

Returns:

out – Output (codomain) vector.

Return type:

struphy.polar.basic.PolarVector

transpose(conjugate=False)

Returns the transposed operator.

struphy.polar.linear_operators.dot_inner_tp_rings(blocks_e1_e2, blocks_e3, v, out)

Maps either

1. “polar rings” of v to polar coeffs of out (blocks[m][:].shape[0] = n_polar[m] polar coeffs),

2. “polar rings” + “first tp ring” of v to “first tp ring” of out (blocks[m][:].shape[0] = n2),

and performs a Kronecker product in eta_3 dirction (k-indices). For notation see Fig. below.

/ / / / /

k = 1 / / / / /

k = 0 / / / / /

——————————- —————–

j = 0 | | | | | j = 1 | | | | |

| | | |

i = 0 | … | i = n_rings | i > n_rings |

| | | | /

| | | | /

j = n2-1 | | | | | /

polar rings | first tp ring | outer tp zone |

Fig. : Indexing of 3d spline tensor-product (tp) basis functions/DOFs/coefficients.

Parameters:

• blocks_e1_e2 (list) – 2D nested list with matrices that map inner tp rings to polar coeffs or “first tp ring” depending on shape.

• blocks_e3 (list) – 2D nested list with matrices that solely act along eta_3 direction.

• v (StencilVector | BlockVector) – Input vector.

• out (PolarVector) – Output vector that is written to.

struphy.polar.linear_operators.dot_parts_of_polar(blocks_e1_e2, blocks_e3, v, out)

Maps either

1. polar coeffs of v to “polar rings” of out (blocks[:][n].shape[1] = n_polar[n] polar coeffs),

2. “first tp ring” of v to “polar rings” + “first tp ring” of out (blocks[:][n].shape[1] = n2),

and performs a Kronecker product in eta_3 dirction (k-indices). For notation see Fig. below.

/ / / / /

k = 1 / / / / /

k = 0 / / / / /

——————————- —————–

j = 0 | | | | | j = 1 | | | | |

| | | |

i = 0 | … | i = n_rings | i > n_rings |

| | | | /

| | | | /

j = n2-1 | | | | | /

polar rings | first tp ring | outer tp zone |

Fig. : Indexing of 3d spline tensor-product (tp) basis functions/DOFs/coefficients.

Parameters:

• blocks_e1_e2 (list) – 2D nested list with matrices that map polar coeffs or “first tp ring” to inner to rings depending on shape.

• blocks_e3 (list) – 2D nested list with matrices that solely act along eta_3 direction.

• v (StencilVector | BlockVector) – Input vector.

• out (PolarVector) – Output vector that is written to.

struphy.polar.linear_operators.dot_pol_pol(blocks_e1_e2, blocks_e3, v, out)

Maps polar coeffs to polar coeffs.

Parameters:

• blocks_e1_e2 (list[list[ndarray]]) – blocks_e1_e2[m][n].shape = (out.pol[m].shape[0], v.pol[n].shape[0])

• v (PolarVector) – Input vector.

• out (PolarVector) – output vector.

struphy.polar.linear_operators.transpose_block_mat(blocks)

Transpose 2D nested list of numpy arrays.

Parameters:

blocks (list[list[ndarray]]) – Block matrix.

Returns:

out – Transposed block matrix.

Return type:

list[list[ndarray]]

struphy.polar.linear_operators.set_blocks(shapes)

Creates zero blocks list of given shapes.

Parameters:

shapes (list) – 2D list of tuples (.,.) holding the shape of individual blocks (shapes[m][n] = (.,.)).

Returns:

blocks – 2D list of 2D matrices (scipy.sparse or nd.ndarray) filled with zeros.

Return type:

list

struphy.polar.linear_operators.check_blocks(blocks, shapes)

Checks if given blocks have correct shape.

Parameters:

• blocks (list) – 2D list of 2D matrices (scipy.sparse or nd.ndarray).

• shapes (list) – 2D list of tuples (.,.) holding the shape of individual blocks (shapes[m][n] = (.,.)).