Symbolic Expressions

Operators

class pystopt.operators.RepresentationInterface(ambient_dim: int, kernel_arguments: Dict[str, Any] | None = None)

Bases: object

A generic interface for a boundary integral-based representation of a PDE. It is usually given as

\[(\mathcal{I} + \mathcal{K})[q] = b,\]

where \(\mathcal{I}\) is usually an algebraic linear operator, \(\mathcal{K}\) is the layer potential (a compact linear operator) and \(b\) is the right-hand side. Solving this equation gives the density and solves the PDE.

ambient_dim

kernel_arguments

__init__(ambient_dim: int, kernel_arguments: Dict[str, Any] | None = None) → None

pystopt.operators.make_sym_density(self: RepresentationInterface, name: str = 'sigma') → Expression | ndarray

pystopt.operators.make_sym_density(self: BeltramiRepresentation, name: str = 'sigma') → Expression | ndarray

pystopt.operators.make_sym_density(self: ScalarNeumannRepresentation, name: str = 'sigma') → Expression | ndarray

pystopt.operators.make_sym_density(self: TwoPhaseStokesRepresentation, name: str = 'q') → ndarray

pystopt.operators.make_sym_density(self: TwoPhaseSingleLayerElectricRepresentation, name='q')

Returns:

a symbolic expression for the density with the given name, depending on the dimensions of the representation.

pystopt.operators.make_sym_operator(self: RepresentationInterface, density: Expression | ndarray, *args: Any, **kwargs: Any) → Expression | ndarray

pystopt.operators.make_sym_operator(self: BeltramiRepresentation, sigma: Expression | ndarray) → Expression | ndarray

pystopt.operators.make_sym_operator(self: ScalarNeumannRepresentation, sigma: Expression | ndarray) → Expression | ndarray

pystopt.operators.make_sym_operator(self: TwoPhaseSingleLayerStokesRepresentation, density, *, viscosity_ratio_name='viscosity_ratio')

pystopt.operators.make_sym_operator(self: TwoPhaseSingleLayerElectricRepresentation, density, *, ratio_name='R')

Returns:

a symbolic expression for a boundary integral operator applied to the given density.

pystopt.operators.prepare_sym_rhs(self: RepresentationInterface, b: Expression | ndarray, *args: Any, **kwargs: Any) → Expression | ndarray

pystopt.operators.prepare_sym_rhs(self: BeltramiRepresentation, b: Expression | ndarray) → Expression | ndarray

pystopt.operators.prepare_sym_rhs(self: ScalarNeumannRepresentation, b: Expression | ndarray) → Expression | ndarray

pystopt.operators.prepare_sym_rhs(self: TwoPhaseSingleLayerStokesRepresentation, b, *, viscosity_ratio_name='viscosity_ratio')

pystopt.operators.prepare_sym_rhs(self: TwoPhaseSingleLayerElectricRepresentation, b, *, ratio_name='R')

Returns:

a modified right-hand side b that conforms to the integral operator representation.

class pystopt.operators.BeltramiRepresentation(op: BeltramiOperator)

Bases: RepresentationInterface

Boundary integral representation of a Beltrami-type surface operator. The theorerical background for this type of representation can be found in [ONeil2017].

[ONeil2017] (1,2)

M. O’Neil, Second-Kind Integral Equations for the Laplace-Beltrami Problem on Surfaces in Three Dimensions, 2017, arXiv.

__init__(op: BeltramiOperator) → None

Parameters:

precond – a string identifier for the type of preconditioning applied to the operator (as described in [ONeil2017]). Allowed values are "left" and "right".

class pystopt.operators.LaplaceBeltramiRepresentation(ambient_dim: int, *, dim: int | None = None, precond: str | None = None)

Bases: BeltramiRepresentation

class pystopt.operators.YukawaBeltramiRepresentation(ambient_dim: int, *, dim: int | None = None, precond: str | None = None, yukawa_lambda_name: str = 'k')

