From 20c0e84cefe7aaca52d1b2043c715248825fe127 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 13:49:16 +1030 Subject: [PATCH 1/2] refactor: StructuredGrid2D extends StructuredGrid2DGeometry Same consolidation as the 3D grid. The 2D support re-implemented the geometry that StructuredGrid2DGeometry already provided, and the geometry copy had the same stub global_node_indices (elements came back as ij pairs). - StructuredGrid2D inherits the geometry (nodes, elements, point location, index conversions) and keeps the interpolation machinery - move the pure-geometry members into the geometry class: bilinear weights, position_to_local_coordinates / dof_coefs, cell corners / vertices, node_indexes_to_position, global_cell_indices and evaluate_value - implement global_node_indices in the geometry class - fix print_geometry, which indexed element [2] of a 2D vector - BoundingBox.structured_grid imports loop_common.supports lazily (supports now import geometry) A reference test, generated from master before the change, pins the support's behaviour (162 arrays across three grids, including exact boundary points). Co-Authored-By: Claude Sonnet 5.5 --- .../src/loop_common/geometry/_bounding_box.py | 4 +- .../geometry/_structured_grid_2d.py | 164 ++++++- .../supports/_2d_structured_grid.py | 407 +----------------- .../data/structured_grid_2d_reference.npz | Bin 0 -> 75292 bytes .../tests/test_structured_grid_2d_geometry.py | 74 ++++ .../test_structured_grid_2d_reference.py | 182 ++++++++ 6 files changed, 427 insertions(+), 404 deletions(-) create mode 100644 packages/loop_common/tests/data/structured_grid_2d_reference.npz create mode 100644 packages/loop_common/tests/test_structured_grid_2d_geometry.py create mode 100644 packages/loop_common/tests/test_structured_grid_2d_reference.py diff --git a/packages/loop_common/src/loop_common/geometry/_bounding_box.py b/packages/loop_common/src/loop_common/geometry/_bounding_box.py index d49781295..ab177be6e 100644 --- a/packages/loop_common/src/loop_common/geometry/_bounding_box.py +++ b/packages/loop_common/src/loop_common/geometry/_bounding_box.py @@ -8,7 +8,6 @@ # from LoopStructural.utils.exceptions import LoopValueError from loop_common.math import rng -from loop_common.supports import StructuredGrid logger = getLogger(__name__) @@ -667,6 +666,9 @@ def structured_grid( name: str = "bounding_box", local_coordinates: bool = False, ): + # imported here because loop_common.supports builds on loop_common.geometry + from loop_common.supports import StructuredGrid + # python is passing a reference to the cell_data, vertex_data dicts so we need to # copy them to make sure that different instances of StructuredGrid are not sharing the same # underlying objects diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py b/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py index e4c8c8b74..7c48836b6 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py @@ -56,12 +56,10 @@ def elements(self) -> np.ndarray: return self.global_node_indices(self.cell_corner_indexes(cell_indexes)) def print_geometry(self): - logger.info(f"Origin: {self.origin[0]:f} {self.origin[1]:f} {self.origin[2]:f}") - logger.info( - f"Cell size: {self.step_vector[0]:f} {self.step_vector[1]:f} {self.step_vector[2]:f}" - ) + logger.info(f"Origin: {self.origin[0]:f} {self.origin[1]:f}") + logger.info(f"Cell size: {self.step_vector[0]:f} {self.step_vector[1]:f}") max = self.origin + self.nsteps_cells * self.step_vector - logger.info(f"Max extent: {max[0]:f} {max[1]:f} {max[2]:f}") + logger.info(f"Max extent: {max[0]:f} {max[1]:f}") def cell_centres(self, global_index: np.ndarray) -> np.ndarray: cell_indexes = self.global_index_to_cell_index(global_index) @@ -142,5 +140,157 @@ def global_index_to_node_index(self, global_index): cell_indexes[:, 1] = global_index // self.nsteps[0, None] % self.nsteps[1, None] return cell_indexes - def global_node_indices(self, node_indexes): - return node_indexes + def global_node_indices(self, node_indexes: np.ndarray) -> np.ndarray: + """Flat node indexes from (..., 2) node ij indexes.""" + return self._global_indices(node_indexes, self.nsteps) + + def bilinear(self, local_coords: np.ndarray) -> np.ndarray: + """ + returns the bilinear interpolation for the local coordinates + Parameters + ---------- + x - double, array of doubles + y - double, array of doubles + z - double, array of doubles + + Returns + ------- + array of interpolation coefficients + + """ + + return np.array( + [ + (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]), + local_coords[:, 0] * (1 - local_coords[:, 1]), + (1 - local_coords[:, 0]) * local_coords[:, 1], + local_coords[:, 0] * local_coords[:, 1], + ] + ).T + + def position_to_local_coordinates(self, pos: np.ndarray) -> np.ndarray: + """ + Convert from global to local coordinates within a cel + Parameters + ---------- + pos - array of positions inside + + Returns + ------- + localx, localy, localz + + """ + # TODO check if inside mesh + + # calculate local coordinates for positions + local_coords = np.zeros(pos.shape) + local_coords[:, 0] = ( + (pos[:, 0] - self.origin[None, 0]) % self.step_vector[None, 0] + ) / self.step_vector[None, 0] + local_coords[:, 1] = ( + (pos[:, 1] - self.origin[None, 1]) % self.step_vector[None, 1] + ) / self.step_vector[None, 1] + + return local_coords + + def position_to_dof_coefs(self, pos: np.ndarray): + """ + global posotion to interpolation coefficients + Parameters + ---------- + pos + + Returns + ------- + + """ + local_coords = self.position_to_local_coordinates(pos) + weights = self.bilinear(local_coords) + return weights + + def _global_indices(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarray: + if len(indexes.shape) == 1: + raise ValueError("Indexes must be a 2D array") + if indexes.shape[-1] != 2: + raise ValueError("Last dimension of cell indexing needs to be ijk indexing") + original_shape = indexes.shape + indexes = indexes.reshape(-1, 2) + gi = indexes[:, 0] + nsteps[0] * indexes[:, 1] + return gi.reshape(original_shape[:-1]) + + def global_cell_indices(self, indexes: np.ndarray) -> np.ndarray: + return self._global_indices(indexes, self.nsteps_cells) + + def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: + + original_shape = node_indexes.shape + node_indexes = node_indexes.reshape((-1, 2)) + xy = np.zeros((node_indexes.shape[0], 2), dtype=float) + xy[:, 0] = self.origin[0] + self.step_vector[0] * node_indexes[:, 0] + xy[:, 1] = self.origin[1] + self.step_vector[1] * node_indexes[:, 1] + xy = xy.reshape(original_shape) + return xy + + def position_to_cell_corners(self, pos): + """Get the global indices of the vertices (corner) nodes of the cell containing each point. + + Parameters + ---------- + pos : np.array + (N, 2) array of xy coordinates representing the positions of N points. + + Returns + ------- + globalidx : np.array + (N, 4) array of global indices corresponding to the 4 corner nodes of the cell + each point lies in. If a point lies outside the support, its corresponding entry + will be set to -1. + inside : np.array + (N,) boolean array indicating whether each point is inside the support domain. + """ + corner_index, inside = self.position_to_cell_index(pos) + corners = self.cell_corner_indexes(corner_index) + globalidx = self.global_node_indices(corners) + # if global index is not inside the support set to -1 + globalidx[~inside] = -1 + return globalidx, inside + + def position_to_cell_vertices(self, pos): + """Get the vertices of the cell a point is in + + Parameters + ---------- + pos : np.array + Nx3 array of xyz locations + + Returns + ------- + np.array((N,3),dtype=float), np.array(N,dtype=int) + vertices, inside + """ + gi, inside = self.position_to_cell_corners(pos) + + node_indexes = self.global_index_to_node_index(gi.flatten()) + return self.node_indexes_to_position(node_indexes), inside + + def evaluate_value(self, evaluation_points: np.ndarray, property_array: np.ndarray): + """ + Evaluate the value of of the property at the locations. + Trilinear interpolation dot corner values + + Parameters + ---------- + evaluation_points np array of locations + property_name string of property name + + Returns + ------- + + """ + idc, inside = self.position_to_cell_corners(evaluation_points) + v = np.zeros(idc.shape) + v[:, :] = np.nan + + v[inside, :] = self.position_to_dof_coefs(evaluation_points[inside, :]) + v[inside, :] *= property_array[idc[inside, :]] + return np.sum(v, axis=1) diff --git a/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py b/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py index 974821b6d..85e977ef9 100644 --- a/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py +++ b/packages/loop_common/src/loop_common/supports/_2d_structured_grid.py @@ -7,6 +7,8 @@ import numpy as np +from loop_common.geometry import StructuredGrid2DGeometry + from ..math.finite_difference_stencil import Operator from . import SupportType from ._base_support import BaseSupport @@ -14,8 +16,14 @@ logger = logging.getLogger(__name__) -class StructuredGrid2D(BaseSupport): - """ """ +class StructuredGrid2D(StructuredGrid2DGeometry, BaseSupport): + """A 2D structured grid used as an interpolation support. + + The grid geometry (node and cell positions, point location, index conversions, bilinear + weights and node-value interpolation) comes from + :class:`loop_common.geometry.StructuredGrid2DGeometry`; this class adds the interpolation + machinery. ``nsteps`` is the number of **nodes** in each direction. + """ dimension = 2 @@ -33,57 +41,15 @@ def __init__( nsteps - 2d list or numpy array of ints step_vector - 2d list or numpy array of int """ - if origin is None: - origin = np.zeros(2) - if nsteps is None: - nsteps = np.array([10, 10]) - if step_vector is None: - step_vector = np.ones(2) - + StructuredGrid2DGeometry.__init__(self, origin=origin, nsteps=nsteps, step_vector=step_vector) self.type = SupportType.StructuredGrid2D - self.nsteps = np.ceil(np.array(nsteps)).astype(int) - self.step_vector = np.array(step_vector) - self.origin = np.array(origin) - self.maximum = origin + self.nsteps * self.step_vector - - self.dim = 2 - self.nsteps_cells = self.nsteps - 1 - self.n_cell_x = self.nsteps[0] - 1 - self.n_cell_y = self.nsteps[1] - 1 self.properties = {} - - # calculate the node positions using numpy (this should probably not - # be stored as it defeats - # the purpose of a structured grid - - # self.barycentre = self.cell_centres(np.arange(self.n_elements)) - self.regions = {} self.regions["everywhere"] = np.ones(self.n_nodes).astype(bool) - @property - def nodes(self): - max = self.origin + self.nsteps_cells * self.step_vector - x = np.linspace(self.origin[0], max[0], self.nsteps[0]) - y = np.linspace(self.origin[1], max[1], self.nsteps[1]) - xx, yy = np.meshgrid(x, y, indexing="ij") - return np.array([xx.flatten(order="F"), yy.flatten(order="F")]).T - - @property - def n_nodes(self): - return self.nsteps[0] * self.nsteps[1] - def set_nelements(self, nelements) -> int: raise NotImplementedError("Cannot set number of elements for 2D structured grid") - @property - def n_elements(self): - return self.nsteps_cells[0] * self.nsteps_cells[1] - - @property - def element_size(self): - return np.prod(self.step_vector) - @property def barycentre(self): return self.cell_centres(np.arange(self.n_elements)) @@ -91,339 +57,6 @@ def barycentre(self): # @property # def barycentre(self): # return self.cell_centres(np.arange(self.n_elements)) - @property - def elements(self) -> np.ndarray: - global_index = np.arange(self.n_elements) - cell_indexes = self.global_index_to_cell_index(global_index) - - return self.global_node_indices(self.cell_corner_indexes(cell_indexes)) - - def print_geometry(self): - print(f"Origin: {self.origin[0]:f} {self.origin[1]:f} {self.origin[2]:f}") - print( - f"Cell size: {self.step_vector[0]:f} {self.step_vector[1]:f} {self.step_vector[2]:f}" - ) - max = self.origin + self.nsteps_cells * self.step_vector - print(f"Max extent: {max[0]:f} {max[1]:f} {max[2]:f}") - - def cell_centres(self, global_index: np.ndarray) -> np.ndarray: - """[summary] - - [extended_summary] - - Parameters - ---------- - global_index : [type] - [description] - - Returns - ------- - [type] - [description] - """ - cell_indexes = self.global_index_to_cell_index(global_index) - cell_centres = np.zeros((cell_indexes.shape[0], 2)) - - cell_centres[:, 0] = ( - self.origin[None, 0] - + self.step_vector[None, 0] * 0.5 - + self.step_vector[None, 0] * cell_indexes[:, 0] - ) - cell_centres[:, 1] = ( - self.origin[None, 1] - + self.step_vector[None, 1] * 0.5 - + self.step_vector[None, 1] * cell_indexes[:, 1] - ) - return cell_centres - - def position_to_cell_index(self, pos: np.ndarray) -> tuple[np.ndarray, np.ndarray]: - """[summary] - - [extended_summary] - - Parameters - ---------- - pos : [type] - [description] - - Returns - ------- - [type] - [description] - """ - inside = self.inside(pos) - cell_indexes = np.zeros((pos.shape[0], 2)) - cell_indexes[:, 0] = pos[:, 0] - self.origin[None, 0] - cell_indexes[:, 1] = pos[:, 1] - self.origin[None, 1] - cell_indexes /= self.step_vector[None, :] - return cell_indexes.astype(int), inside - - def inside(self, pos: np.ndarray) -> np.ndarray: - # check whether point is inside box - inside = np.ones(pos.shape[0]).astype(bool) - for i in range(self.dim): - inside *= pos[:, i] > self.origin[None, i] - inside *= ( - pos[:, i] - < self.origin[None, i] + self.step_vector[None, i] * self.nsteps_cells[None, i] - ) - return inside - - def check_position(self, pos: np.ndarray) -> np.ndarray: - """[summary] - - [extended_summary] - - Parameters - ---------- - pos : [type] - [description] - - Returns - ------- - [type] - [description] - """ - - if len(pos.shape) == 1: - pos = np.array([pos]) - if len(pos.shape) != 2: - raise ValueError("Position array needs to be a list of points or a point") - - return pos - - def bilinear(self, local_coords: np.ndarray) -> np.ndarray: - """ - returns the bilinear interpolation for the local coordinates - Parameters - ---------- - x - double, array of doubles - y - double, array of doubles - z - double, array of doubles - - Returns - ------- - array of interpolation coefficients - - """ - - return np.array( - [ - (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]), - local_coords[:, 0] * (1 - local_coords[:, 1]), - (1 - local_coords[:, 0]) * local_coords[:, 1], - local_coords[:, 0] * local_coords[:, 1], - ] - ).T - - def position_to_local_coordinates(self, pos: np.ndarray) -> np.ndarray: - """ - Convert from global to local coordinates within a cel - Parameters - ---------- - pos - array of positions inside - - Returns - ------- - localx, localy, localz - - """ - # TODO check if inside mesh - - # calculate local coordinates for positions - local_coords = np.zeros(pos.shape) - local_coords[:, 0] = ( - (pos[:, 0] - self.origin[None, 0]) % self.step_vector[None, 0] - ) / self.step_vector[None, 0] - local_coords[:, 1] = ( - (pos[:, 1] - self.origin[None, 1]) % self.step_vector[None, 1] - ) / self.step_vector[None, 1] - - return local_coords - - def position_to_dof_coefs(self, pos: np.ndarray): - """ - global posotion to interpolation coefficients - Parameters - ---------- - pos - - Returns - ------- - - """ - local_coords = self.position_to_local_coordinates(pos) - weights = self.bilinear(local_coords) - return weights - - def neighbour_global_indexes(self, mask=None, **kwargs): - """ - Get neighbour indexes - - Parameters - ---------- - kwargs - indexes array specifying the cells to return neighbours - - Returns - ------- - - """ - indexes = None - if "indexes" in kwargs: - indexes = kwargs["indexes"] - if "indexes" not in kwargs: - gi = np.arange(self.n_nodes) - indexes = self.global_index_to_node_index(gi) - edge_mask = ( - (indexes[:, 0] > 0) - & (indexes[:, 0] < self.nsteps[0] - 1) - & (indexes[:, 1] > 0) - & (indexes[:, 1] < self.nsteps[1] - 1) - ) - indexes = indexes[edge_mask, :].T - # ii = [] - # jj = [] - # for i in range(1, self.nsteps[0] - 1): - # for j in range(1, self.nsteps[1] - 1): - # ii.append(i) - # jj.append(j) - # indexes = np.array([ii, jj]) - # indexes = np.array(indexes).T - if indexes.ndim != 2: - print(indexes.ndim) - return - # determine which neighbours to return default is diagonals included. - if mask is None: - mask = np.array([[-1, 0, 1, -1, 0, 1, -1, 0, 1], [1, 1, 1, 0, 0, 0, -1, -1, -1]]) - neighbours = indexes[:, None, :] + mask[:, :, None] - return (neighbours[0, :, :] + self.nsteps[0, None, None] * neighbours[1, :, :]).astype( - np.int64 - ) - - def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: - """ - Returns the indexes of the corners of a cell given its location xi, - yi, zi - - Parameters - ---------- - x_cell_index - y_cell_index - z_cell_index - - Returns - ------- - - """ - corner_indexes = np.zeros((cell_indexes.shape[0], 4, 2), dtype=np.int64) - xcorner = np.array([0, 1, 0, 1]) - ycorner = np.array([0, 0, 1, 1]) - corner_indexes[:, :, 0] = ( - cell_indexes[:, None, 0] + corner_indexes[:, :, 0] + xcorner[None, :] - ) - corner_indexes[:, :, 1] = ( - cell_indexes[:, None, 1] + corner_indexes[:, :, 1] + ycorner[None, :] - ) - return corner_indexes - - def global_index_to_cell_index(self, global_index): - """ - Convert from global indexes to xi,yi,zi - - Parameters - ---------- - global_index - - Returns - ------- - - """ - # determine the ijk indices for the global index. - # remainder when dividing by nx = i - # remained when dividing modulus of nx by ny is j - cell_indexes = np.zeros((global_index.shape[0], 2), dtype=np.int64) - cell_indexes[:, 0] = global_index % self.nsteps_cells[0, None] - cell_indexes[:, 1] = global_index // self.nsteps_cells[0, None] % self.nsteps_cells[1, None] - return cell_indexes - - def global_index_to_node_index(self, global_index): - cell_indexes = np.zeros((global_index.shape[0], 2), dtype=np.int64) - cell_indexes[:, 0] = global_index % self.nsteps[0, None] - cell_indexes[:, 1] = global_index // self.nsteps[0, None] % self.nsteps[1, None] - return cell_indexes - - def _global_indices(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarray: - if len(indexes.shape) == 1: - raise ValueError("Indexes must be a 2D array") - if indexes.shape[-1] != 2: - raise ValueError("Last dimension of cell indexing needs to be ijk indexing") - original_shape = indexes.shape - indexes = indexes.reshape(-1, 2) - gi = indexes[:, 0] + nsteps[0] * indexes[:, 1] - return gi.reshape(original_shape[:-1]) - - def global_cell_indices(self, indexes: np.ndarray) -> np.ndarray: - return self._global_indices(indexes, self.nsteps_cells) - - def global_node_indices(self, indexes: np.ndarray) -> np.ndarray: - return self._global_indices(indexes, self.nsteps) - - def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: - - original_shape = node_indexes.shape - node_indexes = node_indexes.reshape((-1, 2)) - xy = np.zeros((node_indexes.shape[0], 2), dtype=float) - xy[:, 0] = self.origin[0] + self.step_vector[0] * node_indexes[:, 0] - xy[:, 1] = self.origin[1] + self.step_vector[1] * node_indexes[:, 1] - xy = xy.reshape(original_shape) - return xy - - def position_to_cell_corners(self, pos): - """Get the global indices of the vertices (corner) nodes of the cell containing each point. - - Parameters - ---------- - pos : np.array - (N, 2) array of xy coordinates representing the positions of N points. - - Returns - ------- - globalidx : np.array - (N, 4) array of global indices corresponding to the 4 corner nodes of the cell - each point lies in. If a point lies outside the support, its corresponding entry - will be set to -1. - inside : np.array - (N,) boolean array indicating whether each point is inside the support domain. - """ - corner_index, inside = self.position_to_cell_index(pos) - corners = self.cell_corner_indexes(corner_index) - globalidx = self.global_node_indices(corners) - # if global index is not inside the support set to -1 - globalidx[~inside] = -1 - return globalidx, inside - - def evaluate_value(self, evaluation_points: np.ndarray, property_array: np.ndarray): - """ - Evaluate the value of of the property at the locations. - Trilinear interpolation dot corner values - - Parameters - ---------- - evaluation_points np array of locations - property_name string of property name - - Returns - ------- - - """ - idc, inside = self.position_to_cell_corners(evaluation_points) - v = np.zeros(idc.shape) - v[:, :] = np.nan - - v[inside, :] = self.position_to_dof_coefs(evaluation_points[inside, :]) - v[inside, :] *= property_array[idc[inside, :]] - return np.sum(v, axis=1) - def evaluate_gradient(self, evaluation_points, property_array): T = np.zeros((evaluation_points.shape[0], 2, 4)) _vertices, T, elements, inside = self.get_element_gradient_for_location(evaluation_points) @@ -478,24 +111,6 @@ def get_element_for_location( a = self.position_to_dof_coefs(pos) return vertices, a, elements, inside - def position_to_cell_vertices(self, pos): - """Get the vertices of the cell a point is in - - Parameters - ---------- - pos : np.array - Nx3 array of xyz locations - - Returns - ------- - np.array((N,3),dtype=float), np.array(N,dtype=int) - vertices, inside - """ - gi, inside = self.position_to_cell_corners(pos) - - node_indexes = self.global_index_to_node_index(gi.flatten()) - return self.node_indexes_to_position(node_indexes), inside - def onGeometryChange(self): pass diff --git a/packages/loop_common/tests/data/structured_grid_2d_reference.npz b/packages/loop_common/tests/data/structured_grid_2d_reference.npz new file mode 100644 index 0000000000000000000000000000000000000000..1a272d143d2cc2d12a7695d650e3ea8a9bdcdff5 GIT binary patch literal 75292 zcmeFa1z6Q<);-Q6JF-QC@dbO}gzN_RlF zSd*JtS?icvk(ru1fP%n-gx);(?MvYM3rPGH@&-D5XNP#3ctKK_VuJDq(~K0;aD&vs z^pvD1HEGF+802KNh%hJwJE0yz(MJKudasZ^JUBz+V~S&fs(^Y9LM96i={^V|{u?BT z1Vl2oNDzOEL`_S_*!bUx2k8&+u#qCW-{OJ4#bcsjZ(w3$a^Ez9e{LGxf50^06z~?m zvyEwNDZwq4J1?188t56A-iN~fFN8t^a?eNBWBEq6IG*3yX}<2czBM zcy?>wg~<(BDJBz7WT~g>n~0KY*xN893_P`smg!|6{`3+6A2A#zXfE z^QqBo;TFx)Uk&pup4IH{Bl(wk=q8-+d+5)AIQ%Dp@PCAVb!!orTZ?EL+@+LTxBSbD z@%s11_!$GjPZ)GetqjcWGRl7fNm8`l{;f4&Z;|{aeBmE_k6Yb858ppXv2cV$B9EZH zPd?RXe3?v8e7{gt8;lJ6`ECWAWMu zNd7I#I~Ly$ZucKRKm0k+zs@s0x=c{7e?hwRfDD{K-D@acJno4|!efzy7a|F$A_-t3 z36De)p82{^Kv*mm!@chFY+;LGAcqC2FE}%7E&2{UVA%5XhhDh(h|KHeqbdqmAWup; z!M6uAzkL+=?Yqc8Gwx;643l@ox5K>4VwxJ34q7)!%LG&H9E)WPDfyPMou8s}~)>+L#lHEh;#5br{E$gG_PjNxQQt3I}{*4ENA*> zC{)9>v;S%(<->>ceJ%Bk%Y5MA@VXDF-Dz;aaD!SP=bIQ|#6S6^yex&+ffNa%5FmJF z%7EbU;GWrL=K7Ux;bPxLp1GNof%Q!aP_s7sWzD5#V5+TSPeyT@Al)IqBPahcki^Im zsAdirpP~X4=YT2E)Id?3KGT3e&+z931b%|VWD8QK%|vvbSQi+l<>}V*_#xn)E#9;8H?rPW-7TGu`CU2S?4T^PZN!1_qHq=Q2tjs7d5M3w8FIuks zP1{{Gs7!PgT;D_1W1H!Ky95AUld3*own?^ti5;NE++I(%AR%87U?5C$f3KObnuHkU-X*2#KnV62dtc)M;cCw( zCe)FH>}$`OqLop9Ko4=7Qw0{H+em|WT`=~?6U31QX-8V*+4G7ZiS3?OCc?N9W5iZJ0ekP-y*-=Mu`3y`Ax(Bjdk{u#I$mGDmeIQ@jf@>1LgaupR zdD#cLAo)`Ck!ENd4nd#q{pm$S2e053E$*!y?<{y5idtrtraG2ZKLsQFcf1Tgg7Kf; zmwn47KV}^52Zh%UL72fFnsXRn67)i6JYLUwPQksc!D$77AGx*4cMuH{f%Iy~`E@o^ z4D(=1pybIQT1q>9&m5O{9HL+urvVwk6KHy_m4kHeJje(nJW2M9W4Pnf)Df3r`$g_R zdJO^*gVXd?dTSuT5xb+SXSI-V=x4Z$WiXPN<^{wV7ROT!3b4bApd?j@FD?-e4)>o& z6}GM~u`5o&(fYi7DwrNU9uX%-a4oLr3S)gnJ+Mv8v6}G{ei&Pt+ulRZ9$i+PDJ?i`%m4EGw z|C4bkjG<fg@5S}HoD&NX0T*u zFjn?@d}wfSK!usbR>3J&=JetnF5E^s9mpYe$cU>cW>)kRTp>4e<(AGC+C8PUU7xAU^26JYCtIe$vl#M_qSBndGmFWg^tc`dlHHu=inBbqnVP<(S(`OJm1M zGwkc<(yl_w5W)C%*M_dTsyjddCHC^Ytb%04YlV6>NUddZi%ZFUBS6SUSNT>rPG{Ey z)ljRe!G^iV8D~=)m(M|L4Hk*%uMMxs8A67QF5O1Xlgr8uue0mthJZ#D*Cn%9palKp zVwbMq5p9ihS3OsauuEi^b8=Ua$=ArgokJdfT-EpoTzyB^fgh>v&7Y;bzh@f%?W^){ZQo~pDjg4aABE_8eW#-~E?M1DtMNWoiP zZ!)bZM7bQ}WY=CK^exI53tkNlm=c@5)wRsyA&i7hjJD?wf@PVJ{G}m-h=k>*@!FJq zn0RdrFsRSJD01;*`b~e0#KyBRzuK68*Oo`s5r`&`N{Vnge!-GVmBJG^!uk$C>+iwa zi5*#lfwVw%v@q&CaGl8(@Znt($w7#f-6C&{5=PDkhi(d1jLsR$GeXE5Vt$bfXA&6l zxv`#z7BUTap_%Z&hJA7u{c!G=$QAyU`E8u`13gcclOOd~q)qoHf*}wp7|j}kybzd! zUdXQ=7@e$347`Mlo&wnrQP#s#pRx0Ln6icE>B3UaS(}@5%uY@a+Zg>@OUszz z?P^G`WQ8^tf|DUvrlDhwPP^In*+{G3PU~_Ib67y^y1s6CnTCE zH?`y>)#+T>N}|}TmNNpP-wrBX17lf0gZ*Te?RR9bD?beqx78z8BGA?p`mClmFkKQh z`g^Vq`quKsY=rd;q6=eTd+AM!U$N#T&H+S&~3&E;pHS#hiN$vtZi$XNa7>f zmr23wf-DNqZjC1uhE0Ovqq*U-?+%+sT|zi3TVA2b=X)>V3Wj#qIT{T!bBZSgYHa9o z8$H=HCwvFtpvdDkrbqrVLCsY2c%IQia}H_M)U@dgf4g?OLXjie5)-ns%eOz(r7nUW zp{g8J#u%>>he)T1+HUyO1?>q8o+)Z;>o)9*h($}UP(wsvQnDb_O=xQ_)weEmpXn1% z;R(5+kBGaxQ&8ev97AcnZR#Yb*A)@ZD>E(y&5*cv=&DX8i5{8G6+6o~a6TS7{aWyr zh48tNC#jh!;ylL=ZC7?6e2?;mMZo;4LO7(P(|3g479Ql|J)yo!^FHNDxX;!|SmwCj zkEEDYLfbT^V_foub(xn@!UM+uTS{Y(<|$p;E{wrClHIUl71}_RShEeYToF~#*Yt)0 zBjGI{UeJC-fS)@p0QR}V2|m~i%WPwi6qq4rKke8f|Lr`%1P5Dxn~5?1C=+X&>E7gF z9o@ULN_VL58F=yW_t(UTGxp~J4_f&pDRaEy)25ztLF3vCi-_W>2!8m0GY_jED2M`u zSo-QIj#h6R-e$~e@6MOClu=C-QJqn*W|UadQ~_%CtCa_rdkOuFgKc`E2`uDXyRW7t zOUfEKEu_xEKcra+DCs~5EQjUA*T!$l7|GEy1FtQp>EVECr#|W1v_LgknfN--awQXp z&gerop#@20e{D}zbf_YTMR`>6J%3w#R@HRSK5~_4(b*ULMck*H+`W7;J(~$xomMF| zBGmZr^S#fZ{3FTXsJ<$kf90Pzu+tCLJyeoz^hz9xuoQ8J2(nWdeJ_;%_AGHTj;9OV zPYR1;y7p0&z+uNcA6=aY*K~!iV=)b+QmrfvXWh9q&6yxO5@K)^CKK(07jRk05Lm(U z@~Q}0vr%3iuAMD{tAggM{p3iE z^J&h`O|W{}=~cWocK z^1k`a2kG;S@9hg~myTjX7IzU0)JPGg{ix`+b@vVl>hM(z>uvPq1qYMgpAqql)19{Kz|t@%1e@iVnAKNBy^Tz>0Z(ghE54$V z(FkSnz+uFTb67@k#(J`Gg|hq>+frf75gB^X^KshXTWjI7$fx_H=X>Mn>zT_2&MdRf zO9e;UkUkehBT+?K%@-Uy7t+upP>uA1N844X=ME){u|$WgYwi<2iH7o;Hm?ye*KeyE zz!Q3704CyK{+^9FfK@tDX_Zd2RB zFJ-Y>6ne7blCKZ>zR=h4NGj25LK0>0C}811`!3O^NQu)Us=9MW!T5aQfiMj&iu%?j z(Dkx8s3?m;9F;+YOvNWwZk?!zi`ZVKSI(-vsn4zelRzP2uVVcY*rCAG7S8f>-M3n_ zL-aXKi<=K6KPKmO19UW+qD5X_cTtz8Pdhv!d>l0694C~4Q#U7@!%TZ%E_WTe+S6w~ z5fX=Fh;`WDjl48ddssr5(4hcC6k}?h4 z{U)qda>Pmaai@Zfye@tW__(Js!5q+v{RViTs&!h*4a3v{L7#xn<4j>h!$gz5)E1!R zVNg$@ZZhniX#1|jA}s^~8G3q@3lr}i-sb5qG1F!q;=tkxuYtj9W8nUAR)=#e!5`#>vxj_?EhzMAO-KBzd{`53jkbN#aSQ^D#W?9J3eM zV&N%PwvD$x4q-zUYFo9|Jbk^&(m+^&q6DMRcvLz1X}+nYhU@4#%NwNSm=_;B(Px>o z`-_@*Gh{WU-V0~9CU#7cvOgT2*|RgMsZGfT9!UpZ@=6oNu=d4md-b=pa1O1|7DLsG zA=xastEeBDd9Hy}tw+2m${^=3eV>Hgkmw$lJGQdjy>GG6%_Zrfiv#O>fiKH-4nksW z)eW6{DjpzehsLlEV@X$m=P?!B0c-q97M5146OMXwKT+BKHL5Nq?_;mZmSdhGluxRW z6G4rb)Ra%JgsV7CKzlu(cE^~t6&@B#$!N61RCCXWW$1BhM@F#OR6Tj!o;b;v#k`m5 zU_?GU;x2$;MONF#=f;_z)8T{cm(oVHhz4EUAdHVkod1-Y>UtKg7^}HM5?$|^E=xns z0sY5V`UOwEaQBUNP|plEn_dghQn%O2*e4-h#gpF8M703Ej*&V}RD1+=Gu;gvvdq#_ zZ^EAL&qD>XoMrbrAT1Us%_KYKSsb!F>E&yc(hrZqQ_DktyMSvI{T}*lC?xJ!zo7Y` zwmz$R#R=)_Q5Uu%3=Qm6c&$EOB}j?I2w&O9T2&0uNDL@#P2Ed&=MlaH+Ct=!l=r1rKqgw*bv`K7`NtK(JR303sLOrN(G}qrec)T2ItZF zEY9*|Tx#-8YQaI=tC!p+jn_44M{ICNd)`3p@*F@ zukC@%Zrg-(7sNVqAcK#JaKgDoQL)l42HAs#!qhg_Hv*`aTh+8MzDb8X<5!>wTcztC za%H&4_xZjMaDvzX#rP_&U$v+ODIlHt_82_f;aL|G&wkCxRyImuoEiidDIVt;cZKBOhzOm8Av1*rBObrhMJ$Ted7ir$d zpU3sepDO4J%4npBxw`XbT(<{DYqbS2bSn!_sAczad*`4jlI=C(TN~_I0u)N=<^gXW zRse+?@nt;KsIXTGs2*aeWsuHz$*omHlWmOj@cK;oXhBC}6+4*c@r7-|kM-{_dG~wC zzt?SIv@dv}DK?iW-s_V(RF3;)G#63zTo#;kkhiK@bY~qaZB*YioL7!Pu;ngeHK3L%?{NDCPzYiM7q>oY5`um=;30{N=Jm0?vUK%2DKIa{qWL`zos&-prU9iORjhRJk9!Z!(QwQqtbC1rz0@_ z=W`ZYWsUHaGWGPQUiI4zokVUkfC;r3CwB7vEgh%+29m^#`snS(0?0#zFTjpjS_1>s zSfkMTt@m-u(?aN5X13mXY;mrmvr|Egq>Eo&fs3^O6BlCpF*e-69)s#}P~ezZREV0r zSt)XD^1ipu29?v`BLIRL3yZcJi>^%{OW)-&V05-M4mILSyVqx>5&^Rcxk$ZYJp>lx z#v`Xg(&5wXV>4r0C5(jgiL2p$C<~SJ*5hdQ3t97-%1O+7$1*?H`R-Z-eksX+FZBJo z2l(OD|09(?QMShLTSnd;&FUGOY2KV6{Zgp9SGM}mzjRx#{L`x2-_yVJ+j=Fva*um4aszpi>e9We+&sTB zjA;wh!Xrm55L!+!Dfn1jK$5NEbJ0`6{-j!{(G?wlC*o4w#vf-tOX5ivm|}K$qU(|O ziFV>G>!V!?qhc^83cG<8YjrrJxYM>-;?I!~EWzB3>2i zz1LOmSfJ?^rqfUh>({NE-KlWQjW*(D;;Sq+b%WDOqaA;K&XO4K zw03FEa$vi$t-Z88b(Qhh+Gvh8R(i7eiO%8-Hs9U1cKA-`U^? zLvOwM%bEKRzrOhqV7Gq#r%O0D-3{M625w4X-#P|_#p6%~%Qy^V33|aZCR1r3B38VH zD8Q;yY2YGOyoa2{`>hgY=! zd<1_msoigIxS7>{ZE(o{>d6}C4n`CW_Kw`Wjh`MO(I*yS(JE!{(O2htHfJ7N>rzVHGX}CWXq~uP^-1M znFs8T?aE|%>;PGL-Co-|6o^CKd=}4EIDu`qj+!usC|k5^*W{A0<2?v1`h|A=QF-&s z0Z&HuUh{`;V57I_x>~ONV=~NV&rPZ}H>hFklrBUP5wDJsIInm6)sj0qj-D2h;go&s za~T)%*umwYrmgArZ-07>&$Ml;v1i>7TE;NbtCh!`vDXO(V+BX-A^>IU6oAcERI(K5 zbioW%*%Z5;HL*lCH9f&Vmz&HwG5^vMHI`{{e1YUKFNd$Pca+hMTEzsKe&T_M(S8!m z*B)4LSzUF_TV)_p(#a2Fhd6zc8Jc|7RZ%W`i#82qPfe<-T=bct~04N9WTPag%^@w5k?Xve*o=0 zc!ic>iL;Nl#H=u@v?N3bkk!)rk$OA4{#;K6w!{Yl%)}Izru& zjDpEsU(=V?;^jnKjD`s0)6XbrZdGaz_msc7&Bohfn{)F+0Qtd9Q*sJ2P{$05Q0X26 z8?)+g6)e6uMm(C@kQAy&62Z|)`V?)<$*q*b3A>@^KyjIQh)k#hYnj(?FD;wjOep}J zF80pu1|7W_ zBYQ3_ND14WA1*4s8d|T11e?u`2V?741d|Z+`e_m~yC8TKtnI(>bX&luWatNrqO-F- zN|)~)UIc}*>S%_5>f0zyiGLADzaPhY!bPjp9&c!Z6L^DMDyhibwv2bhSM@Xh>3;M5lvG#{f=C&xG@I~B-PKD zqfb2xYASNW@08SOG~x9{HBV4Z<(m z*_4+6NI5Eq(gzUi6O4Vh7!Ls^Yi|Atg%7P^MK!Vu6#=Zdr_Lj*qIp7iBG_+TrJR%p zm_Ent9`aEmnj@uJgCU>uf506)Em2v;&2|Xc;HlRI8?V6|Sl&h0OANEO)x>HH3l2766$vBh&0F?E8hLR-xOT9)b@isfX) z#zmibYOY{6t-%%e;cOmCSWO3&y_3c~4OH+?O9VhUiXEu4`e6C5zu;f)Z4-8Ce{5gn zNxWy5$ZgG4SKVC%oxK|fIwMyOnabn($@vIp==GB8V4s$!_WKo!6>#W>o+?m0%IFM) z`K+^Xl+aQW)4~R)Q~`m^j&CoWL5-=P^+nKCk_N>@LSmdhOrgq(kE@hk6s3MrSb$Io#0z?m5ggk8KI zPcpwC%z)lz#u!hjAM)`D97BF`xNQDFGTzg)oF#JeXk?SFSffAh3k$3(%uKGHjBMP1 z_!BdB333P2E`GcX8#3%vZI3yLK_q#j1IF^SlQp^aXzA!5iFaqobz6d@kS@8(S*gZ=McS8MUS$ z9j6+bEfOz5z(Atyn_A^)bJULu#HoG3P~3Y}Cny{PDrH|CD?{cB9X`9Hf~Mioj;+%| zOdpj(k$=Xc&JMORE|`xil=(Jk#u?i=3QP!>4eunRw;4r;;KENGu}+mwatL|JVjrV1 z;`3Jt>2cv~2r05smI{+X@_4B0dp&GFx0?LACDXCKZGBhMHM3OvR^|WJRruddlm8!g zn*3Zr{A*9&HJbdljQHO&;(yDC|1Bf_w~YATGU9*Bi2p4k{qe^EBOH()hna0Z4l@%&OJF@X;CVPXmO;+EH$|U832EkU)6gNwdUc?R1jV#WpfomP zAuWnhXH3;?_3g~`HPF`{*#ey400wANFA@~=%UGf_2UaF(YV(R@Hqn3uMH@nTKu5b= zvyc5*%XKk`=NZ>Fbp}#M|qy z@h7bp5WlS!{M@|#*QTMqHO+sfcNrTT$aAQ(z`)Vy?kf8;)*pCc zdBB1(x1?d6=v%ZqH$2&QB|2q3xBF>#WVdH`sMxGZtIBACl7;5gwTizLw#D)1osQj! zVx1~Hmg$}8Vm*5u2K}#j&IIhsi#GdJ+C|QzMdt@KPMF$adxwC;`8r<0YC!F)^uI$H6k4@^sV*N0OEs>288eMKKepPMP<$Awr zMvc9Gw(?~D$k7Cq-XI?sEdq1pgI8atN1jlTGX64N)pM2MctdO0C>Y>@g=F6wzh3*v z>G=8`+{wA?{8rnZH=C+!8M+Bzq9~*n!o?W;a7U}GXXVE z5+T=|MCEMZpxGFW)F9^EiczQIP5%`jp{NaCs`bc5Y6at}Ab@Cy|M^Y;y11}dqlVl^ z6XG-&R6V9XrJ672uVr}Hg=J{FgQmJCeA4+s>Fy%6^ zWSnJ#b(ToaAGaTTEO?M-Kfc`i;X-trfSb>~pa>9yj4J(}!EUKz%11Ew&9Sl8E|yBC zgJ>lAluI-xDwXS+n&Q_NrkMc`r8uvjRlIvn-)hhjnYm&8?&ZB_cxZwi$p?uWuQ-1aI0PDUF`>6PZ^S zCA$Jxd&Kxj%2zZDbHb=pxui~w=JL{_h%SO5bnkXTxe=d=T%@Ti3|AfKOfEUvJ=5v% zBWAQ@&$PmpnJ{50Js)UB$km@~)-BMgBsss{2CkAcLUo1CkL`aw=QkS_Z=GQ-qjEHB z5Jr>*UsX}j(<`3W8MJ6DiHyKb*U7fv8iqBH-dhkoI6>7<3qVAEFf~>|nuYzr#z9OA zR;Odac1>f}ZwycR{fl9X1T_(CnFzIVMFg6HFf7vg;oS9Xn66NbsU_oy_tZzy8VtQb z>x{%XnynuD?)5E_Yt;Zj0qr@Wu22ajPdx0oXi4w-%K%xi?fsnQT30R}ksyCpvUbh3 z+Ke0XWw-&`uSPNhi1JBH*x!hXqz0%wu6`=lkfSX*fLslIt|RUfR7YuzT2mE9H0>e_ zPUb}INtD~)##5tAH9VYmnLB;qa^uaf$Y$)HHo$V_)uEgnuaX^eAt8f-6qO3t&&Ut+ z^&9iOsI+;Qt2LX(A4xx5WK2I5RD3MMfX{HHr!0&Em}OuC0=EpiLI(k9*ki8!$k*IENjKbuyp>E||4$h8rUs_bH(L~Xz29SA^T+#<7 zZ&{W1xVKQU^=o~7m?)Esl=r-bFuHGAwK4S-XQqh3;)NP2)m2}q{^|s=H&NUN23Eke z1ftA@bqs?Ez)Sd?H7?=|Fnz6C6Sg@pIwu!x#1NsRy?GgyIL z9+i;WRM}To=Lu|XYQSJGa+3b=0)h-rb$fZEp8( zl&5pM8IdJ((o~wTR7GBFoNZDC8Y$`Sq`;uAYY4HD(CdoxjwB;fqhZ_-+G=v&r2B>N zl5@-3*Xd?AOgwG1>ZW5UNDn9^d)k%_7b3yCkfj!b?o5OVPgmfpGj|gf63&&I2TKdX z(i&S0#HncpFJyRn$t^~(viX`a0$U3Z?jsPG*DuK%kFW~cPChE>y?4OkJDgZp%Depqbj&mNCexW< zAYYp1rUX#GC;o!96)%8@S<+p{*Pzn;eu~uRRF(7t3~ERvY4T)EgAdVGdi=xHn%Cmq z+n;vDV3P{I#ul8R2b*==JmWJmzKzi^A%u9T>}EmWteY9&wHxZmI&`Ie97*d|KU{8c zz1tP~UU~52MU54sOb_rAgSn_UtPBl9-6-|3Bx*rgxR@ziPdW^0(ix{Go$GsU#GU!a z8G5gyrM!WIu1$?5)RY>Bq4?#+KSa#bmmY@L%9urP)gWSfgO&A>djl`GZ6u6Yn$wR)~*#>b>uS2pN| zD<@Ra+tv1thcEJOW;p#6@8yZc-0THYh=b9HjaRwFOqZ!Arj=9gRs#b{+M$JR*!t-? za#j>|S5C9xMOdigI7V(&lEcTvUO?CS!KhEq-ZA4%5^wLAl$n7;$yrDMaQL21DCh1@ zn9kO8WhUTC&Gow4WBXER$G`%5)GifX+)Ou4mHhjC;nEuBUkpwX?TF+Kf&1zEV6Sfs zE$*USws(LlO=)qFE+vu~o0Fz~um>DIA0J$4_n!Mg!o-V8h1N4xqOeUakK!$Cwz$ShmoO4$^bhr;kGmXfFus+Yx-rqj0CeD`imUh;)Pb?K5Mg`xDQNYd++{n zPops6Y!f*j?@@Lt6nnScGt3vp)0L|0)?2LojW>X12Sj->(}%M)XBr{!&U$DiT1T9W zd=)Ty=RYGQBav;Yh@w(q*7UItQZxh+F+tRo4G%3M$!L%96xk?93cvxSD&7`Q-EeWv zw5Wg@P|JXL*Ik7XV;u9(qFwwA!!S_##|SzcrMDz%Tt-7HMx=u3gHFhpqqk0q2` zx)K9(nyZxZ+~Cf8yw56*x=btN_JOWsfuMWj1vYrJ%*#mF^3zS3Mc zDyWoCh>IWAQAgG3uxDU3Dl*YTlku9!ns#RC<3UV{qC~|zB0N+;D{1iNMbzmy#ro|Q z^yl6AKd{UDe|DSx57_8GykDm+KPSH365h36e%r5`nd?|;+-S2|k!joiB<0fa<6uSN z&%4F`y2-7L399Y3TMWGJrpGP90_0r?=)I1%pJmDZ1m4{a)$a$af2%rKC|C)?U)s#! z*Ru7k&6NIyHY=naJ2E&x!kFp@c_iMz3I^*=1~+_PMvc-!_T=vX`aA3K{ZE1Z;@Ib) zH(h|{_w|lG6hYjnmi#S`I{bn0#Q$MTzr8AO$G02Ff&UrbZqWVCV*ggY1?D~4xpm{6 ze4?$jk(%jm6)tbd|K`t4_1B%8KHocNZxlEGKq(Lgx%*Bj5b5@W&*~;)nBJE*oBn~r z>QC`JGt|!iCj@ei!4{lfl7VW9z3h9m@Zg(Njx+%l&noAr`vf5e3{@6nz(Fd_y=-EZ zu;6@a$Zxuj8m0zTW`9Gr^e;3h|JsOrcSf`Yhg5%iU%hAbq1ePNmTzYK8Ov|o8v4V0 z|3XA`H;8=hBMLkOpxz?-c5~-DqQ4%l_3xH3 zMg6@lq_Jz4sQ{CSq`B_Hv$(>F97qv{kiB^7h%ngYD zyqM|J$kg^*&wNvW{K1i*#fUk6h=l&dU9>+~b^mq|^i8qt`^}PHMDV_4x?FL^CxUt0 z22ThI9%ocw_(ZVN96NbF;T zZsvgZnSWQg*HENeOB>2p45?|7H6btGsVV#Xb4THXUt=@7u~)3FI$G9QOqL^i{ETDm z$!DoLY2)d^TyN=ap$}X+qW&)OF3hj!=-AcP^Qd!b%NpVjll#)d`#i-PmK#^08Q?xb zD&ZtytG-Uc+j-Xiv?*x;Hru?W{`ug17U*6oC z%6?Wp{ndg03(}{K2snzjO!{?(|4(S2f`ju>cp!))yrHI|a>w@5^~&`i(&3S82Z{C& z!0G0I(cmeT*aEqK%jmk`8u1n{&RwA1UIV*#i}Y7{yqi?~?MwIvHxeW21;IFNH81=MM|MYe2$>vmEybKI* zH;#VBe0()9B03~+`h#=hrrP03@00FGI<;gyU~!+IlgG- zJUTCf&*7-OXi751u-g(EOQHpf-n<=wtHBO6c`FbF`|4yvP+ zS`FKypJh!>opoW5{Q%K0-OgrDn)T6XVG@)ra{sbe>`|1@>i#^7vh(>W@FjV3pNr+? z;_SmUM)G6tD=_<@#aY}F8OH`QY_taL%EkaJ4ZE}2%hw_Jz<%d5x9TdlOY|KV=Iey4 z^Csu4VhG$?m*-GxjBtJa7ei|=Cg@k%S*qc-$WJ^0aRAL`5agzF26n)gu?N=X zmzF1Ib5-ZDtN6b?mHL@af9apQE!^$J!#~Zfe~+f!Z`5`FIfnf#jf#H9u>W(~sDBgM z{~dAEC@r6?+gSZ3X8C6QmKX3uY6P^&L;KM~r&%;^&mMryG`!Si&sX+o z1#@)ZGCFx;$C2&IF9uj!j~%+;>~I#GfdiwkSzCT&!}Z>ogTUn(5&OU z{-Te{Yc^P$IKHi=MIKbI++%FIr8!zhH@dGe+Os)#EjMbVRCHY_HY!x#)@f~-6fY1L z4~Kgpn$`#?-)aY>VfBMBpIITyh$`43uSI3bjjd`|K7&F5mBbTyo@FG1m@wzM|)D_l|t0 zvZ9gVdg8QXvtmDbxyZU+_WJ#Pa>3lznZgyonVrDbvf%RON6oNjTzb~K6<=s<7Rg?R zufQp}aj!NiUU(moTVHTDUQb+g&t6w-WF21GQ(PgBxgB4nwYt?@mz@zWiLIU65CHU| z*d>9Jz(CoH$*bzdYqxu=qMy6!{vHqC(RASZ^bYc;r*41G{PCxnsOg_}VQ=#p{vQ(X zFGW!a?yj|-u}3@xX`2%#uK7SB8#~djDJNLoc~oZ@-m$mN9W3&sO~SZPvt%niB0Y?h zG!;$Pj{OCiWrXIDq8x8rN7d_h8p^|=655knTaCs2X8i|?2bZe0d*-$Ws$8GpAA{$A zT5~4=6Y_a$6YBkzvgIWK%LY%x+#c?{ME<91Tc^h2HI6gP7YoQw2>?DwO7Bpof`U-} zgRapLHwH`dU8o-CW;@{>`@5tdvWF$MvD5)qKebA5V`q9k9iqPiTD^gguoVa-}ckDbqv|FHRp)gpPJwFFOBt~y znSpo#17A2YpNG2aV6U4DKTdr$J+m?l&$Un02yi2Rw3ivXVy;+x{8g9_4I6a6d|fwP z3OWA))uW14AM6D6NtsQtl#0{BKFT8r_-CY*vAG#9Uz8j6pUHWwxu8>D`UMJYYMtSe z@HefsRRJ_*2 z6ovk7wsldD=>=ykaiRzHQbG%{Bi;csbcCeHOdJFKb018EA!w898CND<;2&?JVKHI4 zXW^1cp4Xt2Tl7w>t%1a7c=vH(q@}2d zU_`~(44nr_&oM*EiWRcAv<&;C>TCVs6R0Nz4_ye_p2zcv+vpPQh`Zo+fDe-Kw%bCY z8$0>r(1~NZmZ-*Es%Moo_HLO%q)wEIXVHMK>|DN)^S?L+=){!KIvq{D2`$dA@xe}tnOk`oLE>TbB;M|PP zv0?tRl;n86(oElZBtWS-9@i|;?=#E{HU*NesP&H_H=!Hn zQT^0s20x!ym_3xXUXX*ja_$i^J8DS6UFFA8oG%u24?AfPQ*K^n>$3I6K!oJ9?13Y3 zvq1kW?eZP~6gZ8qF3SOesRJR}!<}BgU*4)$d3h3(9Y|U=dI43=U+C)9Necf|moQVn z4nE4bjvPR5;@&QYA(p+1%de*J z22WrG8Cg({SBOx_wjVfX1$-Hx@o)EN;p59%cC!#ULhfhb_Fg|@^;&)idGK+i_((y% zcxcf9y#7(ESOwwwQj3U?r1i8B%1CZiaC@{Pwh6JFZkQ*0AJOZ}37iMRM_l}oh7j4k zy<~(IET=`JIgkPHZ|m{8tiAIRD$lKJCjq$XBP-E42Cymi!yDB!?dF}znc8rkW#H@^j~Bn$IEMTw>Z9G>!$ zjh%20+ccy;<`V`;wW!R;*KP=G&YloTH!3#TUJV7(?&VXjsj_mg#YvJ|vku&v@|4#< zVTOZ3TCr9DG6e>Sw(O@>Xu#*g99^CY1(j2?U(6_d@pHTPTTZm{;DNUZ80!xc@ZZ!z zUA@&pjl9u9g#>FO6r|C4W+?TPmM6^ERIHUrkT-2j+%n3CpkP|>xz2mbRr&B@lyK!m zlptW71$z6-Fj&Jd0uw;vE?Xn(99DAwfkRSE_zJ!-iGu@+eVyM&+a0+XKeR9&QYYFp zk5_}Ut=mykh+!3h@dY#2r>K=vwvENLi->ZyYO=B-0-Tb(qMSqg``1xlGv4t=+EI8#A$6{NvE8Cmnk4 zd~kG|D-LU(qYJf%jdxenwGG_Atc!20PHPcct3~C`ZjL-R_4{%@l$e1L5Q02lY5zQK;ht5%! zK&0=6k|sqDXatyeUBFhncj=+4i$NKK_tQ?Rz<1RYtU<5bpVS1RYUkU|t0d(T^49q) zz0H0Jy?Goa$wRi1YW3mue)lR$q+YZ29hsF z{zF|K2oep~B1s`zbH5W>vyoKI*IWfk2}(AoBeIXXGoz#MvLwMTLJFF@8Y?PlxP?1H zHrC}8BNpp>dl_c;K$V%7H6@-;!ngn)m=)qY>yL;-ji0%0(ydQ<7^aFQNIfn~wG+9x ze(aKC;c4P0P!(j}&O%ZoGiHMl^f{s@SZBSN?QrjyB4ZD>#ZwJ3!a)?Y&U1S!sWWT!_3z8%nyCzYQ5XE zN!?-16FH7`f>DY=z1$Yc8GRyPo7|uZmyznt?#kB+G9zanL@z~z@pZm#NEm(osjxLxcdx1=CP;lG|Nm?Vj+J^dKJVA@8bIhT*>(6`OMs*n)6AOb#hbY>^qPb^EhlA-aA3* zd!zW;38jUcZcBaBsc_j_uvk}gL?wOo#0Mzcemi~ z?(Xg!^qt%Z>2#*M(=+!?-}~09BD~Q)#gD!ojXM0gk%VmARv$0dg-=X6jLXS`bhp>yuMtFxFUCuE&X>TkTjl-c~ z5b1HXf`f>HOZF1n@6E2hXtEnxCY_FgUmOYs`p)QvBvG);BsAf`4D8*ew*^Bv)4nP- zKREA^*3kRpLw(6SmNi7W4sID)E7nh+G?X}%uNAEEq5DjZW1ms9>2D@XY@NXCJW-_U z^2@YCdOKv%;f$?K&tzGbKFB*Lt(ql!7eXU-aB%sHZv$du9+c zD5qg;i65o0DV5~mas#7W`gLYwyKz_NWD&^@YSxjZYBblhX3T^r;#0@65iu2(qvsnU zS1o2TzGkC=I;OIWCzv0jVxPt*H@mQ_Cr$nZAvjf$<`_$(#sOpyp-63|zpksj zpakEzty*`U6lK`X0GWI@* z;X$Ly*Bb2hkGh{M?Q3RWe69_x?s&V^qWabeSlx?g?9ufs$7~K{pShDlgq$~OOGV0d$r6OKC%G!zyQYImu|;~QiLKMlgXy%21u0{N zlDp5m^F^F*;!QH`(6K(@g|b88256RmW24_b{&iAXQ|G!85FK*W($*1$6+{)dGTL*klXNRHT|u{}M+-Z@ul!Y>6Xse1&U&HG916 z&C9FL3s34OL#m$;mU++&;%u)nUYqM!JNLq%!-HVPr<1AZ`JyL7XcFjGQ%Bv{R zo8qj@~oX{?_Z{OK3)9Vtmp3=Kp5%r!Q>QXx6W1jA{on_ahR0VUE$N4{;3i$La5leZG_Zg0=H z(~dNoW&zGMFEDGRca=VI3|b|TSl)2?wmnz-MeNRxXM*o!!C#6EehZ2E|Aa8;dsUKs zAF973YQLuj{-AMv#Zlfj&h~32?BCCs|4Q}zMoK+Y8o0>a0QnT{_m!7QVc;asz){UQ z*vt2Yn);9b{*cep<6YGHd1HQII^JKCy_ebpzcwbmukH&o)pL=-)mw(sl7)?o z3T7O^gg19lsJEXT(p;hS#t}D4!|F0~GxK`MteK6MFk5SPMOIv~I!i7ss<%ii%Ka>J-dk=l&L<7)7HwI}B`cX~v@&ms!lUKne1Kdx;4$G;$g7S)&spN$p1 zr)2YVwJeA+fWnw)`hm5KM833h-$tx&XVO1_r#|?a0_?bN91x!R1ql1`<>rEdjvqI2 zYqEDq?9$Q#K`TFE3lIpqq;SzZhTH8S;%*DRUN_B%c{#Z=HeH7(aaf!e)UJWjMnzp_ z-Ih>&0M@%Z;4yZ$ma~TL%o()4w&#c~VWhv-*F)4-+pCTS?%ndRXEIhA(n$eUoV1f2 z+;IhwZ%=3fKf~ayL=!f18s0g<&Y*I>!Wyk@pV}IVBcbLR%9St8v7fV86rk1RDaAXl z+in$!9`a-03v#HB{n+_{AJ-H0Tx&Hvy8VGJb(ACZt?cy|xBLN5RpaZd>vQjxo+=MLdWoFrAWmhwBd<9LxgLpXqSHsUVy+dv zXhAM7a#!8a!~tb$feSqRVJTUATDthp-Tk#)zr|FA{n{?8{|f{2 zziO_&)Xx4(oY)^h#&JZ_H}_5WMWX#>6#g+;)#Q)c@f)k-kIc?rkyUB_QR(-*^w(7t zWdj#9R~$P6WiV|D?#Dk&`PbotmA}z|pD

