implement general _normal_vector_norms_for_common_dofs with linear combinations of normal vectors

LDDsimulation/domainPatch.py

@@ -391,7 +391,7 @@ class DomainPatch(df.SubDomain):
         return form_and_rhs
 
 
-    def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
+    def calc_gl0_term(self, debug: bool = True) -> None: # tp.Dict[str, fl.Form]
         """calculate the gl0 terms for all interfaces of the subdomain.
 
         assemble the gl0 terms for all interfaces of the subdomain and save them
@@ -503,15 +503,21 @@ class DomainPatch(df.SubDomain):
                     # 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).
+                    # We take the mean of the pressure on both interfaces, so
+                    # we multiply these by 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]:
+                            n_norm = self._normal_vector_norms_for_common_dofs[interf_ind][phase][common_dof]
                             # from the standpoint of the subdomain we are currently on,
                             # each dof can belong maximally to two interfaces.
                             # because we have not
-                            gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
+                            neumann_flux_assembled[common_dof] = neumann_flux_assembled[common_dof]/n_norm
+                            robin_pressure_assembled[common_dof] = robin_pressure_assembled[common_dof]/2
+
+                    gl0_assembled = neumann_flux_assembled - robin_pressure_assembled
                     # now, add the just assembled and corrected term to the Global
                     # gli_assembled_tmp vector
                     gli_assembled_tmp[phase] += gl0_assembled
@@ -652,12 +658,16 @@ class DomainPatch(df.SubDomain):
                     phi_a  = self.testfunction[phase]
                     gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
                     gli_pressure_part_array = gli_pressure_part.get_local()
-                    gli = gli_pressure_part_array - gli_prev_of_neighbour
                     # check if there are shared dofs with another interface.
                     if dofs_in_common_with_other_interfaces[phase].size != 0:
                         for common_dof in dofs_in_common_with_other_interfaces[phase]:
+                            # n_norm = self._normal_vector_norms_for_common_dofs[ind][phase][common_dof]
+                            gli_prev_of_neighbour[common_dof] = gli_prev_of_neighbour[common_dof]/2
+
                             # see previous case.
-                            gli[common_dof] = 0.5*gli[common_dof]
+                            gli_pressure_part_array[common_dof] = 0.5*gli_pressure_part_array[common_dof]
+
+                    gli = gli_pressure_part_array - gli_prev_of_neighbour
                     # now, add the just assembled and corrected term to the Global
                     # gli_assembled_tmp vector
                     gli_assembled_tmp[phase] += gli
@@ -897,6 +907,7 @@ class DomainPatch(df.SubDomain):
         interfaces. This is used by the functions calc_gli_term and calc_gl0_term
         """
         self._normal_vector_norms_for_common_dofs = dict()
+        self.mesh.init(0)
         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
@@ -908,10 +919,27 @@ class DomainPatch(df.SubDomain):
                         {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.
+                        common_vertex_index = self.dof2vertex[phase][dof]
+                        for index, vertex in enumerate(df.vertices(self.mesh)):
+                            if index == common_vertex_index:
+                                common_vertex = vertex
+                                # print(f"\non subdomain{self.subdomain_index}, interface{interface_ind} and phase {phase}")
+                                # print(f"vertex_{index} = {vertex.point()[:]}, type(vertex) = {type(vertex)}")
+                        linear_combination_of_normals = np.zeros(np.size(vertex.point()[:]), dtype=float)
+                        # normals = []
+                        # print(f"common_vertex = {common_vertex}",
+                        #       f"common_vertex_index = {common_vertex_index}")
+                        for facet_containing_vertex in df.facets(common_vertex):#df.cells(common_vertex):
+                            if facet_containing_vertex.exterior():
+                                # normals.append(facet_containing_vertex.normal()[:])
+                                # print(f"type(facet_normal) = {type(facet_containing_vertex.normal()[:])}")
+                                linear_combination_of_normals += facet_containing_vertex.normal()[:]
+
+                        # print(f"normals of common_vertex: {common_vertex} are ", normals)
+                        # print(f"linear combination is: {linear_combination_of_normals[:]}.")
+                        # print(f"with length {np.linalg.norm(linear_combination_of_normals, ord=2, axis=None)}. For comparison, np.sqrt(2) = {np.sqrt(2)}.")
                         self._normal_vector_norms_for_common_dofs[interface_ind][phase].update(
-                            {dof: np.sqrt(2)}
+                            {dof: np.linalg.norm(linear_combination_of_normals)}
                         )