From e9dadcd1fdb71eb4f28eef45f8a44fa0b2e52e52 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 14 Jun 2019 16:06:02 +0200
Subject: [PATCH] rewrite init_interface_values
---
LDDsimulation/LDDsimulation.py | 17 +-
LDDsimulation/boundary_and_interface.py | 273 ++++++++++++++++----
LDDsimulation/domainPatch.py | 327 +++++++++++++++++++++---
RR-2-patch-test-case/RR-2-patch-test.py | 2 +-
4 files changed, 526 insertions(+), 93 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 84c441f..dc8eb64 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -417,7 +417,7 @@ class LDDsimulation(object):
if iteration > 1:
subsequent_iter_error = dict()
for phase in subdomain.has_phases:
- error = df.Function(subdomain.function_space[phase])
+ error = df.Function(subdomain.function_space["pressure"][phase])
error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
error_calculated = df.norm(error, 'L2')
subsequent_iter_error.update(
@@ -678,7 +678,7 @@ class LDDsimulation(object):
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
pa_exact.t = t
- tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space[phase])
+ tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])
tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
file.write(tmp_exact_pressure, t)
t += self.timestep_size
@@ -829,8 +829,8 @@ class LDDsimulation(object):
def _init_interfaces(self, interface_def_points: tp.List[tp.List[df.Point]] = None,#
adjacent_subdomains: tp.List[np.ndarray] = None) -> None:
- """ generate interfaces and interface marker functions for suddomains and
- set the internal lists used by the LDDsimulation class.
+ """ generate interface objects and interface marker functions for suddomains
+ used by the LDDsimulation class.
This initialises a list of interfaces and global interface marker functions
used by the LDDsimulation class.
@@ -887,13 +887,14 @@ class LDDsimulation(object):
# The reason is that we want to mark outer boundaries with the value 0.
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)
+ self.interface[glob_interface_num].init_facet_dictionaries(self.interface_marker, marker_value, subdomain_index=0)
# print(self.interface[glob_interface_num])
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+"interface_marker_global.pvd"
df.File(filepath) << self.interface_marker
+
def _init_subdomains(self) -> None:
""" generate dictionary of opjects of class DomainPatch.
@@ -991,7 +992,7 @@ class LDDsimulation(object):
interpolation_degree = self.interpolation_degree
for num, subdomain in self.subdomain.items():
- V = subdomain.function_space
+ V = subdomain.function_space["pressure"]
# p0 contains both pw_0 and pnw_0 for subdomain[num]
p0 = self.initial_conditions[num]
for phase in subdomain.has_phases:
@@ -1050,7 +1051,7 @@ class LDDsimulation(object):
else:
self.outerBC.update({subdom_index: dict()})
self.dirichletBC_dfExpression.update({subdom_index: dict()})
- V = subdomain.function_space
+ V = subdomain.function_space["pressure"]
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.
@@ -1099,7 +1100,7 @@ class LDDsimulation(object):
if np.fabs(time - self.starttime) < self.tol:
pass
else:
- V = subdomain.function_space
+ V = subdomain.function_space["pressure"]
boundary_marker = subdomain.outer_boundary_marker
for boundary_index, boundary in enumerate(subdomain.outer_boundary):
boundary_marker_value = boundary.marker_value
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index c8bca3c..eb5033a 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -13,6 +13,7 @@ dofs between subdomains.
import dolfin as df
import numpy as np
import typing as tp
+import copy as cp
class BoundaryPart(df.SubDomain):
@@ -240,15 +241,28 @@ class Interface(BoundaryPart):
# global index ob the subdomain
self.global_index = global_index
self.marker_value = self.global_index+1
+ # 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
+ # self.init_facet_dictionaries()
+ self.facet_to_vertex_coordinates = dict()
+ # dictionary containing maps mapping subdomain mesh indices of facets
+ # to the index of the facet in the global mesh and vice versa.
+ self.local_to_global_facet_map = dict()
+ self.global_to_local_facet_map = dict()
+ # dictionary saving per subdomain a list of outer normals per facet defined
+ # on a mesh passed to the method self.init_facet_dictionaries() which
+ # populates this dictionary.
+ self.outer_normals = dict()
# 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
+ self.pressure_values = dict()
+ self.pressure_values_prev = dict()
# dictionaries holding the forms for the decoupling terms
- self.gli_term = None
- self.gli_term_prev = None
+ self.gli_term = dict()
+ self.gli_term_prev = dict()
# dictionary holding the current iteration number i for each of the adjacent
# subdomains with index in adjacent_subdomains
self.current_iteration = {#
@@ -311,58 +325,231 @@ class Interface(BoundaryPart):
# 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,#
+ interface_marker: df.MeshFunction, #
+ interface_marker_value: int,
+ debug: bool = True) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
+ """ for a given function space function_space, calculate the dof indices
+ and coordinates along 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
+
+ # Output #type tp.Dict[int, tp.Dict[int, np.ndarray]]
+ the first key is the index of
+ the facet relative to the numbering
+ of the mesh that interface_marker
+ is defined on. The value to
+ that key is a dictionary
+ containing global (not cell)
+ dof indix and corresponding
+ coordinates of the dof as
+ key: value pairs.
+
+ """
+ dofmap = function_space.dofmap()
+ element = function_space.element()
+ mesh_entity_iterator = df.SubsetIterator(interface_marker, interface_marker_value)
+ dofs_and_coords = dict()
+ for mesh_entity in mesh_entity_iterator:
+ # 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.
+ for cell in df.cells(mesh_entity):
+ # we actually only have one facet contained in this mesh_entity
+ # but we still need to loop over „all“ facets since df.facets(mesh_entity)
+ # returns an iterator object.
+ facet_cell = cell
+ facet_cell_index = facet_cell.index()
+
+ # the cell is a triangle and for each dof numbered locally for that
+ # cell, we have the coordinates in a numpy array.
+ facet_cell_dof_coordinates = element.tabulate_dof_coordinates(facet_cell)
+ # the cell dof map gives to each local (to the cell) numbering of the
+ # dofs the index of the dof realitve to the global dof numbering of
+ # space, i.e. the index to access the dof with function.vector()[index].
+ facet_cell_dofmap = dofmap.cell_dofs(facet_cell_index)
+
+ for local_dof_ind, dof_coord in enumerate(facet_cell_dof_coordinates):
+ # we only want to store the dofs that are on the boundary. The
+ # cell is a triangle and only has one facet that belongs to
+ # the interface. we want to store only dofs that lie on that particular
+ # facet
+ if self._is_on_boundary_part(dof_coord):
+ global_dof_ind = facet_cell_dofmap[local_dof_ind]
+ dofs_and_coords[facet_index].update({global_dof_ind: dof_coord})
+
+
+ if debug:
+ print(f"newly defined facet_dof_to_coordinates_dictionary of interface[{interface_marker_value-1}]")
+ print(dofs_and_coords.items())
+ return cp.deepcopy(dofs_and_coords)
+
+
+ def init_facet_dictionaries(self,#
+ interface_marker: df.MeshFunction,#
+ interface_marker_value: int,
+ subdomain_index: int,
+ debug: bool = True) -> None:
+ """ create dictionary for each index subdomain_index containing the
+ index of each facet belonging to the interface (with respect to the mesh of
+ the subdomain subdomain_index) as key and stores a
+ list of the coordinates of each vertex (numpy arrays) belonging to the facet as
+ values.
+ Additionally, create dictionary storing per facet the outer normals, i.e.
+ pairs {facet_index: facet.normal()}.
+ This is needed by SubDomain.calc_gl0_term() for accessing the normals.
+
+ The MeshFunction interface_marker is supposed to be defined on the
+ subdomain mesh belonging to subdomain subdomain_index.
+
+ # 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
+ """
+ self.facet_to_vertex_coordinates.update({subdomain_index: dict()})
+ self.outer_normals.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)
+ for mesh_entity in mesh_entity_iterator:
+ # because mesh_entity doesn't have a method to access coordinates we
+ # need to cast it into an acutal facet object.
+ for f in df.facets(mesh_entity):
+ # we actually only have one facet contained in this mesh_entity
+ # but we still need to loop over „all“ facets since df.facets(mesh_entity)
+ # returns an iterator object.
+ facet = f
+
+ facet_vertices = []
+ for vertex in df.vertices(facet):
+ # we only need the x,y coordinates. vertex.point()[:] returns
+ # the array and cannot be sliced. the result can be indexed, hence,
+ # vertex.point()[:][0:2]
+ facet_vertices.append(vertex.point()[:][0:2])
+
+ self.facet_to_vertex_coordinates[subdomain_index].update(
+ {facet.index(): facet_vertices}
+ )
+ if subdomain_index == 0:
+ self.outer_normals[subdomain_index].update(
+ {facet.index(): None}
+ )
+ else:
+ self.outer_normals[subdomain_index].update(
+ {facet.index(): facet.normal()[:][0:2]}
+ )
+ 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:
+ """ 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)"""
+ if subdomain_index == 0:
+ raise RuntimeError("expeced a subdomain index != 0 for calculation",
+ "of local_to_global_facet_map.")
+
+ self.local_to_global_facet_map.update({subdomain_index: dict()})
+ self.global_to_local_facet_map.update({subdomain_index: dict()})
+ global_facet2coord = self.facet_to_vertex_coordinates[0]
+ subdom_facet2coord = self.facet_to_vertex_coordinates[subdomain_index]
+ for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
+ for global_facet_index, global_facet_vertex in global_facet2coord.items():
+ # vertices are not necessarily in the same order
+ vertices_are_equal = np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
+ 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
+ if both_vertices_are_equal:
+ self.local_to_global_facet_map[subdomain_index].update(
+ {local_facet_index: global_facet_index}
+ )
+ self.global_to_local_facet_map[subdomain_index].update(
+ {global_facet_index: local_facet_index}
+ )
+ if debug:
+ print("local_to_global_facet_map: ", self.local_to_global_facet_map[subdomain_index].items())
+ print("global_to_local_facet_map: ", self.global_to_local_facet_map[subdomain_index].items())
+
def init_interface_values(self,#
interface_marker: df.MeshFunction,#
- interface_marker_value: int) -> None:
+ interface_marker_value: int,
+ function_space: tp.Dict[str, tp.Dict[str, df.FunctionSpace]],
+ subdomain_index: int,
+ has_phases: tp.List[str],
+ ) -> None:
""" allocate dictionaries that are used to store the interface dof values
of the pressure, previous pressure, gl decoupling term as well as previous
gl update term respectively as pairs(parent_mesh_index: dof_value)
- The MeshFunction interface_marker is supposed to be one defined on the
- global mesh, so that the indices used in the dictionary are the global
- vertex indices and not indices of some submesh.
-
# Parameters
interface_marker: #type df.MeshFunction marker function marking
edges of i_mesh
according to self.inside
interface_marker_value: #type int marker value of interface_marker
"""
- vertex_indices = self._vertex_indices(interface_marker, #
- interface_marker_value,#
- print_vertex_indices = False)
- self.pressure_values = dict()
- #self.pressure_values is a dictionay which contains for each index in
- # adjacent_subdomains a dictionary for wetting and nonwetting phase with
- # as a value yet another dictionary holding the interface values. phew.
- self.pressure_values.update({self.adjacent_subdomains[0]: dict()})
- self.pressure_values.update({self.adjacent_subdomains[1]: dict()})
- self.pressure_values_prev = dict()
- self.pressure_values_prev.update({self.adjacent_subdomains[0]: dict()})
- self.pressure_values_prev.update({self.adjacent_subdomains[1]: dict()})
- self.gli_term = dict()
- self.gli_term.update({self.adjacent_subdomains[0]: dict()})
- self.gli_term.update({self.adjacent_subdomains[1]: dict()})
- self.gli_term_prev = dict()
- self.gli_term_prev.update({self.adjacent_subdomains[0]: dict()})
- self.gli_term_prev.update({self.adjacent_subdomains[1]: dict()})
- # we initialise all dictionaries to accomodate two phases, to be more
- # flexible with coupling Richards and Two-phase models.
- for subdom_ind in self.adjacent_subdomains:
- for phase in ['wetting', 'nonwetting']:
- self.pressure_values[subdom_ind].update({phase: dict()})
- self.pressure_values_prev[subdom_ind].update({phase: dict()})
- # the gli_term will hold dofs of the gli terms for communications
- # they are initiated with 0
- self.gli_term[subdom_ind].update({phase: dict()})
- self.gli_term_prev[subdom_ind].update({phase: dict()})
- for vert_ind in vertex_indices:
- self.pressure_values[subdom_ind][phase].update({vert_ind: 0.0})
- self.pressure_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
- self.gli_term[subdom_ind][phase].update({vert_ind: 0.0})
- self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
- # end init_interface_values
+ function_name = ["pressure", "gli"]
+ for function in function_name:
+ if function == "pressure":
+ self.pressure_values.update({subdomain_index: dict()})
+ self.pressure_values_prev.update({subdomain_index: dict()})
+ else: #function == "gli":
+ self.gli_term.update({subdomain_index: dict()})
+ self.gli_term_prev.update({subdomain_index: dict()})
+
+ for phase in has_phases:
+ if function == "pressure":
+ self.pressure_values[subdomain_index].update({phase: dict()})
+ self.pressure_values_prev[subdomain_index].update({phase: dict()})
+ else: #function == "gli":
+ self.gli_term[subdomain_index].update({phase: dict()})
+ self.gli_term_prev[subdomain_index].update({phase: dict()})
+
+ space = function_space[function][phase]
+ dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
+ # the dof dictionaries are grouped by the global (with respect to the
+ # global mesh) facet index of each facet along the interface.
+ # each dof is then saved with its koordinates as tuples as keys.
+ for subdom_mesh_facet_index, facet_dof_coord_dict in dof_coordinates.items():
+ global_mesh_facet_index = self.local_to_global_facet_map[subdomain_index][subdom_mesh_facet_index]
+ if function == "pressure":
+ self.pressure_values[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+ self.pressure_values_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+ else: #function == "gli":
+ self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+ self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+
+ for dof_coord in facet_dof_coord_dict.values():
+ # we want to use the coordinates as key in a dictionary.
+ # Because dof_coord is a numpy array and numpy array are
+ # not hashable, we convert the numpy array into a tuple,
+ # which in turn can then be used as a dict key.
+ coord_tupel = (dof_coord[0], dof_coord[1])
+ if function == "pressure":
+ self.pressure_values[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ self.pressure_values_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ else: #function == "gli":
+ self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ # end init_interface_values
+
def read_pressure_dofs(self, from_function: df.Function, #
interface_dofs: np.array,#
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 1d01e0d..bbf83b8 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -197,7 +197,7 @@ class DomainPatch(df.SubDomain):
# Dictionary holding the exact solution expression if we have a case with exact solution.
# Is set by LDDsimualtion._init_exact_solution_expression()
self.pressure_exact = None
- # Dictionary of FEM function of self.function_space holding the interface
+ # Dictionary of FEM function of self.function_space["pressure"] 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
@@ -209,6 +209,8 @@ class DomainPatch(df.SubDomain):
self._init_dof_maps()
self._init_markers()
self._init_measures()
+ self._calc_interface_dof_indices_and_coordinates()
+ self._init_interface_values()
self._calc_dof_indices_of_interfaces()
self._calc_interface_has_common_dof_indices()
self._calc_normal_vectors_norms_for_common_dofs()
@@ -230,7 +232,7 @@ class DomainPatch(df.SubDomain):
self.gli_function = dict()
for phase in self.has_phases:
self.gli_function.update(
- {phase: df.Function(self.function_space[phase])}
+ {phase: df.Function(self.function_space["pressure"][phase])}
)
# END constructor
@@ -424,7 +426,7 @@ class DomainPatch(df.SubDomain):
pnw_prev = self.pressure_prev_timestep['nonwetting']
pc_prev = pnw_prev - pw_prev
- n = df.FacetNormal(self.mesh)
+ # n = df.FacetNormal(self.mesh)
S = self.saturation
for phase in self.has_phases:
# viscosity
@@ -467,28 +469,197 @@ class DomainPatch(df.SubDomain):
# 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()
- # normal_flux = df.dot(flux, n).get_local()
- # print("type(normal_flux)", type(normal_flux))
- V = self.function_space[phase]
- degree = V.ufl_element().degree()
- interface_mesh = df.restriction(self.interface_marker, interface.marker_value)
- # W = VectorFunctionSpace(interface_mesh, 'Q', 0)
- W = FunctionSpace(interface_mesh, 'Q', 0)
+ 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
- normal_flux_Q0 = df.project(df.dot(flux,n), W)
- 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}
- )
+ # 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}
+ # )
### END calc_gl0_term
@@ -510,7 +681,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[phase])
+ gli = df.Function(self.function_space["pressure"][phase])
# neighbour index
neighbour = interface.neighbour[subdomain]
# update the current_iteration number of subdomain stored in the interface
@@ -618,22 +789,32 @@ class DomainPatch(df.SubDomain):
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.function_space["pressure"]
self.trialfunction
self.pressure
self.testfunction
"""
+ function_name = ["pressure", "gli", "flux"]
+
self.function_space = dict()
- self.trialfunction = dict()
- self.testfunction = dict()
- # self.neighbouring_interface_pressure = dict()
- for phase in self.has_phases:
- self.function_space.update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
- self.trialfunction.update({phase: df.TrialFunction(self.function_space[phase])})
- self.testfunction.update({phase: df.TestFunction(self.function_space[phase])})
- # self.neighbouring_interface_pressure.update(
- # {phase: df.interpolate(df.Constant(0), self.function_space[phase])}
- # )
+ for function in function_name:
+ self.function_space.update({function: dict()})
+ if function == "pressure":
+ self.trialfunction = dict()
+ self.testfunction = dict()
+ # self.neighbouring_interface_pressure = dict()
+ for phase in self.has_phases:
+ if function == "pressure":
+ self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
+ self.trialfunction.update({phase: df.TrialFunction(self.function_space["pressure"][phase])})
+ self.testfunction.update({phase: df.TestFunction(self.function_space["pressure"][phase])})
+ elif function == "gli":
+ degree = self.function_space["pressure"][phase].ufl_element().degree()
+ self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'DG', degree)})
+ else:
+ degree = self.function_space["pressure"][phase].ufl_element().degree()
+ self.function_space[function].update({phase: df.VectorFunctionSpace(self.mesh, 'DG', degree)})
+
def _init_dof_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
@@ -649,15 +830,28 @@ class DomainPatch(df.SubDomain):
# print(f"on subdomain{self.subdomain_index} parent_vertex_indices: {self.parent_mesh_index}")
self.dof2vertex = dict()
self.vertex2dof = dict()
+ self.dofmap = dict()
+ self.dofmap.update({"pressure": dict()})
+ self.dofmap.update({"gli": dict()})
+ self.dofmap.update({"flux": dict()})
for phase in self.has_phases:
self.dof2vertex.update(#
- {phase :df.dof_to_vertex_map(self.function_space[phase])}#
+ {phase :df.dof_to_vertex_map(self.function_space["pressure"][phase])}#
)
# print(f"on subdomain{self.subdomain_index} dof2vertex: {self.dof2vertex[phase]}")
self.vertex2dof.update(#
- {phase :df.vertex_to_dof_map(self.function_space[phase])}#
+ {phase :df.vertex_to_dof_map(self.function_space["pressure"][phase])}#
)
# print(f"on subdomain{self.subdomain_index} vertex2dof: {self.vertex2dof[phase]}")
+ self.dofmap["pressure"].update(
+ {phase: self.function_space["pressure"][phase].dofmap()}
+ )
+ self.dofmap["gli"].update(
+ {phase: self.function_space["gli"][phase].dofmap()}
+ )
+ self.dofmap["flux"].update(
+ {phase: self.function_space["flux"][phase].dofmap()}
+ )
def _init_markers(self) -> None:
""" define boundary markers
@@ -677,6 +871,12 @@ class DomainPatch(df.SubDomain):
marker_value = self.interface[glob_index].marker_value
# each interface gets marked with the global interface index.
self.interface[glob_index].mark(self.interface_marker, marker_value)
+ # init facet_to_vertex_coordinates_dictionray for subdoaain mesh
+ self.interface[glob_index].init_facet_dictionaries(
+ 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.
@@ -712,9 +912,54 @@ class DomainPatch(df.SubDomain):
# measure over the interfaces
self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.interface_marker)
+
+ def _calc_interface_dof_indices_and_coordinates(self):
+ """ calculate dictionaries containing for each local facet index of
+ interfaces a dictionary containing {interface_dof_index: dof_coordinates}
+ """
+ function_name = ["pressure", "gli", "flux"]
+ marker = self.interface_marker
+ self.interface_dof_indices_and_coordinates = dict()
+ for function in function_name:
+ self.interface_dof_indices_and_coordinates.update(
+ {function: dict()}
+ )
+ function_space = self.function_space[function]
+ for interface_index in self.has_interface:
+ interface = self.interface[interface_index]
+ marker_value = interface.marker_value
+ self.interface_dof_indices_and_coordinates[function].update(
+ {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)}
+ )
+
+
+ def _init_interface_values(self):
+ """ allocate the dictionaries for each interface used to store the
+ interface values for communication.
+ """
+ # we don't need to communicate flux values directly.
+ function_name = ["pressure", "gli"]
+ marker = self.interface_marker
+
+ for interf in self.has_interface:
+ marker_value = self.interface[interf].marker_value
+ space = self.function_space
+ self.interface[interf].init_interface_values(
+ interface_marker=marker,
+ interface_marker_value=marker_value,
+ function_space=space,
+ subdomain_index=self.subdomain_index,
+ has_phases=self.has_phases,
+ )
+
+
def _calc_dof_indices_of_interfaces(self):
""" calculate dof indices of each interface """
- V = self.function_space
+ V = self.function_space["pressure"]
marker = self.interface_marker
for ind in self.has_interface:
# the marker value for the interfaces is always global interface index+1
@@ -850,7 +1095,7 @@ class DomainPatch(df.SubDomain):
file.write(self.pressure[phase], t)
# if self.pressure_exact:
# self.pressure_exact[phase].t = t
- # tmp_exact_pressure = df.interpolate(self.pressure_exact[phase], self.function_space[phase])
+ # tmp_exact_pressure = df.interpolate(self.pressure_exact[phase], self.function_space["pressure"][phase])
# tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
# file.write(tmp_exact_pressure, t)
else:
@@ -860,7 +1105,7 @@ class DomainPatch(df.SubDomain):
)
file.write(self.pressure[phase], i)
# create function to copy the subsequent errors to
- error = df.Function(self.function_space[phase])
+ error = df.Function(self.function_space["pressure"][phase])
error.assign(self.pressure[phase] - self.pressure_prev_iter[phase])
error.rename("pi_"+"{}".format(phase)+"-pim1_"+"{}".format(phase), #
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
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 14c9635..aa8aff2 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -86,7 +86,7 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 10
+mesh_resolution = 4
##################### TIME #####################################################
timestep_size = 0.01
number_of_timesteps = 100
--
GitLab