From 73807ab0626b3e0a446b1864bfeddcabb854ddf0 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 31 May 2019 15:35:46 +0200
Subject: [PATCH] implement general _normal_vector_norms_for_common_dofs with
linear combinations of normal vectors
---
LDDsimulation/domainPatch.py | 42 ++++++++++++++++++++++++++++++------
1 file changed, 35 insertions(+), 7 deletions(-)
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 6f873c3..d2de50b 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -391,7 +391,7 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
- def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
+ def calc_gl0_term(self, debug: bool = True) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl0 terms for all interfaces of the subdomain.
assemble the gl0 terms for all interfaces of the subdomain and save them
@@ -503,15 +503,21 @@ class DomainPatch(df.SubDomain):
# since we are looping over the interfaces we correct each dof
# that is common with another interfaces with the corresponding
# factor norm(n_1 + 1_2).
+ # We take the mean of the pressure on both interfaces, so
+ # we multiply these by 1/2
if dofs_in_common_with_other_interfaces[phase].size != 0:
if debug:
print(f"the following dofs (phase = {phase}) of interface{interf_ind} are in common with \n",
f" other interfaces: {dofs_in_common_with_other_interfaces[phase]}")
for common_dof in dofs_in_common_with_other_interfaces[phase]:
+ n_norm = self._normal_vector_norms_for_common_dofs[interf_ind][phase][common_dof]
# from the standpoint of the subdomain we are currently on,
# each dof can belong maximally to two interfaces.
# because we have not
- gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
+ neumann_flux_assembled[common_dof] = neumann_flux_assembled[common_dof]/n_norm
+ robin_pressure_assembled[common_dof] = robin_pressure_assembled[common_dof]/2
+
+ gl0_assembled = neumann_flux_assembled - robin_pressure_assembled
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
gli_assembled_tmp[phase] += gl0_assembled
@@ -652,12 +658,16 @@ class DomainPatch(df.SubDomain):
phi_a = self.testfunction[phase]
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
gli_pressure_part_array = gli_pressure_part.get_local()
- gli = gli_pressure_part_array - gli_prev_of_neighbour
# check if there are shared dofs with another interface.
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
+ # n_norm = self._normal_vector_norms_for_common_dofs[ind][phase][common_dof]
+ gli_prev_of_neighbour[common_dof] = gli_prev_of_neighbour[common_dof]/2
+
# see previous case.
- gli[common_dof] = 0.5*gli[common_dof]
+ gli_pressure_part_array[common_dof] = 0.5*gli_pressure_part_array[common_dof]
+
+ gli = gli_pressure_part_array - gli_prev_of_neighbour
# now, add the just assembled and corrected term to the Global
# gli_assembled_tmp vector
gli_assembled_tmp[phase] += gli
@@ -897,6 +907,7 @@ class DomainPatch(df.SubDomain):
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
@@ -908,10 +919,27 @@ class DomainPatch(df.SubDomain):
{phase: dict()}
)
for dof in self._interface_has_common_dof_indices[interface_ind][phase]:
- # at the moment we only determine this correctly if the patch
- # is rectangular and parallel to the x and y axis.
+ 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.sqrt(2)}
+ {dof: np.linalg.norm(linear_combination_of_normals)}
)
--
GitLab