diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 9d3a9f3988d2dbddb9e245416f905be2ace2ce7a..507e2b116741cc140f076d52dc732d3ec5c9f2fe 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -555,216 +555,216 @@ class Interface(BoundaryPart):
# 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(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\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[:])
+ # 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(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\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(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\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 ####
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 07f5f0e32daf3e158f99b1128736d486443ff59e..f70f5a5ec4257dd834bff8f590df5dc5b310fcf1 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -211,28 +211,22 @@ class DomainPatch(df.SubDomain):
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()
+ # self._calc_dof_indices_of_interfaces()
+ # self._calc_interface_has_common_dof_indices()
+ # self._calc_normal_vectors_norms_for_common_dofs()
if self.include_gravity:
self.gravity_term = dict()
self._calc_gravity_expressions()
else:
self.gravity_term = 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 holding the same dofs gli_assembled above
# but as a df.Function. This is needed for writing out the gli terms when
# analyising the convergence of a timestep.
self.gli_function = dict()
for phase in self.has_phases:
self.gli_function.update(
- {phase: df.Function(self.function_space["pressure"][phase])}
+ {phase: df.Function(self.function_space["gli"][phase])}
)
# END constructor
@@ -249,10 +243,10 @@ class DomainPatch(df.SubDomain):
for glob_if_index in self.has_interface:
interface = self.interface[glob_if_index]
interface_string.append(interface.__str__())
- common_int_dof_str = ".\n Dofs common with other interfaces:"\
- +"{}".format([(ind, dof) for ind, dof in self._interface_has_common_dof_indices[glob_if_index].items()])
- interface_string.append(common_int_dof_str)
-
+ # common_int_dof_str = ".\n Dofs common with other interfaces:"\
+ # +"{}".format([(ind, dof) for ind, dof in self._interface_has_common_dof_indices[glob_if_index].items()])
+ # interface_string.append(common_int_dof_str)
+ #
interface_string = " ".join(interface_string)
if self.outer_boundary is None:
@@ -278,16 +272,22 @@ class DomainPatch(df.SubDomain):
to interface.pressure_values or to
interface.pressure_values_prev
"""
- from_func = self.pressure
+ subdomain = self.subdomain_index
for ind in self.has_interface:
+ interface = self.interface[ind]
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,
- previous_iter = previous_iter)
+ 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]
+ if previous_iter:
+ save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_ind]
+ else:
+ save_dict = interface.pressure_values[subdomain][phase][global_facet_ind]
+ for dof_index, dof_coord in p_dof_coordinates[facet_ind].items():
+ dof_coord_tupel = (dof_coord[0], dof_coord[1])
+ save_dict.update({dof_coord_tupel: self.pressure[phase].vector()[dof_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.
@@ -830,102 +830,102 @@ class DomainPatch(df.SubDomain):
)
- def _calc_dof_indices_of_interfaces(self):
- """ calculate dof indices of each interface """
- 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
- marker_value = self.interface[ind].marker_value
- 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])
-
- 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]}
- )
- # for phase in self.has_phases:
- # print(f"self._interface_has_common_dof_indices[{interf_ind}][{phase}]", self._interface_has_common_dof_indices[interf_ind][phase])
- # print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
-
-
- def _calc_normal_vectors_norms_for_common_dofs(self):
- """ calculate norm(n1+n2) where n1 and n2 are normal vectors of adjacent
- interfaces. This is used by the functions calc_gli_term and calc_gl0_term
- """
- self._normal_vector_norms_for_common_dofs = dict()
- self.mesh.init(0)
- for interface_ind in self.has_interface:
- # if one phase has common dofs to two interfaces, the both phases
- # have, so it suffices to check for the wetting phase as it is always
- # present.
- if self._interface_has_common_dof_indices[interface_ind]['wetting'].size != 0:
- self._normal_vector_norms_for_common_dofs.update({interface_ind: dict()})
- for phase in self.has_phases:
- self._normal_vector_norms_for_common_dofs[interface_ind].update(
- {phase: dict()}
- )
- for dof in self._interface_has_common_dof_indices[interface_ind][phase]:
- common_vertex_index = self.dof2vertex[phase][dof]
- for index, vertex in enumerate(df.vertices(self.mesh)):
- if index == common_vertex_index:
- common_vertex = vertex
- # print(f"\non subdomain{self.subdomain_index}, interface{interface_ind} and phase {phase}")
- # print(f"vertex_{index} = {vertex.point()[:]}, type(vertex) = {type(vertex)}")
- linear_combination_of_normals = np.zeros(np.size(vertex.point()[:]), dtype=float)
- # normals = []
- # print(f"common_vertex = {common_vertex}",
- # f"common_vertex_index = {common_vertex_index}")
- for facet_containing_vertex in df.facets(common_vertex):#df.cells(common_vertex):
- if facet_containing_vertex.exterior():
- # normals.append(facet_containing_vertex.normal()[:])
- # print(f"type(facet_normal) = {type(facet_containing_vertex.normal()[:])}")
- linear_combination_of_normals += facet_containing_vertex.normal()[:]
-
- # print(f"normals of common_vertex: {common_vertex} are ", normals)
- # print(f"linear combination is: {linear_combination_of_normals[:]}.")
- # print(f"with length {np.linalg.norm(linear_combination_of_normals, ord=2, axis=None)}. For comparison, np.sqrt(2) = {np.sqrt(2)}.")
- self._normal_vector_norms_for_common_dofs[interface_ind][phase].update(
- {dof: np.linalg.norm(linear_combination_of_normals)}
- )
+ # def _calc_dof_indices_of_interfaces(self):
+ # """ calculate dof indices of each interface """
+ # 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
+ # marker_value = self.interface[ind].marker_value
+ # 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])
+ #
+ # 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]}
+ # )
+ # # for phase in self.has_phases:
+ # # print(f"self._interface_has_common_dof_indices[{interf_ind}][{phase}]", self._interface_has_common_dof_indices[interf_ind][phase])
+ # # print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
+
+
+ # def _calc_normal_vectors_norms_for_common_dofs(self):
+ # """ calculate norm(n1+n2) where n1 and n2 are normal vectors of adjacent
+ # interfaces. This is used by the functions calc_gli_term and calc_gl0_term
+ # """
+ # self._normal_vector_norms_for_common_dofs = dict()
+ # self.mesh.init(0)
+ # for interface_ind in self.has_interface:
+ # # if one phase has common dofs to two interfaces, the both phases
+ # # have, so it suffices to check for the wetting phase as it is always
+ # # present.
+ # if self._interface_has_common_dof_indices[interface_ind]['wetting'].size != 0:
+ # self._normal_vector_norms_for_common_dofs.update({interface_ind: dict()})
+ # for phase in self.has_phases:
+ # self._normal_vector_norms_for_common_dofs[interface_ind].update(
+ # {phase: dict()}
+ # )
+ # for dof in self._interface_has_common_dof_indices[interface_ind][phase]:
+ # common_vertex_index = self.dof2vertex[phase][dof]
+ # for index, vertex in enumerate(df.vertices(self.mesh)):
+ # if index == common_vertex_index:
+ # common_vertex = vertex
+ # # print(f"\non subdomain{self.subdomain_index}, interface{interface_ind} and phase {phase}")
+ # # print(f"vertex_{index} = {vertex.point()[:]}, type(vertex) = {type(vertex)}")
+ # linear_combination_of_normals = np.zeros(np.size(vertex.point()[:]), dtype=float)
+ # # normals = []
+ # # print(f"common_vertex = {common_vertex}",
+ # # f"common_vertex_index = {common_vertex_index}")
+ # for facet_containing_vertex in df.facets(common_vertex):#df.cells(common_vertex):
+ # if facet_containing_vertex.exterior():
+ # # normals.append(facet_containing_vertex.normal()[:])
+ # # print(f"type(facet_normal) = {type(facet_containing_vertex.normal()[:])}")
+ # linear_combination_of_normals += facet_containing_vertex.normal()[:]
+ #
+ # # print(f"normals of common_vertex: {common_vertex} are ", normals)
+ # # print(f"linear combination is: {linear_combination_of_normals[:]}.")
+ # # print(f"with length {np.linalg.norm(linear_combination_of_normals, ord=2, axis=None)}. For comparison, np.sqrt(2) = {np.sqrt(2)}.")
+ # self._normal_vector_norms_for_common_dofs[interface_ind][phase].update(
+ # {dof: np.linalg.norm(linear_combination_of_normals)}
+ # )
def _calc_gravity_expressions(self):
@@ -984,9 +984,9 @@ class DomainPatch(df.SubDomain):
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
)
file.write(error, i)
- self.gli_function[phase].vector().set_local(self.gli_assembled[phase])
- self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
- "all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
- file.write(self.gli_function[phase], i)
+ # self.gli_function[phase].vector().set_local(self.gli_assembled[phase])
+ # self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
+ # "all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
+ # file.write(self.gli_function[phase], i)
# END OF CLASS DomainPatch