fix bug in _calc_corresponding_dof_indices(self, debug=True):

LDDsimulation/domainPatch.py

@@ -851,7 +851,7 @@ class DomainPatch(df.SubDomain):
                         )
 
 
-    def _calc_corresponding_dof_indices(self):
+    def _calc_corresponding_dof_indices(self, debug=True):
         """ calculate dictionary which for each interface and each phase holds
         for each facet index, the dof indices of the pressures and flux components
         corresponding to a given dof index of the gli function.
@@ -868,11 +868,11 @@ class DomainPatch(df.SubDomain):
                 self.interface_corresponding_dof_index[interf_ind].update(
                     {phase: dict()}
                     )
-                # get dictionaries of global dof indices and coordinates along
-                # the interface. These dictionaries have the facet indices of
+                # get dictionaries of dof indices and coordinates along
+                # the interface. These dictionaries have the indices of
                 # facets belonging to the interface as keys (with respect to
                 # the submesh numbering) and the dictionary
-                # containing pairs {global_dof_index: dof_coordinate} for all
+                # containing pairs {dof_index: dof_coordinate} for all
                 # dofs along the facet as values.
                 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]
@@ -882,7 +882,6 @@ class DomainPatch(df.SubDomain):
                 # with respect to the numbering of the subdomain submesh.
                 # Global facet index refers to the facet numbering
                 # of the mesh of the whole computational domain.
-                # local2global_facet = interface.local_to_global_facet_map[subdomain]
                 for local_facet_ind in interface.outer_normals[subdomain].keys():
                     self.interface_corresponding_dof_index[interf_ind][phase].update(
                         {local_facet_ind: dict()}
@@ -899,9 +898,10 @@ class DomainPatch(df.SubDomain):
                     for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
                         self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind].update(
                             {gli_dof_index: {"flux_x": None,
-                                             "flux_x": None,
+                                             "flux_y": None,
                                              "pressure": None}}
                         )
+                        # flux_x dofs
                         for flux_x_dof_ind, dof_coord in flux_x_dofs_along_facet.items():
                             a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
                             b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
@@ -909,10 +909,17 @@ 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:
                             raise RuntimeError(f"didn't find flux_x_dof_ind corresponding to dict key ({gli_dof_coord})",
                                                 "something is wrong")
 
+                        # flux_y dofs
                         for flux_y_dof_ind, dof_coord in flux_y_dofs_along_facet.items():
                             a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
                             b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
@@ -920,10 +927,16 @@ 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})",
                                                 "something is wrong")
 
+                        # pressure dofs
                         for p_dof_ind, dof_coord in p_dofs_along_facet.items():
                             a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
                             b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
@@ -931,6 +944,11 @@ 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")