From cd22b41bf7d7239fc3cf9495127c07ad6377f991 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 5 Jun 2020 19:51:56 +0200
Subject: [PATCH] remove if condition with gravity in calculation of gl0

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

diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index af91ef7..217e4ec 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -504,6 +504,7 @@ class DomainPatch(df.SubDomain):
 
                 flux_function = df.project(flux, self.function_space["flux"][phase])
                 flux_x, flux_y = flux_function.split()
+                # functions get initialised with values zero.
                 gl0 = df.Function(self.function_space["gli"][phase])
                 # # get dictionaries of global dof indices and coordinates along
                 # # the interface. These dictionaries have the facet indices of
@@ -528,7 +529,9 @@ class DomainPatch(df.SubDomain):
                         flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
                         flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
 
-                        if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
+                        if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+                            # note that even if we include gravity, the correct
+                            # flux has been set above (flux=0 without gravity).
                             gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
                                 + flux_y.vector()[flux_y_dof_index]*normal[1]
                         else:
-- 
GitLab