From 87cc2cd4a5db93cb63cd229a4dc5802e8ddd0cec Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 5 Apr 2019 20:07:05 +0200
Subject: [PATCH] fix further bugs

---
 LDDsimulation/LDDsimulation.py          |  5 +-
 LDDsimulation/domainPatch.py            | 94 ++++++++-----------------
 RR-2-patch-test-case/RR-2-patch-test.py |  2 +-
 3 files changed, 34 insertions(+), 67 deletions(-)

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 0fa2541..ba5dc8f 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -491,8 +491,9 @@ class LDDsimulation(object):
                         relative_permeability = self.relative_permeability[subdom_num],#
                         saturation = self.saturation[subdom_num],#
                         timestep_size = self.timestep_size[subdom_num],#
-                        tol = self.tol
-                    )})
+                        tol = self.tol)
+                    }
+                )
 
 
 
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index f318f15..80774e5 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -274,9 +274,11 @@ class DomainPatch(df.SubDomain):
         if phase == "wetting":
             # gather interface forms
             interface_forms = []
+            # the marker values of the interfacese are global interfaces index + 1
             for interface in self.has_interface:
+                marker_value = interface + 1
                 # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
-                interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
+                interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
             # form2 is different for wetting and nonwetting
             form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
@@ -284,15 +286,16 @@ class DomainPatch(df.SubDomain):
             # gather interface forms
             interface_forms = []
             for interface in self.has_interface:
+                marker_value = interface + 1
                 if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
                     # if the neighbour of our subdomain (with index self.subdomain_index)
                     # assumes a Richards model, we assume constant normalised atmospheric pressure
                     # and no interface term is practically appearing.
                     # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
-                    interface_forms.append((df.Constant(0)*phi_a)*ds(interface))
+                    interface_forms.append((df.Constant(0)*phi_a)*ds(marker_value))
                 else:
                     # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
-                    interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
+                    interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
             form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             # the non wetting phase has a + before instead of a - before the bracket
             rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
@@ -329,22 +332,19 @@ class DomainPatch(df.SubDomain):
         # vector (numpy array) that we will use to add all dofs of the gli terms
         # into one array.
         gli_assembled_tmp = dict()
-        if self.isRichards:
+        for phase in self.has_phases:
             # print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
-            gli_assembled_tmp = {
-                'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float)
-                }
-        else:
-            gli_assembled_tmp = {
-                'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float),
-                'nonwetting': np.zeros(self.pressure['nonwetting'].vector()[:].size, dtype=float)
-                }
+            gli_assembled_tmp.update(
+                {phase: np.zeros(self.pressure[phase].vector()[:].size, dtype=float)}
+                )
+
         for ind in self.has_interface:
+            marker_value = ind + 1
             interface = self.interface[ind]
             # neighbour index
             neighbour = interface.neighbour[subdomain]
             neighbour_iter_num = interface.current_iteration[neighbour]
-            ds = self.ds(ind)
+            ds = self.ds(marker_value)
             # needed for the read_gli_dofs() functions
             interface_dofs = self._dof_indices_of_interface[ind]
             # for each interface, this is a list of dofs indices belonging to another
@@ -540,40 +540,22 @@ class DomainPatch(df.SubDomain):
         ### END for ind in self.has_interface:
         # save the result of the calculation to the subdomain and to the interface
         # dictionaries.
-        if self.isRichards:
+        for phase in self.has_phases:
             self.gli_assembled.update(
-                {'wetting': gli_assembled_tmp['wetting']}
+                {phase: gli_assembled_tmp[phase]}
                 )
             for ind in self.has_interface:
                 # self._calc_gli_term_on(interface, iteration)
                 interface = self.interface[ind]
                 # write gli dofs to interface dicts for communiction
