From 931f4dbcafa51650557213d0a8950e04059945cb Mon Sep 17 00:00:00 2001 From: markriegler Date: Fri, 16 Aug 2024 17:14:24 +0200 Subject: [PATCH] Getting it to the state I wanted it to be without all the rebase issues --- ...lerkin_laplace_problem_field_integrator.py | 22 +- splinepy/helpme/integrate.py | 200 ++++++++++++++---- tests/helpme/test_integrate.py | 114 +++++----- 3 files changed, 229 insertions(+), 107 deletions(-) diff --git a/examples/iga/galerkin_laplace_problem_field_integrator.py b/examples/iga/galerkin_laplace_problem_field_integrator.py index acad6cb68..e6dff58b0 100644 --- a/examples/iga/galerkin_laplace_problem_field_integrator.py +++ b/examples/iga/galerkin_laplace_problem_field_integrator.py @@ -5,15 +5,15 @@ -\Delta u = f with source term f=1. -First, homogeneous Dirichlet boundary conditions are applied. Then, inhomogeneous ones -are applied. +First, homogeneous Dirichlet boundary conditions are applied. Then, homogeneous +and inhomogeneous Dirichlet boundary conditions are applied on 3 of 4 boundaries. """ import numpy as np import splinepy as sp -# Test Case +# Number of refinements for the solution field n_refine = 15 @@ -130,7 +130,7 @@ def show_solution(geometry, solution_field): fi.assemble_vector(poisson_rhs) # Homogeneous Dirichlet boundary conditions - fi.assign_homogeneous_dirichlet_boundary_conditions() + fi.apply_homogeneous_dirichlet_boundary_conditions() fi.solve_linear_system() # Plot geometry and field show_solution(geometry, solution_field) @@ -138,10 +138,18 @@ def show_solution(geometry, solution_field): # Inhomogeneous Dirichlet boundary conditions def dirichlet_function(points): """ - On the boundary apply: g(x,y) = x/10 + On the boundary apply: g(x,y) = x/4 """ - return points[:, 0] / 10 + return points[:, 0] / 4 - fi.apply_dirichlet_boundary_conditions(dirichlet_function) + # Assemble again to override previous boundary conditions + fi.assemble_matrix(poisson_lhs) + fi.assemble_vector(poisson_rhs) + + # Apply boundary conditions on 3 boundaries + fi.apply_homogeneous_dirichlet_boundary_conditions(west=True) + fi.apply_dirichlet_boundary_conditions( + dirichlet_function, south=True, north=True + ) fi.solve_linear_system() show_solution(geometry, solution_field) diff --git a/splinepy/helpme/integrate.py b/splinepy/helpme/integrate.py index f029b4637..b884362ca 100644 --- a/splinepy/helpme/integrate.py +++ b/splinepy/helpme/integrate.py @@ -2,7 +2,6 @@ import numpy as _np -from splinepy._base import SplinepyBase as _SplinepyBase from splinepy.utils.data import cartesian_product as _cartesian_product from splinepy.utils.data import has_scipy as _has_scipy @@ -21,7 +20,7 @@ def _get_integral_measure(spline): .. math:: \\mathcal{J}_S = det(\\mathbf(J)) - If the physical dimension is bigger then the parametric dimension it will + If the physical dimension is bigger than the parametric dimension it will return .. math:: @@ -30,12 +29,13 @@ def _get_integral_measure(spline): Parameters ---------- spline : Spline / Multipatch - For parametric and physical dimension + The spline object for which the measure of one of its patches is determined Returns ------- measure : Callable - single patch only + A function which computes the jacobian's determinant for a single patch. + It takes a patch and query positions in the parametric domain as input """ # Check dimensionality if spline.dim == spline.para_dim: @@ -60,6 +60,7 @@ def measure(spline_patch, positions): def _default_quadrature_orders(spline): + # Expected degree for quadrature based on spline's degrees and parametric dimension expected_degree = spline.degrees * spline.para_dim + 1 if spline.is_rational: expected_degree += 2 @@ -71,7 +72,10 @@ def _default_quadrature_orders(spline): def _get_quadrature_information(spline, orders=None): """ - Select appropriate integration order (gauss-legendre) + Select appropriate integration order (gauss-legendre). + + Determines the integration points and weights for numerical integration + based on the given spline and the optionally given quadrature orders. Determinant of a polynomial spline with para_dim==dim has degree @@ -89,16 +93,17 @@ def _get_quadrature_information(spline, orders=None): Parameters ---------- spline : Spline - Spline for integration + The spline object for which the quadrature information is determined. orders : array-like (optional) - Orders along every parametric dimension + Orders along every parametric dimension. If not provided, default quadrature + orders will be used. Returns ------- positions : np.ndarray - quadrature position in unit-square + Quadrature points in unit-square weights : np.ndarray - quadrature weights + Quadrature weights """ # Determine integration points @@ -133,22 +138,24 @@ def _get_quadrature_information(spline, orders=None): def volume(spline, orders=None): r"""Compute volume of a given spline + Uses quadrature to compute volume. Can handle patches with multiple elements. + Parameters ---------- spline : Spline (self if called via integrator) - splinepy - spline type + The spline object for which the volume is computed. orders : array-like (optional) order for gauss quadrature Returns ------- volume : float - Integral of dim-dimensional object + The computed volume of the spline """ from splinepy.spline import Spline as _Spline - # Check i_nput type + # Check input type if not isinstance(spline, _Spline): raise NotImplementedError("integration only works for splines") @@ -179,9 +186,12 @@ def parametric_function( Parameters ---------- spline : Spline - (self if called via integrator) + The geometry over which the function is integrated function : Callable + The user-defined function to integrate. It takes points in the parametric + dimension as input and outputs a scalar or an array of scalars. orders : optional + Quadrature orders for numerical integration Returns ------- @@ -288,6 +298,19 @@ class Transformation: ) def __init__(self, spline, solution_field=None, orders=None): + """ + Setup transformation class with a given geometry. + + Parameters + ------------ + spline: spline + The geometry + solution_field: None or spline + Solution field as spline function. If not given, supports and quadrature + will be calculated for geometry + orders: None or list + Quadrature orders. If not given, default quadrature orders will be used + """ self._spline = spline self._solution_field = solution_field if solution_field is not None: @@ -389,17 +412,20 @@ def quadrature_weights(self): return self._quad_weights def get_element_grid_id(self, element_id): - """Compute element ID in grid + """Compute element ID in grid. + + The grid follows splinepy's ordering: first it goes into the x-direction, + after that into the y-direction. Parameters ---------- element_id: int - ID of element of spline + ID of spline's element in element grid Returns --------- element_grid_id: list - ID of element in grid + The grid ID of the element. """ if self._para_dim == 3: raise NotImplementedError( @@ -409,12 +435,12 @@ def get_element_grid_id(self, element_id): return self._grid_ids[element_id, :] def get_element_quad_points(self, element_id): - """For given element computes quad points + """Compute the quadrature points for a given element. Parameters ----------- element_id: int - ID of element in spline's element + ID of spline's element in element grid Returns ----------- @@ -450,7 +476,7 @@ def get_element_support(self, element_id): Parameters ------------ element_id: int - ID of element + ID of spline's element in element grid Returns --------- @@ -474,7 +500,7 @@ def jacobian(self, element_id): Parameters ------------ element_id: list - ID of element in grid + ID of spline's element in element grid Returns --------- @@ -494,7 +520,7 @@ def jacobian_inverse(self, element_id): Parameters ------------ element_id: list - ID of element in grid + ID of spline's element in element grid Returns --------- @@ -544,8 +570,8 @@ def compute_all_element_quad_points(self, recompute=False): Parameters ---------- - recompute: bool - Recompute quadrature points + recompute: bool (optional) + If True, recompute the qudrature points. Default is no recomputation """ if self._all_element_quad_points is not None and not recompute: return @@ -574,8 +600,8 @@ def compute_all_supports(self, recompute=False): Parameters -------------- - recompute: bool - Recompute quadrature points + recompute: bool (optional) + If True, recomputes the supports. The default is no recomputation """ if self._all_supports is not None and not recompute: return @@ -596,8 +622,8 @@ def compute_all_element_jacobians(self, recompute=False): Parameters ---------- - recompute: bool - Recompute Jacobians + recompute: bool (optional) + If True, recomputes the Jacobians. Default option is no recomputation """ if self._all_jacobians is not None and not recompute: return @@ -615,8 +641,8 @@ def compute_all_element_jacobian_inverses(self, recompute=False): Parameters ---------- - recompute: bool - Recompute Jacobians' inverses + recompute: bool (optional) + If True, recompute Jacobians' inverses. Default is no recomputation """ if self._all_jacobian_inverses is not None and not recompute: return @@ -673,7 +699,7 @@ def compute_all_element_measures(self, recompute=False): ) -class FieldIntegrator(_SplinepyBase): +class FieldIntegrator: __slots__ = ( "_helpee", "_solution_field", @@ -686,7 +712,21 @@ class FieldIntegrator(_SplinepyBase): ) def __init__(self, geometry, solution_field=None, orders=None): - """ """ + """ + Sets up solution field, its mapper, precomputes the transformation and + calculate the number of DoFs. + + Parameters + ---------- + geometry: spline + The geometry + solution_field: None or spline + Solution field. If not given, quadrature and supports will be calculated + using the geometry + orders: None or list + Quadrature order in each dimension. If not given, default quadrature + will be used + """ self._helpee = geometry if solution_field is None: self._solution_field = geometry.copy() @@ -711,7 +751,14 @@ def __init__(self, geometry, solution_field=None, orders=None): self.reset(orders) def reset(self, orders=None): - """ """ + """Sets up the transformation and resets the lhs and rhs. + + Parameters + ------------ + orders: None or list + If given, these orders will be used for quadrature. Otherwise, default + quadrature orders will be used. + """ self._trafo = Transformation( self._helpee, self._solution_field, orders ) @@ -729,7 +776,14 @@ def precompute_transformation(self): @property def supports(self): - """ """ + """ + Get the quadrature points' supports. + + Returns + ------- + supports: np.ndarray + The supports + """ if self._supports is None: self._supports = self._helpee.supports(self._trafo.all_quad_points) @@ -995,9 +1049,31 @@ def check_if_assembled(self): if self._global_system_matrix is None or self._global_rhs is None: raise ValueError("System is not yet fully assembled") - def get_boundary_dofs(self): + def get_boundary_dofs( + self, + all_boundaries=True, + north=False, + east=False, + south=False, + west=False, + ): """ - Get indices of boundary dofs. + Get indices of boundary dofs. Assumes that control points are arranged + in the following way: first one is the southwest corner, then the first + dimension goes from west to east and second dimension goes from south to + north. + + Parameters + ------------------ + all_boundaries: bool + If True, get indices of all boundary dofs and ignores values for + north, south, east and west. On the other, if any boundary is selected, + this value will be ignored + north: bool + If True, return indices for north boundary of geometry + east: bool + south: bool + west: bool Returns ------------- @@ -1012,26 +1088,60 @@ def get_boundary_dofs(self): multi_index = relevant_spline.multi_index + multi_indices = [ + multi_index[0, :], + multi_index[-1, :], + multi_index[:, 0], + multi_index[:, -1], + ] + + # If at least one boundary is True, set all_boundaries to False + if north or east or south or west: + all_boundaries = False + + boundaries = ( + [True] * 4 if all_boundaries else [west, east, south, north] + ) + indices = _np.unique( _np.hstack( - ( - multi_index[0, :], - multi_index[-1, :], - multi_index[:, 0], - multi_index[:, -1], - ) + [ + index + for index, boundary in zip(multi_indices, boundaries) + if boundary + ] ) ) return indices - def assign_homogeneous_dirichlet_boundary_conditions(self): + def apply_homogeneous_dirichlet_boundary_conditions( + self, + all_boundaries=True, + north=False, + east=False, + south=False, + west=False, + ): """ Assembles homogeneous Dirichlet boundary conditions + + Parameters + ------------- + all_boundaries: bool + If True, get indices of all boundary dofs and ignores values for + north, south, east and west + north: bool + If True, sets homogeneous Dirichlet condition on north boundary + east: bool + south: bool + west: bool """ self.check_if_assembled() - indices = self.get_boundary_dofs() + indices = self.get_boundary_dofs( + all_boundaries, north, east, south, west + ) self._global_system_matrix[indices, :] = 0 self._global_system_matrix[indices, indices] = 1 @@ -1087,7 +1197,9 @@ def apply_dirichlet_boundary_conditions( dof_vector = self.L2_projection(function) # Get relevant dofs - indices = self.get_boundary_dofs() + indices = self.get_boundary_dofs( + all_boundaries, north, east, south, west + ) if return_values: return dof_vector[indices], indices diff --git a/tests/helpme/test_integrate.py b/tests/helpme/test_integrate.py index 8232f65e1..3d96193d5 100644 --- a/tests/helpme/test_integrate.py +++ b/tests/helpme/test_integrate.py @@ -4,59 +4,6 @@ import splinepy -def test_transformation_class(np_rng): - """Test element transformation of single patch""" - # Create quadratic spline - spline = splinepy.BSpline( - degrees=[2, 2], - knot_vectors=[[0, 0, 0, 1, 1, 1], [0, 0, 0, 1, 1, 1]], - control_points=splinepy.utils.data.cartesian_product( - [np.linspace(0, 1, 3), np.linspace(0, 1, 3)] - ), - ) - - # Randomly insert knots along both parametric dimensions - spline.insert_knots(0, np_rng.random(2)) - spline.insert_knots(1, np_rng.random(2)) - - # Create transformation class for spline - trafo = splinepy.helpme.integrate.Transformation(spline) - - # Check whether element quadrature points lie inside element - ukvs = spline.unique_knots - trafo.compute_all_element_quad_points() - # Check quadrature points of each element - for element_id, quad_points in enumerate(trafo.all_quad_points): - grid_id = trafo.get_element_grid_id(element_id) - # Extract the corners of the current element - element_corners = [ - ukv[grid_dim_id : grid_dim_id + 2] - for ukv, grid_dim_id in zip(ukvs, grid_id) - ] - # Check if quadrature points lie within element - for dim, corners in enumerate(element_corners): - assert np.all( - (quad_points[:, dim] > corners[0]) - & (quad_points[:, dim] < corners[1]) - ), f"Quadrature points do not lie within element for dimension {dim}" - - # For given spline, all Jacobians should be identity matrix - eye = np.eye(spline.para_dim) - trafo.compute_all_element_jacobian_inverses() - for element_jacobians in trafo.all_jacobians: - for jacobian_at_quad_point in element_jacobians: - assert np.allclose( - eye, jacobian_at_quad_point - ), "All Jacobians should be identity matrix" - - # For created spline, all determinants should equal one - trafo.compute_all_element_jacobian_determinants() - assert np.allclose( - trafo.all_jacobian_determinants, - np.ones_like(trafo.all_jacobian_determinants), - ), "All Jacobians' determinants should be equal to one" - - def test_volume_integration_1D(np_rng): """ Test volume integration for splines using numerical integration of the @@ -87,20 +34,20 @@ def test_volume_integration_embedded(np_rng): Test volume integration for splines using numerical integration of the Jacobi-Determinant """ - # Test 1D -> 2D + # Test 1D -> 2D volume integration for Bezier expected_result = 2.0**1.5 bezier = splinepy.Bezier(degrees=[1], control_points=[[0, 0], [2, 2]]) assert np.allclose(bezier.integrate.volume(), expected_result) - # test for other types same spline + # For same curve, test other spline types assert np.allclose(bezier.bspline.integrate.volume(), expected_result) assert np.allclose( bezier.rationalbezier.integrate.volume(), expected_result ) assert np.allclose(bezier.nurbs.integrate.volume(), expected_result) - # Check if equal after refinement + # Check if volume is equal after degree elevation bezier.elevate_degrees([0, 0, 0]) assert np.allclose(bezier.integrate.volume(), expected_result) @@ -230,6 +177,7 @@ def test_assertions(np_rng): def test_function_integration(np_rng): col1_factor = 2 + # Define vector-valued constant function def volume_function(x): vf = np.ones((x.shape[0], 2)) # scale it with a factor to get a different value @@ -239,7 +187,61 @@ def volume_function(x): bezier = splinepy.Bezier( degrees=[1, 2], control_points=np_rng.random((6, 2)) ) + # Test function integration for constant functions assert np.allclose( [bezier.integrate.volume(), col1_factor * bezier.integrate.volume()], bezier.integrate.parametric_function(volume_function), ) + + +def test_transformation_class(np_rng): + """Test element transformation of single patch""" + # Create quadratic spline + spline = splinepy.BSpline( + degrees=[2, 2], + knot_vectors=[[0, 0, 0, 1, 1, 1], [0, 0, 0, 1, 1, 1]], + control_points=splinepy.utils.data.cartesian_product( + [np.linspace(0, 1, 3), np.linspace(0, 1, 3)] + ), + ) + + # Randomly insert knots along both parametric dimensions + spline.insert_knots(0, np_rng.random(2)) + spline.insert_knots(1, np_rng.random(2)) + + # Create transformation class for spline + trafo = splinepy.helpme.integrate.Transformation(spline) + + # Check whether element quadrature points lie inside element + ukvs = spline.unique_knots + trafo.compute_all_element_quad_points() + # Check quadrature points of each element + for element_id, quad_points in enumerate(trafo.all_quad_points): + grid_id = trafo.get_element_grid_id(element_id) + # Extract the corners of the current element + element_corners = [ + ukv[grid_dim_id : grid_dim_id + 2] + for ukv, grid_dim_id in zip(ukvs, grid_id) + ] + # Check if quadrature points lie within element + for dim, corners in enumerate(element_corners): + assert np.all( + (quad_points[:, dim] > corners[0]) + & (quad_points[:, dim] < corners[1]) + ), f"Quadrature points do not lie within element for dimension {dim}" + + # For given spline, all Jacobians should be identity matrix + eye = np.eye(spline.para_dim) + trafo.compute_all_element_jacobian_inverses() + for element_jacobians in trafo.all_jacobians: + for jacobian_at_quad_point in element_jacobians: + assert np.allclose( + eye, jacobian_at_quad_point + ), "All Jacobians should be identity matrix" + + # For created spline, all determinants should equal one + trafo.compute_all_element_jacobian_determinants() + assert np.allclose( + trafo.all_jacobian_determinants, + np.ones_like(trafo.all_jacobian_determinants), + ), "All Jacobians' determinants should be equal to one"