diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index eb5033a72a7fba2650452cc497e05ded5fda163b..9d3a9f3988d2dbddb9e245416f905be2ace2ce7a 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -241,6 +241,8 @@ class Interface(BoundaryPart):
# global index ob the subdomain
self.global_index = global_index
self.marker_value = self.global_index+1
+ # dictionary holding the facet objects
+ self.facets = dict()
# dictionary containing for each facet belonging to the interface (with
# respect to a subdomain mesh) the coordinates of the the vertices
# belonging to that facet. This dictionary is created by
@@ -420,6 +422,7 @@ class Interface(BoundaryPart):
"""
self.facet_to_vertex_coordinates.update({subdomain_index: dict()})
self.outer_normals.update({subdomain_index: dict()})
+ self.facets.update({subdomain_index: dict()})
# Get a subset iterator for mesh entities along the facets marked by
# interface_marker and value interface_marker_value.
mesh_entity_iterator = df.SubsetIterator(interface_marker, interface_marker_value)
@@ -442,6 +445,7 @@ class Interface(BoundaryPart):
self.facet_to_vertex_coordinates[subdomain_index].update(
{facet.index(): facet_vertices}
)
+ self.facets[subdomain_index].update({facet.index(): facet})
if subdomain_index == 0:
self.outer_normals[subdomain_index].update(
{facet.index(): None}
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index bbf83b8dae8c6387b8cf497f4c756106c5150f50..07f5f0e32daf3e158f99b1128736d486443ff59e 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -452,7 +452,6 @@ class DomainPatch(df.SubDomain):
# the wetting phase is always present and we always need
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
- # TODO: add gravity term!!!!!
if self.include_gravity:
# z_a included the minus sign already, see self._calc_gravity_expressions
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
@@ -465,202 +464,53 @@ 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)
- # we need to treat the flux part of gl0 and the robin (pressure)
- # part of gl0 seperately at the intersection points of two
- # interfaces (more details in comment below). Therefor we save
- # the robin pressure part and the neumann flux part separately
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()
- flux_x_coeficients = flux_x.vector()
- flux_y_coeficients = flux_y.vector()
- # print(f"flux_x coefficients: ", flux_x_coef[:])
- flux_x_dof_map = flux_x.function_space().dofmap()
- flux_y_dof_map = flux_y.function_space().dofmap()
- flux_dofmap = self.dofmap["flux"][phase]
- # flux_dofmap1 = flux_dofmap[0]
- # flux_dofmap2 = flux_dofmap[0]
-
- gli_dofmap = self.dofmap["gli"][phase]
gl0 = df.Function(self.function_space["gli"][phase])
- boundary_facets_and_cells = dict()
- interface_entity_iterator = df.SubsetIterator(self.interface_marker, interface.marker_value)
- for mesh_entity in interface_entity_iterator:
- # facet_index = mesh_entity.global_index()
- # global_facet_index = mesh_entity.global_index()
- # print(f"facet_index = {facet_index}\n",
- # f"global_facet_index = {global_facet_index}\n")
- for f in df.facets(mesh_entity):
- facet = f
- # facet_index = f.global_index()
- # global_facet_index = f.global_index()
- # print(f"facet_index = {facet_index}\n",
- # f"global_facet_index = {global_facet_index}\n")
-
- facet_cell_iterator = df.cells(facet)
- # only one cell of the mesh contains facet, but we still
- # have to loop, because df.cells(facet) returns an iterator
- # for cells.
- for facet_cell in facet_cell_iterator:
- boundary_facets_and_cells.update({facet: facet_cell})
- # for facet in df.facets(self.mesh):
- # if facet.exterior():
- # boundary_facets.update({facet: facet.normal()})
-
- flux_x_element = flux_x.function_space().element()
- flux_y_element = flux_y.function_space().element()
- gli_element = gl0.function_space().element()
- p_element = self.function_space["pressure"][phase].element()
- for facet, facet_cell in boundary_facets_and_cells.items():
- facet_cell_index = facet_cell.index()
- facet_index = facet.index()
- global_facet_cell_index = facet_cell.global_index()
-
- print(f"facet_cell_index = {facet_cell_index}\n",
- f"global_facet_cell_index = {global_facet_cell_index}")
- # since the flux is constant on each cell, we can take no
- # matter which dof to calculate F*n.
- n = facet.normal()[:]
- print(flux_x_element.tabulate_dof_coordinates(facet_cell))
- flux_x_facet_cell_dof_coordinates = flux_x_element.tabulate_dof_coordinates(facet_cell)
- flux_y_facet_cell_dof_coordinates = flux_y_element.tabulate_dof_coordinates(facet_cell)
- gli_facet_cell_dof_coordinates = gli_element.tabulate_dof_coordinates(facet_cell)
- p_facet_cell_dof_coordinates = p_element.tabulate_dof_coordinates(facet_cell)
-
- gli_cell_dofmap = gli_dofmap.cell_dofs(facet_cell_index)
- flux_x_cell_dofmap = flux_x_dof_map.cell_dofs(facet_cell_index)
- flux_dofmap1_cell = flux_dofmap.cell_dofs(facet_cell_index)
-
- flux_y_cell_dofmap = flux_y_dof_map.cell_dofs(facet_cell_index)
- print("flux_dofmap: ", flux_dofmap1_cell)
- print("flux_x_dofmap: ", flux_x_cell_dofmap)
- print("flux_y_dofmap: ", flux_y_cell_dofmap)
- p_cell_dofmap = p_dofmap.cell_dofs(facet_cell_index)
- print(f"gli_cell_dofmap = {gli_cell_dofmap}")
-
- facet_vertex_coordinates = []
- for vertex in df.vertices(facet):
- # we don't need the 3d coordinates only x and y.
- vertex_coord = vertex.point()[:][0:2]
- facet_vertex_coordinates.append(vertex_coord)
- print(f"type of vertex.point()[:]", type(vertex_coord))
- print('vertex is on interface:', interface._is_on_boundary_part(vertex_coord))
- print(vertex_coord)
-
- p_dofs_along_facet = dict()
- for local_dof_ind, dof_coord in enumerate(p_facet_cell_dof_coordinates):
- if interface._is_on_boundary_part(dof_coord):
- p_dofs_along_facet.update({local_dof_ind: dof_coord})
- print(p_dofs_along_facet.items())
-
- flux_x_dofs_along_facet = dict()
- for local_dof_ind, dof_coord in enumerate(flux_x_facet_cell_dof_coordinates):
- if interface._is_on_boundary_part(dof_coord):
- flux_x_dofs_along_facet.update({local_dof_ind: dof_coord})
- print(flux_x_dofs_along_facet.items())
-
- flux_y_dofs_along_facet = dict()
- for local_dof_ind, dof_coord in enumerate(flux_y_facet_cell_dof_coordinates):
- if interface._is_on_boundary_part(dof_coord):
- flux_y_dofs_along_facet.update({local_dof_ind: dof_coord})
- print(flux_y_dofs_along_facet.items())
-
- gli_dofs_along_facet = dict()
- for local_dof_ind, dof_coord in enumerate(gli_facet_cell_dof_coordinates):
- if interface._is_on_boundary_part(dof_coord):
- gli_dofs_along_facet.update({local_dof_ind: dof_coord})
- print(gli_dofs_along_facet.items())
-
- for local_dof_ind, gli_dof_coord in gli_dofs_along_facet.items():
- flux_x_dof_ind = -1
- for local_ind, coord in flux_x_dofs_along_facet.items():
- if np.array_equal(coord, gli_dof_coord):
- flux_x_dof_ind = flux_x_cell_dofmap[local_ind]
-
- flux_y_dof_ind = -1
- for local_ind, coord in flux_y_dofs_along_facet.items():
- if np.array_equal(coord, gli_dof_coord):
- flux_y_dof_ind = flux_y_cell_dofmap[local_ind]
-
- p_dof_ind = -1
- for local_ind, coord in p_dofs_along_facet.items():
- if np.array_equal(coord, gli_dof_coord):
- p_dof_ind = p_cell_dofmap[local_ind]
-
-
- print(f"local_dof_ind = {local_dof_ind} and type is {type(local_dof_ind)}")
- print(f"cell_dof_map = {gli_cell_dofmap}")
- gli_dof_ind = gli_cell_dofmap[local_dof_ind]
- print(f"local_dof_ind = {local_dof_ind} and global ind is {gli_dof_ind}")
- print(f"nomal vector components = [{n[0]}, {n[1]}]")
- print(f"p_dof_ind = {p_dof_ind}")
- gl0.vector()[gli_dof_ind] = flux_x.vector()[flux_x_dof_ind]*n[0] + flux_y.vector()[flux_y_dof_ind]*n[1] - robin_pressure[p_dof_ind]
- print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_ind]}")
-
-
-
- # print(cell_dofs_flux_x)
- # cell_dofs_flux_y = flux_y_dof_map.cell_dofs(facet_cell_index)
- # cell_dofs_gli = gli_dofmap.cell_dofs(facet_cell_index)
- # flux_x_coef = flux_x.vector()
- # # print(f"flux_x coefficients: ", flux_x_coef[:])
- # flux_space_dof_map = flux_x.function_space().dofmap()
- # print("local to global dofs:", flux_space.dofmap().tabulate_local_to_global_dofs())
- # print("dofmap index map:", flux_space.dofmap().index_map())
- # print("entity dofs:", flux_space.dofmap().entity_dofs(self.mesh, 2))
- # print("block size of dofmap:", flux_space.dofmap().block_size())
- # interface_facets = df.SubsetIterator(self.interface_marker, interface.marker_value)
- # for facet in interface_facets:
- # print(facet)
- # print(facet.entities(0))
- # how_many_cells = 0
- # for cell in df.cells(facet):
- # how_many_cells += 1
- # print(f"Cell containing facet {facet} has dofs", flux_space_dof_map.cell_dofs(cell.index()))
- # print("entity dofs of facet:", flux_space_dof_map.entity_dofs(self.mesh, 2, [facet.entities(0)[1]]))
- # # print("entity dofs of facet:", flux_space_dof_map.entity_dofs(self.mesh, 2, [facet.entities(0)[1]]))
- # print(f"the facet {facet.index()} has {how_many_cells} cells")
- # for vertex in df.vertices(facet):
- # print(vertex.point()[:])
- # print(f"dofmap flux_space:", flux_space.dofmap().dofs())
- # print(f"cell dof flux_space:", flux_space.dofmap().cell_dofs(0))
-
- # print(flux_space.tabulate_dof_coordinates())
- # print("Mesh Cells that contain a boundary facet (alt)")
- # for mesh_cell in df.cells(self.mesh):
- # for facet in df.facets(mesh_cell):
- # if facet.exterior():
- # print(facet.entities(0))
- # for vertex in df.vertices(facet):
- # print(vertex.point()[:])
- # print(mesh_cell, mesh_cell.get_vertex_coordinates())
-
-
-
-
- # flux_function = df.project(flux_x, W)
- # Project the normal flux onto V so that we can
- # gl0 = df.project(df.dot(flux_function,n) - robin_pressure, gli_space)
- ds = self.ds
- marker_value = interface.marker_value
- phi_a = self.testfunction[phase]
- gl0_form = gl0*phi_a*ds(marker_value)
- gl0_assembled = df.assemble(gl0_form)
- # # flux_x, flux_y = flux_FEM.split(deepcopy=True)
- # normal_flux = normal_flux_FEM.vector()
- # interface_dofs = self._dof_indices_of_interface[interf_ind][phase]
- # parent = self.parent_mesh_index
- # d2v = self.dof2vertex[phase]
- # for dof_index in interface_dofs:
- # gli = normal_flux[dof_index] - robin_pressure[dof_index]
- # # 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.
- # interface.gli_term_prev[subdomain_ind][phase].update(
- # {parent[d2v[dof_index]]: gli}
- # )
-
+ gli_dof_coordinates = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+ p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+ flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ for local_facet_ind, normal in interface.outer_normals[subdomain].items():
+ gli_dofs_along_facet = gli_dof_coordinates[local_facet_ind]
+ p_dofs_along_facet = p_dof_coordinates[local_facet_ind]
+ flux_x_dofs_along_facet = flux_dof_coordinates["x"][local_facet_ind]
+ flux_y_dofs_along_facet = flux_dof_coordinates["y"][local_facet_ind]
+
+ # for each gli dof with index gli_dof_ind, we need the
+ # dof indices for the flux and the pressure corresponding
+ # to the coordinates of the dof, in order to calculate gl0
+ # dof wise.
+ for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
+ flux_x_dof_index = -1
+ for flux_x_dof_ind, dof_coord in flux_x_dofs_along_facet.items():
+ if np.array_equal(dof_coord, gli_dof_coord):
+ flux_x_dof_index = flux_x_dof_ind
+
+ flux_y_dof_index = -1
+ for flux_y_dof_ind, dof_coord in flux_y_dofs_along_facet.items():
+ if np.array_equal(dof_coord, gli_dof_coord):
+ flux_y_dof_index = flux_y_dof_ind
+
+ p_dof_index = -1
+ for p_dof_ind, dof_coord in p_dofs_along_facet.items():
+ if np.array_equal(dof_coord, gli_dof_coord):
+ 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]}")
+
+ # save this newly calculated dof to the interface
+ # dictionary gli_term_prev
+ global_facet_ind = local2global_facet[local_facet_ind]
+ dof_coord_tupel = (gli_dof_coord[0], gli_dof_coord[1])
+ interface.gli_term_prev[subdomain][phase][global_facet_ind].update(
+ {dof_coord_tupel: gl0.vector()[gli_dof_index]}
+ )
### END calc_gl0_term
@@ -681,7 +531,7 @@ class DomainPatch(df.SubDomain):
subdomain = self.subdomain_index
interface = self.interface[interface_index]
# gli should be initialised by zero.
- gli = df.Function(self.function_space["pressure"][phase])
+ 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
@@ -689,19 +539,6 @@ class DomainPatch(df.SubDomain):
interface.current_iteration[subdomain] = iteration
# save the iteration number of the neighbour
neighbour_iter_num = interface.current_iteration[neighbour]
- # needed for the read_gli_dofs() functions
- interface_dofs = self._dof_indices_of_interface[interface_index][phase]
- # for each interface, the following 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[interface_index]
- if debug:
- print(f"On subdomain{subdomain}, we have:\n",
- f"Interface{interface_index}, Neighbour index = {neighbour}\n",
- f"dofs on interface:\n{interface_dofs['wetting']}\n",
- f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
- f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
# read previous neighbouring pressure values from interface.
@@ -745,32 +582,51 @@ class DomainPatch(df.SubDomain):
"and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
+ "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
- parent = self.parent_mesh_index
- d2v = self.dof2vertex[phase]
- # read gli_prev term from the neighbour.
- gli_prev = interface.gli_term_prev[neighbour][phase]
dt = self.timestep_size
Lambda = self.lambda_param[phase]
- for dof_index in interface_dofs:
- glob_vert_ind = parent[d2v[dof_index]]
- gli_value = -2*Lambda*pressure_prev[glob_vert_ind] - gli_prev[glob_vert_ind]
- gli.vector()[dof_index] = gli_value
- gli_save_dictionary[glob_vert_ind] = 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,
- # if we don't save it to gli_term_prev of the neighbour.
- # similarly for the pressure values.
- if iteration == neighbour_iter_num:
- # 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].update(#
- {phase: interface.gli_term[neighbour][phase].copy()}
+
+ 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]
+ # 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():
+ # find the right previous pressure dof of the interface
+ # dictionary to use for the calculation.
+ p_previous_dof = -9999999.9999
+ for p_dof_coord_tupel, p_prev_dof in p_prev.items():
+ p_coord = np.array([p_dof_coord_tupel[0], p_dof_coord_tupel[1]])
+ if np.array_equal(gli_dof_coord, p_coord):
+ p_previous_dof = p_prev_dof
+ # same with the gli_previous_dof
+ gli_previous_dof = -9999999.9999
+ for gli_prev_coord_tupel, gli_prev_dof in gli_prev.items():
+ gli_prev_coord = np.array([gli_prev_coord_tupel[0], gli_prev_coord_tupel[1]])
+ if np.array_equal(gli_dof_coord, gli_prev_coord):
+ gli_previous_dof = gli_prev_dof
+
+ 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
+
+ # In the following case, the gli_term of the neighbour is saved to
+ # gli_term currently and will be overwritten in the next step,
+ # if we don't save it to gli_term_prev of the neighbour.
+ # similarly for the pressure values.
+ if iteration == neighbour_iter_num:
+ # 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()}
+ )
+ interface.pressure_values_prev[neighbour][phase].update(
+ {global_facet_ind: interface.pressure_values[neighbour][phase][global_facet_ind].copy()}
)
- interface.pressure_values_prev[neighbour][phase].update(
- {phase: interface.pressure_values[neighbour][phase].copy()}
- )
return gli
### END calc_gli_term
@@ -850,7 +706,13 @@ class DomainPatch(df.SubDomain):
{phase: self.function_space["gli"][phase].dofmap()}
)
self.dofmap["flux"].update(
- {phase: self.function_space["flux"][phase].dofmap()}
+ {phase: dict()}
+ )
+ self.dofmap["flux"][phase].update(
+ {"x": self.function_space["flux"][phase].sub(0).dofmap()}
+ )
+ self.dofmap["flux"][phase].update(
+ {"y": self.function_space["flux"][phase].sub(1).dofmap()}
)
def _init_markers(self) -> None:
@@ -932,9 +794,20 @@ class DomainPatch(df.SubDomain):
{interface_index: dict()}
)
for phase in self.has_phases:
- self.interface_dof_indices_and_coordinates[function][interface_index].update(#
- {phase: interface.dofs_and_coordinates(function_space[phase], marker, marker_value)}
- )
+ if function == "flux":
+ self.interface_dof_indices_and_coordinates[function][interface_index].update(#
+ {phase: dict()}
+ )
+ self.interface_dof_indices_and_coordinates[function][interface_index][phase].update(#
+ {"x": interface.dofs_and_coordinates(function_space[phase].sub(0), marker, marker_value)}
+ )
+ self.interface_dof_indices_and_coordinates[function][interface_index][phase].update(#
+ {"y": interface.dofs_and_coordinates(function_space[phase].sub(1), marker, marker_value)}
+ )
+ else:
+ self.interface_dof_indices_and_coordinates[function][interface_index].update(#
+ {phase: interface.dofs_and_coordinates(function_space[phase], marker, marker_value)}
+ )
def _init_interface_values(self):