diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 364ff748e23cc1287cfd3c53a21f29ab7484046d..d5f6413a4b28286fc099a98ec09db244cb48de11 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -422,6 +422,10 @@ class LDDsimulation(object):
for num in range(1, self._number_of_subdomains+1):
self.mesh_subdomain.append(df.SubMesh(mesh, self.domain_marker, num))
+ if self.write2file['meshes_and_markers']:
+ filepath = self.output_dir+"domain_marker.pvd"
+ df.File(filepath) << self.domain_marker
+
def _init_interfaces(self, interface_def_points: tp.List[tp.List[df.Point]] = None,#
adjacent_subdomains: tp.List[np.ndarray] = None) -> None:
@@ -484,6 +488,10 @@ class LDDsimulation(object):
self.interface[num].mark(self.interface_marker, num+1)
self.interface[num].init_interface_values(self.interface_marker, num+1)
+ if self.write2file['meshes_and_markers']:
+ filepath = self.output_dir+"interface_marker_global.pvd"
+ df.File(filepath) << self.interface_marker
+
def _init_subdomains(self) -> None:
""" generate list of opjects of class DomainPatch.
"""
@@ -549,41 +557,24 @@ class LDDsimulation(object):
V = subdomain.function_space
# p0 contains both pw_0 and pnw_0 for subdomain[num]
p0 = self.initial_conditions[num]
- if subdomain.isRichards:
+ for phase in subdomain.has_phases:
# note that the default interpolation degree is 2
- pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
- subdomain.pressure = {
- 'wetting': df.Function(V['wetting'])
- }
- subdomain.pressure_prev_timestep = {
- 'wetting': df.Function(V['wetting'])
- }
- subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
- # print("vector()",subdomain.pressure_prev_iter['wetting'].vector().get_local())
- pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
- subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
- subdomain.pressure['wetting'].vector().set_local(pw_prev)
+ pa0 = df.Expression(p0[phase], domain = mesh, degree = interpolation_degree)
+ pa0_interpolated = df.interpolate(pa0, V[phase])
+ # alocate memory
+ subdomain.pressure.update(
+ {phase: df.Function(V[phase])}
+ )
+ subdomain.pressure_prev_iter.update(
+ {phase: df.Function(V[phase])}
+ )
+ subdomain.pressure_prev_timestep.update(
+ {phase: df.Function(V[phase])}
+ )
+ subdomain.pressure[phase].assign(pa0_interpolated)
+ subdomain.pressure_prev_iter[phase].assign(pa0_interpolated)
+ subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
- else:
- pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
- pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
- subdomain.pressure = {
- 'wetting': df.Function(V['wetting']),#
- 'nonwetting': df.Function(V['nonwetting'])
- }
- subdomain.pressure_prev_timestep = {
- 'wetting': df.Function(V['wetting']),
- 'nonwetting': df.Function(V['nonwetting'])
- }
- subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
- 'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
-
- pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
- pnw_prev = subdomain.pressure_prev_iter['nonwetting'].vector().get_local()
- subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
- subdomain.pressure['wetting'].vector().set_local(pw_prev)
- subdomain.pressure_prev_timestep['nonwetting'].vector().set_local(pnw_prev)
- subdomain.pressure['nonwetting'].vector().set_local(pnw_prev)
# populate interface dictionaries with pressure values.
subdomain.write_pressure_to_interfaces(initial_values = True)
@@ -611,33 +602,18 @@ class LDDsimulation(object):
pDC = self.dirichletBC_Expressions[num]
boundary_marker = subdomain.outer_boundary_marker
- if subdomain.isRichards:
+ for phase in subdomain.has_phases:
# note that the default interpolation degree is 2
if np.abs(time) < self.tol :
# time = 0 so the Dolfin Expression has to be created first.
- pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
- self.dirichletBC_dfExpression[num].update({'wetting': pDCw})
- else:
- pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
-
- self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
- else:
- if np.abs(time) < self.calc_tol :
- # time = 0 so the Dolfin Expression has to be created first.
- 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.dirichletBC_dfExpression[num].update(
- {'wetting': pDCw},#
- {'nonwetting': pDCnw},#
- )
+ pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
+ self.dirichletBC_dfExpression[num].update({phase: pDCa})
else:
- pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
- pDCnw = self.dirichletBC_dfExpression[num]['nonwetting'].t = time
+ pDCw = self.dirichletBC_dfExpression[num][phase].t = time
- self.outerBC[num].update(#
- {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
- {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
- )
+ self.outerBC[num].update(
+ {phase: df.DirichletBC(V[phase], pDCa, boundary_marker, 1)}
+ )
else:
# case subdomain.outer_boundary is None
print("subdomain is an inner subdomain and has no outer boundary.\n",
@@ -657,35 +633,12 @@ class LDDsimulation(object):
# p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
f = self.sources[subdomain_index][phase]
# note that the default interpolation degree is 2
- f_expression = df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
+ f_expression = {
+ phase: df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
+ }
# print(f_expression.__dir__())
# subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
- subdomain.source.update({phase : f_expression})
+ subdomain.source = f_expression.copy()
else:
# note that the default interpolation degree is 2
subdomain.source[phase].t = time
-
- # if np.abs(time) < self.calc_tol:
- # # here t = 0 and we have to initialise the sources.
- # mesh = subdomain.mesh
- # # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
- # f = self.sources[subdomain_index]
- # if subdomain.isRichards:
- # # note that the default interpolation degree is 2
- # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
- # # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
- # subdomain.source = {'wetting' : fw}
- # else:
- # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
- # fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
- # subdomain.source = {'wetting' : fw,#
- # 'nonwetting' : fnw}
- # else:
- # # here the sources are already initiated and stored in the subdomain.
- # # they must be updated in time.
- # if subdomain.isRichards:
- # # note that the default interpolation degree is 2
- # subdomain.source['wetting'].t = time
- # else:
- # subdomain.source['wetting'].t = time
- # subdomain.source['nonwetting'].t = time
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 32ecc6e3794b6f67936f0ca69ef3401b34af530e..1ba9132dfec044c7a4b6e68a7f897b87419e9cb0 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -159,14 +159,14 @@ class DomainPatch(df.SubDomain):
self.iteration_number = 0
self.dx = None
self.ds = None
- # solution function(s) for the form set by _init_function_space
- self.pressure = None
+ # solution function(s) for the form set by LDDsimulation._init_initial_values()
+ self.pressure = dict()
# this is being populated by LDDsimulation._eval_sources()
self.source = dict()
- # FEM function holding the previous iteration.
- self.pressure_prev_iter = None
- # FEM function holding the pressures of the previous timestep t_n.
- self.pressure_prev_timestep = None
+ # dictionary of FEM function for each phase holding the previous iteration.
+ self.pressure_prev_iter = dict()
+ # dictionary of FEM function for each phase holding the pressures of the previous timestep t_n.
+ self.pressure_prev_timestep = dict()
# test function(s) for the form set by _init_function_space
self.testfunction = None
# Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 2b8eeaedfec4cc34f8109bc46746955e06314c0f..03678bd0ccea5ed98ec41702cd42d8ca8148d452 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -213,9 +213,11 @@ simulation.set_parameters(output_dir = "./",#
simulation.initialise()
# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
# subdomain marker functions
+output_dir = simulation.output_dir
domain_marker = simulation.domain_marker
mesh_subdomain = simulation.mesh_subdomain
+
# mesh = mesh_subdomain[0]
# interface_marker = df.MeshFunction('int', mesh, mesh.topology().dim()-1)
# interface_marker.set_all(0)
@@ -232,13 +234,6 @@ interface_marker = simulation.interface_marker
subdoms = simulation.subdomain
-# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-df.File('./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('./domain_marker.pvd') << domain_marker
-# df.File('./test_domain_layered_soil_solution.pvd') << u
df.File('./subdomain1_interface_marker.pvd') << simulation.subdomain[1].interface_marker
df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1