From d076a8a3dd7b242055f83806a4ce3083ba074a62 Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Thu, 6 Apr 2023 09:19:23 -0500 Subject: [PATCH 1/5] Make it possible to query FEM mappings from Python --- yt/utilities/lib/element_mappings.pyx | 28 +++++++++++++++++++++++++++ yt/utilities/mesh_code_generation.py | 6 +++--- 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/yt/utilities/lib/element_mappings.pyx b/yt/utilities/lib/element_mappings.pyx index c94654c7071..bd5c54c9825 100644 --- a/yt/utilities/lib/element_mappings.pyx +++ b/yt/utilities/lib/element_mappings.pyx @@ -115,6 +115,34 @@ cdef class ElementSampler: val = self.sample_at_unit_point(mapped_coord, field_values) return val + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def map_reals_to_unit(self, + np.float64_t[:,::1] vertices, + np.float64_t[:,::1] positions): + cdef double mapped_x[3] + cdef int i, n + # We have N vertices, which each have three components. + cdef np.ndarray[np.float64_t, ndim=2] output_coords + output_coords = np.zeros((positions.shape[0], positions.shape[1]), dtype="float64") + # Now for each position, we map + for n in range(positions.shape[0]): + self.map_real_to_unit(mapped_x, &vertices[0,0], &positions[n, 0]) + for i in range(3): + output_coords[n, i] = mapped_x[i] + return output_coords + + def sample_at_real_points(self, + np.float64_t[:,::1] vertices, + np.float64_t[::1] field_values, + np.float64_t[:,::1] positions): + cdef np.ndarray[np.float64_t, ndim=1] output_values + output_values = np.zeros(positions.shape[0], dtype="float64") + for n in range(positions.shape[0]): + output_values[n] = self.sample_at_real_point( + &vertices[0,0], &field_values[0], &positions[n,0]) + return output_values cdef class P1Sampler1D(ElementSampler): ''' diff --git a/yt/utilities/mesh_code_generation.py b/yt/utilities/mesh_code_generation.py index 9f34523298a..7e9fc11a46e 100644 --- a/yt/utilities/mesh_code_generation.py +++ b/yt/utilities/mesh_code_generation.py @@ -83,11 +83,11 @@ def _compute_jacobian(self): assert self.num_vertices == len(self.N) assert self.num_dim == self.num_mapped_coords - X = MatrixSymbol("vertices", self.num_vertices, self.num_dim) + self.X = MatrixSymbol("vertices", self.num_vertices, self.num_dim) self.fx = MatrixSymbol("fx", self.num_dim, 1) - physical_position = MatrixSymbol("phys_x", self.num_dim, 1) + self.physical_position = MatrixSymbol("phys_x", self.num_dim, 1) - self.f = (self.N.T * Matrix(X)).T - physical_position + self.f = (self.N.T * Matrix(self.X)).T - self.physical_position self.J = symarray("J", (self.num_dim, self.num_dim)) for i in range(self.num_dim): From 804bbd153bb316b7d49a58d6e1846108978afd4c Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Thu, 7 Sep 2023 10:04:52 -0500 Subject: [PATCH 2/5] Add check_contains and change position shape loop --- yt/utilities/lib/element_mappings.pyx | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/yt/utilities/lib/element_mappings.pyx b/yt/utilities/lib/element_mappings.pyx index bd5c54c9825..8e5ce4ffc7e 100644 --- a/yt/utilities/lib/element_mappings.pyx +++ b/yt/utilities/lib/element_mappings.pyx @@ -95,6 +95,24 @@ cdef class ElementSampler: cdef int check_inside(self, double* mapped_coord) nogil: pass + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def check_contains(self, np.float64_t[:,::1] vertices, + np.float64_t[:,::1] positions): + cdef np.ndarray[np.float64_t, ndim=2] mapped_coords + cdef np.ndarray[np.uint8_t, ndim=1] mask + mapped_coords = self.map_reals_to_unit(vertices, positions) + mask = np.zeros(mapped_coords.shape[0]) + cdef double[3] mapped_coord + cdef int i, j + for i in range(mapped_coords.shape[0]): + for j in range(mapped_coords.shape[1]): + mapped_coord[j] = mapped_coords[i, j] + mask[i] = self.check_inside(mapped_coord) + return mask + + @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) @@ -129,7 +147,7 @@ cdef class ElementSampler: # Now for each position, we map for n in range(positions.shape[0]): self.map_real_to_unit(mapped_x, &vertices[0,0], &positions[n, 0]) - for i in range(3): + for i in positions.shape[1]: output_coords[n, i] = mapped_x[i] return output_coords From 175cefa16d47d2edc467ed0a042b569764ff61e2 Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Tue, 12 Sep 2023 11:43:20 -0500 Subject: [PATCH 3/5] Update yt/utilities/lib/element_mappings.pyx Co-authored-by: Chris Havlin --- yt/utilities/lib/element_mappings.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/utilities/lib/element_mappings.pyx b/yt/utilities/lib/element_mappings.pyx index 8e5ce4ffc7e..ac55f448df3 100644 --- a/yt/utilities/lib/element_mappings.pyx +++ b/yt/utilities/lib/element_mappings.pyx @@ -147,7 +147,7 @@ cdef class ElementSampler: # Now for each position, we map for n in range(positions.shape[0]): self.map_real_to_unit(mapped_x, &vertices[0,0], &positions[n, 0]) - for i in positions.shape[1]: + for i in range(positions.shape[1]): output_coords[n, i] = mapped_x[i] return output_coords From 4c0f874b273ab8b0a5c9364987bb5baeb9a1a24c Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Tue, 12 Sep 2023 11:43:26 -0500 Subject: [PATCH 4/5] Update yt/utilities/lib/element_mappings.pyx Co-authored-by: Chris Havlin --- yt/utilities/lib/element_mappings.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/utilities/lib/element_mappings.pyx b/yt/utilities/lib/element_mappings.pyx index ac55f448df3..742a80f4372 100644 --- a/yt/utilities/lib/element_mappings.pyx +++ b/yt/utilities/lib/element_mappings.pyx @@ -103,7 +103,7 @@ cdef class ElementSampler: cdef np.ndarray[np.float64_t, ndim=2] mapped_coords cdef np.ndarray[np.uint8_t, ndim=1] mask mapped_coords = self.map_reals_to_unit(vertices, positions) - mask = np.zeros(mapped_coords.shape[0]) + mask = np.zeros(mapped_coords.shape[0], dtype=np.uint8) cdef double[3] mapped_coord cdef int i, j for i in range(mapped_coords.shape[0]): From 59dc296db2d1716f67cb6462715e21b031d6ac13 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 12 Sep 2023 20:54:51 +0000 Subject: [PATCH 5/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/utilities/lib/element_mappings.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/utilities/lib/element_mappings.pyx b/yt/utilities/lib/element_mappings.pyx index f228e6b677e..2a1fb70b113 100644 --- a/yt/utilities/lib/element_mappings.pyx +++ b/yt/utilities/lib/element_mappings.pyx @@ -103,7 +103,7 @@ cdef class ElementSampler: cdef np.ndarray[np.float64_t, ndim=2] mapped_coords cdef np.ndarray[np.uint8_t, ndim=1] mask mapped_coords = self.map_reals_to_unit(vertices, positions) - mask = np.zeros(mapped_coords.shape[0], dtype=np.uint8) + mask = np.zeros(mapped_coords.shape[0], dtype=np.uint8) cdef double[3] mapped_coord cdef int i, j for i in range(mapped_coords.shape[0]):