diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py index 4cfcc6f578f5e5527657d8a15bfebc2ec06f565e..fadc059ed801675a744c85fa5f8ec307be17c377 100644 --- a/LDDsimulation/domainPatch.py +++ b/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")