zhVgr?gMOq83G6`_FT-Z*mhPI17eE@ z#9ODti;4*Rxx?5 znT}ZF#JqAn7qRxSstpp~N?$%VCj#iOgbF_v8&&f+X{*)$^1kNBp8J8e8YI3?DZvpFb!5j`ZD+8yl50XfH?=$jnnK7L zP44W<***~DXIj%zCvZN(`e`&}*lRDX%#)72aNIC(Z`ywO;ucK=kq$$Vy0Mp>)Q@sn zx#(%1x-Cn;w(!nlrukF~l9c?7)7xADA$d0#oE_J!`pDCP8&z|F_8hUO&7OTN$GbPpTm?Nc77&Jplx0X^0CvB#s1j-|E-Sdi8NeVjK9Wp0fRBkwNLv8UEj z`^_Ij)}H62t5c09hGMxeHJ_UKewaQCLS2g)!Mr?q!^(Zl3)fr+l5=H#+p^~ABafq| zu~BnSLHgG7jxTPL*dk_ccXm$y+=xl!c)Sfh(ZxG$-oo1&OVwr4rHp7;hB00L%JycX zGQaMcU?u^2ZZ_|fjdk(bO9ASxe$sAu?#^;F+%XawLG-)^axZJFh^?@I{gl3cY5n*H z0{iHhgbh#XHbl6s^?JPhtKx>Rg2UF0>VymjyC16wp*-_P?qlGWYJhJsa1VL7UsUhO zs*?8-RL2$Pup{PY3z;Y*6d2~@DXw10wGFz|{u6WPkhm4;Wa<2%@QDp>!%h~~EzG9o zGViW5CMiUzeQI&jc)H%TP_5W*ioWQRDa0RKwfSWYj4vw6J_VbAa|RoOqk!MBVg-(! z^u8enue?SsjZ7l~;S@lmsG@wBe?II}8XF_KJ%vq6uL@;kT|x{es=u(1!`o16&t~Iw zUEM?z0V?6SM@;~7?mj>^a}w#JbTY$U!6A_z*tr+K63y53PHb&o*AXe{oi;Kq#!Nt# za|Fk68aXpbnc0H0vN^F9x@;Uf6&Wus<;vii=u_G4DvsCYF1~9OSY&I<=89Gd5+WzJ z%$3v|aP>?zA127ML!Rmk1K5r#~xT)473jAGc+9G*^IOfP0AOEir5g>Xs~xVoaB z2px(MvY8STNFf~H3WRDdnJ_;Q_3I|2l2+XkUGjQdLpGYFS?!uAn>XA{!EQ-Fby5H zZxdE)Pm!@x?5cb%$G{eYtoQX56_A3>c9Gu5AAkpdV`Al{mTA=Sv@`BD5`1D-UnST> zOkPSL$WTCGN&+;V_B>3$l#0-GR6ow~-~`08^pcTgJD1Q03^G0x0jFbm8;Kjf>D^v^ zija9Sb94zu20B#$n(tlz!mBX1Hzaw9%8qwFxP)Bb$wSbCEYH{Uu;r{34B>#My0UB+ zw60v9WZ4y#f`MAz;En6n0i9CiPLV)!Gpd>=VyoA>828*412C@APUm%*&IDuUlY75YRj$AB37R z^ht_^@Xw8Fc0|#aDpC#=xE$wtu41<-^~uRA6?(1AIaE*bDJmtn0XM-+BgbKGl-Emg z{)yg7G=q7WmfYh+Ufjrdlga>`#cI{*G1FZ&)F5<+nNSAsL{d+;7ke>AOli)g*V4~%HvKKBIQw_rT}&C*x=n*NOgDn; zLt_{71nm3pKWZTo5KU7t0$qw?xw(>s+mX&2gVJ~_HAmT$2`GM88W$BxWs3f^V(a~K zq=Q%TkfC{OM>z})l`a>65AsL4>i5Me-<{*WgYf(@UiI%xs9{}v?w79{@`B>xs9|34BW(@TW|?o&5^V)OV<%$EN99A*eD zj2jKC2f6QavIWi!c`;N>9}An|N*X3KycaR3w5^q>t{$$6{?3o7H4h`G>*i)1>*^UN z9cs3{z3j5N6+nR{50w2`!l^6~Yo3@O4kafc0V_&( zACOdx_KO>A0jr{a+1B9Zst8HnOAw1=QDMJm6J>uqx3lNKlGuYy-l=f3i((H@xd=Ii zG2=;OxhtmN@Y`87RH8U)*k9hM>tUSl9ybm*bFuZh8G&-hy0W-a4Jgs6iE)c`Q*bMF zrLe=BfH#Jae|r8T=n00~i^l3s%QssIvDM>5K}VG|Y^8fH&pFCDyIh~!p=M^}zi^d1Llm`Y7-+6Ox3S5<|uR53DbDr|8_l%_gk}u`=SsNWc3M`w7VFCnh z1~~vCz(3+~h7#9ek3QaxNv!5^4BQEEE@hBgactZ^cQ$a@bPC~QXehdsxnpoGx@oU> zwQ||K%>bmH(?)Ria?Vnxasni7Lph3$k*+vCa>yp`p=>4g5g*?DaA(wBw5s#Pi-4&vHOqhR8+^A#{n+&>nvM^ zI`4;%Ox+R-P2w`f)8*=yDUg

