From 9d3d817df3d5955bcee6ad0efb200ed69528fca3 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 14 Mar 2019 18:05:09 +0100
Subject: [PATCH] update Interface.write_dofs() and Interface.read_dofs()

---
 RR-2-patch-test.py |  21 +++++--
 domainPatch.py     | 139 ++++++++++++++++++++++++++++++++++++---------
 2 files changed, 127 insertions(+), 33 deletions(-)

diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 41057ba..ddc6a82 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -271,15 +271,21 @@ print("vertices of subdomain1: \n", coordinates1)
 # dom1_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker1, 1, print_coordinates = True)
 dom1_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker1, 1, print_vertex_indices = True)
 
-interface[0].write_dofs(fem_function = f1test, #
+interface[0].write_dofs(from_function = f1test, #
                         interface_dofs = dofs_on_interface1,#
                         dof_to_vert_map = dom1_d2v,#
-                        local_to_parent_vertex_map = dom1_parent_vertex_indices)
+                        local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+                        phase = 'wetting',#
+                        pressure = True,#
+                        subdomain_ind = 1)
 
-interface[0].read_dofs(fem_function = f2test, #
+interface[0].read_dofs(to_function = f2test, #
                         interface_dofs = dofs_on_interface2,#
                         dof_to_vert_map = dom2_d2v,#
-                        local_to_parent_vertex_map = dom2_parent_vertex_indices)
+                        local_to_parent_vertex_map = dom2_parent_vertex_indices,#
+                        phase = 'wetting',#
+                        pressure = True,#
+                        subdomain_ind = 1)
 
 i=0
 for x in coordinates2:
@@ -365,10 +371,13 @@ for iterations in range(1):
     df.solve(A1,U1,b1)
     print('solution:\n', U1)
     p1_i.vector()[:] = U1
-    interface[0].write_dofs(fem_function = p1_i, #
+    interface[0].write_dofs(from_function = p1_i, #
                             interface_dofs = dofs_on_interface1,#
                             dof_to_vert_map = dom1_d2v,#
-                            local_to_parent_vertex_map = dom1_parent_vertex_indices)
+                            local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+                            phase = 'wetting',#
+                            pressure = True,#
+                            subdomain_ind = 1)
 #
 # Save mesh to file
 df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
diff --git a/domainPatch.py b/domainPatch.py
index b48d556..6ee9480 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -596,22 +596,42 @@ class Interface(BoundaryPart):
                     self.gl_values[subdom_ind][phase].update({vert_ind: 0.0})
                     self.gl_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
 
-    def write_dofs(self, fem_function: df.Function, #
+    def write_dofs(self, from_function: df.Function, #
                    interface_dofs: np.array,#
                    dof_to_vert_map: np.array,#
-                   local_to_parent_vertex_map: np.array) -> None:
-        """ write dofs of fem_function with indices interface_dofs to self.pressure_values
-
-        Write the degrees of freedom of the function fem_function with indices
-        interface_dofs, i.e. the dofs of the interface, to the dictionary of the
-        interface self.pressure_values.
+                   local_to_parent_vertex_map: np.array,#
+                   phase: str,#
+                   subdomain_ind: int,#
+                   pressure: bool = False,#
+                   gl: bool = False,#
+                   previous_iter: bool = False,#
+                   ) -> None:
+        """ write dofs of from_function with indices interface_dofs to self.pressure_values
+
+        Write the degrees of freedom of the function from_function with indices
+        interface_dofs, i.e. the dofs of the interface, to one of the dictionaries
+        of the interface
+            self.pressure_values[subdomain_ind][phase], (if pressure == True)
+            self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
+            self.gl_values[subdomain_ind][phase], (if gl == True)
+            self.gl_values_prev[subdomain_ind][phase] (if previous_iter == True)
 
         # Parameters
-        fem_function:               #type: df.Function, function to save the dofs to
+        from_function:              #type: df.Function, function to save the dofs to
         interface_dofs:             #type: np.array,    index array of dofs to be saved
         dof_to_vert_map:            #type: np.array,    dof to vert map of the function
                                                         space on the submesh.
         local_to_parent_vertex_map: #type: np.array,
+        phase                       #type: str          is either 'wetting' or
+                                                        'nonwetting' and determins
+                                                        the dict to be written to.
+        subdomain_ind               #type: int          subdomain index
+        pressure                    #type: bool         determins wether pressure
+                                                        values are being written or not
+        gl                          #type: bool         determins wether gl
+                                                        values are being written or not
+        previous_iter               #type: bool         determins wether a previous
+                                                        iterate is being written or not
         """
         if self.pressure_values == None:
             raise RuntimeError("self.pressure_values not initiated. You forgot to \
@@ -619,33 +639,74 @@ class Interface(BoundaryPart):
 
         parent = local_to_parent_vertex_map
         d2v = dof_to_vert_map
-        for dof_index in interface_dofs:
-            # fem_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-            # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-            self.pressure_values.update({parent[d2v[dof_index]]: fem_function.vector()[dof_index]})
-
-        print("interface values\n", self.pressure_values)
+        # chose where to write the dofs to
+        if pressure:
+            if not previous_iter:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+                print("have written pressure interface values\n", self.pressure_values[subdomain_ind][phase])
+            else:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+                print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
+        if gl:
+            if not previous_iter:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    self.gl_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+                print("have written gl interface values\n", self.gl_values[subdomain_ind][phase])
+            else:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    self.gl_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+                print("have written previous gl interface values\n", self.gl_values_prev[subdomain_ind][phase])
 
-    def read_dofs(self, fem_function: df.Function, #
+    def read_dofs(self, to_function: df.Function, #
                   interface_dofs: np.array,#
                   dof_to_vert_map: np.array,#
-                  local_to_parent_vertex_map: np.array) -> None:
-        """ read dofs off self.pressure_values and overwrite the dofs of fem_function
+                  local_to_parent_vertex_map: np.array,#
+                  phase: str,#
+                  subdomain_ind: int,#
+                  pressure: bool = False,#
+                  gl: bool = False,#
+                  previous_iter: bool = False,#
+                  ) -> None:
+        """ read dofs off self.pressure_values and overwrite the dofs of from_function
             with indices interface_dofs
 
-        Read the degrees of freedom of the interface self.pressure_values and
-        overwrite the dofs of the function fem_function with indices
+        Read the degrees of freedom of one of the dictionaries of the interface
+            self.pressure_values[subdomain_ind][phase], (if pressure == True)
+            self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
+            self.gl_values[subdomain_ind][phase], (if gl == True)
+            self.gl_values_prev[subdomain_ind][phase] (if previous_iter == True)
+        and overwrite the dofs of the function to_function with indices
         interface_dofs, i.e. update the interface dofs of the function.
         .
 
         # Parameters
-        fem_function:               #type: df.Function, function to overwrite the
+        to_function:               #type: df.Function, function to overwrite the
                                                         dofs of.
         interface_dofs:             #type: np.array,    index array of interface
-                                                        dofs of fem_function
+                                                        dofs of from_function
         dof_to_vert_map:            #type: np.array,    dof to vert map of the function
                                                         space on the submesh.
         local_to_parent_vertex_map: #type: np.array,
+        phase                       #type: str          is either 'wetting' or
+                                                        'nonwetting' and determins
+                                                        the dict to be read of.
+        subdomain_ind               #type: int          subdomain index
+        pressure                    #type: bool         determins wether pressure
+                                                        values are being read or not
+        gl                          #type: bool         determins wether gl
+                                                        values are being read or not
+        previous_iter               #type: bool         determins wether a previous
+                                                        iterate is being read or not
         """
         if self.pressure_values == None:
             raise RuntimeError("self.pressure_values not initiated. You forgot to \
@@ -653,12 +714,36 @@ class Interface(BoundaryPart):
 
         parent = local_to_parent_vertex_map
         d2v = dof_to_vert_map
-        for dof_index in interface_dofs:
-            # fem_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-            # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-            fem_function.vector()[dof_index] = self.pressure_values[parent[d2v[dof_index]]]
-
-        print("overwritten function \n", fem_function.vector()[:])
+        if pressure:
+            if not previous_iter:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    to_function.vector()[dof_index] = self.pressure_values[subdomain_ind][phase][parent[d2v[dof_index]]]
+                print("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
+                print("wrote these dofs to function \n", to_function.vector()[:])
+            else:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    to_function.vector()[dof_index] = self.pressure_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
+                print("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
+                print("wrote these dofs to function \n", to_function.vector()[:])
+        if gl:
+            if not previous_iter:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    to_function.vector()[dof_index] = self.gl_values[subdomain_ind][phase][parent[d2v[dof_index]]]
+                print("read gl interface values\n", self.gl_values[subdomain_ind][phase])
+                print("wrote these dofs to function \n", to_function.vector()[:])
+            else:
+                for dof_index in interface_dofs:
+                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                    to_function.vector()[dof_index] = self.gl_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
+                print("read gl interface values\n", self.gl_values_prev[subdomain_ind][phase])
+                print("wrote these dofs to function \n", to_function.vector()[:])
 
     #### Private Methods ####
 
-- 
GitLab