add DomainPatch._calc_normal_vectors_norms_for_common_dofs method

LDDsimulation/domainPatch.py

@@ -211,6 +211,7 @@ class DomainPatch(df.SubDomain):
         self._init_measures()
         self._calc_dof_indices_of_interfaces()
         self._calc_interface_has_common_dof_indices()
+        self._calc_normal_vectors_norms_for_common_dofs()
         if self.include_gravity:
             self.gravity_term = dict()
             self._calc_gravity_expressions()
@@ -481,26 +482,35 @@ class DomainPatch(df.SubDomain):
                             flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
                         else:
                             flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
-
-                    gl0 = df.dot(flux, n) - Lambda[phase]*pa_prev
-                    gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
+                    # we need to treat the flux part of gl0 and the robin (pressure)
+                    # part of gl0 seperately at the intersection points of two
+                    # interfaces (more details in comment below). Therefor we save
+                    # the robin pressure part and the neumann flux part separately
+                    robin_pressure = Lambda[phase]*pa_prev
+                    robin_pressure_assembled = df.assemble(dt*robin_pressure*phi_a*ds).get_local()
+                    neumann_flux_assembled = df.assemble(dt*df.dot(flux, n)*phi_a*ds).get_local()
+                    # gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
 
                     # if dofs_in_common_with_other_interfaces[phase].size == 0
                     # there are no dofs that lie on another interface and we can
                     # skip the following step.
+                    # If we have dofs in common with another interface we get two
+                    # different gl0 terms from both interfaces apriori, so in
+                    # order to get one single dof for the intersection we do:
+                    # we replace the neuman flux by the flux in dircetion of the
+                    # linear combination of both adjacent normal vectors:
+                    # flux_neumann = df.dot(flux, n_1 + 1_2)/norm(n_1 + 1_2)
+                    # since we are looping over the interfaces we correct each dof
+                    # that is common with another interfaces with the corresponding
+                    # factor norm(n_1 + 1_2).
                     if dofs_in_common_with_other_interfaces[phase].size != 0:
                         if debug:
                             print(f"the following dofs (phase = {phase}) of interface{interf_ind} are in common with \n",
                                 f" other interfaces: {dofs_in_common_with_other_interfaces[phase]}")
                         for common_dof in dofs_in_common_with_other_interfaces[phase]:
                             # from the standpoint of the subdomain we are currently on,
-                            # each dof can belong maximally to two interfaces. On these
-                            # vertices, common to two interfaces, the dofs of the
-                            # corresponding gli terms of both interfaces are averaged
-                            # to give one value for this vertex. Because we achieve this
-                            # by adding the assembled gli terms together we have to
-                            # multiply the dofs that appear in more than one interface
-                            # by 1/2.
+                            # each dof can belong maximally to two interfaces.
+                            # because we have not
                             gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
                     # now, add the just assembled and corrected term to the Global
                     # gli_assembled_tmp vector
@@ -882,6 +892,29 @@ class DomainPatch(df.SubDomain):
                 #     print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
 
 
+    def _calc_normal_vectors_norms_for_common_dofs(self):
+        """ calculate norm(n1+n2) where n1 and n2 are normal vectors of adjacent
+        interfaces. This is used by the functions calc_gli_term and calc_gl0_term
+        """
+        self._normal_vector_norms_for_common_dofs = dict()
+        for interface_ind in self.has_interface:
+            # if one phase has common dofs to two interfaces, the both phases
+            # have, so it suffices to check for the wetting phase as it is always
+            # present.
+            if self._interface_has_common_dof_indices[interface_ind]['wetting'].size != 0:
+                self._normal_vector_norms_for_common_dofs.update({interface_ind: dict()})
+                for phase in self.has_phases:
+                    self._normal_vector_norms_for_common_dofs[interface_ind].update(
+                        {phase: dict()}
+                    )
+                    for dof in self._interface_has_common_dof_indices[interface_ind][phase]:
+                        # at the moment we only determine this correctly if the patch
+                        # is rectangular and parallel to the x and y axis.
+                        self._normal_vector_norms_for_common_dofs[interface_ind][phase].update(
+                            {dof: np.sqrt(2)}
+                        )
+
+
     def _calc_gravity_expressions(self):
         """ calculate the gravity term if it is included"""
         g = df.Constant(self.gravity_acceleration) #