Skip to content
Snippets Groups Projects
Select Git revision
  • 8cee0973b56adca3c1c3992048bbc03b19973fd7
  • master default protected
  • parallelistation_of_subdomain_loop
  • parallelistation_test
  • nonzero_gli_term_for_nonwetting_in_TPR
  • testing_updating_pwi_as_pwiminus1_into_calculation_of_pnwi
  • v2.2
  • v2.1
  • v2.0
  • v1.0
10 results

boundary_and_interface.py

Blame
  • boundary_and_interface.py 38.97 KiB
    #!/bin/python3
    """ BoundaryPart and Interface class
    
    This module defines the classes BoundaryPart and its child, the class Interface
    used by the class LDDsimulation. The main purpose of the class BoundaryPart is
    to define the method „inside“ which is used to mark boundary parts of each subdomain.
    The class Interface is derived from BoundaryPart but additionally defines dictionaries
    to hold interface pressuers and the gli terms of the L-scheme as well as read and
    write methods to those dictionaries that are used in the communication of the
    dofs between subdomains.
    
    """
    import dolfin as df
    import numpy as np
    import typing as tp
    import copy as cp
    import numpy.linalg as la
    
    
    class BoundaryPart(df.SubDomain):
        """BoundaryPart class defines method self.inside to mark 1D polygonal parts
           of the boundary.
    
        The class BoundaryPart defines the method self.inside, which returns boolean
        True if a point x on the mesh of a given 2D polygonal subdomain lies on the
        part of the 1D polygonal subdomain boundary defined by the dolfin points in
        the array vertices.
    
        # Parameters
        vertices           #type List[df.Points], dolfin points defining a polyongal
                            part of the domain boundary.
        tol                #type float, the calculation tolerance
        internal           #type bool,  set weather the 1D polygone defined by vertices
                            is considered internal (an interface) or really a part
                            of the boundary. This influences how the public method
                            inside is defined.
        marker_value        #type int, optional     marker value assotiated
                                                        with the boundary part
        """
    
        def __init__(self, #
                     vertices: tp.List[df.Point], #
                     tol = None, #
                     internal: bool = False,
                     marker_value: int = None):
            # because we declare our own __init__ method for the derived class BoundaryPart,
            # we overwrite the __init__-method of the parent class df.SubDomain. However,
            # we need the methods from the parent class, so we call the __init__-method
            # from df.SubDomain to access it.
            super().__init__()
    
            self.vertices = vertices
            # if tol is not None, set it to passed value, else keep the default
            if tol:
                self.tol = tol
            else:
                self.tol = df.DOLFIN_EPS
            self.internal = internal
            if marker_value is None:
                self.marker_value = 10^10
            else:
                self.marker_value = marker_value
    
        ###### Public interface
        def inside(self, x: np.ndarray, on_boundary) -> bool:
            """inside returns true if x is on_boundary and on the part of the boundary
               defined by vertices
    
            # Parameters
            self         #type class BoundaryPart, reference to the class object
            x            #type np.ndarray,   grid point, used by parent class method df.SubDomain.mark.
            """
            if self.internal: # don't include test for on_boundary
                return self._is_on_boundary_part(x) #on_boundary and
    
            # include test for on_boundary in the default internal=False case
            return on_boundary and self._is_on_boundary_part(x)
        # END inside
    
        def coordinates(self, #
                        boundary_marker: df.MeshFunction, #
                        boundary_marker_value: int,#
                        print_coordinates: bool = False) -> np.array:
            """ return coordinates and numbering indices of nodes on mesh marked by
                boundary_marker with value boundary_marker.
    
            coordinates() returns coordinates and numbering indices of nodes on
            the mesh marked by boundary_marker with value
            boundary_marker. boundary_marker is a dolfin.MeshFunction of dimension
            boundary_marker.mesh().topology().dim()-1 that is supposed to mark edges
            of the mesh according to the method self.inside().
    
            # Parameters
            boundary_marker:        #type df.MeshFunction   marker function marking edges of i_mesh
                                                            according to self.inside
            boundary_marker_value:  #type int               marker value of boundary_marker
            print_coordinates       #type boolean           specify wether to print out the
                                                            coordinates or not.
            """
            mesh_coordinates = boundary_marker.mesh().coordinates()
            mesh_dimension = boundary_marker.mesh().topology().dim()
            # boundary_marker is a function defined on edges. Therefore
            # sum(boundary_marker.array() == boundary_marker_value) gives the number
            # of edges on the boundary. We need the number of nodes, however.
            number_of_boundary_nodes = sum(boundary_marker.array() == boundary_marker_value) + 1
            # we need one mesh_dimension + 1 columns to store the the index of the node.
            node_coordinates = np.zeros(shape = (number_of_boundary_nodes, mesh_dimension + 1)) #, dtype=float
            node_num = 0
            mesh_vertex_index = 0
            for x in mesh_coordinates:
                if self._is_on_boundary_part(x):
                    node_coordinates[node_num,:] = [x[0], x[1], mesh_vertex_index]
                    node_num += 1
    
                mesh_vertex_index += 1
    
            if print_coordinates:
                print("Interface Coordinates: \n", node_coordinates)
    
            return node_coordinates
    
    
        ###### Private methods
        def _is_on_boundary_part(self, point: np.ndarray) -> bool:
            """_is_on_boundary_part returns true if dist(point,boundary_part) < tol
    
            The function _is_on_boundary_part checks if point (a 2D/3D point) is
            close or on the part of the boundary defined by the array of dolfin points in vertices.
    
            # Parameters
            point          #type np.ndarray,    grid point to be testet.
            """
            #print("Annotations:", self._is_on_boundary_part.__annotations__)
            #print("Documentation String:", self._is_on_boundary_part.__doc__)
            tol = self.tol
    
            p = point
            # print("length of vertices:", len(self.vertices))
            # print("range(len(self.vertices))", range(len(self.vertices)))
            # check if p lies between p1 und p2 for all subsequent points in vertices
            for i in range(len(self.vertices) - 1):
                p1 = self.vertices[i]
                p2 = self.vertices[i+1]
                if self._is_on_line_segment(p1, p2, p):
                    return True
    
            # if _is_on_boundary_part reaches to here, p didn't lie on any line segment.
            return False
    
        def _is_on_line_segment(self, point1: df.Point, #
                               point2: df.Point, #
                               point: df.Point, #
                               dimension: int = 2) -> bool:
            """_is_on_line_segment returns true if dist(point,(point2 - point1)) < tol
    
            The function _is_on_line_segment checks if point (a 2D/3D point) is
            close or on the line segment defined by the two points point2 and point1.
            """
            #print("Annotations:", is_on_line_segment.__annotations__)
            #print("Documentation String:", is_on_line_segment.__doc__)
            tol = self.tol
            # point comes from the mesh and is inserted by the df.SubDomain.mark method.
            # It is already and np.ndarray and in the desired format.
            p =  point
            # point1 and point2 are df.Points. numpy functions expect points as
            # coordinate arrays and not as dolfin.Points.
            # Therefore we convert dolfin.Points to coordinate arrays.
            p1 = point1.array()
            p2 = point2.array()
            # dolfin Points are by default 3D. Zeros get appended to the 2d points we
            # used to define our subdomains with when using point1.array().
            # In 2d we cut this off.
            if dimension == 2:
                p1 = np.array([p1[0], p1[1]])
                p2 = np.array([p2[0], p2[1]])
    
            # check if p is either p1 or p2
            xmin = min(p1[0], p2[0])
            xmax = max(p1[0], p2[0])
            ymin = min(p1[1], p2[1])
            ymax = max(p1[1], p2[1])
    
            # check there holds p1[0] < p[0] < p2[0]. If not, p  cannot be on the line segment
            # same needs to be done for p1[1] < p[1] < p2[1]
    
            # the test ((p[0] - xmax)*(p[0] - xmin) < 0) might fail if p[0] and one
            # of p1[0] or p2[0] are equal. In this case the line pp1 or pp2 is
            # vertical. We need to still calculate the distance to the line segment
            # in this case. Same is true for ((p[1] - ymax)*(p[1] - ymin) < 0)
            x_is_in_range = ((p[0] - xmax)*(p[0] - xmin) < tol) or np.fabs((p[0] - xmax)*(p[0] - xmin)) < tol
            y_is_in_range = ((p[1] - ymax)*(p[1] - ymin) < tol) or np.fabs((p[1] - ymax)*(p[1] - ymin)) < tol
            if x_is_in_range and y_is_in_range:
                # here we have a chance to actually lie on the line segment.
                segment_direction = (p2 - p1)/np.linalg.norm(p2 - p1)
                distance = (p - p1) - np.dot(segment_direction, p - p1)*segment_direction
                # check again for equality with p1 and p2 just to be sure. =)
                if (np.linalg.norm(p - p1) < tol) or (np.linalg.norm(p - p2) < tol):
                    # print(f"grid node at {p} was close to either {p1} or {p2}")
                    return True
                elif np.linalg.norm(distance) < tol:
                    # print(f"point to be tested: p = {p}")
                    # print(f"the distance vector is: {distance}")
                    # print(f" Distance of grid node at {p} was close to segment {p1} - {p2}")
                    return True
                else:
                    return False
            else:
                return False
    
        def _print_annotations(self):
            print("Annotations:", self.__annotations__)
            print("Documentation String:", self.__doc__)
    
    ### End Class BoundaryPart
    
    class Interface(BoundaryPart):
        """Interface class adds interface specific methods like dof mapping from
            neighbouring domains to the class BoundaryPart.
    
        The class Interface builds on the class BoundaryPart, but adds mesh dependent
        functionality like vertex number mappings of the interface for the two neighbouring
        subdomains.
    
        # Parameters
        domain_index            #type int,   (global) index of the patch we are currently on
        adjacent_subdomains     #type int,   (global) index of the neighbour patch
                                of whitch this interface is an interface to.
        """
        def __init__(self, #
                    global_index: int,#
                    adjacent_subdomains: np.array = None,#
                    # this is an array of bools
                    isRichards: np.array = None,#
                    **kwds):
    
    
            # make the parent class methods available.
            super().__init__(**kwds)
            #
            self.isRichards = isRichards
            self.adjacent_subdomains = adjacent_subdomains
            # global index ob the subdomain
            self.global_index = global_index
            self.marker_value = self.global_index+1
            # dictionary holding the facet objects
            self.facets = dict()
            # dictionary containing for each facet belonging to the interface (with
            # respect to a subdomain mesh) the coordinates of the the vertices
            # belonging to that facet. This dictionary is created by
            # self.init_facet_dictionaries()
            self.facet_to_vertex_coordinates = dict()
            # dictionary containing maps mapping subdomain mesh indices of facets
            # to the index of the facet in the global mesh and vice versa.
            self.local_to_global_facet_map = dict()
            self.global_to_local_facet_map = dict()
            # dictionary saving per subdomain a list of outer normals per facet defined
            # on a mesh passed to the method self.init_facet_dictionaries() which
            # populates this dictionary.
            self.outer_normals = dict()
            # pressure_values is a dictionary of dictionaries, one for each of the
            # adjacent subdomains, holding the dof values over the interface of the
            # pressure (solution) values for communication. It is initialised by the
            # method init_interface_values().
            self.pressure_values = dict()
            self.pressure_values_prev = dict()
            # dictionaries holding the forms for the decoupling terms
            self.gli_term = dict()
            self.gli_term_prev = dict()
            # dictionary holding the current iteration number i for each of the adjacent
            # subdomains with index in adjacent_subdomains
            self.current_iteration = {#
                self.adjacent_subdomains[0]: 0,#
                self.adjacent_subdomains[1]: 0
            }
            # dictionary holding a bool variable for each subdomain index in self.adjacent_subdomains
            # indicating wether the neighbouring subdomain w.r.t. that index assumes
            # Richards model or not.
            self.neighbour_isRichards = {
                self.adjacent_subdomains[0]: self.isRichards[self.adjacent_subdomains[1]],#
                self.adjacent_subdomains[1]: self.isRichards[self.adjacent_subdomains[0]]
            }
            # dictionary holding the neighbouring domain index relative to the key index
            # this is used for calculating the gl_terms.
            self.neighbour = {
                self.adjacent_subdomains[0]: self.adjacent_subdomains[1],#
                self.adjacent_subdomains[1]: self.adjacent_subdomains[0]#
            }
    
        # MAGIC METHODS
        def __str__(self):
            """print out a few usefull informations of the interface"""
            string = f"\nInterface{self.global_index} with "\
                  +f" global marker_value = {self.marker_value}.\n"\
                  +f" Neighbouring subdomains are subdomain{self.adjacent_subdomains[0]}"\
                  +f" and subdomain{self.adjacent_subdomains[1]}"
            return string
    
        #### Public Methods #####
        # def set_domain_index(self, domain_index: int):
        #     self.domain_index = domain_index
    
        # def dofs_on_interface(self,#
        #                     function_space: df.FunctionSpace,#
        #                     interface_marker: df.MeshFunction, #
        #                     interface_marker_value: int) -> np.array:
        #     """ return the indices of the degrees of freedom on the interface.
        #
        #     dofs_on_interface() returns the indices of the degrees of freedom on the
        #     interface with respect to the dof numbering on function_space.mesh().
        #     In essence, it is the restriction of the dof_to_vertex_map to the interface.
        #
        #     # Parameters
        #     function_space:         #type df.functionSpace  the function space of whitch
        #                                                     to calculate the dof_to_vertex_map
        #                                                     restriction from.
        #     interface_marker:       #type df.MeshFunction   marker function marking edges of i_mesh
        #                                                     according to self.inside
        #     interface_marker_value: #type int               marker value of interface_marker
        #     """
        #     interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value,
        #                                                     print_vertex_indices = False)
        #     # print(f"local interface vertices: {interface_vertex_indices}")
        #     dof2vertex = df.dof_to_vertex_map(function_space)
        #     # print(f"dof2vertex map: {dof2vertex}\n")
        #     # check which dofs actually lie on the interface using the dof2vertex map.
        #     is_a_dof_on_interface = np.isin(dof2vertex, interface_vertex_indices, assume_unique=True)
        #     # print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
        #     # print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
        #     return np.nonzero(is_a_dof_on_interface)[0]
    
    
        def dofs_and_coordinates(self, function_space: df.FunctionSpace,#
                                       interface_marker: df.MeshFunction, #
                                       interface_marker_value: int,
                                       debug: bool = False) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
            """ for a given function space function_space, calculate the dof indices
            and coresponding coordinates along the interface. The pair {global_dof_index: dof_coordinate}
            is saved as value in a dictionary which has global facet indices as keys.
    
            # Parameters
            function_space:         #type df.functionSpace  the function space of whitch
                                                            to calculate the dof_to_vertex_map
                                                            restriction from.
            interface_marker:       #type df.MeshFunction   marker function marking edges of i_mesh
                                                            according to self.inside
            interface_marker_value: #type int               marker value of interface_marker
    
            # Output                #type tp.Dict[int, tp.Dict[int, np.ndarray]]
                                                        the first key is the index of
                                                        the facet relative to the numbering
                                                        of the mesh that interface_marker
                                                        is defined on. The value to
                                                        that key is a dictionary
                                                        containing global (not cell)
                                                        dof index and corresponding
                                                        coordinates of the dof as
                                                        key: value pairs.
    
            """
            dofmap = function_space.dofmap()
            element = function_space.element()
            facet_mesh_entity_iterator = df.SubsetIterator(interface_marker, interface_marker_value)
            dofs_and_coords = dict()
            for facet_mesh_entity in facet_mesh_entity_iterator:
                # this is the facet index relative to the Mesh
                # that interface_marker is defined on.
                interface_facet_index = facet_mesh_entity.index()
                # if debug:
                #     print("\n")
                # for interface_facet in df.facets(facet_mesh_entity):
                #     if debug:
                #         print(f"interface facet with mesh entity index: {facet_mesh_entity.index()}")
                #         print(f"interface facet with facet index: {interface_facet.index()}")
                #     # we only need the x,y coordinates. vertex.point()[:] returns
                #     # the array and cannot be sliced. the result can be indexed, hence,
                #     # vertex.point()[:][0:2]
                #     # facet_midpoint = interface_facet.midpoint()[:][0:2]
                #     # print(f"facet: {interface_facet}")
                #     # print("facet coordinates: ")
                facet_points = []
                for vertex in df.vertices(facet_mesh_entity):
                    # if debug:
                    #     print(f"vertex: {vertex} with coordinates {vertex.point()[:][0:2]}")
                    facet_points.append(vertex.point())
    
    
                dofs_and_coords.update({interface_facet_index: dict()})
                # because facet_mesh_entity doesn't have a method to access coordinates we
                # need to cast it into an acutal facet object.
                for cell in df.cells(facet_mesh_entity):
                    # we actually only have one facet contained in this facet_mesh_entity
                    # but we still need to loop over „all“ facets since df.facets(facet_mesh_entity)
                    # returns an iterator object.
                    facet_cell = cell
                    facet_cell_index = facet_cell.index()
                    # if debug:
                    #     print(f"facet_cell_index on the cell belonging to facet_mesh_entity {interface_facet_index} is {facet_cell.index()}")
                    # cell_facets_entity = df.facets(cell)
                    # cell_coordinates = cell.get_vertex_coordinates()
    
                # the cell is a triangle and for each dof numbered locally for that
                # cell, we have the coordinates in a numpy array. The numbering is the
                # numbering local to the cell.
                facet_cell_dof_coordinates = element.tabulate_dof_coordinates(facet_cell)
                # the cell dof map gives to each local (to the cell) numbering of the
                # dofs the index of the dofs realitve to the global dof numbering of
                # space, i.e. the index to access the dof with function.vector()[index].
                facet_cell_dofmap = dofmap.cell_dofs(facet_cell_index)
    
                for local_dof_ind, dof_coord in enumerate(facet_cell_dof_coordinates):
                    # we only want to store the dofs that are on the boundary. The
                    # cell is a triangle and only has one facet that belongs to
                    # the interface. we want to store only dofs that lie on that particular
                    # facet
                    if self._is_on_line_segment(facet_points[0], facet_points[1], dof_coord):
                        global_dof_ind = facet_cell_dofmap[local_dof_ind]
                        dofs_and_coords[interface_facet_index].update({global_dof_ind: dof_coord})
    
    
            if debug:
                print(f"The newly defined facet_dof_to_coordinates_dictionary of interface[{interface_marker_value-1}]")
                number_of_interface_facets = 0
                for facets in dofs_and_coords.items():
                    number_of_interface_facets += 1
                print(f"contains {number_of_interface_facets} interface facets, so {number_of_interface_facets + 1} dofs")
                for facet_index, facet_dofs in dofs_and_coords.items():
                    print(f"on facet {facet_index} we have the dofs", facet_dofs)
            return cp.deepcopy(dofs_and_coords)
    
    
        def init_facet_dictionaries(self,#
                                interface_marker: df.MeshFunction,#
                                interface_marker_value: int,
                                subdomain_index: int,
                                debug: bool = False) -> None:
            """ create dictionary for each index subdomain_index containing the
            index of each facet belonging to the interface (with respect to the mesh of
            the subdomain subdomain_index) as key and stores a
            list of the coordinates of each vertex (numpy arrays) belonging to the facet as
            values.
            Additionally, create dictionary storing per facet the outer normals, i.e.
            pairs {facet_index: facet.normal()}.
            This is needed by SubDomain.calc_gl0_term() for accessing the normals.
    
            The MeshFunction interface_marker is supposed to be defined on the
            subdomain mesh belonging to subdomain subdomain_index.
    
            # Parameters
            interface_marker:       #type df.MeshFunction   marker function marking
                                                            edges of i_mesh
                                                            according to self.inside
            interface_marker_value: #type int               marker value of interface_marker
            """
            self.facet_to_vertex_coordinates.update({subdomain_index: dict()})
            self.outer_normals.update({subdomain_index: dict()})
            self.facets.update({subdomain_index: dict()})
            # Get a subset iterator for mesh entities corresponding to the facets marked by
            # interface_marker and value interface_marker_value.
            interface_facet_mesh_entities = df.SubsetIterator(interface_marker, interface_marker_value)
            # we actually want to access things like the outer normal and coordinates
            # of the vertices belonging to the facet. The entity iterator provides a
            # way to iterate over facets but as instances of the more abstract class
            # dolfin.cpp.mesh.MeshEntity. We need to find for each facet_entity the
            # actual instance of class dolfin.cpp.mesh.Facet, which then allows us
            # to access the outer normal.
            if debug:
                print(f"\non subdomain {subdomain_index}")
            for facet_mesh_entity in interface_facet_mesh_entities:
                # df.facets(facet_mesh_entity) gives us for each facet_entity the
                # an instance of class dolfin.cpp.mesh.Facet. Yes, it's annoying,
                # since we actually only have one facet (of class dolfin.cpp.mesh.Facet)
                # contained in  facet_mesh_entity but we still need to loop over
                # „all“ facets since df.facets(facet_mesh_entity) returns an iterator object.
                # if debug:
                #     print(f"\nfacet_mesh_entity: {facet_mesh_entity} of type {type(facet_mesh_entity)}")
                for facet_instance in df.facets(facet_mesh_entity):
                    facet = facet_instance
                    # if debug:
                    #     print(f"facet (type: {type(facet)}) corresponding to facet_mesh_entity: {facet}")
    
                facet_vertices = []
                for vertex in df.vertices(facet):
                    # we only need the x,y coordinates. vertex.point()[:] returns
                    # the array and cannot be sliced. the result can be indexed, hence,
                    # vertex.point()[:][0:2]
                    facet_vertices.append(vertex.point()[:][0:2])
    
                self.facet_to_vertex_coordinates[subdomain_index].update(
                    {facet.index(): facet_vertices}
                )
                # if debug:
                #     print(f"saved facet vertices as {self.facet_to_vertex_coordinates[subdomain_index].items()}")
                self.facets[subdomain_index].update({facet.index(): facet})
                if subdomain_index == 0:
                    # we do not have outer normals if this interface is part of the
                    # global domain (subdomain_index = 0).
                    self.outer_normals[subdomain_index].update(
                        {facet.index(): None}
                    )
                else:
                    self.outer_normals[subdomain_index].update(
                        {facet.index(): facet.normal()[:][0:2]}
                    )
                    if debug:
                        print(f"normal has norm: {la.norm(facet.normal()[:][0:2])}")
                    #     print(f"saved outer normall as {self.outer_normals[subdomain_index].items()}")
                    # print(f"outer normal  {facet.normal()[:][0:2]}: belonging to facet {facet}")
                    # problem_point = np.array([0.7, -0.2])
                    # vertex1ispoint = np.allclose(problem_point, self.facet_to_vertex_coordinates[subdomain_index][facet.index()][0])
                    # vertex2ispoint = np.allclose(problem_point, self.facet_to_vertex_coordinates[subdomain_index][facet.index()][1])
                    # if vertex1ispoint or vertex2ispoint:
                    #     print(f"on subdomain {subdomain_index}:")
                    #     print(f"the normal of the facet belonging to {self.facet_to_vertex_coordinates[subdomain_index][facet.index()][0]} and point {self.facet_to_vertex_coordinates[subdomain_index][facet.index()][1]}, problem point = {problem_point} is", facet.normal()[:][0:2])
                    #     # raise RuntimeError('stopped for debugging')
            if subdomain_index != 0:
                # if subdomain_index is not equal to 0, we are on an actual SubDomain
                # and self.init_facet_dictionaries() has been run for the global domain
                # already. This means we actually can compute a local2global facet map.
                self.init_facet_maps(subdomain_index)
            if debug:
                print(f"newly defined facet_to_vertex_coordinates_dictionary[{subdomain_index}]")
                print(self.facet_to_vertex_coordinates[subdomain_index].items())
                print(f"outer normals per facet are: ", self.outer_normals[subdomain_index].items())
    
    
    
        def init_facet_maps(self, subdomain_index: int, debug: bool = False) -> None:
            """ calculate the mapping that maps the local numbering of interface
            facets (w.r.t. the mesh of subdomain subdomain_index) to the global
            index with respect to the global domain (subdomain 0)"""
            if subdomain_index == 0:
                raise RuntimeError("expeced a subdomain index != 0 for calculation",
                                         "of local_to_global_facet_map.")
    
            self.local_to_global_facet_map.update({subdomain_index: dict()})
            self.global_to_local_facet_map.update({subdomain_index: dict()})
            global_facet2coord = self.facet_to_vertex_coordinates[0]
            subdom_facet2coord = self.facet_to_vertex_coordinates[subdomain_index]
            for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
                for global_facet_index, global_facet_vertex in global_facet2coord.items():
                    # vertices are not necessarily in the same order
                    # np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
                    vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    both_vertices_are_equal = (vertices_are_equal == 2)
                    if both_vertices_are_equal:
                        self.local_to_global_facet_map[subdomain_index].update(
                            {local_facet_index: global_facet_index}
                        )
                        self.global_to_local_facet_map[subdomain_index].update(
                            {global_facet_index: local_facet_index}
                        )
            if debug:
                print(f"\n on subdomain {subdomain_index}")
                print("local_to_global_facet_map: ", self.local_to_global_facet_map[subdomain_index].items())
                print("global_to_local_facet_map: ", self.global_to_local_facet_map[subdomain_index].items())
                print(f"\nlocal_to_global_facet_map sanity check:")
                for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
                    global_index = self.local_to_global_facet_map[subdomain_index][local_facet_index]
                    global_facet_vertex = global_facet2coord[global_index]
                    vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    both_vertices_are_equal = (vertices_are_equal == 2)
                    if both_vertices_are_equal:
                        print(f"facet with local index {local_facet_index} has correctly been identyfied to global number {global_index}")
                    else:
                        print("FALSE FALSE FALSE")
                        print(f"facet with local index {local_facet_index} has INCORRECTLY been identyfied to global number {global_index}")
    
                print(f"\nglobal_to_local_facet_map sanity check:")
                for global_facet_index, global_facet_vertex in global_facet2coord.items():
                    local_index = self.global_to_local_facet_map[subdomain_index][global_facet_index]
                    local_facet_vertex = subdom_facet2coord[local_index]
                    vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
                    vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
                    both_vertices_are_equal = (vertices_are_equal == 2)
                    if both_vertices_are_equal:
                        print(f"facet with global number {global_facet_index} has correctly been identyfied to local index {local_index}")
                    else:
                        print("FALSE FALSE FALSE")
                        print(f"facet with global number {global_facet_index} has INCORRECTLY been identyfied to local index {local_index}")
    
    
        def init_interface_values(self,#
                                interface_marker: df.MeshFunction,#
                                interface_marker_value: int,
                                function_space: tp.Dict[str, tp.Dict[str, df.FunctionSpace]],
                                subdomain_index: int,
                                has_phases: tp.List[str],
                                ) -> None:
            """ allocate dictionaries that are used to store the interface dof values
                of the pressure, previous pressure, gl decoupling term as well as previous
                gl update term respectively as pairs(parent_mesh_index: dof_value)
    
            # Parameters
            interface_marker:       #type df.MeshFunction   marker function marking
                                                            edges of i_mesh
                                                            according to self.inside
            interface_marker_value: #type int               marker value of interface_marker
            """
    
            function_name = ["pressure", "gli"]
            for function in function_name:
                if function == "pressure":
                    self.pressure_values.update({subdomain_index: dict()})
                    self.pressure_values_prev.update({subdomain_index: dict()})
                else:  #function == "gli":
                    self.gli_term.update({subdomain_index: dict()})
                    self.gli_term_prev.update({subdomain_index: dict()})
    
                for phase in has_phases:
                    if function == "pressure":
                        self.pressure_values[subdomain_index].update({phase: dict()})
                        self.pressure_values_prev[subdomain_index].update({phase: dict()})
                    else:  #function == "gli":
                        self.gli_term[subdomain_index].update({phase: dict()})
                        self.gli_term_prev[subdomain_index].update({phase: dict()})
    
                    space = function_space[function][phase]
                    dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
                    # the dof dictionaries are grouped by the global (with respect to the
                    # global mesh) facet index of each facet along the interface.
                    # each dof is then saved with its koordinates as tuples as keys.
                    for subdom_mesh_facet_index, facet_dof_coord_dict in dof_coordinates.items():
                        global_mesh_facet_index = self.local_to_global_facet_map[subdomain_index][subdom_mesh_facet_index]
                        if function == "pressure":
                            self.pressure_values[subdomain_index][phase].update({global_mesh_facet_index: dict()})
                            self.pressure_values_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
                        else:  #function == "gli":
                            self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
                            self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
    
                        for dof_coord in facet_dof_coord_dict.values():
                            # we want to use the coordinates as key in a dictionary.
                            # Because dof_coord is a numpy array and numpy array are
                            # not hashable, we convert the numpy array into a tuple,
                            # which in turn can then be used as a dict key.
                            coord_tupel = (dof_coord[0], dof_coord[1])
                            if function == "pressure":
                                self.pressure_values[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                                self.pressure_values_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                            else:  #function == "gli":
                                self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                                self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                        # end init_interface_values
    
    
    
        #### Private Methods ####
    
        def _vertex_indices(self, #
                        interface_marker: df.MeshFunction, #
                        interface_marker_value: int,#
                        print_vertex_indices: bool = False) -> np.array:
            """ return indices of nodes on the mesh marked by
                interface_marker with value interface_marker.
    
            _vertex_indices() returns indices of nodes on
            the mesh marked by interface_marker with value
            interface_marker. interface_marker is a dolfin.MeshFunction of dimension
            interface_marker.mesh().topology().dim()-1 that is supposed to mark edges of i_mesh
            according to the method super().inside().
    
            # Parameters
            interface_marker:       #type df.MeshFunction   marker function marking edges of i_mesh
                                                            according to self.inside
            interface_marker_value: #type int               marker value of interface_marker
            print_vertex_indices      #type boolean           specify wether to print out the
                                                            node or not.
            """
            mesh_coordinates = interface_marker.mesh().coordinates()
            mesh_dimension = interface_marker.mesh().topology().dim()
            # interface_marker is a function defined on edges. Therefore
            # sum(interface_marker.array() == interface_marker_value) gives the number
            # of edges on the boundary. We need the number of nodes, however.
            number_of_interface_vertices = sum(interface_marker.array() == interface_marker_value) + 1
    
            # print(f"interface marker array",interface_marker.array() == interface_marker_value)
            # print(f"facets marked by interface marker", interface_marker.array())
            # print(f"interface{self.global_index} has coordinates {self.coordinates(interface_marker, interface_marker_value)}")
    
            # print(f"\nDetermined number of interface vertices as {number_of_interface_vertices}")
            # we need one mesh_dimension + 1 columns to store the the index of the node.
            vertex_indices = np.zeros(shape = number_of_interface_vertices, dtype=int)
            # print(f"allocated array for vertex_indices\n", vertex_indices) #
            interface_vertex_number = 0
            # mesh_vertex_index = 0
            # print(f"\n we are one interface{self.global_index}",
            #       f" we determined {number_of_interface_vertices} interface vertices.")
            for vert_num, x in enumerate(mesh_coordinates):
                if self._is_on_boundary_part(x):
                    # print(f"Vertex {x} with index {vert_num} is on interface")
                    # print(f"interface_vertex_number = {interface_vertex_number}")
                    # print(f"dfPoint = ({x}) is on interface{self.global_index}")
                    vertex_indices[interface_vertex_number] = vert_num
                    interface_vertex_number += 1
    
            if print_vertex_indices:
                print(f"mesh node indices: \n", vertex_indices)
    
            return vertex_indices
    ### End Class Interface