Skip to content
Snippets Groups Projects
Commit afcd64e6 authored by David Seus's avatar David Seus
Browse files

lauffähige neue Version. depricate a lot of methods :(

parent 48b646f1
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
LDDsimulation/boundary_and_interface.py
+ 210
− 210
View file @ afcd64e6
...@@ -555,216 +555,216 @@ class Interface(BoundaryPart): ...@@ -555,216 +555,216 @@ class Interface(BoundaryPart):
# end init_interface_values # end init_interface_values
def read_pressure_dofs(self, from_function: df.Function, # # def read_pressure_dofs(self, from_function: df.Function, #
interface_dofs: np.array,# # interface_dofs: np.array,#
dof_to_vert_map: np.array,# # dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,# # local_to_parent_vertex_map: np.array,#
phase: str,# # phase: str,#
subdomain_ind: int,# # subdomain_ind: int,#
previous_iter: bool = False,# # previous_iter: bool = False,#
) -> None: # ) -> None:
""" read dofs of from_function with indices interface_dofs and write to self.pressure_values # """ 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 # 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 # interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
of the interface # of the interface
self.pressure_values[subdomain_ind][phase], (if pressure == True) # self.pressure_values[subdomain_ind][phase], (if pressure == True)
self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True) # self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
#
# Parameters # # Parameters
from_function: #type: df.Function, function to save the dofs of # from_function: #type: df.Function, function to save the dofs of
interface_dofs: #type: np.array, index array of dofs to be saved # 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 # dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh. # space on the submesh.
local_to_parent_vertex_map: #type: np.array, # local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or # phase #type: str is either 'wetting' or
'nonwetting' and determins # 'nonwetting' and determins
the dict to be written to. # the dict to be written to.
subdomain_ind #type: int subdomain index # subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous # previous_iter #type: bool determins wether a previous
iterate is being written or not # iterate is being written or not
""" # """
if self.pressure_values == None: # if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \ # raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") # run self.init_interface_values")
#
parent = local_to_parent_vertex_map # parent = local_to_parent_vertex_map
d2v = dof_to_vert_map # d2v = dof_to_vert_map
if not previous_iter: # if not previous_iter:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]}) # 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]) # # print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
else: # else:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]}) # 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]) # # 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, # # def read_gli_dofs(self, from_array: np.array, #
interface_dofs: np.array,# # interface_dofs: np.array,#
dof_to_vert_map: np.array,# # dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,# # local_to_parent_vertex_map: np.array,#
phase: str,# # phase: str,#
subdomain_ind: int,# # subdomain_ind: int,#
previous_iter: bool = False,# # previous_iter: bool = False,#
) -> None: # ) -> None:
""" Read dofs of from_array with indices interface_dofs and write to self.gli_term # """ 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 # 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 # interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
of the interface # of the interface
self.gli_term[subdomain_ind][phase], (if gl == True) # self.gli_term[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True) # self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
#
# Parameters # # Parameters
from_array: #type: np.array, array to save the dofs of # from_array: #type: np.array, array to save the dofs of
interface_dofs: #type: np.array, index array of dofs to be saved # 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 # dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh. # space on the submesh.
local_to_parent_vertex_map: #type: np.array, # local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or # phase #type: str is either 'wetting' or
'nonwetting' and determins # 'nonwetting' and determins
the dict to be written to. # the dict to be written to.
subdomain_ind #type: int subdomain index # subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous # previous_iter #type: bool determins wether a previous
iterate is being written or not # iterate is being written or not
""" # """
if self.pressure_values == None: # if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \ # raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") # run self.init_interface_values")
#
parent = local_to_parent_vertex_map # parent = local_to_parent_vertex_map
d2v = dof_to_vert_map # d2v = dof_to_vert_map
if not previous_iter: # if not previous_iter:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_array is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]}) # 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]) # # 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: # else:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_array is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]}) # 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]) # # 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, # # def write_pressure_dofs(self, to_function: df.Function, #
interface_dofs: np.array,# # interface_dofs: np.array,#
dof_to_vert_map: np.array,# # dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,# # local_to_parent_vertex_map: np.array,#
phase: str,# # phase: str,#
subdomain_ind: int,# # subdomain_ind: int,#
# pressure: bool = False,# # # pressure: bool = False,#
# gl: bool = False,# # # gl: bool = False,#
previous_iter: bool = False,# # previous_iter: bool = False,#
) -> None: # ) -> None:
""" read dofs off self.pressure_values and overwrite the dofs of to_function # """ read dofs off self.pressure_values and overwrite the dofs of to_function
with indices interface_dofs # with indices interface_dofs
#
Read the degrees of freedom of one of the dictionaries of the interface # 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[subdomain_ind][phase], (if pressure == True)
self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True) # self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
and overwrite the dofs of the function to_function with indices # and overwrite the dofs of the function to_function with indices
interface_dofs, i.e. update the interface dofs of the function. # interface_dofs, i.e. update the interface dofs of the function.
. # .
#
# Parameters # # Parameters
to_function: #type: df.Function, function to overwrite the # to_function: #type: df.Function, function to overwrite the
dofs of. # dofs of.
interface_dofs: #type: np.array, index array of interface # interface_dofs: #type: np.array, index array of interface
dofs of from_function # dofs of from_function
dof_to_vert_map: #type: np.array, dof to vert map of the function # dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh. # space on the submesh.
local_to_parent_vertex_map: #type: np.array, # local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or # phase #type: str is either 'wetting' or
'nonwetting' and determins # 'nonwetting' and determins
the dict to be read of. # the dict to be read of.
subdomain_ind #type: int subdomain index # subdomain_ind #type: int subdomain index
previous_iter #type: bool determins wether a previous # previous_iter #type: bool determins wether a previous
iterate is being read or not # iterate is being read or not
""" # """
if self.pressure_values == None: # if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \ # raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") # run self.init_interface_values")
#
parent = local_to_parent_vertex_map # parent = local_to_parent_vertex_map
d2v = dof_to_vert_map # d2v = dof_to_vert_map
#
if not previous_iter: # if not previous_iter:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]]] # 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("read pressure interface values\n", self.pressure_values[subdomain_ind][phase])
# print("wrote these dofs to function \n", to_function.vector()[:]) # # print("wrote these dofs to function \n", to_function.vector()[:])
else: # else:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# from_function.vector() is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]]] # 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("read previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
# print("wrote these dofs to function \n", to_function.vector()[:]) # # print("wrote these dofs to function \n", to_function.vector()[:])
#
#
def write_gli_dofs(self, to_array: np.array, # # def write_gli_dofs(self, to_array: np.array, #
interface_dofs: np.array,# # interface_dofs: np.array,#
dof_to_vert_map: np.array,# # dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,# # local_to_parent_vertex_map: np.array,#
phase: str,# # phase: str,#
subdomain_ind: int,# # subdomain_ind: int,#
previous_iter: bool = False,# # previous_iter: bool = False,#
) -> None: # ) -> None:
""" read dofs off self.gli_term and overwrite the dofs stored in the array # """ read dofs off self.gli_term and overwrite the dofs stored in the array
to_array with indices interface_dofs # to_array with indices interface_dofs
#
Read the degrees of freedom of one of the dictionaries of the interface # 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[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == 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 # 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. # interface_dofs, i.e. update the interface dofs of the gli terms.
. # .
#
# Parameters # # Parameters
to_array: #type: np.array, array to overwrite the # to_array: #type: np.array, array to overwrite the
dofs of. # dofs of.
interface_dofs: #type: np.array, index array of interface # interface_dofs: #type: np.array, index array of interface
dofs of from_function # dofs of from_function
dof_to_vert_map: #type: np.array, dof to vert map of the function # dof_to_vert_map: #type: np.array, dof to vert map of the function
space on the submesh. # space on the submesh.
local_to_parent_vertex_map: #type: np.array, # local_to_parent_vertex_map: #type: np.array,
phase #type: str is either 'wetting' or # phase #type: str is either 'wetting' or
'nonwetting' and determins # 'nonwetting' and determins
the dict to be read of. # the dict to be read of.
subdomain_ind #type: int subdomain index # subdomain_ind #type: int subdomain index
#
previous_iter #type: bool determins wether a previous # previous_iter #type: bool determins wether a previous
iterate is being read or not # iterate is being read or not
""" # """
if self.pressure_values == None: # if self.pressure_values == None:
raise RuntimeError("self.pressure_values not initiated. You forgot to \ # raise RuntimeError("self.pressure_values not initiated. You forgot to \
run self.init_interface_values") # run self.init_interface_values")
#
parent = local_to_parent_vertex_map # parent = local_to_parent_vertex_map
d2v = dof_to_vert_map # d2v = dof_to_vert_map
if not previous_iter: # if not previous_iter:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# to_array is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]]] # 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("read gl interface values\n", self.gli_term[subdomain_ind][phase])
# print("wrote these dofs to array \n", to_array[:]) # # print("wrote these dofs to array \n", to_array[:])
else: # else:
for dof_index in interface_dofs: # for dof_index in interface_dofs:
# to_array is orderd by dof and not by vertex numbering of the mesh. # # 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. # # 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]]] # 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("read gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
# print("wrote these dofs to array \n", to_array[:]) # # print("wrote these dofs to array \n", to_array[:])
#### Private Methods #### #### Private Methods ####
... ...
......
LDDsimulation/domainPatch.py
+ 123
− 123
View file @ afcd64e6
...@@ -211,28 +211,22 @@ class DomainPatch(df.SubDomain): ...@@ -211,28 +211,22 @@ class DomainPatch(df.SubDomain):
self._init_measures() self._init_measures()
self._calc_interface_dof_indices_and_coordinates() self._calc_interface_dof_indices_and_coordinates()
self._init_interface_values() self._init_interface_values()
self._calc_dof_indices_of_interfaces() # self._calc_dof_indices_of_interfaces()
self._calc_interface_has_common_dof_indices() # self._calc_interface_has_common_dof_indices()
self._calc_normal_vectors_norms_for_common_dofs() # self._calc_normal_vectors_norms_for_common_dofs()
if self.include_gravity: if self.include_gravity:
self.gravity_term = dict() self.gravity_term = dict()
self._calc_gravity_expressions() self._calc_gravity_expressions()
else: else:
self.gravity_term = None 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 # 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 # but as a df.Function. This is needed for writing out the gli terms when
# analyising the convergence of a timestep. # analyising the convergence of a timestep.
self.gli_function = dict() self.gli_function = dict()
for phase in self.has_phases: for phase in self.has_phases:
self.gli_function.update( self.gli_function.update(
{phase: df.Function(self.function_space["pressure"][phase])} {phase: df.Function(self.function_space["gli"][phase])}
) )
# END constructor # END constructor
...@@ -249,10 +243,10 @@ class DomainPatch(df.SubDomain): ...@@ -249,10 +243,10 @@ class DomainPatch(df.SubDomain):
for glob_if_index in self.has_interface: for glob_if_index in self.has_interface:
interface = self.interface[glob_if_index] interface = self.interface[glob_if_index]
interface_string.append(interface.__str__()) interface_string.append(interface.__str__())
common_int_dof_str = ".\n Dofs common with other interfaces:"\ # 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()]) # +"{}".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.append(common_int_dof_str)
#
interface_string = " ".join(interface_string) interface_string = " ".join(interface_string)
if self.outer_boundary is None: if self.outer_boundary is None:
...@@ -278,16 +272,22 @@ class DomainPatch(df.SubDomain): ...@@ -278,16 +272,22 @@ class DomainPatch(df.SubDomain):
to interface.pressure_values or to to interface.pressure_values or to
interface.pressure_values_prev interface.pressure_values_prev
""" """
from_func = self.pressure subdomain = self.subdomain_index
for ind in self.has_interface: for ind in self.has_interface:
interface = self.interface[ind]
for phase in self.has_phases: for phase in self.has_phases:
self.interface[ind].read_pressure_dofs(from_function = from_func[phase], # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][ind][phase]
interface_dofs = self._dof_indices_of_interface[ind][phase],# local2global_facet = interface.local_to_global_facet_map[subdomain]
dof_to_vert_map = self.dof2vertex[phase],# for facet_ind in interface.facets[subdomain].keys():
local_to_parent_vertex_map = self.parent_mesh_index,# global_facet_ind = local2global_facet[facet_ind]
phase = phase,# if previous_iter:
subdomain_ind = self.subdomain_index, save_dict = interface.pressure_values_prev[subdomain][phase][global_facet_ind]
previous_iter = previous_iter) 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]: 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. """ return the governing form and right hand side for phase phase as a dictionary.
...@@ -830,102 +830,102 @@ class DomainPatch(df.SubDomain): ...@@ -830,102 +830,102 @@ class DomainPatch(df.SubDomain):
) )
def _calc_dof_indices_of_interfaces(self): # def _calc_dof_indices_of_interfaces(self):
""" calculate dof indices of each interface """ # """ calculate dof indices of each interface """
V = self.function_space["pressure"] # V = self.function_space["pressure"]
marker = self.interface_marker # marker = self.interface_marker
for ind in self.has_interface: # for ind in self.has_interface:
# the marker value for the interfaces is always global interface index+1 # # the marker value for the interfaces is always global interface index+1
marker_value = self.interface[ind].marker_value # marker_value = self.interface[ind].marker_value
self._dof_indices_of_interface.update(# # self._dof_indices_of_interface.update(#
{ind: dict()} # {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: # 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]) # self._dof_indices_of_interface[ind].update(#
# print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase])) # {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_normal_vectors_norms_for_common_dofs(self): #
""" calculate norm(n1+n2) where n1 and n2 are normal vectors of adjacent # def _calc_interface_has_common_dof_indices(self):
interfaces. This is used by the functions calc_gli_term and calc_gl0_term # """ populate dictionary self._interface_has_common_dof_indices
""" #
self._normal_vector_norms_for_common_dofs = dict() # Determin for each interface in self.has_interface the dof indices that
self.mesh.init(0) # also belong to another interface and store them in the dictionary
for interface_ind in self.has_interface: # self._interface_has_common_dof_indices
# if one phase has common dofs to two interfaces, the both phases # This method populates the dictionary
# have, so it suffices to check for the wetting phase as it is always # self._interface_has_common_dof_indices
# present. # """
if self._interface_has_common_dof_indices[interface_ind]['wetting'].size != 0: # for interf_ind, interf_dofs in self._dof_indices_of_interface.items():
self._normal_vector_norms_for_common_dofs.update({interface_ind: dict()}) # # initialise bool numpy arry for the search.
for phase in self.has_phases: # if self.isRichards:
self._normal_vector_norms_for_common_dofs[interface_ind].update( # is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
{phase: dict()} # else:
) # is_part_of_neighbour_interface_w = np.zeros(interf_dofs['wetting'].size, dtype=bool)
for dof in self._interface_has_common_dof_indices[interface_ind][phase]: # is_part_of_neighbour_interface_nw = np.zeros(interf_dofs['nonwetting'].size, dtype=bool)
common_vertex_index = self.dof2vertex[phase][dof] #
for index, vertex in enumerate(df.vertices(self.mesh)): # for other_interf_ind, other_interf_dofs in self._dof_indices_of_interface.items():
if index == common_vertex_index: # if other_interf_ind is not interf_ind:
common_vertex = vertex # if self.isRichards:
# print(f"\non subdomain{self.subdomain_index}, interface{interface_ind} and phase {phase}") # is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
# print(f"vertex_{index} = {vertex.point()[:]}, type(vertex) = {type(vertex)}") # else:
linear_combination_of_normals = np.zeros(np.size(vertex.point()[:]), dtype=float) # is_part_of_neighbour_interface_w += np.isin(interf_dofs['wetting'], other_interf_dofs['wetting'], assume_unique=True)
# normals = [] # is_part_of_neighbour_interface_nw += np.isin(interf_dofs['nonwetting'], other_interf_dofs['nonwetting'], assume_unique=True)
# print(f"common_vertex = {common_vertex}", # # this carries a numpyarray of the indices which belong to other
# f"common_vertex_index = {common_vertex_index}") # # interfaces, too
for facet_containing_vertex in df.facets(common_vertex):#df.cells(common_vertex): # self._interface_has_common_dof_indices.update({interf_ind: dict()})
if facet_containing_vertex.exterior(): # if self.isRichards:
# normals.append(facet_containing_vertex.normal()[:]) # self._interface_has_common_dof_indices[interf_ind].update(
# print(f"type(facet_normal) = {type(facet_containing_vertex.normal()[:])}") # {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
linear_combination_of_normals += facet_containing_vertex.normal()[:] # )
# else:
# print(f"normals of common_vertex: {common_vertex} are ", normals) # self._interface_has_common_dof_indices[interf_ind].update(
# print(f"linear combination is: {linear_combination_of_normals[:]}.") # {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w],
# print(f"with length {np.linalg.norm(linear_combination_of_normals, ord=2, axis=None)}. For comparison, np.sqrt(2) = {np.sqrt(2)}.") # 'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
self._normal_vector_norms_for_common_dofs[interface_ind][phase].update( # )
{dof: np.linalg.norm(linear_combination_of_normals)} # # 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): def _calc_gravity_expressions(self):
...@@ -984,9 +984,9 @@ class DomainPatch(df.SubDomain): ...@@ -984,9 +984,9 @@ class DomainPatch(df.SubDomain):
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase) "pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
) )
file.write(error, i) file.write(error, i)
self.gli_function[phase].vector().set_local(self.gli_assembled[phase]) # self.gli_function[phase].vector().set_local(self.gli_assembled[phase])
self.gli_function[phase].rename("gli_function_"+"{}".format(phase),# # self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
"all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase)) # "all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
file.write(self.gli_function[phase], i) # file.write(self.gli_function[phase], i)
# END OF CLASS DomainPatch # END OF CLASS DomainPatch
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment