Select Git revision
domainPatch.py
domainPatch.py 70.09 KiB
"""Class definitions for subdomains representing soil patches for LDD-TPR simulation
This file contains class definitions for the class DomainPatch which defines the
subdomains used for the LDD simulation. Each subdomain of this class holds
information about permeability, porosity, the model/physics used on the patch
as well as information about which parts of the boundary are outer boundary of
the global domain and which parts are interfaces to neighboring soil patches of
the same class.
Additionally, a class BoundaryPart is defined, which accepts a List of dolfin
Points vertices, representing the vertices of a 1D polygone. The class contains
the method BoundaryPart.inside(meshpoint,on_boundary) which is used to mark
faces of a given subdomain mesh that lie on the 1D-polygone defined by vertices.
Assumed is the use of python 3.6
"""
import dolfin as df
import mshr
import numpy as np
# import domainPatch as dp
import typing as tp
import ufl as fl
# Interface = tp.NewType('Interface', tp.any)
#import domainPatch as domainPatch
class DomainPatch(df.SubDomain):
""" DomainPatch Class defines subdomains used in the LDD simulation
The class DomainPatch defines subdomains for the LDD simulation. Each subdomain
of this class holds information about the phyisics involved like viscosity,
porosity, the model/physics that is assumed on the patch
as well as information about which parts of the boundary are outer boundary of
the global domain and which parts are interfaces to neighboring soil patches of
the same class.
The constructor accepts
# Parameters
isRichards #type: bool return true if Richards model is
assumed on patch.
mesh #type df::Mesh & Dolfin mesh object, submesh of the
subdomain.
porosity #type float porosity of the the soil in the patch
assumed to be constant in space.
viscosity #type tp.List[float] viscosity(ies) of the the
fluid phases in the patch viscosity[0] is assumed
to be the viscosity of the wetting phase and if
present (isRichards == False) viscosity[1] is
assumed to be the viscosity of the non-wetting phase.
interfaces: #type tp.List[dp.Interface] List of Interface class
objects from the LDDsimulation
class.
has_interface #type tp.List[int]: list of indices (w.r.t. the global
numbering of interfaces in
LDDsimulation.interfaces) indicating
which interfaces are part of the
subdomain.
L #type tp.List[float] L parameters for the L-scheme. In
case isRichards == False (i.e
two-phase flow is assumed on the
patch) we have two parameters, one
for each equation.
lambda_param #tp.List[float] L parameters for the L-scheme. In
case isRichards == False (i.e
two-phase flow is assumed on the
patch) we have two parameters, one
for each equation.
## Public Variables
isRichards #type bool set by input, see above
mesh #type df.Mesh set by input, see above
porosity #type float set by input, see above
viscosity #type tp.List[float] set by input, see above
interfaces: #type tp.List[dp.Interface] set by input, see above
has_interface #type tp.List[int] set by input, see above
L #type tp.List[float] set by input, see above
lambda_param #type tp.List[float] set by input, see above
function_space #type tp.Dict[str: df.Function] function space
used for wetting
and nonwetting
pressures. This is
set by
self._init_function_space()
parent_mesh_index #type np.array parent vertex map of the submesh.
boundary_marker #type df.MeshFunction marker function to mark all
boundaries of the subdomain.
The marker value marking
interface i with global index
has_interface[i] will be
has_interface[i] aswell.
## Public Methods
"""
### constructor
def __init__(self, #
subdomain_index: int,#
isRichards: bool,#
mesh: df.Mesh,#
porosity: float,#
viscosity: tp.Dict[str, float],#
outer_boundary_def_points: tp.Dict[str, tp.List[df.Point]],#
# tp.List[Interface] doesn't work, because it hasen't been defined yet.
interfaces: tp.List[object],#
has_interface: tp.List[int],#
L: tp.Dict[str, float],#
lambda_param: tp.Dict[str, float],#
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
timestep_size: float,#
tol = 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.
df.SubDomain.__init__(self)
if tol:
self.tol = tol
else:
self.tol = df.DOLFIN_EPS
### Class Variables/Objects set by the input to the constructor
self.subdomain_index = subdomain_index
self.isRichards = isRichards
self.mesh = mesh
self.porosity = porosity
self.viscosity = viscosity
self.outer_boundary_def_points = outer_boundary_def_points
self.interface = interfaces
self.has_interface = has_interface
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
# timestep size, tau in the paper
self.timestep_size = timestep_size
### Class variables set by methods
# dictionary holding the dof indices corresponding to an interface of
# given interface. see self._calc_dof_indices_of_interfaces()
self._dof_indices_of_interface = dict()
# dictionary holding for each interface in self.has_interface the dof
# indices that are shared with another interface (possibly none). This info is needed for
# the calculation of the gli terms.
self._interface_has_common_dof_indices = dict()
# List of objects of clas outerBoundary initialised by self._init_markers()
# can be NONE if no outer boundary is present.
self.outer_boundary = []
# measures are set by self._init_measures()
self.iteration_number = 0
self.dx = None
self.ds = None
# solution function(s) for the form set by _init_function_space
self.pressure = None
self.source = None
# FEM function holding the previous iteration.
self.pressure_prev_iter = None
# FEM function holding the pressures of the previous timestep t_n.
self.pressure_prev_timestep = None
# test function(s) for the form set by _init_function_space
self.testfunction = None
# Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
# corresponding to the dofs of an FEM function which is initiated
# by calc_gli_term() holding the dofs of the assembled dt*gli*v*ds terms as
# a sparse vector (here, dt = self.timestep_size). This vector is read by the solver and added to the rhs
# which is given by self.govering_problem().
#
self.gli_assembled = dict()
# Dictionary of FEM function of self.function_space holding the interface
# dofs of the previous iteration of the neighbouring pressure. This is used
# to assemble the gli terms in the method self.calc_gli_term().
self.neighbouring_interface_pressure = None
# sometimes we need to loop over the phases and the length of the loop is
# determined by the model we are assuming
if self.isRichards:
self.has_phases =['wetting']
else:
self.has_phases =['wetting', 'nonwetting']
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_markers()
self._init_measures()
self._calc_dof_indices_of_interfaces()
self._calc_interface_has_common_dof_indices()
# END constructor
#### PUBLIC METHODS
def write_pressure_to_interfaces(self, initial_values: bool = False):
""" save the interface values of self.pressure to all neighbouring interfaces
This method is used by the LDDsimulation.LDDsolver method.
# parameters
initial_values bool Flag to toggle wether to write self.pressure
or self.pressure_prev_iter to the interface.
This is needed for the first ever iteration
at time = 0.
"""
if initial_values:
# this assumes, that self.pressure_prev_iter has been set by
# LDDsimulation._init_initial_values().
from_func = self.pressure_prev_iter
else:
from_func = self.pressure
for ind in self.has_interface:
for phase in self.has_phases:
self.interface[ind].read_pressure_dofs(from_function = from_func[phase], #
interface_dofs = self._dof_indices_of_interface[ind][phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = phase,#
subdomain_ind = self.subdomain_index)
def governing_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
""" return the governing form and right hand side for phase phase as a dictionary.
return the governing form ant right hand side for phase phase (without the gli terms) as a dictionary.
CAUTION: the gli terms are not assembled in this method and need to be
calculated by the solver separately and subtracted from the rhs!
"""
# define measures
dx = self.dx
ds = self.ds
dt = self.timestep_size
porosity = self.porosity
if self.isRichards:
# this assumes that atmospheric pressure has been normalised to 0.
pc_prev_iter = self.pressure_prev_iter['wetting']
pc_prev_timestep = self.pressure_prev_timestep['wetting']
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
"Parameter phase = 'nonwetting'. But isRichards = True for",
"current subdomain. \nResetting phase = 'wetting'.\n")
phase = 'wetting'
else:
pw_prev_iter = self.pressure_prev_iter['wetting']
pnw_prev_iter = self.pressure_prev_iter['nonwetting']
pw_prev_timestep = self.pressure_prev_timestep['wetting']
pnw_prev_timestep = self.pressure_prev_timestep['nonwetting']
pc_prev_iter = pnw_prev_iter - pw_prev_iter
pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
La = self.L[phase]
pa = self.trialfunction[phase]
# this is pw_{i-1} in the paper
pa_prev_iter = self.pressure_prev_iter[phase]
phi_a = self.testfunction[phase]
Lambda_a = self.lambda_param[phase]
ka = self.relative_permeability[phase]
S = self.saturation
mu_a = self.viscosity[phase]
source_a = self.source[phase]
# the first part of the form and rhs are the same for both wetting and nonwetting
form1 = (La*pa*phi_a)*dx
rhs1 = (La*pa_prev_iter*phi_a)*dx
# the forms of wetting and nonwetting phase differ by a minus sign
# and by interface_form terms if the neighbour of the current patch
# assumes Richards model
if phase == "wetting":
# gather interface forms
interface_forms = []
for interface in self.has_interface:
# interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
# form2 is different for wetting and nonwetting
form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
elif phase == "nonwetting":
# gather interface forms
interface_forms = []
for interface in self.has_interface:
if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we assume constant normalised atmospheric pressure
# and no interface term is practically appearing.
# interface_forms += (df.Constant(0)*phi_a)*ds(interface)
interface_forms.append((df.Constant(0)*phi_a)*ds(interface))
else:
# interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(interface))
form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
else:
raise RuntimeError('missmatch for input parameter phase.',
'Expected either phase == wetting or phase == nonwetting ')
return None
form = form1 + form2 + sum(interface_forms)
rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
'form': form,#
'rhs_without_gli' : rhs_without_gli
}
return form_and_rhs
def calc_gli_term(self) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl terms for the iteration number of the subdomain
calculate the gl terms for the iteration number of the subdomain and save them
to self.gli_assembled. The term gets saved as an already assembled sparse
ds form to be added to the rhs by the solver. This is due to the need of
having to communicate dofs over the interface. Additionally, the dofs are
being saved to an interface dictionary exactly as the pressures.
Caution: The array self.gli_assembled needs to be added to the rhs with
a minus sign by the solver!
"""
iteration = self.iteration_number
subdomain = self.subdomain_index
dt = self.timestep_size
Lambda = self.lambda_param
# since we loop over all interfaces in self.has_interface we need a temporary
# vector (numpy array) that we will use to add all dofs of the gli terms
# into one array.
gli_assembled_tmp = dict()
if self.isRichards:
# print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
gli_assembled_tmp = {
'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float)
}
else:
gli_assembled_tmp = {
'wetting': np.zeros(self.pressure['wetting'].vector()[:].size, dtype=float),
'nonwetting': np.zeros(self.pressure['nonwetting'].vector()[:].size, dtype=float)
}
for ind in self.has_interface:
interface = self.interface[ind]
# neighbour index
neighbour = interface.neighbour[subdomain]
neighbour_iter_num = interface.current_iteration[neighbour]
ds = self.ds(ind)
# needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind]
# for each interface, this is a list of dofs indices belonging to another
# interface. The dofs corresponding to these indices must be multiplied
# by 1/2 to to get an average of the gli terms for dofs belonging to
# two different interfaces.
dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[ind]
if iteration == 0:
n = df.FacetNormal(self.mesh)
pw_prev_iter = self.pressure_prev_timestep
k = self.relative_permeability
S = self.saturation
mu = self.viscosity
if not neighbour_iter_num == 0:
# save the gl term from the neighbour first to use it in the
# next iteration
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
# in case we have two-phase and the neighbour is two-phase, two,
# we need to Additionally save the nonwetting values.
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
if self.isRichards:
# assuming that we normalised atmospheric pressure to zero,
# pc = pw in the Richards case.
pc_prev = pw_prev_iter['wetting']
else:
pw_prev = pw_prev_iter['wetting']
pnw_prev = pw_prev_iter['nonwetting']
pc_prev = pnw_prev - pw_prev
for phase in self.has_phases:
# viscosity
mu_a = mu[phase]
# relative permeability
ka = k[phase]
# previous pressure iteration
pa_prev = pw_prev_iter[phase]
# testfunction
phi_a = self.testfunction[phase]
if phase == 'wetting':
# the wetting phase is always present and we always need
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
# TODO: add gravitiy term!!!!!
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
gl0 = (df.dot(flux, n) - Lambda[phase]*pa_prev)
gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
# if dofs_in_common_with_other_interfaces[phase].size == 0
# there are no dofs that lie on another interface and we can
# skip the following step.
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
# from the standpoint of the subdomain we are currently on,
# each dof can belong maximally to two interfaces. On these
# vertices, common to two interfaces, the dofs of the
# corresponding gli terms of both interfaces are averaged
# to give one value for this vertex. Because we achieve this
# by adding the assembled gli terms together we have to
# multiply the dofs that appear in more than one interface
# by 1/2.
gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
gli_assembled_tmp[phase] += gl0_assembled
else:
# phase == 'nonwetting'
# if we are in the two-phase case but our neighbour assumes
# Richards equation, we hase a zero communication term.
if interface.neighbour_isRichards[subdomain]:
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we don't have gli terms for the
# nonwetting phase and assemble the terms as zero form.
gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
else:
# in this case the neighbour assumes two-phase flow
# and we need to calculte a gl0 torm for the nonwetting
# phase.
# we need anothre flux in this case
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
gl0 = (df.dot(flux, n) - Lambda[phase]*pa_prev)
gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
# see previous case.
gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
gli_assembled_tmp[phase] += gl0_assembled
# writing out glo to the current iteration dictionary of the
# interface takes place after the loop, see below.
# END for loop over phases
interface.current_iteration[subdomain] += 1
else: # iteration > 0
# the wetting phase update is always to be calculated even if
# only Richards is assumed on the current patch or the neighbour
# only assumes the Richards equation.
# in case neighbour_iter_num == iteration we need to first save
# the gli_term of the neighbour to gli_term_prev.
if neighbour_iter_num == iteration:
# in this case, first save the gli_terms of the neighbour to
# gli_prev_term of the current subdomain. This is because
# we assume, that the neighbouring patch has done a previous
# iteration and we need the value to calculate the gli term with
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
for phase in self.has_phases:
# in case phase == 'nonwetting' only proceed if the neighbouring
# subdomain does not assume a Richards model.
if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
# then read gli_prev term from interface and write it to
# the gli_assembled array of the subdomain
# to be able to sum the gli terms of all interfaces together,
# we need a dummy vector to hold the dofs read from the interface.
gli_readout_tmp = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
interface.write_gli_dofs(to_array = gli_readout_tmp, #
interface_dofs = interface_dofs[phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = phase,#
subdomain_ind = subdomain,#
# this is set to True because we need gli_prev
previous_iter = True
)
# read neighbouring pressure values from interface and write them to
# the function self.neighbouring_interface_pressure
interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure[phase], #
interface_dofs = interface_dofs[phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = phase,#
subdomain_ind = neighbour,#
previous_iter = False)
# use the above to calculate the gli update with
# in the paper p^{i-1}_{3-l}
pa_prev = self.neighbouring_interface_pressure[phase]
phi_a = self.testfunction[phase]
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
gli_p_part = gli_pressure_part.get_local()
gli_prev_neighbour = gli_readout_tmp
gli = gli_p_part + gli_prev_neighbour
# check if there are shared dofs with another interface.
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
# see previous case.
gli[common_dof] = 0.5*gli[common_dof]
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
gli_assembled_tmp[phase] += gli
if neighbour_iter_num -1 == iteration:
# gli_term_prev has been used by the above and can now safely
# be overwritten with the gli_term of the neighbour.
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
else:
print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
f"and iteration = {iteration} on subdomain {subdomain}\n",
"It is expected that this case should only be happening if the \n",
f"neighbour (index {neighbour}) has already finished its iteration, i.e. reached\n",
f"the error criterion. In this case we should have {iteration} > {neighbour_iter_num}.\n",
"If this is not the case revisite DomainPatch.calc_gli_term().")
# finally we are done and can step up the current iteration count.
interface.current_iteration[subdomain] += 1
### END for ind in self.has_interface:
# save the result of the calculation to the subdomain and to the interface
# dictionaries.
if self.isRichards:
self.gli_assembled.update(
{'wetting': gli_assembled_tmp['wetting']}
)
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
# write gli dofs to interface dicts for communiction
interface.read_gli_dofs(from_array = self.gli_assembled['wetting'], #
interface_dofs = interface_dofs['wetting'],#
dof_to_vert_map = self.dof2vertex['wetting'],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = 'wetting',#
subdomain_ind = subdomain,#
previous_iter = False
)
else:
self.gli_assembled.update({
'wetting': gli_assembled_tmp['wetting'],
'nonwetting': gli_assembled_tmp['nonwetting']
})
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
for phase in ['wetting', 'nonwetting']:
# write gli dofs to interface dicts for communiction
interface.read_gli_dofs(from_array = self.gli_assembled[phase], #
interface_dofs = interface_dofs[phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = phase,#
subdomain_ind = subdomain,#
previous_iter = False
)
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None:
""" reset interface.current_iteration[subdomain] = 0 for all interfaces
of the subdomain.
"""
subdomain = self.subdomain_index
for index in self.has_interface:
self.interface[index].current_iteration[subdomain] = 0
#### PRIVATE METHODS
def _init_function_space(self) -> None:
""" create function space for solution and trial functions
Note that P1 FEM is hard coded here, as it is assumed in other methods
aswell. This method sets the class variable
self.function_space
self.pressure
self.testfunction
"""
if self.isRichards:
self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
self.trialfunction = {
'wetting': df.TrialFunction(self.function_space['wetting'])
}
self.testfunction = {
'wetting': df.TestFunction(self.function_space['wetting'])
}
self.neighbouring_interface_pressure = {
'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
}
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
}
self.trialfunction = {
'wetting' : df.TrialFunction(self.function_space['wetting']),
'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
}
self.testfunction = {
'wetting' : df.TestFunction(self.function_space['wetting']),
'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
}
self.neighbouring_interface_pressure = {#
'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
}
def _init_dof_and_vertex_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
This method sets the class Variables
self.parent_mesh_index
self.dof2vertex
self.vertex2dof
"""
mesh_data = self.mesh.data()
# local to parent vertex index map
self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
if self.isRichards:
self.dof2vertex = {#
'wetting' :df.dof_to_vertex_map(self.function_space['wetting'])#
}
self.vertex2dof = {#
'wetting' :df.vertex_to_dof_map(self.function_space['wetting'])#
}
else:
self.dof2vertex = {#
'wetting' :df.dof_to_vertex_map(self.function_space['wetting']),#
'nonwetting':df.dof_to_vertex_map(self.function_space['nonwetting'])#
}
self.vertex2dof = {#
'wetting' :df.vertex_to_dof_map(self.function_space['wetting']),#
'nonwetting':df.vertex_to_dof_map(self.function_space['nonwetting'])#
}
def _init_markers(self) -> None:
""" define boundary markers
This method sets the class Variables
self.interface_marker
self.outer_boundary_marker
"""
self.outer_boundary_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
self.outer_boundary_marker.set_all(0)
self.interface_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
self.interface_marker.set_all(0)
for glob_index in self.has_interface:
# each interface gets marked with the global interface index.
self.interface[glob_index].mark(self.interface_marker, glob_index+1)
# create outerBoundary objects and mark outer boundaries with 1
for index, boundary_points in self.outer_boundary_def_points.items():
# if our domain patch has no outer boundary, this should be reflected
# by the value None in the dictionary
if boundary_points is None:
#
self.outer_boundary = None
else:
self.outer_boundary.append(#
BoundaryPart(vertices = boundary_points,#
internal = False,
tol= self.tol) #self.mesh.hmin()/100
)
self.outer_boundary[index].mark(self.outer_boundary_marker, 1)
def _init_measures(self):
""" define measures for the governing form
This method sets the class Variables
self.dx
self.ds
"""
# domain measure
self.dx = df.Measure('dx', domain = self.mesh)
# measure over the interfaces
self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.interface_marker)
def _calc_dof_indices_of_interfaces(self):
""" calculate dof indices of each interface """
V = self.function_space
marker = self.interface_marker
for ind in self.has_interface:
# the marker value for the interfaces is always global interface index+1
marker_value = ind + 1
self._dof_indices_of_interface.update(#
{ind: dict()}
)
for phase in self.has_phases:
self._dof_indices_of_interface[ind].update(#
{phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value)}
)
print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind][phase])
# if self.isRichards:
# self._dof_indices_of_interface[ind].update(#
# {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value)}
# )
# print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind]['wetting'])
# else:
# self._dof_indices_of_interface[ind].update(#
# {'wetting': self.interface[ind].dofs_on_interface(V['wetting'], marker, marker_value),#
# 'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, marker_value)}
# )
def _calc_interface_has_common_dof_indices(self):
""" populate dictionary self._interface_has_common_dof_indices
Determin for each interface in self.has_interface the dof indices that
also belong to another interface and store them in the dictionary
self._interface_has_common_dof_indices
This method populates the dictionary
self._interface_has_common_dof_indices
"""
for interf_ind, interf_dofs in self._dof_indices_of_interface.items():
# initialise bool numpy arry for the search.
if self.isRichards:
is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
else:
is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
is_part_of_neighbour_interface_nw = np.zeros(interf_dofs['nonwetting'].size, dtype=bool)
for other_interf_ind, other_interf_dofs in self._dof_indices_of_interface.items():
if other_interf_ind is not interf_ind:
if self.isRichards:
is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
else:
is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
is_part_of_neighbour_interface_nw += np.isin(interf_dofs['nonwetting'], other_interf_dofs['nonwetting'], assume_unique=True)
# this carries a numpyarray of the indices which belong to other
# interfaces, too
self._interface_has_common_dof_indices.update({interf_ind: dict()})
if self.isRichards:
self._interface_has_common_dof_indices[interf_ind].update(
{'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
)
else:
self._interface_has_common_dof_indices[interf_ind].update(
{'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]},
{'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
)
# END OF CLASS DomainPatch
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.
"""
def __init__(self, #
vertices: tp.List[df.Point], #
tol = None, #
internal: bool = False):
# 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
###### 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])
# print(f"test if point {p} is on line segment between {p1} or {p2}")
# equality check
if np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol and np.absolute((p[1] - ymax)*(p[1] - ymin)) < tol:
#print(f"point {p} is close to either {p1} or {p2}")
return True
# check there holds p1[0] < p[0] < p2[0]. If not, p cannot be on the line segment
# 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 pp1 or pp2 is vertical. We need to still calculate
# the distance to the line segment in this case
if ((p[0] - xmax)*(p[0] - xmin) < tol) or np.absolute((p[0] - xmax)*(p[0] - xmin)) < tol:
# 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, #
domain_index: int = None,#
adjacent_subdomains: np.array = None,#
# this is an array of bools
isRichards: np.array = None,#
**kwds):
#
self.isRichards = isRichards
self.adjacent_subdomains = adjacent_subdomains
# 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]#
}
# make the parent class methods available.
super().__init__(**kwds)
#### 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)
dof2vertex = df.dof_to_vertex_map(function_space)
# 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)
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 indicis 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("have written pressure interface values\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("have written previous pressure interface values\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("have written gl interface values\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("have written previous gl interface values\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"\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
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}")
vertex_indices[interface_vertex_number] = vert_num
interface_vertex_number += 1
if print_vertex_indices:
print("mesh node indices: \n", vertex_indices)
return vertex_indices
### End Class Interface