diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index dc8eb64ee13a46358a858bac9d289c1a827c3f08..3d1716be7acb17fa0b930e6e3133e0d5e2e4505b 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -1131,6 +1131,7 @@ class LDDsimulation(object):
else:
pass
+
def _init_sources(self, interpolation_degree: int = None):
""" create source expressions at self.starttime and populate subdomain
source dictionaries.
@@ -1148,6 +1149,7 @@ class LDDsimulation(object):
t = self.starttime)}
)
+
def _add_parameters_to_output_dirname(self):
""" function to create a new output directory for given set of input
parameters and adjust self.output_dir accordingly
@@ -1163,6 +1165,7 @@ class LDDsimulation(object):
os.mkdir(self.output_dir)
print(f"\nwill save output of this simulation to path {self.output_dir}")
+
def _init_analyse_timesteps(self):
""" determin timesteps to anayise
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 507e2b116741cc140f076d52dc732d3ec5c9f2fe..46af28ac033f73133cbb9195ec9a63cff7fe3abc 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -298,34 +298,34 @@ class Interface(BoundaryPart):
# def set_domain_index(self, domain_index: int):
# self.domain_index = domain_index
- def dofs_on_interface(self,#
- function_space: df.FunctionSpace,#
- interface_marker: df.MeshFunction, #
- interface_marker_value: int) -> np.array:
- """ return the indices of the degrees of freedom on the interface.
-
- dofs_on_interface() returns the indices of the degrees of freedom on the
- interface with respect to the dof numbering on function_space.mesh().
- In essence, it is the restriction of the dof_to_vertex_map to the interface.
-
- # Parameters
- function_space: #type df.functionSpace the function space of whitch
- to calculate the dof_to_vertex_map
- restriction from.
- interface_marker: #type df.MeshFunction marker function marking edges of i_mesh
- according to self.inside
- interface_marker_value: #type int marker value of interface_marker
- """
- interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value,
- print_vertex_indices = False)
- # print(f"local interface vertices: {interface_vertex_indices}")
- dof2vertex = df.dof_to_vertex_map(function_space)
- # print(f"dof2vertex map: {dof2vertex}\n")
- # check which dofs actually lie on the interface using the dof2vertex map.
- is_a_dof_on_interface = np.isin(dof2vertex, interface_vertex_indices, assume_unique=True)
- # print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
- # print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
- return np.nonzero(is_a_dof_on_interface)[0]
+ # def dofs_on_interface(self,#
+ # function_space: df.FunctionSpace,#
+ # interface_marker: df.MeshFunction, #
+ # interface_marker_value: int) -> np.array:
+ # """ return the indices of the degrees of freedom on the interface.
+ #
+ # dofs_on_interface() returns the indices of the degrees of freedom on the
+ # interface with respect to the dof numbering on function_space.mesh().
+ # In essence, it is the restriction of the dof_to_vertex_map to the interface.
+ #
+ # # Parameters
+ # function_space: #type df.functionSpace the function space of whitch
+ # to calculate the dof_to_vertex_map
+ # restriction from.
+ # interface_marker: #type df.MeshFunction marker function marking edges of i_mesh
+ # according to self.inside
+ # interface_marker_value: #type int marker value of interface_marker
+ # """
+ # interface_vertex_indices = self._vertex_indices(interface_marker, interface_marker_value,
+ # print_vertex_indices = False)
+ # # print(f"local interface vertices: {interface_vertex_indices}")
+ # dof2vertex = df.dof_to_vertex_map(function_space)
+ # # print(f"dof2vertex map: {dof2vertex}\n")
+ # # check which dofs actually lie on the interface using the dof2vertex map.
+ # is_a_dof_on_interface = np.isin(dof2vertex, interface_vertex_indices, assume_unique=True)
+ # # print(f"is_a_dof_on_interface: {is_a_dof_on_interface}")
+ # # print(f"this means the following dofs are interface dofs:", np.nonzero(is_a_dof_on_interface)[0])
+ # return np.nonzero(is_a_dof_on_interface)[0]
def dofs_and_coordinates(self, function_space: df.FunctionSpace,#
@@ -363,6 +363,7 @@ class Interface(BoundaryPart):
# this is the facet index relative to the Mesh
# that interface_marker is defined on.
facet_index = mesh_entity.index()
+
dofs_and_coords.update({facet_index: dict()})
# because mesh_entity doesn't have a method to access coordinates we
# need to cast it into an acutal facet object.
@@ -447,6 +448,8 @@ class Interface(BoundaryPart):
)
self.facets[subdomain_index].update({facet.index(): facet})
if subdomain_index == 0:
+ # we do not have outer normals if this interface is part of the
+ # global domain (subdomain_index = 0).
self.outer_normals[subdomain_index].update(
{facet.index(): None}
)
@@ -454,13 +457,17 @@ class Interface(BoundaryPart):
self.outer_normals[subdomain_index].update(
{facet.index(): facet.normal()[:][0:2]}
)
+ # if subdomain_index is not equal to 0, we are on an actual SubDomain
+ # and self.init_facet_dictionaries() has been run for the global domain
+ # already. This means we actually can compute a local2global facet map.
+ self.init_facet_maps(subdomain_index)
if debug:
print(f"newly defined facet_to_vertex_coordinates_dictionary[{subdomain_index}]")
print(self.facet_to_vertex_coordinates[subdomain_index].items())
print(f"outer normals per facet are: ", self.outer_normals[subdomain_index].items())
- def init_facet_maps(self, subdomain_index: int, debug: bool = False) -> None:
+ def init_facet_maps(self, subdomain_index: int, debug: bool = True) -> None:
""" calculate the mapping that maps the local numbering of interface
facets (w.r.t. the mesh of subdomain subdomain_index) to the global
index with respect to the global domain (subdomain 0)"""
@@ -479,7 +486,7 @@ class Interface(BoundaryPart):
vertices_are_equal += np.array_equal(local_facet_vertex[0], global_facet_vertex[1])
vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[0])
vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[1])
- both_vertices_are_equal = vertices_are_equal == 2
+ both_vertices_are_equal = (vertices_are_equal == 2)
if both_vertices_are_equal:
self.local_to_global_facet_map[subdomain_index].update(
{local_facet_index: global_facet_index}
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index f70f5a5ec4257dd834bff8f590df5dc5b310fcf1..bff23f732bfa5cee2b10b135de22d94cbe1fbe9f 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -278,13 +278,16 @@ class DomainPatch(df.SubDomain):
for phase in self.has_phases:
p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][ind][phase]
local2global_facet = interface.local_to_global_facet_map[subdomain]
- for facet_ind in interface.facets[subdomain].keys():
- global_facet_ind = local2global_facet[facet_ind]
+ # for facet_ind in interface.facets[subdomain].keys():
+ for subdomain_facet_index, facet_dof_coord_dict in p_dof_coordinates.items():
+ global_facet_index = local2global_facet[subdomain_facet_index]
+ # set dictionary to save the values to.
if previous_iter:
- save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_ind]
+ save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_index]
else:
- save_dict = interface.pressure_values[subdomain][phase][global_facet_ind]
- for dof_index, dof_coord in p_dof_coordinates[facet_ind].items():
+ save_dict = interface.pressure_values[subdomain][phase][global_facet_index]
+
+ for dof_index, dof_coord in facet_dof_coord_dict.items():
dof_coord_tupel = (dof_coord[0], dof_coord[1])
save_dict.update({dof_coord_tupel: self.pressure[phase].vector()[dof_index]})
@@ -464,7 +467,7 @@ class DomainPatch(df.SubDomain):
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
else:
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
- robin_pressure = Lambda[phase]*pa_prev.vector() #.vector()
+
flux_function = df.project(flux, self.function_space["flux"][phase])
flux_x, flux_y = flux_function.split()
gl0 = df.Function(self.function_space["gli"][phase])
@@ -501,8 +504,10 @@ class DomainPatch(df.SubDomain):
p_dof_index = p_dof_ind
gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- + flux_y.vector()[flux_y_dof_index]*normal[1] - robin_pressure[p_dof_index]
- print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
+ + flux_y.vector()[flux_y_dof_index]*normal[1] \
+ - Lambda[phase]*pa_prev.vector()[p_dof_index]
+ if debug:
+ print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
# save this newly calculated dof to the interface
# dictionary gli_term_prev
@@ -532,6 +537,7 @@ class DomainPatch(df.SubDomain):
interface = self.interface[interface_index]
# gli should be initialised by zero.
gli = df.Function(self.function_space["gli"][phase])
+
# neighbour index
neighbour = interface.neighbour[subdomain]
# update the current_iteration number of subdomain stored in the interface
@@ -588,13 +594,13 @@ class DomainPatch(df.SubDomain):
gli_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
p_interface_dofs = self.interface_dof_indices_and_coordinates["pressure"][interface_index][phase]
local2global_facet = interface.local_to_global_facet_map[subdomain]
- for facet_ind in interface.facets[subdomain].keys():
- global_facet_ind = local2global_facet[facet_ind]
+ for local_facet_index, gli_facet_dofs2coordinates in gli_interface_dofs.items():
+ global_facet_index = local2global_facet[local_facet_index]
# read gli_prev term from the neighbour.
- gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_ind]
- p_prev = pressure_prev[global_facet_ind]
- gli_save_dict = gli_save_dictionary[global_facet_ind]
- for gli_dof_index, gli_dof_coord in gli_interface_dofs[facet_ind].items():
+ gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_index]
+ p_prev = pressure_prev[global_facet_index]
+ gli_save_dict = gli_save_dictionary[global_facet_index]
+ for gli_dof_index, gli_dof_coord in gli_facet_dofs2coordinates.items():
# find the right previous pressure dof of the interface
# dictionary to use for the calculation.
p_previous_dof = -9999999.9999
@@ -612,7 +618,7 @@ class DomainPatch(df.SubDomain):
gli_value = -2*Lambda*p_previous_dof - gli_previous_dof
gli.vector()[gli_dof_index] = gli_value
gli_dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
- gli_save_dictionary[gli_dof_coord_tupel] = gli_value
+ gli_save_dict.update({gli_dof_coord_tupel: gli_value})
# In the following case, the gli_term of the neighbour is saved to
# gli_term currently and will be overwritten in the next step,
@@ -622,10 +628,10 @@ class DomainPatch(df.SubDomain):
# print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
# f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
interface.gli_term_prev[neighbour][phase].update(#
- {global_facet_ind: interface.gli_term[neighbour][phase][global_facet_ind].copy()}
+ {global_facet_index: interface.gli_term[neighbour][phase][global_facet_index].copy()}
)
interface.pressure_values_prev[neighbour][phase].update(
- {global_facet_ind: interface.pressure_values[neighbour][phase][global_facet_ind].copy()}
+ {global_facet_index: interface.pressure_values[neighbour][phase][global_facet_index].copy()}
)
return gli
@@ -738,7 +744,6 @@ class DomainPatch(df.SubDomain):
interface_marker=self.interface_marker,
interface_marker_value=marker_value,
subdomain_index=self.subdomain_index)
- self.interface[glob_index].init_facet_maps(self.subdomain_index)
# create outerBoundary objects and mark outer boundaries with 1
# in case there is no outer boundary, self.outer_boundary_def_points is
# None.
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index aa8aff2e62bd67427fe0c260b9192b28e2a9d514..3a1af669afbced9982d570056c8d6fb7bd1a884d 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -89,10 +89,10 @@ isRichards = {
mesh_resolution = 4
##################### TIME #####################################################
timestep_size = 0.01
-number_of_timesteps = 100
+number_of_timesteps = 20
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 11
+number_of_timesteps_to_analyse = 0
starttime = 0