From d010e53d14f7bc838a8b27988892dedb75975cfb Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 29 Aug 2019 10:50:35 +0200
Subject: [PATCH] write testing routine and validate that
_calc_corresponding_dof_indices works properly
---
LDDsimulation/domainPatch.py | 65 +++++++++++++++++++++++++++++-------
1 file changed, 53 insertions(+), 12 deletions(-)
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index fadc059..1de1635 100644
--- a/LDDsimulation/domainPatch.py
+++ b/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
--
GitLab