From 17f3c1bbfa88ade09032c920893b621dd09414c8 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 21 Mar 2019 09:54:30 +0100
Subject: [PATCH] fix domainPatch.write_pressure_to_interfaces()

---
 LDDsimulation.py   | 21 +++++++++++++++++++++
 RR-2-patch-test.py |  7 +++++++
 domainPatch.py     |  8 ++++----
 3 files changed, 32 insertions(+), 4 deletions(-)

diff --git a/LDDsimulation.py b/LDDsimulation.py
index b86e73a..b2e2d11 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -328,6 +328,27 @@ class LDDsimulation(object):
                 subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
                                                     'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
 
+    def _init_DirichletBC(self, interpolation_degree: int = 2):
+        """ set initial values
+        """
+        for num, subdomain in self.subdomain.items():
+            mesh = subdomain.mesh
+            V = subdomain.function_space
+            # p0 contains both pw_0 and pnw_0 for subdomain[num]
+            p0 = self.initial_conditions[num]
+            if subdomain.isRichards:
+                # note that the default interpolation degree is 2
+                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
+                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
+            else:
+                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+                pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
+                                                'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
+                                                    'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+
     def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
         """ evaluate time dependent source terms or initialise them if time == 0
         """
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index bb90d7d..bf9c3f7 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -125,6 +125,13 @@ initial_condition = {
     2: {'wetting': '-x[1]*x[1]'}
 }
 
+exact_solution = {
+    1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
+    2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+}
+
+dirichletBC = exact_solution
+
 # def saturation(pressure, subdomain_index):
 #     # inverse capillary pressure-saturation-relationship
 #     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
diff --git a/domainPatch.py b/domainPatch.py
index ad563c4..5e298c6 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -139,7 +139,7 @@ class DomainPatch(df.SubDomain):
         ### Class variables set by methods
         # dictionary holding the dof indices corresponding to an interface of
         # given interface. see self._calc_dof_indices_of_interfaces()
-        self.dof_indices_of_interface = dict()
+        self._dof_indices_of_interface = dict()
         # measures are set by self._init_measures()
         self.iteration_number = 0
         self.dx = None
@@ -165,7 +165,7 @@ class DomainPatch(df.SubDomain):
     def write_pressure_to_interfaces(self):
         """ save the interface values of self.pressure to all neighbouring interfaces"""
         if self.isRichards:
-            for ind in has_interface:
+            for ind in self.has_interface:
                 interface[ind].write_dofs(from_function = self.pressure['wetting'], #
                                         interface_dofs = self._dof_indices_of_interface[ind]['wetting'],#
                                         dof_to_vert_map = self.dof2vertex['wetting'],#
@@ -173,7 +173,7 @@ class DomainPatch(df.SubDomain):
                                         phase = 'wetting',#
                                         subdomain_ind = self.subdomain_index)
         else:
-            for ind in has_interface:
+            for ind in self.has_interface:
                 for phase in ['wetting', 'nonwetting']:
                     interface[ind].write_dofs(from_function = self.pressure[phase], #
                                             interface_dofs = self._dof_indices_of_interface[ind][phase],#
@@ -311,7 +311,7 @@ class DomainPatch(df.SubDomain):
         """ calculate dof indices of each interface """
         V = self.function_space
         marker = self.boundary_marker
-        for ind in has_interface:
+        for ind in self.has_interface:
             self._dof_indices_of_interface.update(#
                 {ind: dict()}
             )
-- 
GitLab