Select Git revision
boundary_and_interface.py
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