From 6b0d09a712d63b2b307f0d2f880ed10cc20d3bf5 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 4 Apr 2019 13:07:15 +0200
Subject: [PATCH] fix implementation of
 DomainPatch._calc_interface_has_common_dof_indices

---
 LDDsimulation/domainPatch.py | 35 ++++++++++++++++++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index bb99744..8d66b56 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -715,8 +715,41 @@ class DomainPatch(df.SubDomain):
     def _calc_interface_has_common_dof_indices(self):
         """ populate dictionary self._interface_has_common_dof_indices
 
-        Determin for each interface in self.has_interface the dof indices
+        Determin for each interface in self.has_interface the dof indices that
+        also belong to another interface and store them in the dictionary
+        self._interface_has_common_dof_indices
+        This method populates the dictionary
+            self._interface_has_common_dof_indices
         """
+        for interf_ind, interf_dofs in self._dof_indices_of_interface.items():
+            # initialise bool numpy arry for the search.
+            if self.isRichards:
+                is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
+            else:
+                is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
+                is_part_of_neighbour_interface_nw = np.zeros(interf_dofs['nonwetting'].size, dtype=bool)
+
+            for other_interf_ind, other_interf_dofs in self._dof_indices_of_interface.items():
+                if other_interf_ind is not interf_ind:
+                    if self.isRichards:
+                        is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
+                    else:
+                        is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
+                        is_part_of_neighbour_interface_nw += np.isin(interf_dofs['nonwetting'], other_interf_dofs['nonwetting'], assume_unique=True)
+            # this carries a numpyarray of the indices which belong to other
+            # interfaces, too
+            self._interface_has_common_dof_indices.update({interf_ind: dict()})
+            if self.isRichards:
+                self._interface_has_common_dof_indices[interf_ind].update(
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
+                )
+            else:
+                self._interface_has_common_dof_indices[interf_ind].update(
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]},
+                    {'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
+                )
+
+
 
 
 # END OF CLASS DomainPatch
-- 
GitLab