From 618bf9fa2f5310c10011913afb19219bcf083cf3 Mon Sep 17 00:00:00 2001 From: carlos-adir Date: Tue, 5 Mar 2024 22:27:24 +0100 Subject: [PATCH] refactor: pre-compute mesh on adaptative integration --- src/compmec/section/integral.py | 101 ++++++++++++++------------------ tests/test_geomprop.py | 13 ++-- 2 files changed, 51 insertions(+), 63 deletions(-) diff --git a/src/compmec/section/integral.py b/src/compmec/section/integral.py index 7b7673b..5b6cb35 100644 --- a/src/compmec/section/integral.py +++ b/src/compmec/section/integral.py @@ -403,7 +403,7 @@ def polygon(vertices: Tuple[Tuple[float]]) -> Tuple[float]: return geomprops @staticmethod - def adaptative(curve) -> Tuple[float]: + def adaptative(curve, tolerance: float = 1e-9) -> Tuple[float]: """ Computes the polynomials integrals over the area defined by the curve @@ -416,15 +416,13 @@ def adaptative(curve) -> Tuple[float]: :type curve: Curve :return: The integral values :rtype: float + """ - tolerance = 1e-6 expoents = Polynomial.expoents integrals = np.zeros(len(expoents), dtype="float64") - integrator = AdaptativePolynomialIntegrator(curve) + integrator = AdaptativePolynomialIntegrator(curve, tolerance=tolerance) for k, (expx, expy) in enumerate(expoents): - integrator.expx = expx - integrator.expy = expy - value = integrator.integrate(tolerance=tolerance) + value = integrator.integrate(expx, expy) integrals[k] = value / (2 + expx + expy) return integrals @@ -437,13 +435,40 @@ class AdaptativePolynomialIntegrator: using open newton cotes quadrature with 3 points """ - def __init__(self, curve, *, max_depth: int = 9): + integ_matrix = ( + np.array([[0, 4, -1], [-4, 0, 4], [1, -4, 0]], dtype="float64") / 3 + ) + + def __init__(self, curve, *, max_depth: int = 9, tolerance=1e-9): self.curve = curve - self.expx = 0 - self.expy = 0 self.max_depth = max_depth - def milne_formula(self, t0: float, t1: float) -> float: + mesh = set() + for t0, t1 in zip(curve.knots, curve.knots[1:]): + subtol = tolerance * (t1 - t0) / (curve.knots[-1] - curve.knots[0]) + mesh |= self.find_mesh(t0, t1, tolerance=subtol) + self.mesh = tuple(sorted(mesh)) + + def find_mesh( + self, t0: float, t1: float, tolerance: float = 1e-9 + ) -> Tuple[float]: + """ + Finds the curve's subdivisions such the area computed + by the adaptative subdivision is lower than the tolerance + """ + subknots = np.linspace(t0, t1, 5) + xvals, yvals = np.transpose(self.curve.eval(subknots)) + area_mid = xvals[::2] @ self.integ_matrix @ yvals[::2] + area_lef = xvals[:3] @ self.integ_matrix @ yvals[:3] + area_rig = xvals[2:] @ self.integ_matrix @ yvals[2:] + diff = abs(area_lef + area_rig - area_mid) + mesh = {t0, t1} + if diff > tolerance: + mesh |= self.find_mesh(t0, subknots[2], tolerance / 2) + mesh |= self.find_mesh(subknots[2], t1, tolerance / 2) + return mesh + + def milne_formula(self, t0: float, t1: float, a: int, b: int) -> float: """ Computes the integral using Milne formula @@ -458,11 +483,14 @@ def milne_formula(self, t0: float, t1: float) -> float: :type t0: float :param t1: The upper interval's value :type t1: float + :param a: The x expoent + :type a: int + :param b: The y expoent + :type b: int :return: The interval's integral value :rtype: float """ - a, b = self.expx, self.expy ts = np.linspace(t0, t1, 5) points = self.curve.eval(ts) xs = points[:, 0] @@ -475,48 +503,7 @@ def milne_formula(self, t0: float, t1: float) -> float: f3 *= xs[3] * (ys[4] - ys[2]) - ys[3] * (xs[4] - xs[2]) return (4 * f1 - 2 * f2 + 4 * f3) / 3 - def integrate_interval( - self, t0: float, t1: float, tolerance: float, depth: int - ) -> float: - """ - Computes the integral - - int_{t0}^{t1} x^expx * y^expy * (x * dy - y * dx) - - with (x(t), y(t)) being the curve - - This function is recursive - - Parameters - ---------- - - :param t0: The lower interval's value - :type t0: float - :param t1: The upper interval's value - :type t1: float - :param tolerance: The stop tolerance - :type tolerance: float - :param depth: The current depth, a key to know when stop - :type depth: int - :return: The interval's integral value - :rtype: float - - """ - tm = 0.5 * (t0 + t1) - value_midl = self.milne_formula(t0, t1) - value_left = self.milne_formula(t0, tm) - value_righ = self.milne_formula(tm, t1) - value_some = value_left + value_righ - if depth < self.max_depth and abs(value_some - value_midl) > tolerance: - value_left = self.integrate_interval( - t0, tm, tolerance / 2, depth + 1 - ) - value_righ = self.integrate_interval( - tm, t1, tolerance / 2, depth + 1 - ) - return value_left + value_righ - - def integrate(self, *, tolerance: float = 1e-9) -> float: + def integrate(self, expx: int, expy: int) -> float: """ Computes the integral over the entire curve @@ -533,8 +520,8 @@ def integrate(self, *, tolerance: float = 1e-9) -> float: :rtype: float """ value = 0 - knots = self.curve.knots - for t0, t1 in zip(knots, knots[1:]): - sub_tolerance = tolerance * (t1 - t0) / (knots[-1] - knots[0]) - value += self.integrate_interval(t0, t1, sub_tolerance, depth=0) + for t0, t1 in zip(self.mesh, self.mesh[1:]): + tm = 0.5 * (t0 + t1) + value += self.milne_formula(t0, tm, expx, expy) + value += self.milne_formula(tm, t1, expx, expy) return value diff --git a/tests/test_geomprop.py b/tests/test_geomprop.py index 570b912..efc48e5 100644 --- a/tests/test_geomprop.py +++ b/tests/test_geomprop.py @@ -115,9 +115,11 @@ def test_shifted_rectangle(self): assert area == width * height assert Qx == area * center[1] assert Qy == area * center[0] - assert Ixx == width * height**3 / 12 + area * center[1] * center[1] + good_Ixx = width * height**3 / 12 + area * center[1] * center[1] + good_Iyy = height * width**3 / 12 + area * center[0] * center[0] + assert abs(Ixx - good_Ixx) < 1e-6 assert Ixy == area * center[0] * center[1] - assert Iyy == height * width**3 / 12 + area * center[0] * center[0] + assert abs(Iyy - good_Iyy) < 1e-6 @pytest.mark.order(3) @pytest.mark.dependency( @@ -157,12 +159,11 @@ def test_read_square(self): assert area == 4 assert Qx == 0 assert Qy == 0 - assert Ixx == 4 / 3 + assert abs(Ixx - 4 / 3) < 1e-9 assert Ixy == 0 - assert Iyy == 4 / 3 + assert abs(Iyy - 4 / 3) < 1e-9 @pytest.mark.order(3) - @pytest.mark.skip() @pytest.mark.timeout(20) @pytest.mark.dependency(depends=["TestToFromJson::test_begin"]) def test_read_circle(self): @@ -176,7 +177,7 @@ def test_read_circle(self): area = square.area() Qx, Qy = square.first_moment() Ixx, Ixy, Iyy = square.second_moment() - assert abs(area - math.pi) < 1e-6 + assert abs(area - math.pi) < 2e-6 assert abs(Qx) < 1e-6 assert abs(Qy) < 1e-6 assert abs(Ixx - math.pi / 4) < 1e-6