Skip to content

Commit

Permalink
Moved functionality to integrate module
Browse files Browse the repository at this point in the history
  • Loading branch information
jzwar committed Jul 7, 2023
1 parent a153f35 commit 1245a1f
Show file tree
Hide file tree
Showing 6 changed files with 191 additions and 99 deletions.
58 changes: 0 additions & 58 deletions splinepy/bezier.py
Original file line number Diff line number Diff line change
Expand Up @@ -248,64 +248,6 @@ def extract_dim(self, dim):

return type(self)(spline=splinepy_core.extract_dim(self, dim))

def volume(self):
r"""Integrates volume
Determinante has degree
.. math::
p_i^{det} = n_{dim} \cdot p_i - 1
cf. [Mantzaflaris et al., 2017,
DOI:http://dx.doi.org/10.1016/j.cma.2016.11.013]
Parameters
----------
None
Returns
-------
volume : float
Integral of dim-dimensional object
"""
from splinepy.utils.data import cartesian_product

# Determine integration points
positions = []
weights = []

# Check dimensionality
if not (self.dim == self.para_dim):
raise ValueError(
"`Volume` of embedded spline depends on projection, "
"integration is aborted"
)

# Determine integration orders
expected_degree = self.degrees * self.para_dim + 1
if self.is_rational:
expected_degree += 2
self._logw("Integration on rational spline is only approximation")

# Gauss-Legendre is exact for polynomials 2*n-1
orders = np.ceil((expected_degree + 1) * 0.5).astype("int")

for order in orders:
# Get legendre quadratuer points
pos, ws = np.polynomial.legendre.leggauss(order)
# from (-1,1) to (0,1)
positions.append(pos * 0.5 + 0.5)
weights.append(ws * 0.5)

# summarize integration points
positions = cartesian_product(positions)
weights = cartesian_product(weights).prod(axis=1)

# Calculate Volume
volume = np.sum(np.linalg.det(self.jacobian(positions)) * weights)
return volume


class Bezier(BezierBase):
r"""
Expand Down
19 changes: 0 additions & 19 deletions splinepy/bspline.py
Original file line number Diff line number Diff line change
Expand Up @@ -305,25 +305,6 @@ def extract_bezier_patches(self):
# correct types
return [settings.NAME_TO_TYPE[p.name](spline=p) for p in patches]

def volume(self):
"""Integrate volume
@todo : Could be more efficient if no bezier extraction would take
place
Check out documentation for BezierBase.volume for more information
Parameters
----------
None
Returns
-------
volume : float
Integral of dim-dimensional object
"""
beziers = self.extract_bezier_patches()
return np.sum([b.volume() for b in beziers])


class BSpline(BSplineBase):
r"""
Expand Down
18 changes: 16 additions & 2 deletions splinepy/helpme/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,17 @@
from splinepy.helpme import create, extract, mapper, multi_index, visualize
from splinepy.helpme import (
create,
extract,
integrate,
mapper,
multi_index,
visualize,
)

__all__ = ["multi_index", "mapper", "create", "extract", "visualize"]
__all__ = [
"create",
"extract",
"integrate",
"mapper",
"multi_index",
"visualize",
]
109 changes: 109 additions & 0 deletions splinepy/helpme/integrate.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
import numpy as np


def volume(spline, orders=None):
r"""Compute volume of a given spline
Determinante has degree
.. math::
p_i^{det} = n_{dim} \cdot p_i - 1
cf. [Mantzaflaris et al., 2017,
DOI:http://dx.doi.org/10.1016/j.cma.2016.11.013]
Parameters
----------
spline : Spline to be integrated
splinepy - spline type
orders : array-like (optional)
order for gauss quadrature
Returns
-------
volume : float
Integral of dim-dimensional object
"""
from splinepy.spline import Spline
from splinepy.utils.data import cartesian_product

# Check input type
if not isinstance(spline, Spline):
raise NotImplementedError("Extrude only works for splines")

# Determine integration points
positions = []
weights = []

# Check dimensionality
if not (spline.dim == spline.para_dim):
raise ValueError(
"`Volume` of embedded spline depends on projection, "
"integration is aborted"
)

# Determine integration orders
if orders is None:
expected_degree = spline.degrees * spline.para_dim + 1
if spline.is_rational:
expected_degree += 2
spline._logw(
"Integration on rational spline is only approximation"
)

# Gauss-Legendre is exact for polynomials 2*n-1
quad_orders = np.ceil((expected_degree + 1) * 0.5).astype("int")
else:
quad_orders = np.ascontiguousarray(orders, dtype=int).flatten()
if quad_orders.size != spline.para_dim:
raise ValueError(
"Integration order must be array of size para_dim"
)

for order in quad_orders:
# Get legendre quadratuer points
pos, ws = np.polynomial.legendre.leggauss(order)
# from (-1,1) to (0,1)
positions.append(pos * 0.5 + 0.5)
weights.append(ws * 0.5)

# summarize integration points
positions = cartesian_product(positions)
weights = cartesian_product(weights).prod(axis=1)

# Calculate Volume
if spline.has_knot_vectors:
volume = np.sum(
[
np.sum(np.linalg.det(b.jacobian(positions)) * weights)
for b in spline.extract.beziers()
]
)
else:
volume = np.sum(np.linalg.det(spline.jacobian(positions)) * weights)
return volume


class Integrator:
"""Helper class to integrate some values on a given spline geometry
Examples
--------
.. code-block:: python
splinepy.helpme.integrate.volumer()
# Equivalent to
spline.integrate.volume()
Parameters
----------
spline : Spline
Spline parend
"""

def __init__(self, spl):
self.spline = spl

def volume(self, *args, **kwargs):
return volume(self.spline, *args, **kwargs)
31 changes: 30 additions & 1 deletion splinepy/spline.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
from splinepy.helpme import visualize
from splinepy.helpme.create import Creator
from splinepy.helpme.extract import Extractor
from splinepy.helpme.integrate import Integrator
from splinepy.utils.data import SplineData


Expand Down Expand Up @@ -652,7 +653,13 @@ class Spline(SplinepyBase, core.PySpline):
None
"""

__slots__ = ("_extractor", "_creator", "_show_options", "_spline_data")
__slots__ = (
"_extractor",
"_creator",
"_integrator",
"_show_options",
"_spline_data",
)

__show_option__ = visualize.SplineShowOption

Expand Down Expand Up @@ -828,6 +835,28 @@ def extract(self):
"""
return _get_helper(self, "_extractor", Extractor)

@property
def integrate(self):
"""Returns spline integrator. Can directly perform integrations
available at `splinepy/helpme/integrate.py`.
Examples
---------
.. code-block :: python
spline_faces = spline.integrate.volume()
Parameters
-----------
None
Returns
--------
integrator: Integrator
"""
return _get_helper(self, "_integrator", Integrator)

@property
def create(self):
"""Returns spline creator Can be used to create new splines using
Expand Down
55 changes: 36 additions & 19 deletions tests/test_volume_integration.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,22 +13,24 @@ def test_volume_integration_1D(self):
# Test 1D
bezier = c.splinepy.Bezier(degrees=[1], control_points=[[0], [2]])

self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

# test for other types same spline
self.assertTrue(c.np.allclose(bezier.bspline.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.rationalbezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.nurbs.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.bspline.integrate.volume(), 2.0))
self.assertTrue(
c.np.allclose(bezier.rationalbezier.integrate.volume(), 2.0)
)
self.assertTrue(c.np.allclose(bezier.nurbs.integrate.volume(), 2.0))

# Check if equal after refinement
bezier.elevate_degrees([0, 0, 0])

self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

# Test knot insertion
bspline = bezier.bspline
bspline.insert_knots(0, c.np.random.rand(10))
self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

def test_volume_integration_2D(self):
"""
Expand All @@ -40,27 +42,29 @@ def test_volume_integration_2D(self):
degrees=[1, 1], control_points=[[0, 0], [2, 0], [0, 1], [2, 1]]
)

self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

# test for other types same spline
self.assertTrue(c.np.allclose(bezier.bspline.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.rationalbezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.nurbs.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.bspline.integrate.volume(), 2.0))
self.assertTrue(
c.np.allclose(bezier.rationalbezier.integrate.volume(), 2.0)
)
self.assertTrue(c.np.allclose(bezier.nurbs.integrate.volume(), 2.0))

# Check if equal after refinement
bezier.elevate_degrees([0, 0, 1])

self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

# Test knot insertion
bspline = bezier.bspline
bspline.insert_knots(0, c.np.random.rand(10))
self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

# Move control points along axis does not change volume
mi = bspline.multi_index
bspline.cps[mi[3, :]][:, 1] += 10
self.assertTrue(c.np.allclose(bezier.volume(), 2.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 2.0))

def test_volume_integration_3D(self):
"""
Expand All @@ -82,22 +86,22 @@ def test_volume_integration_3D(self):
],
)

self.assertTrue(c.np.allclose(bezier.volume(), 3.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 3.0))

# Check if equal after refinement
bezier.elevate_degrees([0, 0, 1, 1, 2, 2])

self.assertTrue(c.np.allclose(bezier.volume(), 3.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 3.0))

# Test knot insertion
bspline = bezier.bspline
bspline.insert_knots(0, c.np.random.rand(10))
self.assertTrue(c.np.allclose(bezier.volume(), 3.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 3.0))

# Move control points along axis does not change volume
mi = bspline.multi_index
bspline.cps[mi[3, :, :]][:, 1] += 10
self.assertTrue(c.np.allclose(bezier.volume(), 3.0))
self.assertTrue(c.np.allclose(bezier.integrate.volume(), 3.0))

def test_complex_geometry(self):
"""Test on a more complex geometry"""
Expand Down Expand Up @@ -125,7 +129,11 @@ def test_complex_geometry(self):
bezier.cps[mi[1:3, 1:3, 1:3]][:] += 0.05 * c.np.random.rand(
*bezier.cps[mi[1:3, 1:3, 1:3]][:].shape
)
self.assertTrue(c.np.allclose(bezier.volume(), bezier_c.volume()))
self.assertTrue(
c.np.allclose(
bezier.integrate.volume(), bezier_c.integrate.volume()
)
)

def test_assertions(self):
"""Test the assertions in volume function"""
Expand All @@ -136,7 +144,16 @@ def test_assertions(self):
ValueError,
"`Volume` of embedded spline depends on projection, "
"integration is aborted",
bezier.volume,
bezier.integrate.volume,
)
bezier = c.splinepy.Bezier(
degrees=[2, 2], control_points=c.np.random.rand(9, 2)
)
self.assertRaisesRegex(
ValueError,
"Integration order must be array of size para_dim",
bezier.integrate.volume,
orders=[2, 4, 5],
)


Expand Down

0 comments on commit 1245a1f

Please sign in to comment.