diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py index ec521c31c..273c77e88 100644 --- a/LoopStructural/datatypes/_bounding_box.py +++ b/LoopStructural/datatypes/_bounding_box.py @@ -114,7 +114,6 @@ def fit(self, locations: np.ndarray): raise LoopValueError( f"locations array is {locations.shape[1]}D but bounding box is {self.dimensions}" ) - print("fitting") self.origin = locations.min(axis=0) self.maximum = locations.max(axis=0) return self diff --git a/LoopStructural/interpolators/__init__.py b/LoopStructural/interpolators/__init__.py index 625c4d43a..edbaeb94c 100644 --- a/LoopStructural/interpolators/__init__.py +++ b/LoopStructural/interpolators/__init__.py @@ -63,8 +63,8 @@ class InterpolatorType(IntEnum): from ..interpolators._finite_difference_interpolator import ( FiniteDifferenceInterpolator, ) -from ..interpolators.piecewiselinear_interpolator import ( - PiecewiseLinearInterpolator, +from ..interpolators._p1interpolator import ( + P1Interpolator as PiecewiseLinearInterpolator, ) from ..interpolators._discrete_fold_interpolator import ( DiscreteFoldInterpolator, diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py index 55fc3d8ec..4610731ab 100644 --- a/LoopStructural/interpolators/supports/_3d_base_structured.py +++ b/LoopStructural/interpolators/supports/_3d_base_structured.py @@ -275,34 +275,8 @@ def _global_indicies(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarra + nsteps[None, 0] * nsteps[None, 1] * indexes[:, 2] ) return gi.reshape(original_shape[:-1]) - # if len(indexes.shape) == 2: - # if indexes.shape[1] != 3 and indexes.shape[0] == 3: - # indexes = indexes.swapaxes(0, 1) - # if indexes.shape[1] != 3: - # logger.error("Indexes shape {}".format(indexes.shape)) - # raise ValueError("Cell indexes needs to be Nx3") - # return ( - # indexes[:, 0] - # + nsteps[None, 0] * indexes[:, 1] - # + nsteps[None, 0] * nsteps[None, 1] * indexes[:, 2] - # ) - # if len(indexes.shape) == 3: - # # if indexes.shape[2] != 3 and indexes.shape[1] == 3: - # # indexes = indexes.swapaxes(1, 2) - # # if indexes.shape[2] != 3 and indexes.shape[0] == 3: - # # indexes = indexes.swapaxes(0, 2) - # if indexes.shape[2] != 3: - # logger.error("Indexes shape {}".format(indexes.shape)) - # raise ValueError("Cell indexes needs to be NxNx3") - # return ( - # indexes[:, :, 0] - # + nsteps[None, None, 0] * indexes[:, :, 1] - # + nsteps[None, None, 0] * nsteps[None, None, 1] * indexes[:, :, 2] - # ) - # else: - # raise ValueError("Cell indexes need to be a 2 or 3d numpy array") - - def cell_corner_indexes(self, cell_indexes): + + 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 @@ -455,5 +429,4 @@ def global_cell_indices(self, indexes) -> np.ndarray: return self._global_indicies(indexes, self.nsteps_cells) def pyvista(self): - pass diff --git a/LoopStructural/interpolators/supports/_3d_structured_tetra.py b/LoopStructural/interpolators/supports/_3d_structured_tetra.py index 5cc9cf5fa..a6ce11f1c 100644 --- a/LoopStructural/interpolators/supports/_3d_structured_tetra.py +++ b/LoopStructural/interpolators/supports/_3d_structured_tetra.py @@ -126,18 +126,21 @@ def _init_face_table(self): el_rel[np.arange(n1.shape[0]), 1] = n2 el_rel = el_rel[el_rel[:, 0] >= 0, :] - el_rel2 = np.zeros((self.neighbours.flatten().shape[0], 2), dtype=int) - el_rel2[:] = -1 + # el_rel2 = np.zeros((self.neighbours.flatten().shape[0], 2), dtype=int) + self.shared_element_relationships[:] = -1 el_pairs = coo_matrix( (np.ones(el_rel.shape[0]), (el_rel[:, 0], el_rel[:, 1])) ).tocsr() i, j = tril(el_pairs).nonzero() - el_rel2[: len(i), 0] = i - el_rel2[: len(i), 1] = j - el_rel2 = el_rel2[el_rel2[:, 0] >= 0, :] + self.shared_element_relationships[: len(i), 0] = i + self.shared_element_relationships[: len(i), 1] = j - faces = element_nodes[el_rel2[:, 0], :].multiply( - element_nodes[el_rel2[:, 1], :] + self.shared_element_relationships = self.shared_element_relationships[ + self.shared_element_relationships[:, 0] >= 0, : + ] + + faces = element_nodes[self.shared_element_relationships[:, 0], :].multiply( + element_nodes[self.shared_element_relationships[:, 1], :] ) shared_faces = faces[np.array(np.sum(faces, axis=1) == 3).flatten(), :] row, col = shared_faces.nonzero() @@ -146,9 +149,14 @@ def _init_face_table(self): shared_face_index = np.zeros((shared_faces.shape[0], 3), dtype=int) shared_face_index[:] = -1 shared_face_index[row.reshape(-1, 3)[:, 0], :] = col.reshape(-1, 3) - self.shared_element_relationships = el_rel2 # [np.arange(n1.shape[0]), 0] = n1 - # self.shared_element_relationships[np.arange(n1.shape[0]), 1] = n2 - # self.shared_elements[np.arange(el_rel2.shape[0]), :] = shared_face_index + + self.shared_elements[ + np.arange(self.shared_element_relationships.shape[0]), : + ] = shared_face_index + # resize + self.shared_elements = self.shared_elements[ + : len(self.shared_element_relationships), : + ] @property def shared_element_norm(self): @@ -187,7 +195,7 @@ def evaluate_value(self, pos: np.ndarray, property_array: np.ndarray) -> np.ndar values[:] = np.nan vertices, c, tetras, inside = self.get_element_for_location(pos) values[inside] = np.sum( - c[inside, :] * property_array[tetras[inside, :]], axis=1 + c[inside, :] * property_array[self.elements[tetras[inside]]], axis=1 ) return values @@ -219,7 +227,8 @@ def evaluate_gradient( ) = self.get_element_gradient_for_location(pos) # grads = np.zeros(tetras.shape) values[inside, :] = ( - element_gradients[inside, :, :] * property_array[tetras[inside, None, :]] + element_gradients[inside, :, :] + * property_array[self.elements[tetras[inside]][:, None, :]] ).sum(2) length = np.sum(values[inside, :], axis=1) # values[inside,:] /= length[:,None] @@ -320,12 +329,12 @@ def get_element_for_location(self, pos: np.ndarray): c_return[inside] = c[mask] tetra_return = np.zeros((pos.shape[0])).astype(int) tetra_return[:] = -1 - print(inside.shape, gi.shape, c.shape, mask.shape) local_tetra_index = np.tile(np.arange(0, 5)[None, :], (mask.shape[0], 1)) local_tetra_index = local_tetra_index[mask] - tetra_global_index = self.tetra_global_index(cell_indexes, local_tetra_index) - print(local_tetra_index) - tetra_return[inside] = tetra_global_index[inside] + tetra_global_index = self.tetra_global_index( + cell_indexes[inside, :], local_tetra_index + ) + tetra_return[inside] = tetra_global_index return vertices_return, c_return, tetra_return, inside def evaluate_shape(self, locations): @@ -350,8 +359,9 @@ def get_elements(self): x = np.arange(0, self.nsteps_cells[0]) y = np.arange(0, self.nsteps_cells[1]) z = np.arange(0, self.nsteps_cells[2]) - - cell_indexes = np.array(np.meshgrid(x, y, z, indexing="ij")).reshape(3, -1).T + ## reverse x and z so that indexing is + zz, yy, xx = np.meshgrid(z, y, x, indexing="ij") + cell_indexes = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T # get cell corners cell_corners = self.cell_corner_indexes(cell_indexes) even_mask = np.sum(cell_indexes, axis=1) % 2 == 0 @@ -375,7 +385,6 @@ def tetra_global_index(self, indices, tetra_index): ------- """ - print("idc", indices.shape) return ( tetra_index + indices[:, 0] * 5 @@ -462,7 +471,7 @@ def evaluate_shape_derivatives(self, pos, elements=None): if elements is None: verts, c, elements, inside = self.get_element_for_location(pos) # np.arange(0, self.n_elements, dtype=int) - print(pos.shape, inside.shape, elements.shape) + return ( self.get_element_gradients(elements), elements, @@ -657,57 +666,57 @@ def get_neighbours(self) -> np.ndarray: odd_mask = ( np.sum(self.global_index_to_cell_index(tetra_index // 5), axis=0) % 2 == 1 ) - odd_mask = ~odd_mask.astype(bool) + odd_mask = odd_mask.astype(bool) # apply masks to masks = [] masks.append( [ np.logical_and(one_mask, odd_mask), - np.array([[-1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3]]), + np.array([[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 4]]), ] ) masks.append( [ np.logical_and(two_mask, odd_mask), - np.array([[1, 0, 0, 2], [0, -1, 0, 1], [0, 0, 1, 4]]), + np.array([[-1, 0, 0, 2], [0, -1, 0, 1], [0, 0, 1, 3]]), ] ) masks.append( [ np.logical_and(three_mask, odd_mask), - np.array([[-1, 0, 0, 4], [0, -1, 0, 3], [0, 0, -1, 2]]), + np.array([[-1, 0, 0, 4], [0, 1, 0, 3], [0, 0, -1, 1]]), ] ) masks.append( [ np.logical_and(four_mask, odd_mask), - np.array([[1, 0, 0, 3], [0, 1, 0, 4], [0, 0, -1, 1]]), + np.array([[1, 0, 0, 3], [0, -1, 0, 4], [0, 0, -1, 2]]), ] ) masks.append( [ np.logical_and(one_mask, ~odd_mask), - np.array([[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 4]]), + np.array([[-1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3]]), ] ) masks.append( [ np.logical_and(two_mask, ~odd_mask), - np.array([[-1, 0, 0, 2], [0, -1, 0, 1], [0, 0, 1, 3]]), + np.array([[1, 0, 0, 2], [0, -1, 0, 1], [0, 0, 1, 4]]), ] ) masks.append( [ np.logical_and(three_mask, ~odd_mask), - np.array([[-1, 0, 0, 4], [0, 1, 0, 3], [0, 0, -1, 1]]), + np.array([[-1, 0, 0, 4], [0, -1, 0, 3], [0, 0, -1, 2]]), ] ) masks.append( [ np.logical_and(four_mask, ~odd_mask), - np.array([[1, 0, 0, 3], [0, -1, 0, 4], [0, 0, -1, 2]]), + np.array([[1, 0, 0, 3], [0, 1, 0, 4], [0, 0, -1, 1]]), ] ) @@ -730,12 +739,10 @@ def get_neighbours(self) -> np.ndarray: global_neighbour_idx = np.zeros((c_xi.shape[0], 4)).astype(int) global_neighbour_idx[:] = -1 global_neighbour_idx = ( - mask[:, 3] - + neigh_cell[:, :, 0] * 5 - + neigh_cell[:, :, 1] * self.nsteps_cells[0] * 5 - + neigh_cell[:, :, 2] * self.nsteps_cells[0] * self.nsteps_cells[1] * 5 - ) - + neigh_cell[:, :, 0] + + neigh_cell[:, :, 1] * self.nsteps_cells[0] + + neigh_cell[:, :, 2] * self.nsteps_cells[0] * self.nsteps_cells[1] + ) * 5 + mask[:, 3] global_neighbour_idx[~inside] = -1 neighbours[logic, 1:] = global_neighbour_idx