fix small bug subm

firedrakeproject · Apr 4, 2024 · 3d1b6bd · 3d1b6bd
1 parent 4a0378f
commit 3d1b6bd
Showing 1 changed file with 10 additions and 5 deletions.
diff --git 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: