polyquad 1.2.6

Frugal generation of quadrature for polytopal domain

polyquad

Frugal generation of quadrature for polytopal domains

Python package associated with this paper.

About

polyquad is a tool that generates quadratures (or cubatures) for non necessarily convex polytopal domains. Important points:

• the distribution of points is shape independent
• the number of integration points depends only on the polynomial order of the integrand
• expensive operations (QR decomposition) are done only once per polynomial order, their results are reused from one shape to the other
• integration points are defined on a bounding box encapsulating the polytope (and can thus fall outside of the polytope), so the integrand should be defined there as well.

Installation

The package is deployed on pypi, so it can be installed simply using pip

pip install polyquad

Note: this package depends on numba. It seems that, at the moment, it's not compatible with the latest python (3.13). Therefore you may need to revert to an older python for a while.

Usage

For a basic usage you have at your disposal 2 functions:

• polyquad.get_quadrature_2d
• polyquad.get_quadrature_3d

These should behave nicely as long as you feed them the right data structures. We will explain these just after. Before we need to talk about mapping.

Mapping

Both polyquad.get_quadrature_2d and polyquad.get_quadrature_3d come with an argument mapping which should be either True or False (False is default).

I you use this library, chances are you already have a bounding box concept in your code, thus your polytope is probably already lying in its bounding box if so put mapping = False (or don't specify anything as False is the default value) to avoid useless mappings.

On the contrary maybe your polytope isn't already in the reference bounding box, in that case you should specify mapping = True.

polyquad.get_quadrature_2d

As the name suggests, the function should be used to get quadratures over a 2d polygon. Let us tackle the case of a simple pentagon:

We first need to declare the coordinates of vertices

verts = np.array(((1,-1), (-1,0), (0,3), (2,3), (3,0)))

Then we need to specify the ordering of vertices, which is quite straightforward in that case

face = np.array((0,1,2,3,4))

Then, getting the quadrature for a polynomial order k is as simple as

points, weights = polyquad.get_quadrature_2d(k, verts, face, mapping = True)

polyquad.get_quadrature_3d

This one is for polyhedra. The only difference lies the declaration of faces. For sake of simplicity we here give the example of a simple pyramid (this library is intended to be used with way more complex polyhedra)

Again we first declare the coordinates of vertices

verts = np.array(((0,   0, 0),
                  (1,   0, 0),
                  (1,   1, 0),
                  (0,   1, 0),
                  (.5, .5, 1)))

Then we need to declare each planar face. A planar face is defined by the list of the index of vertices comprising it. For example the square base of the pyramid contains the four first vertices defined above, so it is defined as : [0,1,2,3] This is done for each face and gathered in a list of list as follows:

faces = [[0,1,2,3],
         [0,1,4],
         [1,2,4],
         [2,3,4],
         [3,0,4]]

Note: if the faces of the polyhedron all have the same number of vertices then faces can be casted to a np.array of shape (number of faces, number of vertices per face). This slightly speeds up the computation.

cite this work

If you use this code for you work, you are kindly invited to cite the associated paper

@article{langlois2024,
	title= {Frugal numerical integration scheme for polytopal domains},
	author= {Langlois,C. and van Putten,T. and B\'eriot,H. and Deckers,E.},
	journal= {Engineering with computers},
	doi={https://doi.org/10.1007/s00366-024-02080-1},
	url = {https://link.springer.com/article/10.1007/s00366-024-02080-1}
	}

---

Feedback and comments can be addressed to corresponding authors of the paper.