-                interface.read_gli_dofs(from_array = self.gli_assembled['wetting'], #
-                               interface_dofs = interface_dofs['wetting'],#
-                               dof_to_vert_map = self.dof2vertex['wetting'],#
+                interface.read_gli_dofs(from_array = self.gli_assembled[phase], #
+                               interface_dofs = interface_dofs[phase],#
+                               dof_to_vert_map = self.dof2vertex[phase],#
                                local_to_parent_vertex_map = self.parent_mesh_index,#
-                               phase = 'wetting',#
+                               phase = phase,#
                                subdomain_ind = subdomain,#
                                previous_iter = False
                                )
-        else:
-            self.gli_assembled.update({
-                'wetting': gli_assembled_tmp['wetting'],
-                'nonwetting': gli_assembled_tmp['nonwetting']
-                })
-            for ind in self.has_interface:
-                # self._calc_gli_term_on(interface, iteration)
-                interface = self.interface[ind]
-                for phase in ['wetting', 'nonwetting']:
-                    # write gli dofs to interface dicts for communiction
-                    interface.read_gli_dofs(from_array = self.gli_assembled[phase], #
-                                   interface_dofs = interface_dofs[phase],#
-                                   dof_to_vert_map = self.dof2vertex[phase],#
-                                   local_to_parent_vertex_map = self.parent_mesh_index,#
-                                   phase = phase,#
-                                   subdomain_ind = subdomain,#
-                                   previous_iter = False
-                                   )
     ### END calc_gli_term
 
     def null_all_interface_iteration_numbers(self) -> None:
@@ -634,22 +616,15 @@ class DomainPatch(df.SubDomain):
         mesh_data = self.mesh.data()
         # local to parent vertex index map
         self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
-        if self.isRichards:
-            self.dof2vertex = {#
-                'wetting' :df.dof_to_vertex_map(self.function_space['wetting'])#
-                }
-            self.vertex2dof = {#
-                'wetting' :df.vertex_to_dof_map(self.function_space['wetting'])#
-                }
-        else:
-            self.dof2vertex = {#
-                'wetting'   :df.dof_to_vertex_map(self.function_space['wetting']),#
-                'nonwetting':df.dof_to_vertex_map(self.function_space['nonwetting'])#
-                }
-            self.vertex2dof = {#
-                'wetting'   :df.vertex_to_dof_map(self.function_space['wetting']),#
-                'nonwetting':df.vertex_to_dof_map(self.function_space['nonwetting'])#
-                }
+        self.dof2vertex = dict()
+        self.vertex2dof = dict()
+        for phase in self.has_phases:
+            self.dof2vertex.update(#
+                {phase :df.dof_to_vertex_map(self.function_space[phase])}#
+                )
+            self.vertex2dof.update(#
+                {phase :df.vertex_to_dof_map(self.function_space[phase])}#
+                )
 
     def _init_markers(self) -> None:
         """ define boundary markers
@@ -664,8 +639,9 @@ class DomainPatch(df.SubDomain):
         self.interface_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
         self.interface_marker.set_all(0)
         for glob_index in self.has_interface:
+            marker_value = glob_index + 1
             # each interface gets marked with the global interface index.
-            self.interface[glob_index].mark(self.interface_marker, glob_index+1)
+            self.interface[glob_index].mark(self.interface_marker, marker_value)
         # create outerBoundary objects and mark outer boundaries with 1
         for index, boundary_points in self.outer_boundary_def_points.items():
             # if our domain patch has no outer boundary, this should be reflected
@@ -709,16 +685,6 @@ class DomainPatch(df.SubDomain):
                     {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value)}
                 )
                 print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind][phase])
-            # if self.isRichards:
-            #     self._dof_indices_of_interface[ind].update(#
-            #         {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value)}
-            #     )
-            #     print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind]['wetting'])
-            # else:
-            #     self._dof_indices_of_interface[ind].update(#
-            #         {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value),#
-            #          'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, marker_value)}
-            #     )
 
     def _calc_interface_has_common_dof_indices(self):
         """ populate dictionary self._interface_has_common_dof_indices
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index fd6c65c..d19ff5e 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -180,7 +180,7 @@ dirichletBC = exact_solution
 #
 # sa
 
-mesh_resolution = 3
+mesh_resolution = 4
 
 # initialise LDD simulation class
 simulation = ldd.LDDsimulation(tol = 1E-12)
-- 
GitLab