write testing routine and validate that _calc_corresponding_dof_indices works properly

LDDsimulation/domainPatch.py

@@ -909,10 +909,6 @@ class DomainPatch(df.SubDomain):
                                 self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
                                     {"flux_x": flux_x_dof_ind}
                                 )
-                            else:
-                                print(f"in flux_x")
-                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({dof_coord}) and ({gli_dof_coord})",
-                                                    "Check on the mesh what is going on.")
 
                         # defensive sanity check if the dof has been set.
                         if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_x"] is None:
@@ -927,10 +923,7 @@ class DomainPatch(df.SubDomain):
                                 self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
                                     {"flux_y": flux_y_dof_ind}
                                 )
-                            else:
-                                print(f"in flux_y")
-                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({dof_coord}) and ({gli_dof_coord})",
-                                                    "Check on the mesh what is going on.")
+
                         # sanity check for flux_y
                         if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_y"] is None:
                             raise RuntimeError(f"didn't find flux_y_dof_ind corresponding to dict key ({gli_dof_coord})",
@@ -944,15 +937,63 @@ class DomainPatch(df.SubDomain):
                                 self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
                                     {"pressure": p_dof_ind}
                                 )
-                            else:
-                                print(f"in pressure")
-                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({dof_coord}) and ({gli_dof_coord})",
-                                                    "Check on the mesh what is going on.")
 
                         if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["pressure"] is None:
                             raise RuntimeError(f"didn't find p_dof_ind corresponding to dict key ({gli_dof_coord})",
                                                 "something is wrong")
 
+        if debug:
+            # routine to check if dofs were correctly matched. Here, we compare
+            # the coordinates corresponding to the different indices. If they
+            # match, we correctly identified the dof indices.
+            for interf_ind in self.has_interface:
+                interface = self.interface[interf_ind]
+                for phase in self.has_phases:
+                    gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+                    p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+                    flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+                    for local_facet_ind in interface.outer_normals[subdomain].keys():
+                        gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
+                        p_dofs_along_facet = p_dof_coordinates[local_facet_ind]
+                        flux_x_dofs_along_facet = flux_dof_coordinates["x"][local_facet_ind]
+                        flux_y_dofs_along_facet = flux_dof_coordinates["y"][local_facet_ind]
+                        for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
+                            corres_flux_x_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_x"]
+                            corres_flux_y_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_y"]
+                            corres_pressure_dof_ind = self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["pressure"]
+                            corres_flux_x_dof_coord = flux_x_dofs_along_facet[corres_flux_x_dof_ind]
+                            corres_flux_y_dof_coord = flux_y_dofs_along_facet[corres_flux_y_dof_ind]
+                            corres_pressure_dof_coord = p_dofs_along_facet[corres_pressure_dof_ind]
+
+                            flux_x_a_matches_b = np.allclose(corres_flux_x_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
+                            flux_x_b_matches_a = np.allclose(gli_dof_coord, corres_flux_x_dof_coord, rtol=1e-10, atol=1e-14)
+                            if flux_x_a_matches_b and flux_x_b_matches_a:
+                                print(f"flux_x: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_flux_x_dof_ind}, {corres_flux_x_dof_coord})")
+                            else:
+                                print(f"in flux_x")
+                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_flux_x_dof_coord}) and ({gli_dof_coord})",
+                                                    "dofs matching might be faulty.")
+
+                            flux_y_a_matches_b = np.allclose(corres_flux_y_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
+                            flux_y_b_matches_a = np.allclose(gli_dof_coord, corres_flux_y_dof_coord, rtol=1e-10, atol=1e-14)
+                            if flux_y_a_matches_b and flux_y_b_matches_a:
+                                print(f"flux_y: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_flux_y_dof_ind}, {corres_flux_y_dof_coord})")
+                            else:
+                                print(f"in flux_y")
+                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_flux_y_dof_coord}) and ({gli_dof_coord})",
+                                                    "dofs matching might be faulty.")
+
+                            pressure_a_matches_b = np.allclose(corres_pressure_dof_coord ,gli_dof_coord, rtol=1e-10, atol=1e-14)
+                            pressure_b_matches_a = np.allclose(gli_dof_coord, corres_pressure_dof_coord, rtol=1e-10, atol=1e-14)
+                            if pressure_a_matches_b and pressure_b_matches_a:
+                                print(f"pressure: gli dof ({gli_dof_index}, {gli_dof_coord}) correctly matched to flux_x dof ({corres_pressure_dof_ind}, {corres_pressure_dof_coord})")
+                            else:
+                                print(f"in pressure")
+                                raise RuntimeError(f"np.allclose(a,b) != np.allclose(b,a) while comparing ({corres_pressure_dof_coord}) and ({gli_dof_coord})",
+                                                    "dofs matching might be faulty.")
+
+
+
 
     def _init_interface_values(self):
         """ allocate the dictionaries for each interface used to store the