diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 9af6a2a580659deb5c6a82860824107d261a8b7e..33901e69e77ee5f0b8abf806063a4c298a5a9c3c 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -116,17 +116,18 @@ class LDDsimulation(object):
self._initialised = False
# dictionary of objects of class DomainPatch initialised by self._init_subdomains()
self.subdomain = dict()
- # dictionary to hold the input Expressions for the external boundary. These are
+ # dictionary to hold the input Expression strings for the external boundary. These are
# turned into df.Expressions objects by the method update_DirichletBC_dictionary()
# this is mostly in place when examples with exact solutions are calculated.
self.dirichletBC_expression_strings = None
# dictionary to hold the df.Expressions for the external boundary created
# from the above self.dirichletBC_expression_strings. These are
# turned into df.dirichletBC objects by the method update_DirichletBC_dictionary()
- self.dirichletBC_dfExpression = dict()
+ # The dictionary is created by self._init_DirichletBC_dictionary
+ self.dirichletBC_dfExpression = None
# dictionary to store the df.dirichletBC objects in. These are used by the
- # solver to set the boundary conditions
- self.outerBC = dict()
+ # solver to set the boundary conditions. Gets initiated by self._init_DirichletBC_dictionary
+ self.outerBC = None
# Dictionary that holdes the exact solution. Gets set by self.set_parameters()
# but can remain None if we don't consider a case with an exact solution.
self.exact_solution = None
@@ -262,6 +263,7 @@ class LDDsimulation(object):
# initialise exact solution expression dictionary if we have an exact
# solution case.
self._init_exact_solution_expression()
+ self._init_DirichletBC_dictionary()
if self.exact_solution:
df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow"))
self.write_exact_solution_to_xdmf()
@@ -355,9 +357,11 @@ class LDDsimulation(object):
# r_cons.close()
# energy.close()
df.list_timings(df.TimingClear.keep, [df.TimingType.wall, df.TimingType.system])
- with open(self.output_dir+"timings.txt", "w") as timefile:
- timefile.write(df.timings)
-
+ # timing_table = df.timings(df.TimingClear.keep, [df.TimingType.wall, df.TimingType.user, df.TimingType.system])
+ # print(timing_table.__get_value__())
+ # with open(self.output_dir+"timings.txt", "w") as timefile:
+ # timefile.write(df.timings)
+
df.dump_timings_to_xml(self.output_dir+"timings.xml",df.TimingClear.keep)
@@ -534,9 +538,10 @@ class LDDsimulation(object):
# reset all interface[has_interface].current_iteration[ind] = 0, for
# index has_interface in subdomain.has_interface
subdomain.null_all_interface_iteration_numbers()
- self.update_DirichletBC_dictionary(time = time, #
- subdomain_index = ind,
- )
+ if subdomain.outer_boundary is not None:
+ self.update_DirichletBC_dictionary(time = time, #
+ subdomain_index = ind,
+ )
# # this is the beginning of the solver, so iteration must be set
# # to 0.
subdomain.iteration_number = 0
@@ -587,7 +592,6 @@ class LDDsimulation(object):
subdomain.calc_gli_term()
for phase in subdomain.has_phases:
- # if not calc_isdone_for_phase[phase]:
# extract L-scheme form and rhs (without gli term) from subdomain.
governing_problem = subdomain.governing_problem(phase = phase)
form = governing_problem['form']
@@ -602,6 +606,8 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36:
print("\nSystem before applying outer boundary conditions:")
print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+ # we need the rhs as a dolfin function, so we overwrite
+ # rhs_without_gli_assembled with the gli_term subtracted from itself.
rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
rhs_assembled = rhs_without_gli_assembled
@@ -886,20 +892,21 @@ class LDDsimulation(object):
# define global interface marker. This is a mesh function on edges of the mesh.
self.interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
self.interface_marker.set_all(0)
- for num, vertices in enumerate(interface_def_points):
+ for glob_interface_num, vertices in enumerate(interface_def_points):
self.interface.append(bi.Interface(vertices=vertices, #
tol=self.tol,#mesh.hmin()/100,
internal=True,#
- adjacent_subdomains = adjacent_subdomains[num],#
+ adjacent_subdomains = adjacent_subdomains[glob_interface_num],#
isRichards = self.isRichards,#
- global_index = num,
+ global_index = glob_interface_num,
)
)#
# marker value gets internally set to global_index+1.
# The reason is that we want to mark outer boundaries with the value 0.
- self.interface[num].mark(self.interface_marker, self.interface[num].marker_value)
- self.interface[num].init_interface_values(self.interface_marker, self.interface[num].marker_value)
- # print(self.interface[num])
+ marker_value = self.interface[glob_interface_num].marker_value
+ self.interface[glob_interface_num].mark(self.interface_marker, marker_value)
+ self.interface[glob_interface_num].init_interface_values(self.interface_marker, marker_value)
+ # print(self.interface[glob_interface_num])
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"interface_marker_global.pvd"
@@ -957,6 +964,12 @@ class LDDsimulation(object):
# for parameter, value in self.linear_solver.parameters.items():
# print(f"parameter: {parameter} = {value}")
# df.info(solver_param, True)
+ if self.write2file['meshes_and_markers']:
+ filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}.pvd"
+ df.File(filepath) << self.subdomain[subdom_num].interface_marker
+ filepath = self.output_dir+f"outer_boundary_marker_subdomain{subdom_num}.pvd"
+ df.File(filepath) << self.subdomain[subdom_num].outer_boundary_marker
+
print(self.subdomain[subdom_num])
@@ -1030,9 +1043,62 @@ class LDDsimulation(object):
write_iter_for_fixed_time = False,
)
+
+ def _init_DirichletBC_dictionary(self, interpolation_degree: int = None):
+ """initialise dictionary that holds the dolfin.DirichletBC objects
+ for each outer boundary of each subdomain.
+
+ We create two dictionaries, self.outerBC and self.dirichletBC_dfExpression.
+ The second one holds the dolfin expressions that is used to create the
+ dirichlet boundary object of class dolfin.dirichletBC as well as update
+ the time of the expression. The second one holds the dolfin.dirichletBC
+ boundary object itself.
+ """
+ if interpolation_degree is None:
+ interpolation_degree = self.interpolation_degree
+
+ # we only need to set Dirichlet boundary conditions if we have an outer
+ # boundary.
+ self.outerBC = dict()
+ self.dirichletBC_dfExpression = dict()
+ for subdom_index, subdomain in self.subdomain.items():
+ if subdomain.outer_boundary is None:
+ self.outerBC.update({subdom_index: None})
+ self.dirichletBC_dfExpression.update({subdom_index: None})
+ else:
+ self.outerBC.update({subdom_index: dict()})
+ self.dirichletBC_dfExpression.update({subdom_index: dict()})
+ V = subdomain.function_space
+ boundary_marker = subdomain.outer_boundary_marker
+ # loop over all outer boundary parts. This is a bit of a brainfuck:
+ # subdomain.outer_boundary gives a dictionary of the outer boundaries of the subdomain.
+ # Since a domain patch can have several disjoint outer boundary parts, the expressions
+ # need to get an enumaration index which starts at 0. So subdomain.outer_boundary[j] is
+ # the dictionary of outer dirichlet conditions of the subdomain and boundary part j for
+ # the wetting and the nonwetting phase, so subdomain.outer_boundary[j]['wetting'] and
+ # subdomain.outer_boundary[j]['nonwetting'] return
+ # the actual expression needed for the dirichlet condition for both phases if present.
+ for boundary_index, boundary in enumerate(subdomain.outer_boundary):
+ # create a dictionary of dictionaries for each outer boundary
+ # part of the subdomain.
+ self.outerBC[subdom_index].update({boundary_index: dict()})
+ self.dirichletBC_dfExpression[subdom_index].update({boundary_index: dict()})
+ # this contains the Expression strings input by the user.
+ pDC = self.dirichletBC_expression_strings[subdom_index][boundary_index]
+ boundary_marker_value = boundary.marker_value
+ for phase in subdomain.has_phases:
+ pDCa = df.Expression(pDC[phase], domain = subdomain.mesh, #
+ degree = interpolation_degree, #
+ t=self.starttime)
+ self.dirichletBC_dfExpression[subdom_index][boundary_index].update({phase: pDCa})
+ self.outerBC[subdom_index][boundary_index].update(
+ {phase: df.DirichletBC(V[phase], pDCa, boundary_marker, boundary_marker_value)}
+ )
+
+
def update_DirichletBC_dictionary(self, subdomain_index: int,#
interpolation_degree: int = None, #
- time: float = 0,#
+ time: float = None,#
):
""" update time of the dirichlet boundary object for subdomain with
index subdomain_index.
@@ -1043,56 +1109,25 @@ class LDDsimulation(object):
if interpolation_degree is None:
interpolation_degree = self.interpolation_degree
- num = subdomain_index
- subdomain = self.subdomain[num]
- # we only need to set Dirichlet boundary conditions if we have an outer
- # boundary.
- if subdomain.outer_boundary is not None:
+ if time is None:
+ time = self.starttime
+ subdomain = self.subdomain[subdomain_index]
+ # in case time == self.starttime, the expression has already been
+ # assembled by self._init_DirichletBC_dictionary.
+ if np.fabs(time - self.starttime) < self.tol:
+ pass
+ else:
V = subdomain.function_space
- # loop over all outer boundary parts. This is a bit of a brainfuck:
- # subdomain.outer_boundary gives a dictionary of the outer boundaries of the subdomain.
- # Since a domain patch can have several disjoint outer boundary parts, the expressions
- # need to get an enumaration index which starts at 0. So subdomain.outer_boundary[j] is
- # the dictionary of outer dirichlet conditions of the subdomain and boundary part j for
- # the wetting and the nonwetting phase, so subdomain.outer_boundary[j]['wetting'] and
- # subdomain.outer_boundary[j]['nonwetting'] return
- # the actual expression needed for the dirichlet condition for both phases if present.
boundary_marker = subdomain.outer_boundary_marker
- for index, boundary in enumerate(subdomain.outer_boundary):
- if np.abs(time-self.starttime) < self.tol:
- # Here the dictionary has to be created first because t = 0.
- # create a dictionary of dictionaries for each subdomain
- self.outerBC.update({num: dict()})
- self.dirichletBC_dfExpression.update({num: dict()})
- # create a dictionary of dictionaries for each outer boundary
- # part of the subdomain.
- self.outerBC[num].update({index: dict()})
- self.dirichletBC_dfExpression[num].update({index: dict()})
- # this contains the Expression strings input by the user.
- pDC = self.dirichletBC_expression_strings[num][index]
- # as with the interfaces, since numbering of the outer boundary
- # parts of each subdomain starts at 0 and subdomain.outer_boundary_marker
- # gets set to 0 on all facets of the mesh, we need to up the value for
- # an actual outer boundary part by 1.
- boundary_marker_value = index+1
+ for boundary_index, boundary in enumerate(subdomain.outer_boundary):
+ boundary_marker_value = boundary.marker_value
for phase in subdomain.has_phases:
- # note that the default interpolation degree is 2
- if np.abs(time - self.starttime) < self.tol :
- # time = 0 so the Dolfin Expression has to be created first.
- pDCa = df.Expression(pDC[phase], domain = subdomain.mesh, #
- degree = interpolation_degree, #
- t=self.starttime)
- self.dirichletBC_dfExpression[num][index].update({phase: pDCa})
- else:
- pDCa = self.dirichletBC_dfExpression[num][index][phase]
- pDCa.t = time
- self.outerBC[num][index].update(
+ pDCa = self.dirichletBC_dfExpression[subdomain_index][boundary_index][phase]
+ pDCa.t = time
+ self.outerBC[subdomain_index][boundary_index].update(
{phase: df.DirichletBC(V[phase], pDCa, boundary_marker, boundary_marker_value)}
)
- else:
- # case subdomain.outer_boundary is None
- print("subdomain is an inner subdomain and has no outer boundary.\n",
- "no Dirichlet boundary has been set.")
+
def _init_exact_solution_expression(self):
""" initialise dolfin expressions for exact solution and populate
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index af44457fafcc05725abb97ee2a92aa559a3b2f55..76b78db1b475e02ccbe6b088112653785c5ca96f 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -32,12 +32,15 @@ class BoundaryPart(df.SubDomain):
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):
+ 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
@@ -51,6 +54,10 @@ class BoundaryPart(df.SubDomain):
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:
@@ -228,6 +235,10 @@ class Interface(BoundaryPart):
# 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
@@ -263,10 +274,6 @@ class Interface(BoundaryPart):
self.adjacent_subdomains[1]: self.adjacent_subdomains[0]#
}
-
- # make the parent class methods available.
- super().__init__(**kwds)
-
# MAGIC METHODS
def __str__(self):
"""print out a few usefull informations of the interface"""
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 94fa2ca690dc12e7adf69f2614886643bf8ef8b8..5bc83173312fabc14e89d3208c218d9424408d28 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -786,24 +786,27 @@ class DomainPatch(df.SubDomain):
# each interface gets marked with the global interface index.
self.interface[glob_index].mark(self.interface_marker, marker_value)
# 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:
+ # in case there is no outer boundary, self.outer_boundary_def_points is
+ # None.
+ if self.outer_boundary_def_points is None:
+ self.outer_boundary = None
+ else:
+ for index, boundary_points in self.outer_boundary_def_points.items():
self.outer_boundary.append(#
- bi.BoundaryPart(vertices = boundary_points,#
- internal = False,
- tol= self.tol) #self.mesh.hmin()/100
+ bi.BoundaryPart(vertices=boundary_points,#
+ internal=False,
+ tol=self.tol,
+ marker_value=index+1
+ ) #self.mesh.hmin()/100
)
# similar to the interfaces self.outer_boundary_def_points is a
# dictionary enumerating all the boundary parts. The enumeration
# starts at 0, so we set the marker value to index+1
- boundary_marker_value = index+1
- self.outer_boundary[index].mark(self.outer_boundary_marker, boundary_marker_value)
-
+ print(f"on subdomain{self.subdomain_index} I have boundary part {index}")
+ boundary_marker_value = self.outer_boundary[index].marker_value
+ self.outer_boundary[index].mark(
+ self.outer_boundary_marker, boundary_marker_value
+ )
def _init_measures(self):
""" define measures for the governing form