Bases: BeltramiRepresentation

class pystopt.operators.ScalarNeumannRepresentation(kernel: Kernel, *, sign: int = 1, kernel_arguments: Dict[str, Any] | None = None, use_l2_weighting: bool = True, use_improved_operator: bool = False)

Bases: RepresentationInterface

class pystopt.operators.LaplaceNeumannRepresentation(ambient_dim: int, *, sign: int = 1, use_l2_weighting: bool = True)

Bases: ScalarNeumannRepresentation

class pystopt.operators.YukawaNeumannRepresentation(ambient_dim: int, *, sign: int = 1, use_l2_weighting: bool = True, yukawa_lambda_name: str = 'k')

Bases: ScalarNeumannRepresentation

Derivatives

This module defined all the derivatives required and supported by the optimization problems. The specific implementations are registered through the functools.singledispatch() mechanism.

The intended use of this functionality, e.g., for a ShapeFunctional is

import pystopt.derivatives as grad
g = grad.shape(cost, ambient_dim, dofdesc=dofdesc)

pystopt.derivatives.shape(self, *args, **kwargs)

pystopt.derivatives.shape(self: TwoPhaseStokesBoundaryCondition, op, density, ad, adjoint_density)

pystopt.derivatives.shape(self: YoungLaplaceJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.shape(self: VariableYoungLaplaceJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.shape(self: GravityJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.shape(self: TwoPhaseStokesRepresentation, op, density, ad, adjoint_density, viscosity_ratio_name='viscosity_ratio')

pystopt.derivatives.shape(self: AXPBYFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: GeometryTrackingFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

pystopt.derivatives.shape(self: CentroidTrackingFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: AreaTrackingFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: VolumeTrackingFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: EnergyFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: TangentialGradientFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: NormalVelocityFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: NormalTestFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: CurvatureTestFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.shape(self: CurvatureVectorTestFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

Computes the shape derivative of self.

pystopt.derivatives.velocity(self, *args, **kwargs)

pystopt.derivatives.velocity(self: TwoPhaseStokesBoundaryCondition, op, density, ad, adjoint_density)

pystopt.derivatives.velocity(self: TractionJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.velocity(self: TwoPhaseStokesRepresentation, op, density, ad, adjoint_density)

pystopt.derivatives.velocity(self: AXPBYFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.velocity(self: GeometryShapeFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.velocity(self: NormalVelocityFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

Computes the derivative with respect to the velocity field of self.

pystopt.derivatives.traction(self, *args, **kwargs)

pystopt.derivatives.traction(self: TwoPhaseStokesBoundaryCondition, op, density, ad, adjoint_density)

pystopt.derivatives.traction(self: TractionJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.traction(self: TwoPhaseStokesRepresentation, op, density, ad, adjoint_density)

pystopt.derivatives.traction(self: AXPBYFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.traction(self: GeometryShapeFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

pystopt.derivatives.traction(self: NormalVelocityFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

Computes the derivative with respect to the traction of self.

pystopt.derivatives.capillary(self, *args, **kwargs)

pystopt.derivatives.capillary(self: TwoPhaseStokesBoundaryCondition, op, density, ad, adjoint_density)

pystopt.derivatives.capillary(self: YoungLaplaceJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.capillary(self: VariableYoungLaplaceJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.capillary(self: GravityJumpCondition, op, density, ad, adjoint_density)

pystopt.derivatives.capillary(self: TwoPhaseStokesRepresentation, op, density, ad, adjoint_density)

pystopt.derivatives.capillary(self: AXPBYFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

pystopt.derivatives.capillary(self: GeometryShapeFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

pystopt.derivatives.capillary(self: NormalVelocityFunctional, ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the derivative with respect to the Capillary number of self.

Shape Functionals

All the cost functionals here implement the derivatives defined in pystopt.derivatives. All the implementations have the signature

pystopt.cost.func(self: ShapeFunctional, ambient_dim, dim, dofdesc)

Specific functionals can be obtained programmatically using

pystopt.cost.cost_functional_from_name(name: str, **kwargs: Any) → ShapeFunctional

Parameters:

name – name of the cost functional. The stringified name is just the name without the Functional part joined by underscores. For example, for GeometryTrackingFunctional is "geometry_tracking".

All cost functionals should inherit from

class pystopt.cost.ShapeFunctional

Bases: object

A generic shape functional.

Shape functionals also support basic arithmetic, i.e. multiplication / division by numeric constants and addition / subtraction between functionals. These operations result in a AXPBYFunctional.

__call__(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Call self as a function.

class pystopt.cost.AXPBYFunctional(func: Tuple[ShapeFunctional, ...], factor: Tuple[Real, ...])

Bases: ShapeFunctional

A linear combination of shape functionals of the form

\[\mathcal{J} = \sum f_i \mathcal{J}_i,\]

where \(f_i\) are given by factor and \(J_i\) are given by func.

func

A tuple of ShapeFunctionals.

factor

A tuple of real numbers multiplying the functionals in func.

class pystopt.cost.GeometryShapeFunctional

Bases: ShapeFunctional

Shape functional that only depends on the geometry.

For example, the Stokes velocity field (or other variables) does not appear in the cost functional. This then sets the relevant derivatives to zero.

Some example geometry-based shape functionals (that have all zero derivatives with respect to the Stokes variables). First, we have a set of geometry tracking cost functionals in various (possibly time-dependent) flows.

class pystopt.cost.GeometryTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None, allow_tangential_gradient: bool)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \int_{\Sigma(t)} \|\mathbf{X}(t) - \mathbf{X}_d(t)\|^2 \,\mathrm{d}S.\]

where \(t\) is a parameter for the interface (in practice, this is time).

ambient_dim

t

A time-like parameter for the geometry, if any.

xd

Reference desired geometry coordinates. To get the actual coordinates for use nodes(), as they can depend on other parameters in the class.

dofdesc

A DOFDescriptor for the geometry to which xd belongs to, if any.

allow_tangential_gradient

If False, only the normal component of the gradient is used. Otherwise, a tangential component is present, that can allow better matching discrete points.

nodes() → ndarray

Returns:

an expression for the nodes based on xd and other parameters describing the desired geometry.

class pystopt.cost.UniformGeometryTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None, allow_tangential_gradient: bool, uinf_d: ndarray)

Bases: GeometryTrackingFunctional

Implements time-dependent geometry tracking in a uniform flow.

The desired geometry \(\mathbf{X}_d\) is parametrized by a time-like variable. In this case, the desired geometry is given by

\[\mathbf{X}_d(t) \triangleq \mathbf{x}_d + t \mathbf{u}^\infty,\]

where \(\mathbf{u}^\infty\) is a uniform farfield velocity field that translates the given geometry.

See also get_uniform_farfield().

uinf_d

Uniform velocity field that translates the geometry xd.

class pystopt.cost.SolidBodyRotationGeometryTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None, allow_tangential_gradient: bool, omega_d: ndarray)

Bases: GeometryTrackingFunctional

Implements time-dependent geometry tracking in a solid body rotation flow.

In this case, the desired geometry is given by

\[\dot{\mathbf{X}}(t) = \omega_d \times \mathbf{X}\]

with the usual restriction of the cross product to the two-dimensional case. In the 2D case, this is just a rotation in the \(xy\) plane, but in 3D it can be any arbitrary rotation around a given axis \(\omega_d\).

See also get_solid_body_rotation_farfield().

omega_d

Axis of the rotation used to define the position of the desired geometry at a given time t.

class pystopt.cost.BezierGeometryTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None, allow_tangential_gradient: bool, points_d: Tuple[ndarray, ...])

Bases: GeometryTrackingFunctional

Implements time-dependent geometry tracking in a uniform flow.

In this case, the desired geometry is given by

\[\dot{\mathbf{X}}(t) = \sum P_i \phi_i'(t)\]

where \(\phi_i\) are the Bernstein polynomials of order order. The actual position of the interface is then just obtained by integration using GeometryTrackingFunctional.xd as an initial condition.

See also get_bezier_uniform_farfield().

npoints

Number of control points in the Bézier curve.

order

Order of the Bézier curve.

points_d

Control points that define the Bézier curve for the velocity field.

class pystopt.cost.HelixGeometryTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None, allow_tangential_gradient: bool, omega_d: float | pytential.symbolic.primitives.var, height_d: float | pytential.symbolic.primitives.var)

Bases: GeometryTrackingFunctional

Implements time-dependent geometry tracking in a helical flow.

In this case, the desired geometry is given by

\[\begin{split}\begin{cases} x(t) = R \cos (\omega_d t), \\ y(t) = R \sin (\omega_d t), \\ z(t) = H_d t \end{cases}\end{split}\]

where \(R\) is the radius, \(\omega_d\) is the omega_d, \(H\) is the height_d and \(\t\) is a time-like variable in \([0, 1]\).

See also get_helical_farfield().

radius

Radius of the helix.

omega_d

Wavenumber that, together with the height_d, determines the number of turns in the helix.

height_d

The height of the helix.

Similarly, we have a set of centroid tracking cost functionals.

class pystopt.cost.CentroidTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \|\mathbf{x}_c(t) - \mathbf{x}_d(t)\|^2,\]

where \(\mathbf{x}_c\) is the centroid and \(\mathbf{x}_d\) is a desired position.

ambient_dim

t

A time-like parameter for the geometry, if any.

xd

Reference desired geometry centroid. To get the actual coordinates for use nodes(), as they can depend on other parameters in the class.

nodes() → ndarray

Returns:

an expression for the centroid based on xd and other parameters describing the desired geometry.

class pystopt.cost.UniformCentroidTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, uinf_d: ndarray)

Bases: CentroidTrackingFunctional

class pystopt.cost.SolidBodyRotationCentroidTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, omega_d: ndarray)

Bases: CentroidTrackingFunctional

class pystopt.cost.ExtensionalCentroidTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, alpha_d: float | pytential.symbolic.primitives.var)

Bases: CentroidTrackingFunctional

class pystopt.cost.HelixCentroidTrackingFunctional(xd: ndarray, t: pytential.symbolic.primitives.var | None, height_d: float | pytential.symbolic.primitives.var, omega_d: float | pytential.symbolic.primitives.var)

Bases: CentroidTrackingFunctional

Other examples of cost functionals that only depend on the geometry are also given.

class pystopt.cost.AreaTrackingFunctional(ad: pytential.symbolic.primitives.var)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \left(\int_\Sigma\,\mathrm{d}S - A_d\right)^2\]

where \(A_d\) is a desired area.

class pystopt.cost.VolumeTrackingFunctional(vd: pytential.symbolic.primitives.var)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \left(\int_\Omega\,\mathrm{d}V - V_d\right)^2\]

where \(V_d\) is a desired volume. The volume in question here is interior to a closed surface \(\Sigma\), so that we can use the Divergence Theorem to easily compute it.

class pystopt.cost.EnergyFunctional(f: pytential.symbolic.primitives.var, gradf: ndarray | None = None)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \int_\Sigma f^2 \,\mathrm{d}S.\]

f

Scalar function used in cost functional.

gradf

If the scalar f is defined in the ambient space and not just on the surface \(\Sigma\), the full gradient is also required to compute the cost gradient.

class pystopt.cost.TangentialGradientFunctional(f: pytential.symbolic.primitives.var, gradf: ndarray | None = None, hessf: ndarray | None = None)

Bases: GeometryShapeFunctional

Implements the evaluation and gradient for the cost functional

\[\mathcal{J} = \frac{1}{2} \int_\Sigma \nabla_\Sigma f \cdot \nabla_\Sigma f \,\mathrm{d}S.\]

f

Scalar function used in cost functional.

gradf

Ambient gradient of f. If not provided, the gradient is computed numerically.

Finally, some example Stokes-based shape functionals are given.

class pystopt.cost.NormalVelocityFunctional(u: ndarray, ud: pytential.symbolic.primitives.var, dofdesc: pytential.symbolic.primitives.DOFDescriptor | None = None, gradu_dot_t: ndarray | None = None, gradud: ndarray | None = None, divu: pytential.symbolic.primitives.var | None = None)

Bases: ShapeFunctional

Implements the evaluation and gradients for the cost functional

\[\mathcal{J} = \frac{1}{2} \int_\Sigma |\mathbf{u} \cdot \mathbf{n} - u_d|^2 \,\mathrm{d}S.\]

u

Velocity field used in the cost functional.

ud

Desired normal component of the velocity field.

dofdesc

A DOFDescriptor for the geometry on which the desired velocity was computed, if any.

gradu_dot_t

Tangential components of the velocity gradient. They are only required to compute the shape gradient. If not provided, this is computed numerically.

gradud

Surface gradient of ud. If ud is constant, this is determined automatically, otherwise it is computed, e.g., using surface_gradient().

pystopt.cost.get_desired_normal_velocity(actx: ArrayContext, places: pytential.GeometryCollection, op: TwoPhaseStokesRepresentation, *, context: Dict[str, Any] | None = None, lambdas: Dict[str, float] | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Dict[str, Any]

Solves the Stokes system on a desired geometry in places described by dofdesc.

Parameters:

• context – a dict of additional variables required to solve the Stokes problem described by op.

• lambdas – viscosity ratios for the desired geometry.

Returns:

a dict which contains the required desired velocity field for NormalVelocityFunctional. The keys in the dictionary match the attribute names.

pystopt.cost.make_normal_velocity_functional_from_stokes(op: TwoPhaseStokesRepresentation, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → NormalVelocityFunctional

The following cost functionals are provided for testing only (used to verify the shape gradients of each quantity).

class pystopt.cost.NormalTestFunctional(f: ndarray, divf: pytential.symbolic.primitives.var)

Bases: GeometryShapeFunctional

The functional is given by

\[\mathcal{J} = \int_\Sigma \mathbf{f} \cdot \mathbf{n} \,\mathrm{d}S.\]

with an (Eulerian) shape gradient given by

\[\nabla_{\mathbf{X}} \mathcal{J} = \nabla \cdot \mathbf{f}\]

f

divf

class pystopt.cost.CurvatureTestFunctional(f: pytential.symbolic.primitives.var, gradf: ndarray, hessf: ndarray)

Bases: GeometryShapeFunctional

The functional is given by

\[\mathcal{J} = \int_\Sigma \kappa f \,\mathrm{d}S.\]

f

divf

hessf

class pystopt.cost.CurvatureVectorTestFunctional(f: ndarray, gradf: pytential.symbolic.primitives.var, hessf: ndarray)

Bases: GeometryShapeFunctional

The functional is given by

\[\mathcal{J} = \int_\Sigma \kappa \mathbf{f} \cdot \mathbf{n} \,\mathrm{d}x.\]

f

divf

hessf

Primitives

All the functions from pytential.symbolic.primitives are exported from here. Additionally, the following primitives are implemented.

pystopt.symbolic.primitives.first_fundamental_form(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the matrix elements of the first fundamental form.

The elements are given by

\[g_{ij} \triangleq \frac{\partial X_k}{\partial \xi_i} \frac{\partial X_k}{\partial \xi_j},\]

where \(\mathbf{X}(\xi)\) is the surface parametrization and \((i, j) \in \{1, d\}\) for \(d\) = dim.

pystopt.symbolic.primitives.inverse_first_fundamental_form(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the matrix elements of the inverse of the first fundamental form given in first_fundamental_form().

pystopt.symbolic.primitives.second_fundamental_form(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the matrix elements of the second fundamental form.

They are given by

\[h_{ij} \triangleq -\frac{\partial n_k}{\partial \xi_i} \frac{\partial X_k}{\partial \xi_j} = n_k \frac{\partial^2 X_k}{\partial \xi_i \partial \xi_j},\]

where \(\mathbf{n}\) is the unit normal vector and \((i, j) \in \{1, d\}\) for \(d\) = dim. Here the first formula is implemented, while pytential.symbolic.primitives.second_fundamental_form() implements the second formula.

pystopt.symbolic.primitives.shape_operator(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the matrix elements of the (intrinsic) shape operator.

They are given by

\[S_{ij} \triangleq -h_{ij} g^{ij}.\]

pystopt.symbolic.primitives.extrinsic_shape_operator(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Computes the matrix elements of the (extrinsic) shape operator, sometimes called the Weingarten map.

Its components form a \(n \times n\) matrix, where \(n\) = ambient_dim. It is given by

\[W_{ij} \triangleq g^{kl} \frac{\partial n_j}{\partial \xi_k} \frac{\partial X_i}{\partial \xi_l}\]

i.e. simply the surface gradient of the normal \(W \triangleq \nabla_\Sigma \mathbf{n}\). It has the property that

\[\frac{\partial X_i}{\partial \xi_k} W_{ij} = \frac{\partial n_j}{\partial \xi_k}.\]

Dotting with the parametrization derivative again retrieves the components of the (intrinsic) shape operator.

pystopt.symbolic.primitives.summed_curvature(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Sum of the principal curvatures or trace of the shape_operator()

\[\kappa \triangleq \operatorname{tr} S = \sum_{i = 0}^{d - 1} \kappa_i.\]

pystopt.symbolic.primitives.mean_curvature(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Average of the principal curvatures from summed_curvature().

\[H = \frac{\kappa}{d} = \frac{1}{d - 1} \operatorname{tr} S.\]

pystopt.symbolic.primitives.gaussian_curvature(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Product of the principal curvatures or determinant of the shape_operator()

\[\kappa_G \triangleq \operatorname{det} S = \prod_{i = 1}^d \kappa_i.\]

For one dimensional surfaces, the determinant is \(0\) and the Gaussian curvature vanishes.

pystopt.symbolic.primitives.squared_curvature(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Sum of squares of the principal curvatures

\[K \triangleq \sum_{i = 1}^d \kappa_i^2.\]

For 2D surfaces in a 3D ambient space, this can also be written in terms of the summed curvature and Gaussian curvature

\[K = \kappa^2 - 2 \kappa_G.\]

pystopt.symbolic.primitives.principal_curvature(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → Expression

Principal curvatures are the eigenvalues of the shape_operator().

Returns:

an ndarray of size dim with the expressions of all principal curvatures.

pystopt.symbolic.primitives.principal_directions(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

Principal directions are the eigenvectors of the shape_operator().

The principal directions returned by this function solve the generalized eigenvalue problem

\[-h q_i = \kappa_i g q_i,\]

where \(\kappa_i\) are the principal curvatures, \(h\) is the second fundamental form and \(g\) is the first fundamental form. To get the actual eigenvectors of the shape operator, it suffices to multiply by the first fundamental form, since

\[-h g^{-1} d_i = S d_i = \kappa_i d_i,\]

where \(d_i \triangleq g q_i\).

pystopt.symbolic.primitives.extrinsic_principal_directions(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None) → ndarray

Principal directions in ambient coordinates, which correspond to the eigenvectors of extrinsic_shape_operator().

They are obtained by

\[D_{i, j} \triangleq \frac{\partial X_j}{\xi_k} d_{i, k}\]

where \(d_i\) are the intrinsic principal directions from principal_directions(). Note that the principal directions form an orthonormal basis of the tangent space, similar to tangential_onb().

pystopt.symbolic.primitives.surface_gradient(ambient_dim: int, operand: Expression, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, is_spectral: bool = False) → Expression

Computes the surface gradient of operand. In terms of the ambient gradient it is given as

\[\nabla_\Sigma f = \nabla f - (\mathbf{n} \cdot \nabla f) \mathbf{n},\]

but it can also be written intrinsically as

\[(\nabla_\Sigma f)_k = g^{ij} \frac{\partial f}{\partial \xi_i} \frac{\partial X_k}{\partial \xi_j}.\]

pystopt.symbolic.primitives.surface_divergence(ambient_dim: int, operand: ndarray, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

Computes the surface divergence of operand. In terms of the ambient gradient it is given as

\[\nabla_\Sigma \cdot \mathbf{f} = \nabla \cdot \mathbf{f} - (\mathbf{n} \cdot \nabla \mathbf{f}) \cdot \mathbf{n},\]

but it can also be written intrinsically as

\[\nabla_\Sigma \cdot \mathbf{f} = g^{ij} \frac{\partial f_k}{\partial \xi_i} \frac{\partial X_k}{\partial \xi_j}.\]

pystopt.symbolic.primitives.surface_laplace_beltrami(ambient_dim: int, operand: Expression, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

Computes the surface Laplace-Beltrami operator. In terms of the ambient gradient it is given as

\[\Delta_\Sigma f = \Delta f - \kappa (\mathbf{n} \cdot \nabla f) - (\mathbf{n} \cdot \nabla \nabla f) \cdot \mathbf{n},\]

but it can also be written intrinsically as

\[\Delta_\Sigma f = g^{ij} \frac{\partial}{\partial \xi_i} \left( g^{kl} \frac{\partial f}{\partial \xi_k} \frac{\partial X_k}{\partial \xi_l} \right) \frac{\partial X_k}{\partial \xi_j}.\]

pystopt.symbolic.primitives.surface_volume(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, groupwise: bool = False)

Volume (or area in 2D) enclosed by a closed surface \(Sigma\)

\[V \triangleq \int_{\Sigma} \frac{1}{d} \mathbf{x} \cdot \mathbf{n} \,\mathrm{d}S.\]

pystopt.symbolic.primitives.surface_centroid(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, groupwise: bool = False) → ndarray

Centroid of a closed surface

\[\mathbf{c} = \frac{1}{|\Omega|} \int_{\Omega} \mathb{x} \,\mathrm{d}S,\]

where \(|\Omega|\) is the volume enclosed by the surface.

pystopt.symbolic.primitives.h_max_from_volume(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

Computes a characteristc (largest) length scale from element areas

\[h_{max} \triangleq \max_{k} \left|\int_{E_k} \,\mathrm{d}S\right|^{1/d}.\]

Note that this metric is mostly suited for meshes where each element is mostly isotropic, i.e. stretches the same in all directions. For a more general metric on simplices, see pytential.symbolic.primitives.h_max().

pystopt.symbolic.primitives.h_min_from_volume(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

Computes a characteristc (smallest) length scale from element areas

\[h_{min} \triangleq \min_{k} \left|\int_{E_k} \,\mathrm{d}S\right|^{1/d}.\]

pystopt.symbolic.primitives.h_min(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

Computes a characteristic (smallest) length scale from element.

Like pytential.symbolic.primitives.h_max(), this based on looking at the stretching of an element.

pystopt.symbolic.primitives.gaussian_bump_function(ambient_dim: int, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, normalized: bool = True, sigma_name: str = 'sigma')

Constructs a simple Gaussian bump function

\[h(\mathbf{x}) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{\mathbf{x} \cdot \mathbf{x}}{2 \sigma^2}\right),\]

where \(\sigma\) is the variance of the Gaussian. If normalized is True, the bump function is normalized with respect to the \(\mathrm{L}^\infty\) norm, i.e. the prefactor is removed.

pystopt.symbolic.primitives.spherical_bump_function(ambient_dim: int, h: Expression, *, normalized=True)

Constructs a bump function

\[\tilde{h}(\mathbf{h}) = \sin \theta h(\mathbf{x}),\]

where h is a given bump function and the angle \(theta\) is defined as

\[\theta = \arccos \frac{z}{r}.\]

This ensures that no issues arise at the poles :math:theta = 0, pi, as the bump function goes to zero.

pystopt.symbolic.primitives.norm(ambient_dim: int, operand: Expression, *, p: Real = 2, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

pystopt.symbolic.primitives.relative_norm(ambient_dim: int, x: Expression, y: Expression, *, p: Real = 2, atol: Real = 1e-14, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None)

pystopt.symbolic.primitives.n_dot(ambient_dim_or_expr: int | Expression, *, dim: int | None = None, dofdesc: pytential.symbolic.primitives.DOFDescriptorLike | None = None, name: str = 'x') → Expression

Geometric Shape Derivatives

pystopt.symbolic.derivatives.adjoint_normal_shape_derivative(ambient_dim, f, *, dim=None, dofdesc=None)

Computes the adjoint of the Eulerian shape derivative of the normal applied to the vector f with respect to the \(L^2\) inner product.

From Lemma 5.5 in [Walker2015], we have that the Eulerian shape derivative of the normal is given by

\[\mathbf{n}'[\mathbf{V}] = -\nabla_\Sigma (\mathbf{V} \cdot \mathbf{n})\]

so, using the surface Divergence Theorem, we get that its adjoint is

\[\mathbf{n}^*[\mathbf{f}] = (\nabla_\Sigma \cdot \mathbf{f} - \kappa \mathbf{f} \cdot \mathbf{n}) \mathbf{n}.\]

pystopt.symbolic.derivatives.adjoint_normal_shape_gradient(ambient_dim, f, *, divf, dim=None, dofdesc=None)

Shape gradient of the functional

\[\int_\Sigma \mathbf{f} \cdot \mathbf{n} \,\mathrm{d}S = \langle \nabla g, \mathbf{V} \cdot \mathbf{n}\rangle.\]

pystopt.symbolic.derivatives.adjoint_curvature_shape_derivative(ambient_dim, f, *, dim=None, dofdesc=None)

Computes the adjoint of the Eulerian shape derivative of the curvature applied to the field f, with respect to the \(L^2\) inner product.

From Lemma 5.6 in [Walker2015], we have that the Eulerian shape derivative of the curvature is given by

\[\kappa'[\mathbf{V}] = -\Delta_\Sigma (\mathbf{V} \cdot \mathbf{n}) - (\kappa^2 - 2 \kappa_G) (\mathbf{V} \cdot \mathbf{n})\]

so, using the surface Divergence Theorem, we get that its adjoint is

\[\kappa^*[f] = -\Delta_\Sigma f - (\kappa^2 - 2 \kappa_G) f\]

pystopt.symbolic.derivatives.adjoint_curvature_shape_gradient(ambient_dim, f, *, gradf_dot_n, lapf, hessf_dot_n, dim=None, dofdesc=None)

Shape gradient of the functional

\[\int_\Sigma \kappa f \,\mathrm{d}S = \langle \nabla g, \mathbf{V} \cdot \mathbf{n}\rangle.\]

pystopt.symbolic.derivatives.adjoint_curvature_vector_shape_derivative(ambient_dim, f, *, dim=None, dofdesc=None)

pystopt.symbolic.derivatives.adjoint_curvature_vector_shape_gradient(ambient_dim, *, f, divf, gradf_dot_n, gradf_dot_t, lapf_dot_n, hessf_dot_n, dim=None, dofdesc=None)