w4RxN!skME=S75ivq7x~}`9=hI~UZO4c6%^*-Q^B{ANA9|1dqZ#)xFv}D5g*D29WYz> z`#`-!tfrzSZo3gvsoht!y9gL`T;d)?Te7J=0vBTG?czr8nb5!|Zaj?UQpB=IW=V@7 z4sRf&LY!8P2%vEk)@slUxylRmhI-kiR}qo#BF&{T0)(D&m<@m~af1zDUcYTeg72ck za~?f9HG+Ccz+MrOc6)Sc2tnLEy2dV00gtlud8Rag5hrwLza@II$i4*xgY39IfG3s{ z7Nx6s^>N8%9UfInK9pzC@#?Y%^L#6Nc+fbyF=X|DoUo+c3ya0`E5DHPj5j_It8PLi z5PB=28%Q^GmIwG+?g54oMW%i;zz4NM&&X0$WId<)kjA6tydEvkwA4Nam+7P!p9mG@|>Qj}RvIh86 zVeBPWP}tXDABG@FU%mu-NV|crO7$R?ZmH>woMLDEPUL@$6$B zR0^n&xcAqwPXUO*_vu_^PUR6Wj0x0e*z`Xs2&-TMni>pBbm=XUW3HE-Gq{fI6XmW( z=Y8>BPM}%>E-lx2pj!B}PG}u7e8*cn5E055Bso8hSorwKH1y-YoRAD;I7!$~6j3o| zICnLHALds9J44f}=PCF|Lu)Lv)2k2cjEC4Y(nmO>@F+KVG38XpxkOYSVCUIASXZys z%E{g|ym!JW1Hnjc1^21WGHLNx>luN2amw7{0c%oQe!5Lg5I?bG%_y$s`q~{TQnR`1 zl-L#yWvFbWL#`9W8-lj=jXM>cAo}Lvb3wUR#)Zks8yIk|Ndc91cUL1Mu`?=5+)mb= zDLR8#Lpw^>3jNltK!X_9AY&ZoE*@K07~g!fj->9|7Em@2c03ia4%>IFfN;E<0zecs zdZO-Hk301!FL5&)HclSJuq<&a)@{TUnFkZ}3xLOiqye9#tx#WzrZ*_4q7rLZLVhgC zY>}*3mZk8%_#J^?vRof49Wbh!slQ+#`xze8`R!+RO9B2G4 z_yc=R0WTbeF{aNmdqbefN)53c&%uGW|6G26cs^T zpl{(*b61Qz33Ormbw<@lR4D>VcJwhZc2EZ}$ULw6V{$$U4s;qfaM5{J~7}i4|qfb1KZOsmq*(jptOUB*&14I^f`SxkJk7*W$*7 znhXHBS8uGaEZ|%n9Ogf-RJ3^5CADv=_nHMaS+@;*+D`;J$bYjEc(Gb0;cg1A$g$A+ z2HpPF2x2ffXLV+0JCzR3H@Y7+$a5qCH z&l67hfDZW7=U|=ZQ%hCm@o9V1XT%6Nlnmx4>x4ts_bqV^8s4+H0v|M#!KW{=HbX>Z z6MxfEV_^ksW@Z;hUj4f87F_RR%I+2B)E5(As(Y-I^J+XQit8DUwfV;SeX>+mR)r|I1(;-zpX!ZkIc@smf^X3x;Cd<0m4bDsa)I0tqL znYUjzRQ@@Z69+6xOi2}GKWE2PPU#Z;6_*n5IZERqV`s`ELMkNM8C2Kf%gQYtOitUW z;3ccpGUuaPF9?rVhR}@_6Rpzw@JCxb+QVAk9P(3I;DX^lHg#r?{V+9xIt~Ij28O;M zjQGv^!I{buo>Suywa5st{pJ=6qvw5&Lc52%@~ad8yS&D>e&u0SPvZ>&s#HQaPYmKY z0xDj2hLXdOqnX-m9D7{+0ARFbQtk|VP@?wnJ8fJ;Hd2Z8Cg11mblsfpJ~mqpQ`KV? z2hN`)TH?HHjx(=R1cjIUCWhBxq+K$p%h2%eTn5pNc{#2aZ$_qgRb~<$ zP4v~Mp)7Fgf~ADoe}fbB2qVQ(M{PxY-2 zbf{Z7O?fc^+dEsQeAW(&G8UjD(_gTSM>S&hEX=~>XdimqZ`bB=C~#Ew=o?dx#j6RZ zeDFH4F}^PHt(4p)mr!4$cHKAXamE7C0gh@C`n)@g;Ubwd3XgE zmlB#TjJQ|aqPr1u%A6bxH`;~C*fYd-qVnt(xT~_IvlnMTC%q1OgC0z=pSZ_SErxS0 z0sOwxkQ}USH*{}^6n(gA{@FXAkitv)>y>f>jm2dWPq5#QpmH4(PCBn$vB_VpAR zol^pP#A9g|dRjoKg|B9G15lTLlC&iIOsSTm1;9vEHadm zt<~rPXC{V@so9N6Hcr2@DHSSFV&Lkz4Y%lnL+0Zbvy%fM+vTvt{VXagqaJs^#`*E_ zW(N=On=*8fInhkR`!gtI1HFzA9J|L;D8p`ItuOQJ1!IWfu^vVru?zsCpEKyzu$T+U zV(}2LYBYzj?G3RS)1ehB3#gaqFrPpo0Zt^DeqPnQD&w#yhmW_9)31kkRlLiGDRDC9 zo}Ut^m-J(<-Y6IqE+R&;)UHIXx3soO-Mr>LeeDfFx^kNp+8!imv2NA$4)-W-4BZ%$ z%?#ms8RfL4HMUaKgz#xNaN?4)z7{{#mQJv&&A*MTbv1x`xi?*2dpxNirm|7vaZ)Og zO-S{Dg|h;}<&dBc)7eERy4Ix>1Ov{i{1hx&O;jar$Tm?^a-3U*#F_G?zEa1z*+p{W zJa!Lo_aqr}HZ^lu!RA`(@s9Z8$;vtsl`KqMt~htX9aWL#*ts|GXQI}I?sO$aU{R3l zZ@ntAfIGECHM_Z4Xuuty&U`XI7C7C)S%~Wz(5D4gnvt^qf(7sN(kIytr!QP^clMdtLPXWWih3DRKil|i8EP{ zM=-$kHB)(X^=rO%56C2um|-@W<~jy7W{t$&3Y?1^7-DCU#`iR7@8WIXd3?xoK@)SjPh^aH28%+N(2J-Zt#b@ z1N&DW_3w@OQsMe*V}6^?`hU?R{@*ZZ72Tus-lQ+Q!h4hc3OsAD`MUwgAMvcF(ZwY9 z27TXQzYS>pZx6#?It+y2$6{10WJ-a1Q@+Igy>I?1I_m<^KR7?$4f@>Rx^!>U_Z{}1 z=CdOIGjKDtDDm69o-fyw?)_u>pC_~eL;1(vV!jIEym!p^gX=G1v@V$cXwdiW`2u(R z#XbLFO6zYp=O?!O=QypuVaxZN)-PPrd%ye-l3Kqkvfr8U*RxuU7Fj^ zH(0)Ljz520o&NQAB>bAzN^}ov{c`OF@$=%lM3~1oJDNXI&6lL@*zaxnK7M~hvL^jG zvwn_W_GEcR^j}-aFA&*p*MWY-z9F{sLfxD7{jBrP)YN;xxx`lke`i(y`TiDY*mrc> z3r(10oQWUJ`ksFJy|?}c_@{c0g6Z!)^!*b0+XPg_Kl09>;fwu>PigfTjQ`0){7`9r zEr#c!d{RGN)Jh6pJGi&)ivWHf=>J4eG&I%G{qFrB?|}3DmtQOY`0GKQToLdp9b_HL zCDLrRxgR6e(oT9M6O+&(*|Q}f36++D121+^(8ocM^7GQ6Yo~;%Vu{_ST4l7Yd6hLp zRy7;xW4(c8*3|U)4^2)55Fmhy-j~)wEpnX@ExN|yqUk(8fvWNl)byL&xn1?(%U-- zv<$vIRy~)hg5KpHbl>lL^-p;2+ll3#L;BbD^Y_HjF9)mpe*Ry_9Q_01{`XKvYgU~& z??(;Mci#M`QS-ap(RW|R&F{>Ce`0d}HFq>c199OYq~vF?wCE%T`Y`OBwIDYsnrHkE zsb~gP`vs}h!K=&)Yg`fifuFu@Gyr_S@%QXs0QuoS^nFbKsdoJjNA-o1&^PS<7h2+d zRR2>y{64VrCmr$g#i)({+htQ7Aj3QNUGu}53s6T{dsJ81yUbguo0Lm|)1WP!QO@k& z`@mg_^M{W3zAov{HSAkY{571)mgdjM_alJxcRc@7QYGE}HNlSnMntd&?KMOKLyez! zWZ&@TreaprC-bS1F(ASD+8Z;YhX$v@ciyN-e8AczD61OcK+>A==&964rEiKtlg(5d zyM?#{>|Oy5ujYtXpG=;GEq5Vtk#qZCCp?VILBB4}!s_0BRjm-nR7^I@R!lad+#U6hKpM0ygo?;wVOabe`8?FDy!eRMD_Ww4 z5Fx@ha>2Ac!?q{^_*uD@70v?~m|+c{?NOWCKI@I_njQk2P2MFryP;5uzr4_Wkd z3>Ns;Yq}N9HY&Bc`62Si*SZ+Ln7OsE7|e<1MT+L70mtTJ92_4LO+oPyK2PN@0o63w zl2#YN;iqlyusSe8c_STec-Bx1EJF9a=@yTmDwirMXlmsk%)k=oY2T- z?ts$9KQ+0O-L`?IV>-vaxJC-fe%snTf=!f-I^$6&`En)gNvSoQm4S=VjC2_0;}mJH z=x2vwutK01&aw8A=wJmFTD(bo@^=UElC%>AqiD=I07Xt|`hxlcv@5@1m=Q~efsE!Y z0~(yBeV!`ylI56()%IMT4PaGjWUFS|)f%VLV?M}JA0$Aa?Kz5WGFf!h>hipOsBvZB zO6M?lvNWr(syahV|FaN#CmtZb$|GLi@hlz|t_$pJ z-p9}wHjo2f5VClIDs@VN6ZxHpGP`qf0hSs;t9_^w&Ylgb>$7u|g-o7R(jjRWehN^R z0oMcE_A!z8J3I(Tgymv3Yzp`tyCS zTgW}=aTOMorKu5lbrQtUU{(~_*)w>!dt}5Acch7mZtMlDtaMS@^pGLqF+dj;cUj5R zr~QJ+hUWE;CWFaB=~qY14XUqm7@(0eARb*96#qvH@vwbO}2fSABojr`qMYhvTd zIabmr*%J@MXK>nYTZ;Y2xZmfc2^_qD<2^i$bq0Gpz@Uf|SLhpLOAdt76ifdH2;7 z!zaJ5OQMcgISKYWd^_9+nDleHhWNESh?%eP?Upv~?@Y$OORZeKr&h`WSP^|6$#Hj3 zrb+8m4H0xB0bhh0Rvsj5VTRZ&IkS^G*C9X4H3&u;suSr<*<(QGA5+vX2^me>!N-f2 zRy)^l&}X)OZPf?pvt6RdqZE%<#zTINENE_Xh!24G9|3DaIfu0x0Yj;Cx^c#0$EtOT zSG~5**-X>T~}=rrdxEXBb46c zdJLT@UNd_mDv&U_h@ia|J(1WbM2)kq_0EuI&LP}i#3l&-LiZGwXvUWjVfmd+b6~D{ zIE*{LY38Ldni}avkgM|P1BdhF(UBh3B<<6~n95Dq5gzHz_=HE5LNQTng^4Zgvck(F zd$tpFxTs~LK_0V4;h;#IiWAue^_ ze5j`cyU{auBEUo`6{&)96$q6MV_M7-^h7r#!lP(fR8m~Q*;@7|w#!{be4Rm+r>PzI zw4cx(rSRW~74&np^$Vo7EHzG;dS5WddHEMw=S;DX1;0H!Y6>hcq*+w!b-3xiqnjV% z+!U!S6pV&E-K1*Oi-%d0zB6Aqsj;rGqCo(KR6sni(OFP^xO&wj1D(Cw$yMZ$26lNy zYU-emrPE?qj-3(Qu=1kk3CkJWBpmmXbcTH; z52FQ~9p{r1^t7tlieUnkEO{)bw?mfFMWXhrgZjw@$GrOJ2L+GCSaUa_2Uq#Dn(MSW z`pDnRy;g@A=dW(ml!rFXXIHXuqkXn51u&JEgxTo7p)8;)cvfxoIDP2>0*6;vknu!E z3&WwYM^(F!%B0N(4@MP=pIWH=)Ah+<({OE5*Eu&sA5TGnA+T${8C{tKr(lzkD~5pA zMd50D9pv>M0+6R4=yMz`@Pv~?*AaaJR%?XQrTNUE0iB<$ZQdvI2a&;Jv=J|sVY%d* zp7a~$TPI+AUL9$P3zD{iY)3gRKq*ZbIQTHzP6KywO<5YuBo4ydfukbngd9Q0s5Z_s{HZq+E>ZLh z$p;+;?GC{V7y`zj{r=VQW|uLuM&)?lUVwTedzY5h!HkK`suk1A%`9IRgd!TMJkzWz z|Hn0~u9)`KD(?u-@}a}XUwCs4%o5CuLw9zc^LjLr(m10Iqsz{|4~B)>-)ViTLmj)m z2i@hCI1@`f8sXldUc%#FosNJkWrQLuGc4I@V?Q{e;-8Pf`AA2Ir|`yZy_~1?K{pc4 zI-`lU_8VK>p7ZwJ0XoA{bvwLQos1a}3|uUt9o5Z^sf>k*E&|VOR`62|Nfs5~BCrkD zDXj4ZYPv)rY%9#b7=5HGrjLE^Fw`NpzxFUV+=sL&C9Uu!hbe*t(z~r)v?#WAp3vK? zrQTGP75=WE0&EO7HplsK8QX11$K_j3WljtlE}7>E4RbKK(To+|OizTid12S4ld<_k z-&kUKSZS5%`#LfdTfenGCk`#LI4md^h7Xw7p1VJcb*FIqy+H+sNXT}}sF z6f7`6AI^>}ps&|RrVZ(aD-flKgKBXGOiShu-!hQ6c{N7P)Jp~ts2~Jv3Cp-kE+A4l zpCw~-IEoU9VoJZF7z_vN@_D5ywD7WA_=Wrv1j-;oj<(p8`h5g>8(djB(hb`*GUEkR zibUbA_s@6$Dm_K;Ve!s2vJuK2i7||-IW_~&wVAJ^LqDR*=lJ3kW%K89 zF!#tRvqfyUK}tpu0SQjZKB_3t^Tyh^sv{d?T5Q8XGqM)QJ=kHnn8LgK>H4(@#wb0e zbp95<1ufw^S%c13d?)JqwL#B5U*ulmY07^XFi@j5ZRc`pSEME$2Ro(rbr9?3op9Wd zt!`B34nq?s4RRz(@wbj$LX!7iv1-lEN^Ef1<#(Xj4tiP1kcMv67J<9lXr&a``}~pM zvog?<%3w76G^b(#FfW|qVtkIg4W!E5Kmygqv%4qt3)prKN_sPJC;z=4;%ck9xRTT$cp9S9D%E4A${7yHWYcD5|cbfG3m= zcPyR)X82g(WAJ>-H;Fi!9ILUl+#p!gZXNVRhQ9m1TWA8r;Rp;sxez~Y19x7k@I)kbdDG!`KjOZ&8loo~GYTYe$w`=mpb8jj7S z^ZkOxl%YG03rz{UrG&TGt2N^0gwUeXjSz-s^)lN*PAYX2M{(@REFkQSdr%g59AQg03}$05XCg3%0r$d;Lx@%-8r8=)Rh*;j(7Hg zvC3FI)=C3on5B=;4EJU`9E7`arZlgPG{&JF_)zA*47!UYaCmZFAY8yaYgZDr?#dGg zferZB)iAwmPAPvKW=4WAgZpr`>LCib-f6}RGo^w4>XF%?u+%d5XpBMwOiTV`sM~1jYiI4jdNH@u@atXUL*7)J`y>4K z8S=k`S$QDD8k@GCP_aMAii@>W+2XPvCo%6TXi`>DlE;yFzFJtft^IhcTPg2MnH}fS zY^zeIq5Q69@0Q+nw_sXVd~VZnsD8gM;|wEBVFzaxKvC$i>{h1|bAqc{mvfdaSJYb2 zBR-*zsb-b7!jpX?Us5@v7g%>@cJ>A&Zo+&N2#|e}oQhVn{`C+xxINLJ@KW1Q&zBB!Pn*0X$9*8PONN!`1xtkd%u1c@D8rzrx?VwJ+k%X3*`@yqm^ zbA+ZyTVp2u2CKJMqqRN6XKnNWF^md*W%KqTT$P5)x>u6&S-uP@#Agk73FdZVI7Jot zL$MzTs#{&^M2^cL+nmoBM=vR-$|NQPa{Fj9?QEv7f5;$xFUx!{zy21r@-G!;2+0Up z?!)f;9OKXVllNx+w@Wa7r62uMnZXub^H^u53Z*0MeNB z)wZ9ggg?)a6McIz^Y4q7Nip@Ie_-El!505ymE+It^Zd@fA8Qq2RY1PIYWI6_>6e)M zZGz&T?BJhgwZHZ57wz|3P8$Go_QNzI?t%W|-hsa|^e3dmcVBz?C$)}0_hJOgGmv|c z`Ul$@#G(!FZA1TVfPBN({U*XN=c{GEqs;u=b-&VIX)dHV2TFNwk5@KlaL!WvB0p5; z|0Am+ew z;Eba}(SrJfl82|py=X2xQ=b?&>tBcIS64F(%mp6fvr|p=zLxswC_6SfxT6eigv5kSTB+(Mn95*y15}mC4GCYGY|9nnNp(qI|FxW z#)FUtro+Ujtm#h&u*gls&8i?OXKa_+t1(G*M5?okx1KA?acEu-KMt!`-8c&vdlXj$ z_2RU(Od|(PSKp-`$t}LSVxDTzJwm$I%c+cIK+-Y@oYBAF??sV%Wpi?eFQx z)ZuW-94I3lgtmx7Q+1?}RJcEeteJ!<(}3j;8|Ui8`{LLAYhV=He7R z1-~PHslii1ses2LmG6?1Z8VOhK}R{?83_KayENlAAMmJ$S_oRjSDU9_zwHWk26E{t zC&dJBdf|G_eO~kt#4BNk-b5Zk>j76@kC`Mn=)o8Z>4E_A#a{56Mz@E!9i6Ts;hHvX z zT4+ztc%CpudxHt#n>8VDP;{>z%mxy9x89M%N+hrJ5TC$KzSxW>3CQ*0UFNvjp>;gb z68)HTrJcccIUyF0RA1eh_<@NL`Y?4^f{3N(tDMoqe~9;JWfNavD`vgPY)n%Api?ZS zcex8TW1PI6Cn>CfBDQRKsA_32AF`LgeM(wtvtv3{AaZ@U(5pM`63y+Za`Wso1%jlG zlc*=l?K6kh8&JS|f*GMVtbGPi0;TLyj9%AZqT}XhtFpa#GmGOn*@WEG(3}?bg>n8E zpqX=|j^Lu*G-LYUor-$PbJz)Z@5LO&dk&B?g4>*OPp!!xc_~Ka_=CV+w4_?|TgQRj z*csByvz}vxdq{`g1rjVjm;+1O+r#T35f0f*h~lHOm^_`2E+NrkA$yn2>Ty||Gef`c z1BY!&%~4{UH4dxE`vy6=$H;-*P*>{&<$N^(0!oL1HF`r&PP7R`bg)`F0hJUgh-rFN z6u%lpg8c)3Ti>JM4EST5j-1|Y@%&BOR8cnFHa)n?fzjKSKae zbl{|Q4%rZ51uvoutPQjw(bm04*VrQ9@L`x&D0O_Kb{Cbqh$snW4(+ueZM8sN6j$h} z&(V5Sw_!wH4e7z%tOv89yFRR|E;3Xhjrg28SWjI;a4aa|;ObBlhH4l#Wo>$*czLBg zbqlNds34b$k#$vK;qpXIU+9^iK|o&DRGo5$)P=d;!y3?c&eu?ty6W8~0^{=+Re;V0$zij^7i_?EZt^Qtr_*3orCRX2*>i!lX{uUwr z79su?A^sL2{uUwr79su?A^sL2{yz~R2D|+L_sNPswRQZbEA66xoUi$L2I;%^aeD$w zWgkQ`jt*N}rF`L~OWQr+=*!iRkpO>O%__#~|_ly4+cwxuu=73#_7 z*O6LbwX4@YeM;W`eDX8=V{MaFL61W1YM zJ!f|-xrGm}RW&=F#0m%Q-dFw?lAZ@WCd2_Cf!akX()p(7;5~fMYNDOX|F-Jb! zJ^whhgHGa!3KwLO@Z>sLoU}A7!GkEigv4X*s(bj{WBe%IWu&^Al@R!tvS%+ux`oI| z`6jMCFCWn(r72>vYA-|;qp<1Zv1pVrdT4ovXzoVnOUH|mzF`!lhUpMaMr>eiuZ7r3 zCO|WHM3n>s}*q1={)FHCLpZJp{1f69s?$FsACt+37MHwWIg)Zt-EgD*~d!-uJ5d{+%a zdch43YiUpUbi?EkfIq`<(}!G!Gfc|&N_99&Bi+7BW77>*Sd)q5g7XYIT!&vbwm}`O z6*eV9?0j{*M=q?%pyT$qK98KFwIicF&J^2KtVEvgZexHQ_+_h9mI`rn@>0Gk^VIXk z%1d0NOIm?usKsY+0sie|Y?j%io~Too{SsBR#A!6oDqa`k6M0GZ48o{IGQRMwn$076 za#{BYRV&USI(NRb{DU+DCfZJtNqwR;8_s?KR-RmoX<1gH*i?cHBrOr>Y@kg6btJxGQPl)YbZ9xe{0=Y0WI0%NKuX zUwhax-jbV^pygdesyTaIn4}!0wY(!=<|}tj>Y-(oC=GLGr$@=c-lJULqvbtFnl?{; z%-=5K>&L8L?mWEBkw&BJk_nh9+L1lpQ98IfwP z*=7`UMt7pYL`{LYE*0p)+=e;A^B!km&B_w+9a}%?mR`(|EN`b?nagcPl%07#+Lj4= zKt~v5&JbyA4vpQg44f2> zF-b^T*7(2xbB$%rhI)9%>fh8l)hHrjS#(7<3`m08++_lPhJ_4W0N zc&aW#8raqPv*{=|5H-_-Ee=v1VDJkrZ zDp6(lMcfK}q$+s&edUTVJx?7eDOKC8PoDF z!=@`8k)-bUY$`D;Y$;2|GoU|ZIut+ViUZa=;*lT?|at9OkI zw7)bel91@OlaqegZCvknm!1R#?q>rl3 zp6;~}ZOee^R7E=r06lF0=YWXHvNBYfqIgzk@Q#aJ!kRJu?GZONa86X}is)t>2F|}l zTIbYI=O0NtV)oZQ=HV1r^^6ajQZhOLBVX{1&&8)%3w7ba&?GplIg3Q$L$cLc@qwoH z8TyGPvrm>UfQI)hYqD7u=e`!Gcuzn`6L&_BI^3O~u z(y;DEGrZEuJb*L$eo1=k$k#{}6ZU1XO_^NvnHgj+HG(VH672NkyxLdgt+EcF}hivChgmY zf39MIrQF-l6pN>a=_*RwLW5~@Wm>4rk+w0LcXW6$Fdg=(s>dX`Xx({l#$=?@P>~%t zKUanh~`7o)PO z7ayf{`Y!A9=8}Cmup`SToJ4>IxfF;Q$LKL&`LU>7jM8c;cB^1=u8*sg-XZ6fW zR~nMG$C+5dd{;!+!mjS@{jR$uDNrZrof=3LxF&CQZCI^(1B`0MQH8Ehcwf@ikws;v zAU&K6M5+JN+_^wyRb_GfyFSR(DwRYLtTALn(A2CPP4iI>rU(l7AWc9-Ar(L5kusoo z#R9`u3234giio9XS|VU1&WzcELT60PNBBUnmk$yF!D06K_|Cm&pYOiDHET6p*Y#Vj z<^J7$_Bm(Y`@fHU+#hLvHavXZ{qL7gnYNgmF~3fp)tuAVJ9Bo*rDH!W95cq#`?s*t zlp#yLDsHL{J+b)BAjl2(I>to(@X_GSfkpL|E;ZLTPw)zy{m7b>U*~sURQP0J?cp@Y zF?$9F4Ug>V?4N$Wa%6dV$zzke?v0q>S2}&u@fVM*T$wcQFW~dH=7yxH1NUTP&c4$z z%iAwuP~Ps{**E8OTeijV$7MsBLK0!PjvM-O{h5OPs{=p%&c$Wauk&5=5}dL_BYs}G z@||^)rg5&Xw1uRD=lbaXp|bG7`@Zcc|$T@bEUb;9)`eH38`-y9l9)JA)2BSF<|t zcws(po>}Laym7OSS8b2BJ*%hm^YmVE;g=&T8mGn|4tBg7cD)O9#e=c0M?N^&%V*H} zMt|p0VR>Ns`&siBB=_^WelhoGwMYF2Z~`xn9dy1+%fCnX9Eb|Jyd=6NGBCZOYr)2E zmW0GCj*A+1VRr`rT^`f6x^%x^u2a+fH;XP7ZSwZI-z{oMOVK^IN|$rP7iR~IffTd$ z_M)OQxvK(Sd&oR&epK>Xf!#fN?n~T|v~|SH(C%sLll>qq49#sRK3nygkK?q)>fnI) zHUyS$d*}J2J8yWHi#HsNY-o5DIyEb?YH&u}i*IO6c2oF(it&CsDvGbI z{jf3=`qatomV5BkgI?ZKTY7A(SvC$%t6w!-$zCyj$~@o0jnRuvG}nLE`5*EBtSoWp zxxvf(l~EtBdbV3}=GE-`Z+-H8&;*yd-Z#tLul0qNwQyR3)s^;XiT|^}w9nb!S!uE* z%h*?rh|-K+G&eRZJ~}pTz?=l;L-8-Hiqiff3Eaw%pAUQZbNZr(!>5;5j>wELtz9p0 z{qEN6?@l$Owuh2tH26Y~{)-<{*{!D=W_9BHV#}6k`SnLDT&C@-O$>Ej=dmM(9TxGo zVRr6+sK7wi{xOP|X~Yd2w29mQTsKK;?2XKJ3Za2x>n8m^1q za%=QZxhoz(P7tK^H!y8|!tax=G7?PSfZUkc8M4^!y+Z3!|%I)wt z))$ORg1g=pUsQ0G(+-a#m0>1=MHi zt@T`Zj^Etg1w%=+AUU`rf8Kaht$C=-3J+xCg@KOHmCEs=K*fmpJ{}sEMLdKA_-Ewd zX;mODK*GC7o-mCA(n;7c@hA*dWfvu)qK^?(FK;B}HOPwG9>~XFnEEzE#Zn!-m%}=B zEF~{$k(ajqgqStX3QY9HInO$B_Ji0>5W_Oq@qBdcWo5D7A@SyHKEPT0E_8A0pq(2P z1*cGRd~9pQHXbo?>8B<65S|=e+t@%&! z5cZ{fp2H>_JhDbTLxZmQ3Mmg{XEFwQ!RJt|V_yxZm~L(*{^zsLNsu7=7#KEE4`{jJ z!h6VOL=0RyN)KqcDZ&GL`)^r0igYL>F9^Q0iGs^$1t9H?P}?Q6d2-UJpMECKo5|?_ zlM7!oOxBGmCb#$*ChLY5liep7CM(97+#(8q=zak^PBmY$WQDs$d|cTpbxeBApYais z6yr#4zL4_owmv;jBU$gm@IDQSnN~j*NULmGw+k9AO&r!ItJ>#D_(XrcrgkGQ?a$WzFvLzmWwYO@om`g#mQ-j zp&_u7V)|X9=C}dj)2IvT3VY4@cJESad}(t>veuEmgypCH$!qv28#B!#qzf-mGTYY(VNZHL4eQ6GUx85rZh zHA%)=T04A7pNurH&w;Zg?*3ad>A{VQZS29P!IYq-wi6%gHG`^cJb}DJsJ*(}ZYgN# z>Pyb_9ZV+B64E9Lb&|D;Kz>^3GMx1W^Ri}9Z&Y}SS(X*tJtXf)ehUhtJi8S)2Gy#8 zlqL7Z@YMXlBw`S3&_2>83fpH>k8D~FMvS9)S-Q0_sUX;j}h1 znWE$n=?n1YtW!n%f#(kr5c3*ko*fE2C*|z167WrIXEJpMF5N_-c{%k@UxGiHy6?2V zbwc7CX-cu7P|>F}9@zYYN7yzvS+-_0_JwPG-N;UBOiYH+qC-^7K4h7=3KA#i0!%!T zp-a@&iI{jfQ1puHLncF zag2H)Umd0B6|2WF>XaJFh#<$PMY~nJqVELoVJ9n(QALi$pu-0gWrrkpMiVYeLg;F8 zvTWCQFn$|$eAN)}p-HMdWG1>aPAPkmmA0roPUb(@*sJSuvN#r%##mRyK-Csn|42ib zvyZ-UGurHx0qWYEtUyH_z_DD**z2fs_VF=pMXv(c2@QQt77wDz7GHuDAb^(0BdasL z44t_OcHAMu_sLG^X>`^TNV?J{3gryNZ0XU^1f_wYPIv__CTN#h3A#ERYvn42U_hzk z`LzHfG(9{Rt7~<#R1mr0Sq3}4+0wuu$d#1h%pj}!5c3?^aR}EOkUh}qN>b@<5?GVD ziHC#xL78A&v12bv7>1T+Crd}by&cPgAw0FbKbE~vSM6jWnkq0SsvTcxTw&~OG^56I za@HeE>~To;NQJ0)j$oBLPV8SYM7y;+Mnx3L&T%MZ^-dNe5lQ4a^N)I+56dWEpt*hZNG@;lCRMAtixCMQ9?u9&c zDN%+R(($BB3VnX=hOoPiDq}U3JgMJFl{}ua7r&HoT3VjX9QY(q^0iFRQ}d*tBR#j{ zu`;<_rts>U%2QZWO5xErsln(EAu59%Gw*6JWx1-vB=r&}TbcMa%GPgWKMiVnQrM6r z-5l64Z1lG>%%G|#RS4sp3sxhfwY-8?Il#mp2?xoSN5dF-&km;2s&Z3{3l0B0+R7)z24mLaC5jXLy zjIvMRvlAXOJ$@&>Pk?|Y=6xFQ&D@8kJ2N-O``J$(*7m? z7V=b>w08_Bee$I^edLbe3#*hJR)*_F7{iNi+XL4PHHJsru?Magb5`%;JkqHa!|hc3 z7Y4>E06!d?42lBsnz@ORypiJt?zy1qXJ=Hh> literal 0 HcmV?d00001 diff --git a/packages/loop_common/tests/test_structured_grid_2d_geometry.py b/packages/loop_common/tests/test_structured_grid_2d_geometry.py new file mode 100644 index 000000000..28b33f398 --- /dev/null +++ b/packages/loop_common/tests/test_structured_grid_2d_geometry.py @@ -0,0 +1,74 @@ +"""StructuredGrid2DGeometry: the geometry and node-value interpolation the 2D support builds on.""" +import logging + +import numpy as np +import pytest +from loop_common.geometry import StructuredGrid2DGeometry +from loop_common.supports import StructuredGrid2D + +GRIDS = { + "unit": dict(origin=[0.0, 0.0], nsteps=[5, 5], step_vector=[1.0, 1.0]), + "aniso": dict(origin=[10.0, -3.0], nsteps=[7, 4], step_vector=[0.5, 2.0]), +} + + +@pytest.fixture(params=list(GRIDS)) +def geometry(request): + return StructuredGrid2DGeometry(**GRIDS[request.param]) + + +def test_the_support_is_a_geometry(): + assert issubclass(StructuredGrid2D, StructuredGrid2DGeometry) + + +def test_elements_are_flat_node_indexes(geometry): + assert geometry.elements.shape == (geometry.n_elements, 4) + nodes = geometry.nodes + # every element's four corners are the corners of one cell + corners = nodes[geometry.elements] + spans = corners.max(axis=1) - corners.min(axis=1) + np.testing.assert_allclose(spans, np.broadcast_to(geometry.step_vector, spans.shape)) + + +def test_index_round_trip(geometry): + cells = np.arange(geometry.n_elements) + ij = geometry.global_index_to_cell_index(cells) + np.testing.assert_array_equal(geometry.global_cell_indices(ij), cells) + nodes = np.arange(geometry.n_nodes) + ij = geometry.global_index_to_node_index(nodes) + np.testing.assert_array_equal(geometry.global_node_indices(ij), nodes) + + +def test_index_helpers_reject_bad_shapes(geometry): + with pytest.raises(ValueError): + geometry.global_cell_indices(np.array([0, 0])) + with pytest.raises(ValueError): + geometry.global_cell_indices(np.zeros((3, 3), dtype=int)) + + +def test_bilinear_reproduces_a_linear_field(geometry): + coefficients = np.array([1.5, -2.0]) + values = geometry.nodes @ coefficients + 3.0 + low = geometry.origin + high = geometry.origin + (geometry.nsteps - 1) * geometry.step_vector + points = np.random.default_rng(0).uniform(low + 1e-6, high - 1e-6, (200, 2)) + np.testing.assert_allclose( + geometry.evaluate_value(points, values), points @ coefficients + 3.0, atol=1e-9 + ) + + +def test_bilinear_returns_node_values_at_interior_nodes(geometry): + values = np.random.default_rng(1).normal(size=geometry.n_nodes) + nodes = geometry.nodes + interior = ( + (nodes > geometry.origin + 1e-9).all(axis=1) + & (nodes < geometry.origin + (geometry.nsteps - 1) * geometry.step_vector - 1e-9).all(axis=1) + ) + np.testing.assert_allclose( + geometry.evaluate_value(nodes[interior], values)[...], values[interior], atol=1e-12 + ) + + +def test_print_geometry_works_on_a_2d_grid(geometry, caplog): + with caplog.at_level(logging.INFO): + geometry.print_geometry() diff --git a/packages/loop_common/tests/test_structured_grid_2d_reference.py b/packages/loop_common/tests/test_structured_grid_2d_reference.py new file mode 100644 index 000000000..993502c14 --- /dev/null +++ b/packages/loop_common/tests/test_structured_grid_2d_reference.py @@ -0,0 +1,182 @@ +"""Pin the behaviour of the 2D structured interpolation support. + +Records a wide set of ``StructuredGrid2D`` outputs (geometry, node/element arrays, point location +including exact boundary points, bilinear weights, value/gradient evaluation, element lookup, +regularisation operators, index conversions, ``set_nelements``) on a few grids and compares them +with ``data/structured_grid_2d_reference.npz``. + +The reference was generated from the code *before* ``StructuredGrid2D`` was changed to extend +``StructuredGrid2DGeometry``. Regenerate after an intended change, and review the diff: + + python packages/loop_common/tests/test_structured_grid_2d_reference.py --update +""" +import sys +import warnings +from pathlib import Path + +import numpy as np +import pytest +from loop_common.supports import StructuredGrid2D + +REFERENCE = Path(__file__).parent / "data" / "structured_grid_2d_reference.npz" + +# (origin, nodes per axis, step vector). StructuredGrid2D takes node counts. +GRIDS = { + "unit": ([0.0, 0.0], [5, 5], [1.0, 1.0]), + "aniso": ([10.0, -3.0], [7, 4], [0.5, 2.0]), + "minimal": ([0.0, 0.0], [2, 2], [3.0, 1.0]), +} +ATTRIBUTES = ( + "nsteps", + "nsteps_cells", + "maximum", + "origin", + "step_vector", + "n_nodes", + "n_elements", + "element_size", + "n_cell_x", + "n_cell_y", + "dim", + "dimension", + "nodes", + "elements", + "barycentre", +) +OPERATOR_WEIGHTS = {"dxx": 1.0, "dyy": 1.0, "dxy": 1.0} + + +def _flatten(value, prefix, out): + if isinstance(value, dict): + for key in sorted(value): + _flatten(value[key], f"{prefix}.{key}", out) + elif isinstance(value, (tuple, list)) and not all(np.isscalar(v) for v in value): + for i, item in enumerate(value): + _flatten(item, f"{prefix}.{i}", out) + else: + out[prefix] = np.asarray(value) + + +def _record(out, key, call): + try: + _flatten(call(), key, out) + except Exception as exc: # noqa: BLE001 - the exception type is part of the behaviour + out[key] = np.asarray(f"EXC:{type(exc).__name__}") + + +def collect() -> dict[str, np.ndarray]: + out: dict[str, np.ndarray] = {} + rng = np.random.default_rng(5) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + for name, (origin, nodes, step) in GRIDS.items(): + grid = StructuredGrid2D(origin=origin, nsteps=nodes, step_vector=step) + _probe(out, name, grid, origin, nodes, step, rng) + return out + + +def _probe(out, prefix, grid, origin, nodes, step, rng): + for attribute in ATTRIBUTES: + _record(out, f"{prefix}/{attribute}", lambda attribute=attribute: getattr(grid, attribute)) + _record(out, f"{prefix}/str_type", lambda: np.asarray(int(grid.type))) + + low = np.array(origin) - 1.0 + high = np.array(origin) + np.array(nodes) * np.array(step) + 1.0 + anywhere = rng.uniform(low, high, (120, 2)) + # points exactly on the boundary of the nodes' extent and of `maximum` + node_max = np.array(origin) + (np.array(nodes) - 1) * np.array(step) + corners = np.array(np.meshgrid(*zip(origin, node_max), indexing="ij")).reshape(2, -1).T + faces = [] + centre = 0.5 * (np.array(origin) + node_max) + for axis in range(2): + for bound in (origin[axis], node_max[axis]): + point = centre.copy() + point[axis] = bound + faces.append(point) + anywhere = np.vstack([anywhere, corners, np.array(faces)]) + interior = rng.uniform(np.array(origin) + 1e-3, node_max - 1e-3, (80, 2)) + + _record(out, f"{prefix}/inside", lambda: grid.inside(anywhere)) + _record(out, f"{prefix}/position_to_cell_index", lambda: grid.position_to_cell_index(anywhere)) + _record(out, f"{prefix}/check_position_1d", lambda: grid.check_position(np.array([1.0, 2.0]))) + _record(out, f"{prefix}/check_position_3d", lambda: grid.check_position(np.zeros((2, 2, 2)))) + _record(out, f"{prefix}/position_to_cell_corners", lambda: grid.position_to_cell_corners(anywhere)) + _record(out, f"{prefix}/position_to_cell_vertices", lambda: grid.position_to_cell_vertices(interior)) + _record(out, f"{prefix}/position_to_local_coordinates", lambda: grid.position_to_local_coordinates(interior)) + _record(out, f"{prefix}/position_to_dof_coefs", lambda: grid.position_to_dof_coefs(interior)) + _record(out, f"{prefix}/bilinear", lambda: grid.bilinear(rng.uniform(0, 1, (10, 2)))) + + cell_index = np.arange(grid.n_elements) + node_index = np.arange(grid.n_nodes) + _record(out, f"{prefix}/cell_centres", lambda: grid.cell_centres(cell_index)) + _record(out, f"{prefix}/global_index_to_cell_index", lambda: grid.global_index_to_cell_index(cell_index)) + node_ijk = grid.global_index_to_node_index(node_index) + _record(out, f"{prefix}/global_index_to_node_index", lambda: node_ijk) + cell_ijk = grid.global_index_to_cell_index(cell_index) + _record(out, f"{prefix}/cell_corner_indexes", lambda: grid.cell_corner_indexes(cell_ijk)) + _record(out, f"{prefix}/global_node_indices", lambda: grid.global_node_indices(node_ijk)) + _record(out, f"{prefix}/global_cell_indices", lambda: grid.global_cell_indices(cell_ijk)) + _record(out, f"{prefix}/node_indexes_to_position", lambda: grid.node_indexes_to_position(node_ijk[:10])) + _record(out, f"{prefix}/neighbour_global_indexes", lambda: grid.neighbour_global_indexes()) + + values = rng.normal(size=grid.n_nodes) + _record(out, f"{prefix}/evaluate_value", lambda: grid.evaluate_value(anywhere, values)) + _record(out, f"{prefix}/evaluate_gradient", lambda: grid.evaluate_gradient(interior, values)) + _record(out, f"{prefix}/get_element_for_location", lambda: grid.get_element_for_location(interior)) + _record( + out, + f"{prefix}/get_element_gradient_for_location", + lambda: grid.get_element_gradient_for_location(interior), + ) + _record( + out, + f"{prefix}/get_operators", + lambda: dict(grid.get_operators(OPERATOR_WEIGHTS)), + ) + _record(out, f"{prefix}/vtk_n_cells", lambda: np.array([grid.vtk(0.0).n_cells])) + _record(out, f"{prefix}/set_nelements", lambda: grid.set_nelements(100)) + + +def _assert_same(key, got, want): + if want.dtype.kind in "US": + assert got.dtype.kind in "US" and str(got) == str(want), f"{key}: {got} != {want}" + return + assert got.dtype.kind not in "US", f"{key}: now raises {got}" + assert got.shape == want.shape, f"{key}: shape {got.shape} != {want.shape}" + assert got.dtype.kind == want.dtype.kind, f"{key}: dtype {got.dtype} != {want.dtype}" + if want.dtype.kind == "f": + np.testing.assert_allclose(got, want, rtol=1e-9, atol=1e-12, err_msg=key, equal_nan=True) + else: + np.testing.assert_array_equal(got, want, err_msg=key) + + +@pytest.fixture(scope="module") +def current(): + return collect() + + +@pytest.fixture(scope="module") +def reference(): + with np.load(REFERENCE) as data: + return {key: data[key] for key in data.files} + + +def test_same_probes_are_recorded(current, reference): + assert sorted(current) == sorted(reference) + + +@pytest.mark.parametrize("grid_name", list(GRIDS)) +def test_structured_grid_2d_behaviour_is_unchanged(grid_name, current, reference): + keys = [k for k in reference if k.startswith(f"{grid_name}/")] + assert keys + for key in keys: + _assert_same(key, current[key], reference[key]) + + +if __name__ == "__main__": + if sys.argv[1:] != ["--update"]: + raise SystemExit(__doc__) + REFERENCE.parent.mkdir(exist_ok=True) + result = collect() + np.savez_compressed(REFERENCE, **result) + print(f"wrote {len(result)} arrays to {REFERENCE}") From b69e3e6cc7cdb95e039b41a8d814b27c8bb4eac2 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sun, 1 Nov 2026 18:40:09 +1030 Subject: [PATCH 2/2] ci: run qgis-plugin-compat on Python 3.12, since plugin requires pyvista>=0.49.2 (needs 3.10+) --- .github/workflows/qgis-compat.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/qgis-compat.yml b/.github/workflows/qgis-compat.yml index 3a5078c17..287f7a92b 100644 --- a/.github/workflows/qgis-compat.yml +++ b/.github/workflows/qgis-compat.yml @@ -36,7 +36,7 @@ jobs: - name: Set up uv and Python uses: astral-sh/setup-uv@v5 with: - python-version: "3.9" + python-version: "3.12" enable-cache: true - name: Install plugin's non-QGIS test requirements