movement package

Submodules

movement.laplacian_smoothing module

class LaplacianSmoother(mesh, timestep, **kwargs)

Bases: PrimeMover

Movement of a mesh is driven by a mesh velocity \(\mathbf{v}\), which is determined by solving a vector Laplace problem

\[\nabla^2_{\boldsymbol{\xi}}\mathbf{v} = \boldsymbol{0}, \quad \boldsymbol{\xi}\in\Omega,\]

under non-zero Dirichlet boundary conditions on a forced boundary section \(\partial\Omega_f\) and zero Dirichlet boundary conditions elsewhere:

\[\begin{split}\mathbf{v} = \left\{\begin{array}{rl} \mathbf{v}_D, & \boldsymbol{\xi}\in\partial\Omega_f\\ \boldsymbol{0}, & \boldsymbol{\xi}\in \partial\Omega\backslash\partial\Omega_f \end{array}\right.\end{split}\]

where the computational coordinates \(\boldsymbol{\xi} := \mathbf{x}(t_0)\) are the physical coordinates at the beginning of the simulation.

move(time, update_boundary_velocity=None, boundary_conditions=None, solver_parameters=None)

Assemble and solve the Laplacian system and update the coordinates.

Parameters:

• time (float) – the current time

• update_boundary_velocity (Callable with a single argument of float type) – function that updates the boundary conditions at the current time

• boundary_conditions (DirichletBC or list thereof) – Dirichlet boundary conditions to be enforced

• solver_parameters (dict) – solver parameters for solving the Laplace equation

movement.math module

equation_of_hyperplane(*points)

Deduce an expression for the equation of a hyperplane passing through a set of points.

Parameters:

points (tuple of tuples) – points the hyperplane passes through

Returns:

a function representing the hyperplane

Return type:

Callable

movement.monge_ampere module

Mesh movement based on solutions of equations of Monge-Ampère type.

MongeAmpereMover(mesh, monitor_function, method='relaxation', **kwargs)

Factory function for generating Monge-Ampère mesh movers.

Movement of a mesh is determined by a monitor_function \(m\) and an equation of Monge-Ampère type,

\[m(\mathbf{x})\det(\mathbf{I} + \mathbf{H}(\phi)) = \theta,\]

for a convex scalar potential \(\phi=\phi(\boldsymbol{\xi})\), which is a function of the coordinates of the computational mesh. Here \(m=m(\mathbf{x})\) is a function of the coordinates of the physical mesh, \(\mathbf{I}\) is the identity matrix, \(\theta\) is a normalisation coefficient and \(\mathbf{H}(\phi)\) denotes the Hessian of \(\phi\) with respect to \(\boldsymbol{\xi}\).

Different implementations solve the Monge-Ampère equation in different ways. If the method argument is set to “relaxation” then it is solved in parabolised form in MongeAmpereMover_Relaxation. If the argument is set to “quasi_newton” then it is solved in its elliptic form using a quasi-Newton method in MongeAmpereMover_QuasiNewton. Descriptions of both methods may be found in [McRae et al., 2018].

The physical mesh coordinates \(\mathbf{x}\) are updated according to

\[\mathbf{x} = \boldsymbol{\xi} + \nabla_{\boldsymbol{\xi}}\phi.\]

Parameters:

• mesh (firedrake.mesh.MeshGeometry) – the physical mesh

• monitor_function (Callable) – a Python function which takes a mesh as input

• method (str) – choose from ‘relaxation’ and ‘quasi_newton’

• phi_init (firedrake.function.Function) – initial guess for the scalar potential

• H_init (firedrake.function.Function) – initial guess for the Hessian

• maxiter (int) – maximum number of iterations for the solver

• rtol (float) – relative tolerance for the residual

• dtol (float) – divergence tolerance for the residual

• pseudo_timestep (float) – pseudo-timestep (only relevant to relaxation method)

• fixed_boundary_segments (list of str or int) – labels corresponding to boundary segments to be fixed with a zero Dirichlet condition. The ‘on_boundary’ label indicates the whole domain boundary

Returns:

the Monge-Ampère Mover object

Return type:

MongeAmpereMover_Relaxation or MongeAmpereMover_QuasiNewton

tangential(v, n)

Return component of v perpendicular to n (assumed normalised).

This is used to project vectors onto the tangent plane of a boundary.

Parameters:

• v (ufl.Expr) – the vector to project

• n (ufl.Expr) – the normal vector

class MongeAmpereMover_Base(mesh, monitor_function, **kwargs)

Bases: PrimeMover

Base class for mesh movers based on the solution of Monge-Ampère type equations.

Currently implemented subclasses: MongeAmpereMover_Relaxation and MongeAmpereMover_QuasiNewton. Descriptions of both methods may be found in [McRae et al., 2018].

apply_initial_guess(phi_init, H_init)

Initialise the approximations to the scalar potential and its Hessian with an initial guess.

By default, both are initialised to zero, which corresponds to the case where the computational and physical meshes coincide.

Parameters:

• phi_init (firedrake.function.Function) – initial guess for the scalar potential

• H_init (firedrake.function.Function) – initial guess for the Hessian

property relative_l2_residual

Returns:

the relative \(L^2\) norm residual.

Return type:

float

property l2_projector

Create a linear solver for obtaining the gradient of the potential using an \(L^2\) projection.

Returns:

the linear solver

Return type:

LinearVariationalSolver

class MongeAmpereMover_Relaxation(mesh, monitor_function, phi_init=None, H_init=None, **kwargs)

Bases: MongeAmpereMover_Base

The standard, elliptic form of the Monge-Ampère equation used for mesh movement is:

\[m(\mathbf{x})\det(\mathbf{I} + \mathbf{H}(\phi)) = \theta,\]

for a convex scalar potential \(\phi=\phi(\boldsymbol{\xi})\), which is a function of the coordinates of the computational mesh. Here \(m=m(\mathbf{x})\) is a user-provided monitor function, which is a function of the coordinates of the physical mesh. \(\mathbf{I}\) is the identity matrix, \(\theta\) is a normalisation coefficient and \(\mathbf{H}(\phi)\) denotes the Hessian of \(\phi\) with respect to \(\boldsymbol{\xi}\).

In this mesh mover, the Monge-Ampère equation is instead solved in a parabolised form using a pseudo-time relaxation,

\[-\frac\partial{\partial\tau}\Delta\phi = m(\mathbf{x})\det(\mathbf{I} + \mathbf{H}(\phi)) - \theta,\]

where \(\tau\) is the pseudo-time variable. Forward Euler is used for the pseudo-time integration (see [McRae et al., 2018] for details).

This approach typically takes tens or hundreds of iterations to converge, but each iteration is relatively cheap.

property pseudotimestepper

Setup the pseudo-timestepper for the relaxation method.

Forward Euler is used for the pseudo-time integration (see [McRae et al., 2018] for details). The pseudo-timestep may be set through the pseudo_timestep keyword argument to the constructor.

Returns:

the pseudo-timestepper

Return type:

LinearVariationalSolver

property equidistributor

Setup the equidistributor for the relaxation method.

The equidistributor solves the following equation:

\[\int_{\Omega} \tau : \mathbf{H} \, \mathrm{d}x = -\int_{\Omega} (\nabla \cdot \tau) \cdot (\nabla \phi) \, \mathrm{d}x + \int_{\partial \Omega} ((\nabla \phi \cdot \widehat{\mathbf{n}}) \cdot \tau) \cdot \widehat{\mathbf{n}} \, \mathrm{d}s, \quad \forall \tau \in \mathbb{P}1^{d \times d}\]

for the Hessian \(\mathbf{H}\), where \(d\) is the spatial dimension and \(\widehat{\mathbf{n}}\) is a normal vector to the boundary.

Returns:

the equidistributor

Return type:

LinearVariationalSolver

move()

Run the relaxation method to convergence and update the mesh.

Returns:

the iteration count

Return type:

int

class MongeAmpereMover_QuasiNewton(mesh, monitor_function, phi_init=None, H_init=None, **kwargs)

Bases: MongeAmpereMover_Base

The standard, elliptic form of the Monge-Ampère equation used for mesh movement is:

\[m(\mathbf{x})\det(\mathbf{I} + \mathbf{H}(\phi)) = \theta,\]

for a convex scalar potential \(\phi=\phi(\boldsymbol{\xi})\), which is a function of the coordinates of the computational mesh. Here \(m=m(\mathbf{x})\) is a user-provided monitor function, which is a function of the coordinates of the physical mesh. \(\mathbf{I}\) is the identity matrix, \(\theta\) is a normalisation coefficient and \(\mathbf{H}(\phi)\) denotes the Hessian of \(\phi\) with respect to \(\boldsymbol{\xi}\).

In this mesh mover, the elliptic Monge-Ampère equation is solved using a quasi-Newton method (see [McRae et al., 2018] for details).

This approach typically takes fewer than ten iterations to converge, but each iteration is relatively expensive.

property equidistributor

Setup the equidistributor for the quasi-newton method.

The equidistributor solves

\[\int_{\Omega} \boldsymbol{\tau} \cdot \mathbf{H} \, \mathrm{d}x + \int_{\Omega} (\nabla \cdot \boldsymbol{\tau}) \cdot (\nabla \phi) \, \mathrm{d}x - \int_{\partial \Omega} (((\nabla \phi) \cdot \widehat{\mathbf{n}}) \cdot \boldsymbol{\tau}) \cdot \widehat{\mathbf{n}} \, \mathrm{d}s - \int_{\Omega} \psi (m \det(\mathbf{I} + \mathbf{H}) - \theta) \, \mathrm{d}x = 0, \quad \forall \boldsymbol{\tau} \in \mathbb{P}1^{d \times d}, \quad \forall \psi \in \mathbb{P}1,\]

for the potential \(\phi\) and its Hessian \(\mathbf{H}\), where \(d\) is the spatial dimension and \(\widehat{\mathbf{n}}\) is a normal vector to the boundary.

Returns:

the equidistributor

Return type:

NonlinearVariationalSolver

move()

Run the quasi-Newton method to convergence and update the mesh.

Returns:

the iteration count

Return type:

int

movement.monitor module

class MonitorBuilder(dim)

Bases: object

Abstract base class for monitor function factories.

get_monitor()

Returns a callable monitor function whose only argument is the mesh.

Returns:

monitor function

Return type:

callable monitor function with a single argument

class ConstantMonitorBuilder(dim)

Bases: MonitorBuilder

Builder class for constant monitor functions.

class BallMonitorBuilder(dim, centre, radius, amplitude, width)

Bases: MonitorBuilder

Builder class for monitor functions focused around ball shapes:

\[m(\mathbf{x}) = 1 + \frac{\alpha} {\cosh^2\left(\beta\left((\mathbf{x}-\mathbf{c})\cdot(\mathbf{x}-\mathbf{c}) -\gamma^2\right)\right)},\]

where \(\mathbf{c}\) is the centre point, \(\alpha\) is the amplitude of the monitor function, \(\beta\) is the width of the transition region, and \(\gamma\) is the radius of the ball.

class RingMonitorBuilder(centre, radius, amplitude, width)

Bases: BallMonitorBuilder

Builder class for monitor functions focused around 2D ring shapes:

\[m(x,y) = 1 + \frac{\alpha} {\cosh^2\left(\beta\left((x-c_0)^2+(y-c_1)^2\right) -\gamma^2\right)},\]

where \((c_0,c_1)\) is the centre point, \(\alpha\) is the amplitude of the monitor function, \(\beta\) is the width of the transition region, and \(\gamma\) is the radius of the ball.

class SolutionBasedMonitorBuilder(dim, solution)

Bases: MonitorBuilder

Abstract base class for monitor factories based on solution data.

projection(mesh)

Project the solution field onto the given mesh.

Parameters:

mesh (firedrake.mesh.MeshGeometry) – the mesh on which the solution is to be projected

Returns:

the projected solution field

Return type:

firedrake.function.Function

class GradientMonitorBuilder(dim, solution, scale_factor)

Bases: SolutionBasedMonitorBuilder

Builder class for monitor functions based on gradients of solutions:

\[m(\mathbf{x}) = 1 + \alpha\frac{\nabla u\cdot\nabla u} {\max_{\mathbf{x}\in\Omega}\nabla u\cdot\nabla u},\]

where \(\alpha\) is a scale factor and \(u\) is the solution field of interest.

recover_gradient(target_space)

Recover the gradient of the solution field projected onto the current mesh.

Parameters:

target_space (firedrake.functionspaceimpl.FunctionSpace) – space to recover gradient in

Returns:

the recovered gradient in vector \(\mathbb{P}1\) space

Return type:

firedrake.function.Function

class HessianMonitorBuilder(dim, solution, scale_factor)

Bases: SolutionBasedMonitorBuilder

Builder class for monitor functions based on Hessians of solutions.

\[m(\mathbf{x}) = 1 + \beta\frac{\mathbf{H}(u):\mathbf{H}(u)} {\max_{\mathbf{x}\in\Omega}\mathbf{H}(u):\mathbf{H}(u)},\]

where \(\beta\) is a scale factor, \(u\) is the solution field of interest, and \(\mathbf{H}(u)\) is the Hessian

recover_hessian(target_space)

Recover the Hessian of the solution field.

Parameters:

target_space (firedrake.functionspaceimpl.FunctionSpace) – space to recover Hessian in

Returns:

the recovered Hessian in tensor \(\mathbb{P}1\) space

Return type:

firedrake.function.Function

class GradientHessianMonitorBuilder(dim, solution, gradient_scale_factor, hessian_scale_factor)

Bases: GradientMonitorBuilder, HessianMonitorBuilder

Builder class for monitor functions based on both gradients and Hessians of solutions.

\[m(\mathbf{x}) = 1 + \alpha\frac{\nabla u\cdot\nabla u} {\max_{\mathbf{x}\in\Omega}\nabla u\cdot\nabla u} + \beta\frac{\mathbf{H}(u):\mathbf{H}(u)} {\max_{\mathbf{x}\in\Omega}\mathbf{H}(u):\mathbf{H}(u)},\]

where \(\alpha\) is a scale factor for the gradient part, \(\beta\) is a scale factor for the Hessian part, \(u\) is the solution field of interest, and \(\mathbf{H}(u)\) is the Hessian

movement.mover module

class PrimeMover(mesh, monitor_function=None, raise_convergence_errors=True, tangling_check=True, quadrature_degree=None)

Bases: ABC

Base class for all mesh movers.

property volume_ratio

Returns:

the ratio of the largest and smallest element volumes.

Return type:

float

property coefficient_of_variation

Returns:

the coefficient of variation (σ/μ) of element volumes.

Return type:

float

abstract move()

Move the mesh according to the method of choice.

to_physical_coordinates()

Switch coordinates to correspond to the physical mesh mathcal{H}_P.

to_computational_coordinates()

Switch coordinates to correspond to the computational mesh mathcal{H}_C.

plural(iterations)

Returns:

‘s’ if iterations should be referred to in the plural sense

Return type:

str

movement.solver_parameters module

movement.spring module

SpringMover(*args, method='lineal', **kwargs)

Movement of a mesh is determined by reinterpreting it as a structure of stiff beams and solving an associated discrete linear elasticity problem. (See [Farhat et al., 1998] for details.)

Parameters:

• mesh (firedrake.mesh.MeshGeometry) – the physical mesh to be moved

• timestep (float) – the timestep length used

• method (str) – flavour of spring-based method to use

class SpringMover_Base(mesh, timestep, **kwargs)

Bases: PrimeMover

Base class for mesh movers based on spring analogies.

property facet_areas

Compute the areas of all facets in the mesh.

In 2D, this corresponds to edge lengths.

property tangents

Compute tangent vectors for all edges in the mesh.

property angles

Compute the argument of each edge in the mesh, i.e. its angle from the \(x\)-axis in the \(x-y\) plane.

assemble_stiffness_matrix(boundary_conditions=None)

Enforce that nodes on certain tagged boundaries do not move.

Parameters:

boundary_conditions (DirichletBC or list thereof) – Dirichlet boundary conditions to be enforced

Returns:

the stiffness matrix with boundary conditions applied

Return type:

numpy.ndarray

class SpringMover_Lineal(mesh, timestep, **kwargs)

Bases: SpringMover_Base

Movement of a mesh is determined by reinterpreting it as a structure of stiff beams and solving an associated discrete linear elasticity problem.

We consider the ‘lineal’ case, as described in [Farhat et al., 1998].

move(time, update_boundary_displacement=None, boundary_conditions=None)

Assemble and solve the lineal spring system and update the coordinates.

Parameters:

• time (float) – the current time

• update_boundary_displacement (Callable with a single argument of float type) – function that updates the boundary conditions at the current time

• boundary_conditions (DirichletBC or list thereof) – Dirichlet boundary conditions to be enforced

class SpringMover_Torsional(*args, **kwargs)

Bases: SpringMover_Lineal

Movement of a mesh is determined by reinterpreting it as a structure of stiff beams and solving an associated discrete linear elasticity problem.

We consider the ‘torsional’ case, as described in [Farhat et al., 1998].

movement.tangling module

class MeshTanglingChecker(mesh, raise_error=True)

Bases: object

A class for tracking whether a mesh has tangled, i.e. whether any of its elements have become inverted.

property jacobian_determinant_ratio

Compute the ratio of the determinant of the Jacobian of the current mesh to the determinant of the Jacobian of the original mesh.

check()

Check whether any element orientations have changed since the tangling checker was created.

Module contents