From 3d1b6bd2c55a6bf52b5da6b519cc7a4415013869 Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Sun, 31 Mar 2024 11:54:48 +0100 Subject: [PATCH] fix small bug subm --- firedrake/mesh.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 0cd0ea4f3e..d32bd58d1c 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -1505,26 +1505,31 @@ def submesh_map_child_parent(self, source_integral_type, source_subset_points, c else: _, target_indices_int, source_indices_int = np.intersect1d(subpoints[target.interior_facets.facets], source_subset_points, return_indices=True) _, target_indices_ext, source_indices_ext = np.intersect1d(subpoints[target.exterior_facets.facets], source_subset_points, return_indices=True) + #print(self.comm.rank, "source_subset_points", source_subset_points, flush=True) + #print(self.comm.rank, "subpoints[target.exterior_facets.facets]", subpoints[target.exterior_facets.facets], flush=True) n_int = len(source_indices_int) n_ext = len(source_indices_ext) n_int_max = self._comm.allreduce(n_int, op=MPI.MAX) n_ext_max = self._comm.allreduce(n_ext, op=MPI.MAX) if n_int_max > 0: assert n_ext_max == 0 - assert n_int == len(source_subset_points) + assert n_int <= len(source_subset_points) target_integral_type = "interior_facet" elif n_ext_max > 0: assert n_int_max == 0 - assert n_ext == len(source_subset_points) + assert n_ext <= len(source_subset_points) target_integral_type = "exterior_facet" else: raise RuntimeError("Can not find a map from source to target.") if not child_parent_map: - target_subset_points = np.empty_like(source_subset_points) if target_integral_type == "interior_facet": - target_subset_points[source_indices_int] = target.interior_facets.facets[target_indices_int] + #target_subset_points = np.empty_like(source_subset_points, shape=(n_int, )) + #target_subset_points[source_indices_int] = target.interior_facets.facets[target_indices_int] + target_subset_points = target.interior_facets.facets[target_indices_int] elif target_integral_type == "exterior_facet": - target_subset_points[source_indices_ext] = target.exterior_facets.facets[target_indices_ext] + #target_subset_points = np.empty_like(source_subset_points, shape=(n_ext, )) + #target_subset_points[source_indices_ext] = target.exterior_facets.facets[target_indices_ext] + target_subset_points = target.exterior_facets.facets[target_indices_ext] else: raise NotImplementedError(f"Not implemented for (source_dim, target_dim, source_integral_type) == ({source_dim}, {target_dim}, {source_integral_type})") if child_parent_map: