From 8a2e7011e8a662234b3126074b5ef051fb85ce08 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Wed, 27 Mar 2019 08:30:42 +0100
Subject: [PATCH] fix Dirichlett boundary implementation
---
LDDsimulation.py | 39 ++++----
RR-2-patch-test.py | 110 +++++++++++++++-------
domainPatch.py | 225 ++++++++++++++++++++++++++++++++++++---------
3 files changed, 277 insertions(+), 97 deletions(-)
diff --git a/LDDsimulation.py b/LDDsimulation.py
index a6848c8..bc3064f 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -95,6 +95,7 @@ class LDDsimulation(object):
# this is actually an array of bools
isRichards: bool,#
interface_def_points: tp.List[tp.List[df.Point]],#
+ outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],#
adjacent_subdomains: tp.List[np.ndarray],#
mesh_resolution: float,#
viscosity: tp.Dict[int, tp.List[float]],#
@@ -103,7 +104,7 @@ class LDDsimulation(object):
lambda_param: tp.Dict[int, tp.Dict[str, float]],#
relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
saturation: tp.Dict[int, tp.Callable[...,None]],#
- timestep: tp.Dict[int, float],#
+ timestep_size: tp.Dict[int, float],#
sources: tp.Dict[int, tp.Dict[str, str]],#
initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
outer_dirichletBC: tp.Dict[int, tp.Dict[str, str]],#
@@ -115,6 +116,7 @@ class LDDsimulation(object):
self.isRichards = isRichards
self.mesh_resolution = mesh_resolution
self.interface_def_points = interface_def_points
+ self.outer_boundary_def_points = outer_boundary_def_points
self.adjacent_subdomains = adjacent_subdomains
self.viscosity = viscosity
self.porosity = porosity
@@ -122,7 +124,7 @@ class LDDsimulation(object):
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
self.saturation = saturation
- self.timestep = timestep
+ self.timestep_size = timestep_size
self.sources = sources
self.initial_conditions = initial_conditions
self.outer_dirichletBC = outer_dirichletBC
@@ -273,18 +275,19 @@ class LDDsimulation(object):
interface_list = self._get_interfaces(subdom_num)
self.subdomain.update(#
{subdom_num : dp.DomainPatch(#
- subdomain_index = subdom_num,#
- isRichards = isR,#
- mesh = self.mesh_subdomain[subdom_num],#
- viscosity = self.viscosity[subdom_num],#
- porosity = self.porosity[subdom_num],#
- interfaces = self.interface,#
- has_interface = interface_list,#
- L = self.L[subdom_num],#
- lambda_param = self.lambda_param[subdom_num],#
- relative_permeability = self.relative_permeability[subdom_num],#
- saturation = self.saturation[subdom_num],#
- timestep = self.timestep[subdom_num],#
+ subdomain_index = subdom_num,#
+ isRichards = isR,#
+ mesh = self.mesh_subdomain[subdom_num],#
+ viscosity = self.viscosity[subdom_num],#
+ porosity = self.porosity[subdom_num],#
+ outer_boundary_def_points = self.outer_boundary_def_points[subdom_num],#
+ interfaces = self.interface,#
+ has_interface = interface_list,#
+ L = self.L[subdom_num],#
+ lambda_param = self.lambda_param[subdom_num],#
+ relative_permeability = self.relative_permeability[subdom_num],#
+ saturation = self.saturation[subdom_num],#
+ timestep_size = self.timestep_size[subdom_num],#
)})
@@ -348,17 +351,17 @@ class LDDsimulation(object):
self.outerBC.update({num: dict()})
pDC = self.outer_dirichletBC[num]
- boundary_marker = subdomain.boundary_marker
+ boundary_marker = subdomain.outer_boundary_marker
if subdomain.isRichards:
# note that the default interpolation degree is 2
pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
- self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 0)})
+ self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
else:
pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
self.outerBC[num].update(#
- {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 0)},#
- {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 0)},#
+ {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
+ {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
)
def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index a765fad..8fbb509 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -23,11 +23,38 @@ sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
df.Point(1.0,1.0),
df.Point(0.0,1.0) ]
+
+# vertex coordinates of the outer boundaries. If it can not be specified as a
+# polygon, use an entry per boundary polygon. This information is used for defining
+# the Dirichlet boundary conditions. If a domain is completely internal, the
+# dictionary entry should be 0: None
+subdomain1_outer_boundary_verts = {
+ 0: [interface12_vertices[0], #
+ df.Point(0.0,1.0), #
+ df.Point(1.0,1.0), #
+ interface12_vertices[1]]
+}
# subdomain2
sub_domain2_vertices = [df.Point(0.0,0.0),
df.Point(1.0,0.0),
interface12_vertices[1],
interface12_vertices[0] ]
+
+subdomain2_outer_boundary_verts = {
+ 0: [interface12_vertices[1], #
+ df.Point(1.0,0.0), #
+ df.Point(0.0,0.0), #
+ interface12_vertices[0]]
+}
+# subdomain2_outer_boundary_verts = {
+# 0: [interface12_vertices[0], df.Point(0.0,0.0)],#
+# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
+# 2: [df.Point(1.0,0.0), interface12_vertices[1]]
+# }
+# subdomain2_outer_boundary_verts = {
+# 0: None
+# }
+
# list of subdomains given by the boundary polygon vertices.
# Subdomains are given as a list of dolfin points forming
# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
@@ -39,6 +66,13 @@ subdomain_def_points = [sub_domain0_vertices,#
sub_domain2_vertices]
# in the below list, index 0 corresponds to the 12 interface which has index 1
interface_def_points = [interface12_vertices]
+
+outer_boundary_def_points = {
+ # subdomain number
+ 1 : subdomain1_outer_boundary_verts,
+ 2 : subdomain2_outer_boundary_verts
+}
+
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
@@ -49,7 +83,7 @@ dt1 = 0.1
dt2 = dt1
# the timestep is taken as a dictionary to be able to run different timesteps on
# different subdomains.
-timestep = {
+timestep_size = {
1: dt1,
2: dt2
}
@@ -145,6 +179,7 @@ simulation.set_parameters(output_dir = "",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
interface_def_points = interface_def_points,#
+ outer_boundary_def_points = outer_boundary_def_points,#
adjacent_subdomains = adjacent_subdomains,#
mesh_resolution = 2,#
viscosity = viscosity,#
@@ -153,7 +188,7 @@ simulation.set_parameters(output_dir = "",#
lambda_param = lambda_param,#
relative_permeability = relative_permeability,#
saturation = sat_pressure_relationship,#
- timestep = timestep,#
+ timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
outer_dirichletBC = dirichletBC,#
@@ -179,32 +214,32 @@ mesh_subdomain = simulation.mesh_subdomain
interface = simulation.interface
interface_marker = simulation.interface_marker
-boundary_marker1 = simulation.subdomain[1].boundary_marker
-boundary_marker2 = simulation.subdomain[2].boundary_marker
+boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
+boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
dx_1 = simulation.subdomain[1].dx
dx_2 = simulation.subdomain[2].dx
ds_1 = simulation.subdomain[1].ds
ds_2 = simulation.subdomain[2].ds
-subdomain_boundary_marker1 = df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
-subdomain_boundary_marker2 = df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
-subdomain_boundary_marker1.set_all(0)
-subdomain_boundary_marker2.set_all(0)
-# print(dir(subdomain_boundary_marker1))
+interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
+interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
+# interface_marker1.set_all(0)
+# interface_marker2.set_all(0)
+# print(dir(interface_marker1))
-interface[0].mark(subdomain_boundary_marker1, 1)
-interface[0].mark(subdomain_boundary_marker2, 1)
-# subdomain_boundary_marker1 = simulation.subdomain[1].boundary_marker
+interface[0].mark(interface_marker1, 1)
+interface[0].mark(interface_marker2, 1)
+# interface_marker1 = simulation.subdomain[1].boundary_marker
# domain measures
# dx1 = simulation.subdomain[1].dx
-dx1 = df.Measure('dx', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1)
-dx2 = df.Measure('dx', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
+dx1 = df.Measure('dx', domain = mesh_subdomain[1])
+dx2 = df.Measure('dx', domain = mesh_subdomain[2])
# boundary measures
-ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1)
-ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
+ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = interface_marker1)
+ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = interface_marker2)
V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1)
p1_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2, t = 0)
@@ -229,18 +264,18 @@ dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
# interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
# interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
print("\non subdomain2: \n")
-# dom2_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker2, 1, print_coordinates = True)
-dom2_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
+# dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
+dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
-dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_marker1, 1)
+dofs_on_interface1 = interface[0].dofs_on_interface(V1, interface_marker1, 1)
print("Dofs on Interface dom1", dofs_on_interface1)
V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
-dofs_on_interface2 = interface[0].dofs_on_interface(V2, subdomain_boundary_marker2, 1)
+dofs_on_interface2 = interface[0].dofs_on_interface(V2, interface_marker2, 1)
print("Dofs on Interface dom1", dofs_on_interface2)
-# markermesh = subdomain_boundary_marker2.mesh()
+# markermesh = interface_marker2.mesh()
# functionSpaceMesh = V2.mesh()
# print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
# print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
@@ -295,8 +330,8 @@ for x in coordinates1:
i += 1
print("\non subdomain1: \n")
print("vertices of subdomain1: \n", coordinates1)
-# dom1_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker1, 1, print_coordinates = True)
-dom1_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker1, 1, print_vertex_indices = True)
+# dom1_interface_coordinates = interface[0].coordinates(interface_marker1, 1, print_coordinates = True)
+dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
interface[0].write_dofs(from_function = f1test, #
interface_dofs = dofs_on_interface1,#
@@ -319,8 +354,8 @@ for x in coordinates2:
print("\non subdomain2: \n")
print("vertices of subdomain2: \n", coordinates2)
-# dom2_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker2, 1, print_coordinates = True)
-dom2_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
+# dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
+dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
print("Parent indices", dom2_parent_vertex_indices)
@@ -359,11 +394,11 @@ p1_0 = df.interpolate(p1_initial, V1)
p2_0 = df.interpolate(p2_initial, V2)
# Initial boundary value for the pressure on interface
-p_gamma1 = df.DirichletBC(V1, df.Constant(0.5), subdomain_boundary_marker1, 1)
-p_gamma2 = df.DirichletBC(V2, df.Constant(0.1), subdomain_boundary_marker2, 1)
+p_gamma1 = df.DirichletBC(V1, df.Constant(0.5), interface_marker1, 1)
+p_gamma2 = df.DirichletBC(V2, df.Constant(0.1), interface_marker2, 1)
-outerBC1 = df.DirichletBC(V1, p1_exact, subdomain_boundary_marker1, 0)
-outerBC2 = df.DirichletBC(V2, p2_exact, subdomain_boundary_marker2, 0)
+outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
+outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
### source terms
@@ -388,16 +423,22 @@ g1_i = -lambda1*p1_i
g2_i = -lambda1*p2_i
for iterations in range(1):
- a1 = (L1*u1*v1)*dx1 + (timestep*df.dot(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx1 + (timestep*lambda1*u1*v1)*ds1
+ # + relative_permeability(saturation(p1_i, 1), 1)* +
+ a1 = (df.inner(df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1
A1 = df.assemble(a1)
+ dsform = (timestep*lambda1*u1*v1)*ds1(1)
+ dsform_vec = df.assemble(dsform)
b1 = df.assemble(rhs1)
p_gamma1.apply(A1,b1)
outerBC1.apply(A1,b1)
u = df.Function(V1)
U1 = u.vector()
df.solve(A1,U1,b1)
- print('solution:\n', U1)
+ print("Form:\n", A1.array())
+ print("RHS:\n", b1.get_local())
+ print("ds term:\n", dsform_vec.array())
+ print('solution:\n', U1.get_local())
p1_i.vector()[:] = U1
interface[0].write_dofs(from_function = p1_i, #
interface_dofs = dofs_on_interface1,#
@@ -408,10 +449,13 @@ for iterations in range(1):
#
# Save mesh to file
df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
df.File('./test_global_interface_marker.pvd') << interface_marker
df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
df.File('./test_domain_marker.pvd') << domain_marker
df.File('./test_domain_layered_soil_solution.pvd') << u
-df.File('./test_subdomain1_boundary_marker.pvd') << subdomain_boundary_marker1
-df.File('./test_subdomain2_boundary_marker.pvd') << subdomain_boundary_marker2
+df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
+df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
+df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
diff --git a/domainPatch.py b/domainPatch.py
index 5e298c6..3855a6d 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -108,6 +108,7 @@ class DomainPatch(df.SubDomain):
mesh: df.Mesh,#
porosity: float,#
viscosity: tp.Dict[str, float],#
+ outer_boundary_def_points: tp.Dict[str, tp.List[df.Point]],#
# tp.List[Interface] doesn't work, because it hasen't been defined yet.
interfaces: tp.List[object],#
has_interface: tp.List[int],#
@@ -115,7 +116,7 @@ class DomainPatch(df.SubDomain):
lambda_param: tp.Dict[str, float],#
relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
saturation: tp.Callable[..., None],#
- timestep: float,#
+ timestep_size: float,#
):
# because we declare our own __init__ method for the derived class BoundaryPart,
# we overwrite the __init__-method of the parent class df.SubDomain. However,
@@ -128,6 +129,7 @@ class DomainPatch(df.SubDomain):
self.mesh = mesh
self.porosity = porosity
self.viscosity = viscosity
+ self.outer_boundary_def_points = outer_boundary_def_points
self.interface = interfaces
self.has_interface = has_interface
self.L = L
@@ -135,11 +137,12 @@ class DomainPatch(df.SubDomain):
self.relative_permeability = relative_permeability
self.saturation = saturation
# timestep size, tau in the paper
- self.timestep = timestep
+ self.timestep_size = timestep_size
### Class variables set by methods
# dictionary holding the dof indices corresponding to an interface of
# given interface. see self._calc_dof_indices_of_interfaces()
self._dof_indices_of_interface = dict()
+ self.outer_boundary = []
# measures are set by self._init_measures()
self.iteration_number = 0
self.dx = None
@@ -153,9 +156,12 @@ class DomainPatch(df.SubDomain):
self.pressure_prev_timestep = None
# test function(s) for the form set by _init_function_space
self.TestFunction = None
+ # # dictionary which is initiated by _init_function_space holding the
+ # # the gli term needed for the form assemply.
+ # self.gli = None
self._init_function_space()
self._init_dof_and_vertex_maps()
- self._init_boundary_markers()
+ self._init_markers()
self._init_measures()
self._calc_dof_indices_of_interfaces()
@@ -188,7 +194,7 @@ class DomainPatch(df.SubDomain):
# define measures
dx = self.dx
ds = self.ds
- dt = self.timestep
+ dt = self.timestep_size
if self.isRichards:
Lw = self.L['wetting']
# this is pw_i in the paper
@@ -201,16 +207,19 @@ class DomainPatch(df.SubDomain):
kw = self.relative_permeability['wetting']
S = self.saturation
mu_w = self.viscosity['wetting']
+ source_w = self.source['wetting']
# we need to have all interfaces in the form
interface_forms = []
for index in self.has_interface:
interface_forms.append((dt*Lambda*pw*v)*ds[index])
form = (Lw*pw*v)*dx + (dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v)))*dx + df.sum(interface_forms)
# assemble rhs for
- gi = self._calc_gi(iteration = iter_num)
rhs_interface_forms = []
+ self._calc_gli_term(iteration = iter_num)
+ # gli = self.gli['wetting']
for index in self.has_interface:
- rhs_interface_forms.append((dt*gi*v)*ds[index])
+ gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
+ rhs_interface_forms.append((dt*gli*v)*ds[index])
rhs = (Lw*pw_prev_iter*v)*dx - ((S(pw_prev_iter) - S(pw_prev_timestep)*v)*dx + dt*source_w*v)*dx - df.sum(rhs_interface_forms)
form_and_rhs = {#
'form': form,#
@@ -240,6 +249,7 @@ class DomainPatch(df.SubDomain):
self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
self.pressure = df.TrialFunction(self.function_space['wetting'])
self.testfunction = df.TestFunction(self.function_space['wetting'])
+ # self.gli = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
@@ -253,6 +263,10 @@ class DomainPatch(df.SubDomain):
'wetting' : df.TestFunction(self.function_space['wetting']),
'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
}
+ # self.gli = {#
+ # 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
+ # 'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
+ # }
def _init_dof_and_vertex_maps(self) -> None:
""" calculate dof maps and the vertex to parent map.
@@ -282,18 +296,36 @@ class DomainPatch(df.SubDomain):
'nonwetting':df.vertex_to_dof_map(self.function_space['nonwetting'])#
}
- def _init_boundary_markers(self) -> None:
+ def _init_markers(self) -> None:
""" define boundary markers
This method sets the class Variables
- self.boundary_marker
+ self.interface_marker
+ self.outer_boundary_marker
"""
- self.boundary_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
- self.boundary_marker.set_all(0)
+ self.outer_boundary_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
+ self.outer_boundary_marker.set_all(0)
+ self.interface_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
+ self.interface_marker.set_all(0)
for glob_index in self.has_interface:
# each interface gets marked with the global interface index.
- self.interface[glob_index].mark(self.boundary_marker, glob_index)
+ self.interface[glob_index].mark(self.interface_marker, glob_index)
+ # create outerBoundary objects and mark outer boundaries with 1
+ for index, boundary_points in self.outer_boundary_def_points.items():
+ # if our domain patch has no outer boundary, this should be reflected
+ # by the value None in the dictionary
+ if boundary_points is None:
+ #
+ self.outer_boundary = None
+ else:
+ self.outer_boundary.append(#
+ BoundaryPart(vertices = boundary_points,#
+ internal = False,
+ tol=self.mesh.hmin()/100 )
+ )
+ self.outer_boundary[index].mark(self.outer_boundary_marker, 1)
+
def _init_measures(self):
""" define measures for the governing form
@@ -305,12 +337,12 @@ class DomainPatch(df.SubDomain):
# domain measure
self.dx = df.Measure('dx', domain = self.mesh)
# measure over the interfaces
- self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.boundary_marker)
+ self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.interface_marker)
def _calc_dof_indices_of_interfaces(self):
""" calculate dof indices of each interface """
V = self.function_space
- marker = self.boundary_marker
+ marker = self.interface_marker
for ind in self.has_interface:
self._dof_indices_of_interface.update(#
{ind: dict()}
@@ -325,8 +357,103 @@ class DomainPatch(df.SubDomain):
'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
)
-
- # END is_Richards
+ # def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
+ # """calculate the gl terms for the iteration number iteration
+ #
+ # calculate the gl terms for the iteration number iteration and save them
+ # to self.gli.
+ # """
+ # for ind in self.has_interface:
+ # # self._calc_gli_term_on(interface, iteration)
+ # interface = self.interface[ind]
+ # subdomain = self.subdomain_index
+ # # neighbour index
+ # neighbour = interface.neighbour[subdomain]
+ # neighbour_iter_num = interface.current_iteration[neighbour]
+ #
+ # Lambda = self.lambda_param
+ # if iteration == 0:
+ # n = df.FacetNormal(self.mesh)
+ # pw_prev_iter = self.pressure_prev_timestep
+ # k = self.relative_permeability
+ # S = self.saturation
+ # mu = self.viscosity
+ #
+ # if not neighbour_iter_num == 0:
+ # # save the gl term from the neighbour first to use it in the
+ # # next iteration
+ # if self.isRichards:
+ # interface.gli_term_prev[subdomain].update(#
+ # {'wetting': interface.gli_term[neighbour]['wetting']}
+ # )
+ # else:
+ # interface.gli_term_prev[subdomain].update(#
+ # {'wetting': interface.gli_term[neighbour]['wetting']},#
+ # {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+ # )
+ #
+ # if self.isRichards:
+ # mu_w = mu['wetting']
+ # kw = k['wetting']
+ # pwn_prev = pw_prev_iter['wetting']
+ # #calc gl0 from the flux term
+ # flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev)
+ # gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
+ # interface.gli_term[subdomain].update({'wetting': gl0})
+ # else:
+ # mu_w = mu['wetting']
+ # mu_nw = mu['nonwetting']
+ # pwn_prev = pw_prev_iter['wetting']
+ # pnwn_prev = pw_prev_iter['nonwetting']
+ # kw = k['wetting']
+ # knw = k['nonwetting']
+ # #calc gl0 from the flux term
+ # fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
+ # fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
+ # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
+ # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
+ # interface.gli_term[subdomain].update({'wetting': gwl0})
+ # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+ #
+ # interface.current_iteration[subdomain] += 1
+ #
+ # else:
+ # if neighbour_iter_num == iteration:
+ # if self.isRichards:
+ # interface.gli_term_prev[subdomain].update(#
+ # {'wetting': interface.gli_term[neighbour]['wetting']}
+ # )
+ #
+ # pwn_prev = pw_prev_iter['wetting']
+ # #calc gl0 from the flux term
+ # flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev)
+ # gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
+ # interface.gli_term[subdomain].update({'wetting': gl0})
+ # else:
+ # mu_w = mu['wetting']
+ # mu_nw = mu['nonwetting']
+ # pwn_prev = pw_prev_iter['wetting']
+ # pnwn_prev = pw_prev_iter['nonwetting']
+ # kw = k['wetting']
+ # knw = k['nonwetting']
+ # #calc gl0 from the flux term
+ # fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
+ # fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
+ # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
+ # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
+ # interface.gli_term[subdomain].update({'wetting': gwl0})
+ # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+ #
+ # else:
+ # print('')
+ # interface.current_iteration[subdomain] += 1
+ # # def _calc_gli_term_on(self, interface: int, iteration: int) -> None:
+ # # """calculate the gl terms for the iteration number iteration
+ # #
+ # # calculate the gl terms for the iteration number iteration and save them
+ # # to self.gli.
+ # # """
+ # # END is_Richards
# END OF CLASS DomainPatch
@@ -550,14 +677,21 @@ class Interface(BoundaryPart):
# method init_interface_values().
self.pressure_values = None
self.pressure_values_prev = None
- # same as for self.pressure but for the interface decoupling terms gl.
- self.gl_values = None
- self.gl_values_prev = None
+ # dictionaries holding the forms for the decoupling terms
+ self.gli_term = None
+ self.gli_term_prev = None
# dictionary holding the current iteration number i for each of the adjacent
# subdomains with index in adjacent_subdomains
- self.current_iteration = dict()
- self.current_iteration.update({self.adjacent_subdomains[0]: 0})
- self.current_iteration.update({self.adjacent_subdomains[1]: 0})
+ self.current_iteration = {#
+ self.adjacent_subdomains[0]: 0,#
+ self.adjacent_subdomains[1]: 0
+ }
+ # dictionary holding the neighbouring domain index relative to the key index
+ # this is used for calculating the gl_terms.
+ self.neighbour = {
+ self.adjacent_subdomains[0]: self.adjacent_subdomains[1],#
+ self.adjacent_subdomains[1]: self.adjacent_subdomains[0]#
+ }
# make the parent class methods available.
@@ -567,9 +701,6 @@ class Interface(BoundaryPart):
# def set_domain_index(self, domain_index: int):
# self.domain_index = domain_index
- def set_adjacent_subdomains(self, adjacent_subdomains: np.array):
- self.adjacent_subdomains = adjacent_subdomains
-
def dofs_on_interface(self,#
function_space: df.FunctionSpace,#
interface_marker: df.MeshFunction, #
@@ -623,25 +754,27 @@ class Interface(BoundaryPart):
self.pressure_values_prev.update({self.adjacent_subdomains[0]: dict()})
self.pressure_values_prev.update({self.adjacent_subdomains[1]: dict()})
- self.gl_values = dict()
- self.gl_values.update({self.adjacent_subdomains[0]: dict()})
- self.gl_values.update({self.adjacent_subdomains[1]: dict()})
- self.gl_values_prev = dict()
- self.gl_values_prev.update({self.adjacent_subdomains[0]: dict()})
- self.gl_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()})
- self.gl_values[subdom_ind].update({phase: dict()})
- self.gl_values_prev[subdom_ind].update({phase: dict()})
+ # the gli_term will hold forms not functions, therefor we leave
+ # the dictionary empty. it will be populated by domainPatch._calc_gli_term()
+ # 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.gl_values[subdom_ind][phase].update({vert_ind: 0.0})
- self.gl_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
def write_dofs(self, from_function: df.Function, #
@@ -661,8 +794,8 @@ class Interface(BoundaryPart):
of the interface
self.pressure_values[subdomain_ind][phase], (if pressure == True)
self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
- self.gl_values[subdomain_ind][phase], (if gl == True)
- self.gl_values_prev[subdomain_ind][phase] (if previous_iter == True)
+ self.gli_term[subdomain_ind][phase], (if gl == True)
+ self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
# Parameters
from_function: #type: df.Function, function to save the dofs to
@@ -706,14 +839,14 @@ class Interface(BoundaryPart):
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.gl_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- print("have written gl interface values\n", self.gl_values[subdomain_ind][phase])
+ self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+ print("have written gl interface values\n", self.gli_term[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.gl_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
- print("have written previous gl interface values\n", self.gl_values_prev[subdomain_ind][phase])
+ self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+ print("have written previous gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
def read_dofs(self, to_function: df.Function, #
interface_dofs: np.array,#
@@ -731,8 +864,8 @@ class Interface(BoundaryPart):
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)
- self.gl_values[subdomain_ind][phase], (if gl == True)
- self.gl_values_prev[subdomain_ind][phase] (if previous_iter == True)
+ self.gli_term[subdomain_ind][phase], (if gl == True)
+ self.gli_term_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.
.
@@ -782,15 +915,15 @@ class Interface(BoundaryPart):
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.gl_values[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read gl interface values\n", self.gl_values[subdomain_ind][phase])
+ to_function.vector()[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 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.gl_values_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
- print("read gl interface values\n", self.gl_values_prev[subdomain_ind][phase])
+ to_function.vector()[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 function \n", to_function.vector()[:])
#### Private Methods ####
--
GitLab