Select Git revision
boundary_and_interface.py 33.31 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
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
# 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 = None
self.pressure_values_prev = None
# dictionaries holding the forms for the decoupling terms
self.gli_term = None
self.gli_term_prev = None
# 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 init_interface_values(self,#
interface_marker: df.MeshFunction,#
interface_marker_value: int) -> 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)
The MeshFunction interface_marker is supposed to be one defined on the
global mesh, so that the indices used in the dictionary are the global
vertex indices and not indices of some submesh.
# 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
"""
vertex_indices = self._vertex_indices(interface_marker, #
interface_marker_value,#
print_vertex_indices = False)
self.pressure_values = dict()
#self.pressure_values is a dictionay which contains for each index in
# adjacent_subdomains a dictionary for wetting and nonwetting phase with
# as a value yet another dictionary holding the interface values. phew.
self.pressure_values.update({self.adjacent_subdomains[0]: dict()})
self.pressure_values.update({self.adjacent_subdomains[1]: dict()})
self.pressure_values_prev = dict()
self.pressure_values_prev.update({self.adjacent_subdomains[0]: dict()})
self.pressure_values_prev.update({self.adjacent_subdomains[1]: dict()})
self.gli_term = dict()
self.gli_term.update({self.adjacent_subdomains[0]: dict()})
self.gli_term.update({self.adjacent_subdomains[1]: dict()})
self.gli_term_prev = dict()
self.gli_term_prev.update({self.adjacent_subdomains[0]: dict()})
self.gli_term_prev.update({self.adjacent_subdomains[1]: dict()})
# we initialise all dictionaries to accomodate two phases, to be more
# flexible with coupling Richards and Two-phase models.
for subdom_ind in self.adjacent_subdomains:
for phase in ['wetting', 'nonwetting']:
self.pressure_values[subdom_ind].update({phase: dict()})
self.pressure_values_prev[subdom_ind].update({phase: dict()})
# the gli_term will hold dofs of the gli terms for communications
# they are initiated with 0
self.gli_term[subdom_ind].update({phase: dict()})
self.gli_term_prev[subdom_ind].update({phase: dict()})
for vert_ind in vertex_indices:
self.pressure_values[subdom_ind][phase].update({vert_ind: 0.0})
self.pressure_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
self.gli_term[subdom_ind][phase].update({vert_ind: 0.0})
self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
# end init_interface_values
def read_pressure_dofs(self, from_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
phase: str,#
subdomain_ind: int,#
previous_iter: bool = False,#
) -> None:
""" read dofs of from_function with indices interface_dofs and write to self.pressure_values
Read the degrees of freedom of the function from_function with indices
interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
of the interface
self.pressure_values[subdomain_ind][phase], (if pressure == True)
self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
# Parameters
from_function: #type: df.Function, function to save the dofs of
interface_dofs: #type: np.array, index array of dofs to be saved
dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or
'nonwetting' and determins
the dict to be written to.
subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous
iterate is being written or not
"""
if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
d2v = dof_to_vert_map
if not previous_iter:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
# print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
else:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
# print(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\n", self.pressure_values_prev[subdomain_ind][phase])
def read_gli_dofs(self, from_array: np.array, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
phase: str,#
subdomain_ind: int,#
previous_iter: bool = False,#
) -> None:
""" Read dofs of from_array with indices interface_dofs and write to self.gli_term
Read the degrees of freedom of the array from_array with indices
interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
of the interface
self.gli_term[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
# Parameters
from_array: #type: np.array, array to save the dofs of
interface_dofs: #type: np.array, index array of dofs to be saved
dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or
'nonwetting' and determins
the dict to be written to.
subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous
iterate is being written or not
"""
if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
d2v = dof_to_vert_map
if not previous_iter:
for dof_index in interface_dofs:
# from_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
# print(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\n", self.gli_term[subdomain_ind][phase])
else:
for dof_index in interface_dofs:
# from_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
# print(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\n", self.gli_term_prev[subdomain_ind][phase])
def write_pressure_dofs(self, to_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
phase: str,#
subdomain_ind: int,#
# pressure: bool = False,#
# gl: bool = False,#
previous_iter: bool = False,#
) -> None:
""" read dofs off self.pressure_values and overwrite the dofs of to_function
with indices interface_dofs
Read the degrees of freedom of one of the dictionaries of the interface
self.pressure_values[subdomain_ind][phase], (if pressure == True)
self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
and overwrite the dofs of the function to_function with indices
interface_dofs, i.e. update the interface dofs of the function.
.
# Parameters
to_function: #type: df.Function, function to overwrite the
dofs of.
interface_dofs: #type: np.array, index array of interface
dofs of from_function
dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or
'nonwetting' and determins
the dict to be read of.
subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous
iterate is being read or not
"""
if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
d2v = dof_to_vert_map
if not previous_iter:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_function.vector()[dof_index] = self.pressure_values[subdomain_ind][phase][parent[d2v[dof_index]]]
# print("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
# print("wrote these dofs to function \n", to_function.vector()[:])
else:
for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_function.vector()[dof_index] = self.pressure_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
# print("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
# print("wrote these dofs to function \n", to_function.vector()[:])
def write_gli_dofs(self, to_array: np.array, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
phase: str,#
subdomain_ind: int,#
previous_iter: bool = False,#
) -> None:
""" read dofs off self.gli_term and overwrite the dofs stored in the array
to_array with indices interface_dofs
Read the degrees of freedom of one of the dictionaries of the interface
self.gli_term[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
and overwrite the values (holding dofs) of the array to_array with indices
interface_dofs, i.e. update the interface dofs of the gli terms.
.
# Parameters
to_array: #type: np.array, array to overwrite the
dofs of.
interface_dofs: #type: np.array, index array of interface
dofs of from_function
dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh.
local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or
'nonwetting' and determins
the dict to be read of.
subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous
iterate is being read or not
"""
if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values")
parent = local_to_parent_vertex_map
d2v = dof_to_vert_map
if not previous_iter:
for dof_index in interface_dofs:
# to_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_array[dof_index] = self.gli_term[subdomain_ind][phase][parent[d2v[dof_index]]]
# print("read gl interface values\n", self.gli_term[subdomain_ind][phase])
# print("wrote these dofs to array \n", to_array[:])
else:
for dof_index in interface_dofs:
# to_array is orderd by dof and not by vertex numbering of the mesh.
# parent needs a mesh vertex index of the submesh, therefore d2v is needed.
to_array[dof_index] = self.gli_term_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
# print("read gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
# print("wrote these dofs to array \n", to_array[:])
#### 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)}")
# for cell in interface_marker[interface_marker.array() == interface_marker_value]:
# print(cell.get_cell_data())
# 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