warp.fem#
The warp.fem
module is designed to facilitate solving physical systems described as differential
equations. For example, it can solve PDEs for diffusion, convection, fluid flow, and elasticity problems
using finite-element-based (FEM) Galerkin methods, and allows users to quickly experiment with various FEM
formulations and discretization schemes.
Integrands#
The core functionality of the FEM toolkit is the ability to integrate constant, linear, and bilinear forms over various domains and using arbitrary interpolation basis.
The main mechanism is the integrand()
decorator, for instance:
@integrand
def linear_form(
s: Sample,
v: Field,
):
return v(s)
@integrand
def diffusion_form(s: Sample, u: Field, v: Field, nu: float):
return nu * wp.dot(
grad(u, s),
grad(v, s),
)
Integrands are normal Warp kernels, meaning any usual Warp function can be used. However, they accept a few special parameters:
Sample
contains information about the current integration sample point, such as the element index and coordinates in element.
Field
designates an abstract field, which will be replaced at call time by the actual field type: for instance aDiscreteField
,field.TestField
orfield.TrialField
defined over an arbitraryFunctionSpace
. A field u can be evaluated at a given sample s using theinner()
operator, i.e,inner(u, s)
, or as a shortcut using the usual call operator,u(s)
. Several other operators are available, such asgrad()
; see the Operators section.
Domain
designates an abstract integration domain. Several operators are also provided for domains, for example in the example below evaluating the normal at the current sample position:@integrand def boundary_form( s: Sample, domain: Domain, u: Field, ): nor = normal(domain, s) return wp.dot(u(s), nor)
Integrands cannot be used directly with warp.launch()
, but must be called through integrate()
or interpolate()
instead.
The root integrand (integrand argument passed to integrate()
or interpolate()
call) will automatically get passed Sample
and Domain
parameters.
Field
parameters must be passed as a dictionary in the fields argument of the launcher function, and all other standard Warp types arguments must be
passed as a dictionary in the values argument of the launcher function, for instance:
integrate(diffusion_form, fields={"u": trial, "v": test}, values={"nu": viscosity})
Basic workflow#
The typical steps for solving a linear PDE are as follow:
Define a
Geometry
(grid, mesh, etc). At the moment, 2D and 3D regular grids, triangle, quadrilateral, tetrahedron and hexahedron meshes are supported.Define one or more
FunctionSpace
, by equipping the geometry elements with shape functions. Seemake_polynomial_space()
. At the moment, continuous/discontinuous Lagrange (\(P_{k[d]}, Q_{k[d]}\)) and Serendipity (\(S_k\)) shape functions of order \(k \leq 3\) are supported.Define an integration
GeometryDomain
, for instance the geometry’s cells (Cells
) or boundary sides (BoundarySides
).Integrate linear forms to build the system’s right-hand-side. Define a test function over the function space using
make_test()
, aQuadrature
formula (or let the module choose one based on the function space degree), and callintegrate()
with the linear form integrand. The result is awarp.array
containing the integration result for each of the function space degrees of freedom.Integrate bilinear forms to build the system’s left-hand-side. Define a trial function over the function space using
make_trial()
, then callintegrate()
with the bilinear form integrand. The result is awarp.sparse.BsrMatrix
containing the integration result for each pair of test and trial function space degrees of freedom. Note that the trial and test functions do not have to be defined over the same function space, so that Mixed FEM is supported.Solve the resulting linear system using the solver of your choice
The following excerpt from the introductory example warp/examples/fem/example_diffusion.py
outlines this procedure:
# Grid geometry
geo = Grid2D(n=50, cell_size=2)
# Domain and function spaces
domain = Cells(geometry=geo)
scalar_space = make_polynomial_space(geo, degree=3)
# Right-hand-side (forcing term)
test = make_test(space=scalar_space, domain=domain)
rhs = integrate(linear_form, fields={"v": test})
# Weakly-imposed boundary conditions on Y sides
boundary = BoundarySides(geo)
bd_test = make_test(space=scalar_space, domain=boundary)
bd_trial = make_trial(space=scalar_space, domain=boundary)
bd_matrix = integrate(y_mass_form, fields={"u": bd_trial, "v": bd_test})
# Diffusion form
trial = make_trial(space=scalar_space, domain=domain)
matrix = integrate(diffusion_form, fields={"u": trial, "v": test}, values={"nu": viscosity})
# Assemble linear system (add diffusion and boundary condition matrices)
bsr_axpy(x=bd_matrix, y=matrix, alpha=boundary_strength, beta=1)
# Solve linear system using Conjugate Gradient
x = wp.zeros_like(rhs)
bsr_cg(matrix, b=rhs, x=x)
Note
The integrate()
function does not check that the passed integrands are actually linear or bilinear forms; it is up to the user to ensure that they are.
To solve non-linear PDEs, one can use an iterative procedure and pass the current value of the studied function DiscreteField
argument to the integrand, on which
arbitrary operations are permitted. However, the result of the form must remain linear in the test and trial fields.
Introductory examples#
warp.fem
ships with a list of examples in the warp/examples/fem
directory illustrating common model problems.
example_diffusion.py
: 2D diffusion with homogeneous Neumann and Dirichlet boundary conditions
example_diffusion_3d.py
: 3D variant of the diffusion problem
example_convection_diffusion.py
: 2D convection-diffusion using semi-Lagrangian advection
example_diffusion_dg0.py
: 2D convection-diffusion using finite-volume and upwind transport
example_diffusion_dg.py
: 2D convection-diffusion using Discontinuous Galerkin with upwind transport and Symmetric Interior Penalty
example_stokes.py
: 2D incompressible Stokes flow using mixed \(P_k/P_{k-1}\) or \(Q_k/P_{(k-1)d}\) elements
example_navier_stokes.py
: 2D Navier-Stokes flow using mixed \(P_k/P_{k-1}\) elements
example_mixed_elasticity.py
: 2D linear elasticity using mixed continuous/discontinuous \(S_k/P_{(k-1)d}\) elements
Advanced usages#
High-order (curved) geometries#
It is possible to convert any Geometry
(grids and explicit meshes) into a curved, high-order variant by deforming them
with an arbitrary-order displacement field using the make_deformed_geometry()
method.
The process looks as follow:
# Define a base geometry
base_geo = fem.Grid3D(res=resolution)
# Define a displacement field on the base geometry
deformation_space = fem.make_polynomial_space(base_geo, degree=deformation_degree, dtype=wp.vec3)
deformation_field = deformation_space.make_field()
# Populate the field value by interpolating an expression
fem.interpolate(deformation_field_expr, dest=deformation_field)
# Construct the deformed geometry from the displacement field
deform_geo = deformation_field.make_deformed_geometry()
# Define new function spaces on the deformed geometry
scalar_space = fem.make_polynomial_space(deformed_geo, degree=scalar_space_degree)
See also example_deformed_geometry.py
for a complete example.
Particle-based quadrature#
The PicQuadrature
provides a way to define Particle-In-Cell quadratures from a set or arbitrary particles,
which can be helpful to develop MPM-like methods.
The particles are automatically bucketed to the geometry cells when the quadrature is initialized.
This is illustrated by the example_stokes_transfer.py
and example_apic_fluid.py
examples.
Partitioning#
The FEM toolkit makes it possible to perform integration on a subset of the domain elements,
possibly re-indexing degrees of freedom so that the linear system contains the local ones only.
This is useful for distributed computation (see warp/examples/fem/example_diffusion_mgpu.py
), or simply to limit the simulation domain to a subset of active cells (see warp/examples/fem/example_stokes_transfer.py
).
A partition of the simulation geometry can be defined using subclasses of GeometryPartition
such as LinearGeometryPartition
or ExplicitGeometryPartition
.
Function spaces can then be partitioned according to the geometry partition using make_space_partition()
.
The resulting SpacePartition
object allows translating between space-wide and partition-wide node indices,
and differentiating interior, frontier and exterior nodes.
Memory management#
Several warp.fem
functions require allocating temporary buffers to perform their computations.
If such functions are called many times in a tight loop, those many allocations and de-allocations may degrade performance.
To overcome this issue, a cache.TemporaryStore
object may be created to persist and reuse temporary allocations across calls,
either globally using set_default_temporary_store()
or at a per-function granularity using the corresponding argument.
Operators#
- warp.fem.normal(domain: Domain, s: Sample)#
Evaluates the element normal at the sample point s. Null for interior points.
- warp.fem.lookup(domain: Domain, x)#
Looks-up the sample point corresponding to a world position x, projecting to the closest point on the domain.
- Arg:
x: world position of the point to look-up in the geometry guess: (optional)
Sample
initial guess, may help perform the query
Notes
Currently this operator is only fully supported for
Grid2D
andGrid3D
geometries. ForTriangleMesh2D
andTetmesh
geometries, the operator requires providing guess.
- warp.fem.measure(domain: Domain, s: Sample)#
Returns the measure (volume, area, or length) determinant of an element at a sample point s
- warp.fem.measure_ratio(domain: Domain, s: Sample)#
Returns the maximum ratio between the measure of this element and that of higher-dimensional neighbours.
- warp.fem.deformation_gradient(domain: Domain, s: Sample)#
Evaluates the gradient of the domain position with respect to the element reference space at the sample point s
- warp.fem.inner(f: Field, s: Sample)#
Evaluates the field at a sample point s. On oriented sides, uses the inner element
- warp.fem.outer(f: Field, s: Sample)#
Evaluates the field at a sample point s. On oriented sides, uses the outer element. On interior points and on domain boundaries, this is equivalent to
inner()
.
- warp.fem.grad(f: Field, s: Sample)#
Evaluates the field gradient at a sample point s. On oriented sides, uses the inner element
- warp.fem.grad_outer(f: Field, s: Sample)#
Evaluates the field gradient at a sample point s. On oriented sides, uses the outer element. On interior points and on domain boundaries, this is equivalent to
grad()
.
- warp.fem.div(f: Field, s: Sample)#
Evaluates the field divergence at a sample point s. On oriented sides, uses the inner element
- warp.fem.div_outer(f: Field, s: Sample)#
Evaluates the field divergence at a sample point s. On oriented sides, uses the outer element. On interior points and on domain boundaries, this is equivalent to
div()
.
- warp.fem.at_node(f: Field, s: Sample)#
For a Test or Trial field, returns a copy of the Sample s moved to the coordinates of the node being evaluated
- warp.fem.curl(f: Field, s: Sample)#
Skew part of the (inner) gradient of the field at s, as a vector such that
wp.cross(curl(u), v) = skew(grad(u)) v
- warp.fem.jump(f: Field, s: Sample)#
Jump between inner and outer element values on an interior side. Zero for interior points or domain boundaries
Integration#
- warp.fem.integrate(integrand, domain=None, quadrature=None, nodal=False, fields=None, values=None, accumulate_dtype=<class 'warp.types.float64'>, output_dtype=None, output=None, device=None, temporary_store=None, kernel_options=None)#
Integrates a constant, linear or bilinear form, and returns a scalar, array, or sparse matrix, respectively.
- Parameters:
integrand (Integrand) – Form to be integrated, must have
integrand()
decoratordomain (GeometryDomain | None) – Integration domain. If None, deduced from fields
quadrature (Quadrature | None) – Quadrature formula. If None, deduced from domain and fields degree.
nodal (bool) – For linear or bilinear form only, use the test function nodes as the quadrature points. Assumes Lagrange interpolation functions are used, and no differential or DG operator is evaluated on the test or trial functions.
fields (Dict[str, FieldLike] | None) – Discrete, test, and trial fields to be passed to the integrand. Keys in the dictionary must match integrand parameter names.
values (Dict[str, Any] | None) – Additional variable values to be passed to the integrand, can be of any type accepted by warp kernel launches. Keys in the dictionary must match integrand parameter names.
temporary_store (TemporaryStore | None) – shared pool from which to allocate temporary arrays
accumulate_dtype (type) – Scalar type to be used for accumulating integration samples
output (array | BsrMatrix | None) – Sparse matrix or warp array into which to store the result of the integration
output_dtype (type | None) – Scalar type for returned results in output is not provided. If None, defaults to accumulate_dtype
device – Device on which to perform the integration
kernel_options (Dict[str, Any] | None) – Overloaded options to be passed to the kernel builder (e.g,
{"enable_backward": True}
)
- warp.fem.interpolate(integrand, dest=None, quadrature=None, fields=None, values=None, device=None, kernel_options=None)#
Interpolates a function at a finite set of sample points and optionally assigns the result to a discrete field or a raw warp array.
- Parameters:
integrand (Integrand) – Function to be interpolated, must have
integrand()
decoratordest (DiscreteField | FieldRestriction | array | None) –
Where to store the interpolation result. Can be either
a
DiscreteField
, or restriction of a discrete field to a domain (frommake_restriction()
). In this case, interpolation will be performed at each node.a normal warp array. In this case, the quadrature argument defining the interpolation locations must be provided and the result of the integrand at each quadrature point will be assigned to the array.
None
. In this case, the quadrature argument must also be provided and the integrand function is responsible for dealing with the interpolation result.
quadrature (Quadrature | None) – Quadrature formula defining the interpolation samples if dest is not a discrete field or field restriction.
fields (Dict[str, FieldLike] | None) – Discrete fields to be passed to the integrand. Keys in the dictionary must match integrand parameters names.
values (Dict[str, Any] | None) – Additional variable values to be passed to the integrand, can be of any type accepted by warp kernel launches. Keys in the dictionary must match integrand parameter names.
device – Device on which to perform the interpolation
kernel_options (Dict[str, Any] | None) – Overloaded options to be passed to the kernel builder (e.g,
{"enable_backward": True}
)
- warp.fem.integrand(func)#
Decorator for functions to be integrated (or interpolated) using warp.fem
- Parameters:
func (Callable) –
- class warp.fem.Sample#
Per-sample point context for evaluating fields and related operators in integrands.
- class warp.fem.Field#
Tag for field-like integrand arguments
- class warp.fem.Domain#
Tag for domain-like integrand arguments
Geometry#
- class warp.fem.Grid2D(res, bounds_lo=None, bounds_hi=None)#
Bases:
Geometry
Two-dimensional regular grid geometry
Constructs a dense 2D grid
- Parameters:
res (vec2i) – Resolution of the grid along each dimension
bounds_lo (vec2f | None) – Position of the lower bound of the axis-aligned grid
bounds_up – Position of the upper bound of the axis-aligned grid
bounds_hi (vec2f | None) –
- class warp.fem.Trimesh2D(tri_vertex_indices, positions, temporary_store=None)#
Bases:
Geometry
Two-dimensional triangular mesh geometry
Constructs a two-dimensional triangular mesh.
- Parameters:
tri_vertex_indices (array) – warp array of shape (num_tris, 3) containing vertex indices for each tri
positions (array) – warp array of shape (num_vertices, 2) containing 2d position for each vertex
temporary_store (TemporaryStore | None) – shared pool from which to allocate temporary arrays
- class warp.fem.Quadmesh2D(quad_vertex_indices, positions, temporary_store=None)#
Bases:
Geometry
Two-dimensional quadrilateral mesh geometry
Constructs a two-dimensional quadrilateral mesh.
- Parameters:
quad_vertex_indices (array) – warp array of shape (num_tris, 4) containing vertex indices for each quad, in counter-clockwise order
positions (array) – warp array of shape (num_vertices, 2) containing 2d position for each vertex
temporary_store (TemporaryStore | None) – shared pool from which to allocate temporary arrays
- class warp.fem.Grid3D(res, bounds_lo=None, bounds_hi=None)#
Bases:
Geometry
Three-dimensional regular grid geometry
Constructs a dense 3D grid
- Parameters:
res (vec3i) – Resolution of the grid along each dimension
bounds_lo (vec3f | None) – Position of the lower bound of the axis-aligned grid
bounds_up – Position of the upper bound of the axis-aligned grid
bounds_hi (vec3f | None) –
- class warp.fem.Tetmesh(tet_vertex_indices, positions, temporary_store=None)#
Bases:
Geometry
Tetrahedral mesh geometry
Constructs a tetrahedral mesh.
- Parameters:
tet_vertex_indices (array) – warp array of shape (num_tets, 4) containing vertex indices for each tet
positions (array) – warp array of shape (num_vertices, 3) containing 3d position for each vertex
temporary_store (TemporaryStore | None) – shared pool from which to allocate temporary arrays
- class warp.fem.Hexmesh(hex_vertex_indices, positions, temporary_store=None)#
Bases:
Geometry
Hexahedral mesh geometry
Constructs a tetrahedral mesh.
- Parameters:
hex_vertex_indices (array) – warp array of shape (num_hexes, 8) containing vertex indices for each hex following standard ordering (bottom face vertices in counter-clockwise order, then similarly for upper face)
positions (array) – warp array of shape (num_vertices, 3) containing 3d position for each vertex
temporary_store (TemporaryStore | None) – shared pool from which to allocate temporary arrays
- class warp.fem.LinearGeometryPartition(geometry, partition_rank, partition_count, device=None, temporary_store=None)#
Creates a geometry partition by uniformly partionning cell indices
- Parameters:
geometry (Geometry) – the geometry to partition
partition_rank (int) – the index of the partition being created
partition_count (int) – the number of partitions that will be created over the geometry
device – Warp device on which to perform and store computations
temporary_store (TemporaryStore) –
- class warp.fem.ExplicitGeometryPartition(geometry, cell_mask, temporary_store=None)#
Creates a geometry partition by uniformly partionning cell indices
- Parameters:
geometry (Geometry) – the geometry to partition
cell_mask (wp.array(dtype=int)) – warp array of length
geometry.cell_count()
indicating which cells are selected. Array values must be either1
(selected) or0
(not selected).temporary_store (TemporaryStore) –
- class warp.fem.Cells(geometry)#
Bases:
GeometryDomain
A Domain containing all cells of the geometry or geometry partition
- Parameters:
geometry (Geometry | GeometryPartition) –
- class warp.fem.Sides(geometry)#
Bases:
GeometryDomain
A Domain containing all (interior and boundary) sides of the geometry or geometry partition
- Parameters:
geometry (Geometry | GeometryPartition) –
- class warp.fem.BoundarySides(geometry)#
Bases:
Sides
A Domain containing boundary sides of the geometry or geometry partition
- Parameters:
geometry (Geometry | GeometryPartition) –
- class warp.fem.FrontierSides(geometry)#
Bases:
Sides
A Domain containing frontier sides of the geometry partition (sides shared with at least another partition)
- Parameters:
geometry (Geometry | GeometryPartition) –
- class warp.fem.Polynomial(value)#
Polynomial family defining interpolation nodes over an interval
- GAUSS_LEGENDRE = 0#
Gauss–Legendre 1D polynomial family (does not include endpoints)
- LOBATTO_GAUSS_LEGENDRE = 1#
Lobatto–Gauss–Legendre 1D polynomial family (includes endpoints)
- EQUISPACED_CLOSED = 2#
Closed 1D polynomial family with uniformly distributed nodes (includes endpoints)
- EQUISPACED_OPEN = 3#
Open 1D polynomial family with uniformly distributed nodes (does not include endpoints)
- class warp.fem.RegularQuadrature(domain, order, family=None)#
Bases:
Quadrature
Regular quadrature formula, using a constant set of quadrature points per element
- Parameters:
domain (GeometryDomain) –
order (int) –
family (Polynomial) –
- class warp.fem.NodalQuadrature(domain, space)#
Bases:
Quadrature
Quadrature using space node points as quadrature points
Note that in contrast to the nodal=True flag for
integrate()
, this quadrature odes not make any assumption about orthogonality of shape functions, and is thus safe to use for arbitrary integrands.- Parameters:
domain (GeometryDomain) –
space (FunctionSpace) –
- class warp.fem.ExplicitQuadrature(domain, points, weights)#
Bases:
Quadrature
Quadrature using explicit per-cell points and weights. The number of quadrature points per cell is assumed to be constant and deduced from the shape of the points and weights arrays.
- Parameters:
domain (GeometryDomain) – Domain of definition of the quadrature formula
points (wp.array2d(dtype=Coords)) – 2d array of shape
(domain.geometry_element-count(), points_per_cell)
containing the coordinates of each quadrature point.weights (wp.array2d(dtype=float)) – 2d array of shape
(domain.geometry_element-count(), points_per_cell)
containing the weight for each quadrature point.
See also:
PicQuadrature
- class warp.fem.PicQuadrature(domain, positions, measures=None, temporary_store=None)#
Bases:
Quadrature
Particle-based quadrature formula, using a global set of points unevenly spread out over geometry elements.
Useful for Particle-In-Cell and derived methods.
- Parameters:
domain (GeometryDomain) – Underlying domain for the quadrature
positions (wp.array(dtype=wp.vecXd) | Tuple[wp.array(dtype=ElementIndex), wp.array(dtype=Coords)]) – Either an array containing the world positions of all particles, or a tuple of arrays containing the cell indices and coordinates for each particle. Note that the former requires the underlying geometry to define a global
Geometry.cell_lookup()
method; currently this is only available forGrid2D
andGrid3D
.measures (wp.array(dtype=float) | None) – Array containing the measure (area/volume) of each particle, used to defined the integration weights. If
None
, defaults to the cell measure divided by the number of particles in the cell.temporary_store (TemporaryStore) – shared pool from which to allocate temporary arrays
Function Spaces#
- warp.fem.make_polynomial_space(geo, dtype=<class 'float'>, dof_mapper=None, degree=1, element_basis=None, discontinuous=False, family=None)#
Equips a geometry with a collocated, polynomial function space. Equivalent to successive calls to
make_polynomial_basis_space()
and make_collocated_function_space.- Parameters:
geo (Geometry) – the Geometry on which to build the space
dtype (type) – value type the function space. If
dof_mapper
is provided, the value type from the DofMapper will be used instead.dof_mapper (DofMapper | None) – mapping from node degrees of freedom to function values, defaults to Identity. Useful for reduced coordinates, e.g.
SymmetricTensorMapper
maps 2x2 (resp 3x3) symmetric tensors to 3 (resp 6) degrees of freedom.degree (int) – polynomial degree of the per-element shape functions
discontinuous (bool) – if True, use Discontinuous Galerkin shape functions. Discontinuous is implied if degree is 0, i.e, piecewise-constant shape functions.
element_basis (ElementBasis | None) – type of basis function for the individual elements
family (Polynomial | None) – Polynomial family used to generate the shape function basis. If not provided, a reasonable basis is chosen.
- Returns:
the constructed function space
- Return type:
CollocatedFunctionSpace
- warp.fem.make_polynomial_basis_space(geo, degree=1, element_basis=None, discontinuous=False, family=None)#
Equips a geometry with a polynomial basis.
- Parameters:
geo (Geometry) – the Geometry on which to build the space
degree (int) – polynomial degree of the per-element shape functions
discontinuous (bool) – if True, use Discontinuous Galerkin shape functions. Discontinuous is implied if degree is 0, i.e, piecewise-constant shape functions.
element_basis (ElementBasis | None) – type of basis function for the individual elements
family (Polynomial | None) – Polynomial family used to generate the shape function basis. If not provided, a reasonable basis is chosen.
- Returns:
the constructed basis space
- Return type:
- warp.fem.make_collocated_function_space(basis_space, dtype=<class 'float'>, dof_mapper=None)#
Constructs a function space from a basis space and a value type, such that all degrees of freedom of the value type are stored at each of the basis nodes.
- Parameters:
geo – the Geometry on which to build the space
dtype (type) – value type the function space. If
dof_mapper
is provided, the value type from the DofMapper will be used instead.dof_mapper (DofMapper | None) – mapping from node degrees of freedom to function values, defaults to Identity. Useful for reduced coordinates, e.g.
SymmetricTensorMapper
maps 2x2 (resp 3x3) symmetric tensors to 3 (resp 6) degrees of freedom.basis_space (BasisSpace) –
- Returns:
the constructed function space
- Return type:
CollocatedFunctionSpace
- warp.fem.make_space_partition(space=None, geometry_partition=None, space_topology=None, with_halo=True, device=None, temporary_store=None)#
Computes the subset of nodes from a function space topology that touch a geometry partition
Either space_topology or space must be provided (and will be considered in that order).
- Parameters:
space (FunctionSpace | None) – (deprecated) the function space defining the topology if space_topology is
None
.geometry_partition (GeometryPartition | None) – The subset of the space geometry. If not provided, use the whole geometry.
space_topology (SpaceTopology | None) – the topology of the function space to consider. If
None
, deduced from space.with_halo (bool) – if True, include the halo nodes (nodes from exterior frontier cells to the partition)
device – Warp device on which to perform and store computations
temporary_store (TemporaryStore | None) –
- Returns:
the resulting space partition
- Return type:
- warp.fem.make_space_restriction(space=None, space_partition=None, domain=None, space_topology=None, device=None, temporary_store=None)#
Restricts a function space partition to a Domain, i.e. a subset of its elements.
One of space_partition, space_topology, or space must be provided (and will be considered in that order).
- Parameters:
space (FunctionSpace | None) – (deprecated) if neither space_partition nor space_topology are provided, the space defining the topology to restrict
space_partition (SpacePartition | None) – the subset of nodes from the space topology to consider
domain (GeometryDomain | None) – the domain to restrict the space to, defaults to all cells of the space geometry or partition.
space_topology (SpaceTopology | None) – the space topology to be restricted, if space_partition is
None
.device – device on which to perform and store computations
temporary_store (Optional[warp.fem.cache.TemporaryStore]) – shared pool from which to allocate temporary arrays
- Return type:
- class warp.fem.ElementBasis(value)#
Choice of basis function to equip individual elements
- LAGRANGE = 0#
Lagrange basis functions \(P_k\) for simplices, tensor products \(Q_k\) for squares and cubes
- SERENDIPITY = 1#
Serendipity elements \(S_k\), corresponding to Lagrange nodes with interior points removed (for degree <= 3)
- NONCONFORMING_POLYNOMIAL = 2#
Simplex Lagrange basis functions \(P_{kd}\) embedded into non conforming reference elements (e.g. squares or cubes). Discontinuous only.
- class warp.fem.SymmetricTensorMapper(dtype, mapping=Mapping.VOIGT)#
Bases:
DofMapper
Orthonormal isomorphism from R^{n (n+1)} to nxn symmetric tensors, using usual L2 norm for vectors and half Frobenius norm, (tau : tau)/2 for tensors.
- Parameters:
dtype (type) –
mapping (Mapping) –
- class warp.fem.SkewSymmetricTensorMapper(dtype)#
Bases:
DofMapper
Orthonormal isomorphism from R^{n (n-1)} to nxn skew-symmetric tensors, using usual L2 norm for vectors and half Frobenius norm, (tau : tau)/2 for tensors.
- Parameters:
dtype (type) –
- class warp.fem.PointBasisSpace(quadrature)#
Bases:
BasisSpace
An unstructured
BasisSpace
that is non-zero at a finite set of points only.The node locations and nodal quadrature weights are defined by a
Quadrature
formula.- Parameters:
quadrature (Quadrature) –
Fields#
- warp.fem.make_test(space, space_restriction=None, space_partition=None, domain=None, device=None)#
Constructs a test field over a function space or its restriction
- Parameters:
space (FunctionSpace) – the function space
space_restriction (SpaceRestriction | None) – restriction of the space topology to a domain
space_partition (SpacePartition | None) – if space_restriction is
None
, the optional subset of node indices to considerdomain (GeometryDomain | None) – if space_restriction is
None
, optional subset of elements to considerdevice – Warp device on which to perform and store computations
- Returns:
the test field
- Return type:
- warp.fem.make_trial(space, space_restriction=None, space_partition=None, domain=None)#
Constructs a trial field over a function space or partition
- Parameters:
space (FunctionSpace) – the function space or function space restriction
space_restriction (SpaceRestriction | None) – restriction of the space topology to a domain
space_partition (SpacePartition | None) – if space_restriction is
None
, the optional subset of node indices to considerdomain (GeometryDomain | None) – if space_restriction is
None
, optional subset of elements to considerdevice – Warp device on which to perform and store computations
- Returns:
the trial field
- Return type:
- warp.fem.make_restriction(field, space_restriction=None, domain=None, device=None)#
Restricts a discrete field to a subset of elements.
- Parameters:
field (DiscreteField) – the discrete field to restrict
space_restriction (SpaceRestriction | None) – the function space restriction defining the subset of elements to consider
domain (GeometryDomain | None) – if
space_restriction
is not provided, theDomain
defining the subset of elements to considerdevice – Warp device on which to perform and store computations
- Returns:
the field restriction
- Return type:
Boundary Conditions#
- warp.fem.normalize_dirichlet_projector(projector_matrix, fixed_value=None)#
Scale projector so that it becomes idempotent, and apply the same scaling to fixed_value if provided
- warp.fem.project_linear_system(system_matrix, system_rhs, projector_matrix, fixed_value=None, normalize_projector=True)#
Projects both the left-hand-side and right-hand-side of a linear system to enforce Dirichlet boundary conditions
If normalize_projector is True, first apply scaling so that the projector_matrix is idempotent
Memory management#
- warp.fem.set_default_temporary_store(temporary_store)#
Globally sets the default
TemporaryStore
instance to use for temporary allocations inwarp.fem
functions.If the default temporary store is set to
None
, temporary allocations are not persisted unless aTemporaryStore
is provided at a per-function granularity.- Parameters:
temporary_store (TemporaryStore | None) –
- warp.fem.borrow_temporary(temporary_store, shape, dtype, pinned=False, requires_grad=False, device=None)#
Borrows and returns a temporary array with specified attributes from a shared pool.
If an array with sufficient capacity and matching desired attributes is already available in the pool, it will be returned. Otherwise, a new allocation will be performed.
- Parameters:
temporary_store (TemporaryStore | None) – the shared pool to borrow the temporary from. If temporary_store is
None
, the global default temporary store, if set, will be used.shape (int | Tuple[int]) – desired dimensions for the temporary array
dtype (type) – desired data type for the temporary array
pinned (bool) – whether a pinned allocation is desired
device – device on which the memory should be allocated; if
None
, the current device will be used.requires_grad (bool) –
- Return type:
- warp.fem.borrow_temporary_like(array, temporary_store)#
Borrows and returns a temporary array with the same attributes as another array or temporary.
- Parameters:
array (array | Temporary) – Warp or temporary array to read the desired attributes from
temporary_store (TemporaryStore | None) – the shared pool to borrow the temporary from. If temporary_store is
None
, the global default temporary store, if set, will be used.
- Return type:
Interfaces#
Interface classes are not meant to be constructed directly, but can be derived from extend the built-in functionality.
- class warp.fem.Geometry#
Interface class for discrete geometries
A geometry is composed of cells and sides. Sides may be boundary or interior (between cells).
- cell_count()#
Number of cells in the geometry
- side_count()#
Number of sides in the geometry
- boundary_side_count()#
Number of boundary sides (sides with a single neighbour cell) in the geometry
- class warp.fem.GeometryPartition(geometry)#
Base class for geometry partitions, i.e. subset of cells and sides
- Parameters:
geometry (Geometry) –
- boundary_side_count()#
Number of geo-boundary sides that are ‘owned’ by this partition
- Return type:
- class warp.fem.GeometryDomain(geometry)#
Interface class for domains, i.e. (partial) views of elements in a Geometry
- Parameters:
geometry (Geometry | GeometryPartition) –
- class ElementKind(value)#
Possible kinds of elements contained in a domain
- property element_kind: ElementKind#
Kind of elements that this domain contains (cells or sides)
- class warp.fem.Quadrature(domain)#
Interface class for quadrature rules
- Parameters:
domain (GeometryDomain) –
- property domain#
Domain over which this quadrature is defined
- total_point_count()#
Total number of quadrature points over the domain
- class warp.fem.FunctionSpace(topology)#
Interface class for function spaces, i.e. geometry + interpolation basis
- Parameters:
topology (SpaceTopology) –
- property topology: SpaceTopology#
Underlying geometry
- trace()#
Trace of the function space over lower-dimensional elements of the geometry
- Return type:
- make_field(space_partition=None)#
Creates a zero-initialized discrete field over the function space holding values for all degrees of freedom of nodes in a space partition
- space_arg:
space_partition: If provided, the subset of nodes to consider
See also:
make_space_partition()
- class warp.fem.SpaceTopology(geometry, nodes_per_element)#
Interface class for defining the topology of a function space.
The topology only considers the indices of the nodes in each element, and as such, the connectivity pattern of the function space. It does not specify the actual location of the nodes within the elements, or the valuation function.
- element_node_indices(out=None)#
Returns a temporary array containing the global index for each node of each element
- trace()#
Trace of the function space over lower-dimensional elements of the geometry
- Return type:
TraceSpaceTopology
- class warp.fem.BasisSpace(topology)#
Interface class for defining a scalar-valued basis over a geometry.
A basis space makes it easy to define multiple function spaces sharing the same basis (and thus nodes) but with different valuation functions; however, it is not a required ingredient of a function space.
See also:
make_polynomial_basis_space()
,make_collocated_function_space()
- Parameters:
topology (SpaceTopology) –
- property topology: SpaceTopology#
Underlying topology of the basis space
- class warp.fem.space.shape.ShapeFunction#
Interface class for defining scalar-valued shape functions over a single element
- class warp.fem.SpacePartition(space_topology, geo_partition)#
- Parameters:
space_topology (SpaceTopology) –
geo_partition (GeometryPartition) –
- node_count()#
Returns number of nodes in this partition
- owned_node_count()#
Returns number of nodes in this partition, excluding exterior halo
- Return type:
- class warp.fem.SpaceRestriction(space_partition, domain, device=None, temporary_store=None)#
Restriction of a space partition to a given GeometryDomain
- Parameters:
space_partition (SpacePartition) –
domain (GeometryDomain) –
temporary_store (TemporaryStore) –
- class warp.fem.DofMapper#
Base class from mapping node degrees of freedom to function values
- class warp.fem.FieldLike#
Base class for integrable fields
- class warp.fem.DiscreteField(space, space_partition)#
Bases:
SpaceField
Explicitly-valued field defined over a partition of a discrete function space
- Parameters:
space (FunctionSpace) –
space_partition (SpacePartition) –
- trace()#
Trace of this field over a lower-dimensional function space
- Return type:
- class warp.fem.field.FieldRestriction(space_restriction, field)#
Restriction of a discrete field to a given GeometryDomain
- Parameters:
space_restriction (SpaceRestriction) –
field (DiscreteField) –
- class warp.fem.field.SpaceField(space, space_partition)#
Bases:
FieldLike
Base class for fields defined over a function space
- Parameters:
space (FunctionSpace) –
space_partition (SpacePartition) –
- class warp.fem.field.TestField(space_restriction, space)#
Bases:
SpaceField
Field defined over a space restriction that can be used as a test function.
In order to reuse computations, it is possible to define the test field using a SpaceRestriction defined for a different value type than the test function value type, as long as the node topology is similar.
- Parameters:
space_restriction (SpaceRestriction) –
space (FunctionSpace) –
- class warp.fem.field.TrialField(space, space_partition, domain)#
Bases:
SpaceField
Field defined over a domain that can be used as a trial function
- Parameters:
space (FunctionSpace) –
space_partition (SpacePartition) –
domain (GeometryDomain) –
- class warp.fem.TemporaryStore#
Shared pool of temporary arrays that will be persisted and reused across invocations of
warp.fem
functions.A
TemporaryStore
instance may either be passed explicitly towarp.fem
functions that accept such an argument, for instanceintegrate.integrate()
, or can be set globally as the default store usingset_default_temporary_store()
.By default, there is no default temporary store, so that temporary allocations are not persisted.
- class warp.fem.cache.Temporary(array, pool=None, shape=None, dtype=None)#
Handle over a temporary array from a
TemporaryStore
.The array will be automatically returned to the temporary pool for reuse upon destruction of this object, unless the temporary is explicitly detached from the pool using
detach()
. The temporary may also be explicitly returned to the pool before destruction usingrelease()
.- Parameters:
array (array) –
pool (TemporaryStore.Pool | None) –
- release()#
Returns the temporary array to the pool