diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 60cbdf8ee4bfb11caaf784878684c34a3dff314b..509d1a29f54e55a81a47a26d0d401875bbf49c2c 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -6,6 +6,7 @@ import sympy as sym
import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
+import functools as ft
#import ufl as ufl
##### Domain and Interface ####
@@ -29,22 +30,83 @@ sub_domain2_vertices = [df.Point(0.0,0.0),
interface12_vertices[0] ]
# 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_vertices[i]) can be used
-# to create the subdomain. subdomain_vertices[0] contains the
-# vertices of the global simulation domain and subdomain_vertices[i] contains the
+# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
+# to create the subdomain. subdomain_def_points[0] contains the
+# vertices of the global simulation domain and subdomain_def_points[i] contains the
# vertices of the subdomain i.
-subdomain_vertices = [sub_domain0_vertices,#
+subdomain_def_points = [sub_domain0_vertices,#
sub_domain1_vertices,#
sub_domain2_vertices]
# in the below list, index 0 corresponds to the 12 interface which has index 1
-interface_vertices = [interface12_vertices]
+interface_def_points = [interface12_vertices]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
-# interface i (i.e. given by interface_vertices[i]).
+# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
+is_Richards = [True, True]
+viscosity = {#
+# subdom_num : viscosity
+ 1 : [1], #
+ 2 : [1]
+}
+
+porosity = {#
+# subdom_num : porosity
+ 1 : 1,#
+ 2 : 1
+}
+
+L = {#
+# subdom_num : subdomain L for L-scheme
+ 1 : [0.25],#
+ 2 : [0.25]
+}
+
+lambda_param = {#
+# subdom_num : lambda parameter for the L-scheme
+ 1 : [4],#
+ 2 : [4]
+}
+
+def rel_perm1(s):
+ # relative permeabilty on subdomain1
+ return s**2
+
+_rel_perm1 = ft.partial(rel_perm1)
+
+def rel_perm2(s):
+ # relative permeabilty on subdomain2
+ return s**3
+_rel_perm2 = ft.partial(rel_perm2)
+
+relative_permeability = {#
+ 1: _rel_perm1,
+ 2: _rel_perm2,
+}
+
+
+# def saturation(pressure, subdomain_index):
+# # inverse capillary pressure-saturation-relationship
+# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+#
+# sa
+
# initialise LDD simulation class
simulation = ldd.LDDsimulation()
-simulation.init_meshes_and_markers(subdomain_vertices, mesh_resolution=2)
+simulation.set_parameters(output_dir = "",#
+ subdomain_def_points = subdomain_def_points,#
+ is_Richards = is_Richards,#
+ interface_def_points = interface_def_points,#
+ adjacent_subdomains = adjacent_subdomains,#
+ mesh_resolution = 2,#
+ viscosity = viscosity,#
+ porosity = porosity,#
+ L = L,#
+ lambda_param = lambda_param,#
+ relative_permeability = relative_permeability)
+
+simulation.initialise()
+# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
# subdomain marker functions
domain_marker = simulation.domain_marker
mesh_subdomain = simulation.mesh_subdomain
@@ -52,232 +114,232 @@ mesh_subdomain = simulation.mesh_subdomain
# mesh = mesh_subdomain[0]
# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
# interface_marker.set_all(0)
-# interface = dp.Interface(vertices=interface_vertices[0], #
+# interface = dp.Interface(vertices=interface_def_points[0], #
# tol=mesh.hmin()/100,#
# internal=True,#
# adjacent_subdomains = adjacent_subdomains[0])
# print("interface vertices are: \n", interface.vertices)
# interface.mark(interface_marker, 1)
-simulation.init_interfaces(interface_vertices, adjacent_subdomains)
+# simulation._init_interfaces(interface_def_points, adjacent_subdomains)
interface = simulation.interface
interface_marker = simulation.interface_marker
-#
-# 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[0].mark(subdomain_boundary_marker1, 1)
-# interface[0].mark(subdomain_boundary_marker2, 1)
-#
-#
-# # domain measures
-# 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)
-#
-# # 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)
-#
-# 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)
-#
-# coordinates = mesh_subdomain[0].coordinates()
-# print("global mesh coordinates of domain: \n", coordinates)
-# coordinates1 = mesh_subdomain[1].coordinates()
+
+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[0].mark(subdomain_boundary_marker1, 1)
+interface[0].mark(subdomain_boundary_marker2, 1)
+
+
+# domain measures
+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)
+
+# 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)
+
+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)
+
+coordinates = mesh_subdomain[0].coordinates()
+print("global mesh coordinates of domain: \n", coordinates)
+coordinates1 = mesh_subdomain[1].coordinates()
+
+coordinates2 = mesh_subdomain[2].coordinates()
+submesh1_data = mesh_subdomain[1].data()
+dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
+print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
#
# coordinates2 = mesh_subdomain[2].coordinates()
-# submesh1_data = mesh_subdomain[1].data()
-# dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
-# print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
-# #
-# # coordinates2 = mesh_subdomain[2].coordinates()
-# # print("\nvertices of subdomain2: \n", coordinates2)
-# submesh2_data = mesh_subdomain[2].data()
-# dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
-# # print("'parent_vertex_indices',1 \n")
-# #
-# #
-# # print("on domain: \n")
-# # 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_vertices = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
-#
-#
-# dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_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)
-# print("Dofs on Interface dom1", dofs_on_interface2)
-#
-# # markermesh = subdomain_boundary_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())
-#
-# dom1_d2v = df.dof_to_vertex_map(V1)
-# dom1_v2d = df.vertex_to_dof_map(V1)
-# print("dof to vertex map for V1:\n", dom1_d2v)
-# print("vertex to dof map for V1:\n", dom1_v2d)
-# dom2_d2v = df.dof_to_vertex_map(V2)
-# dom2_v2d = df.vertex_to_dof_map(V2)
-# print("dof to vertex map for V2:\n", dom2_d2v)
-# print("vertex to dof map for V2:\n", dom2_v2d)
-#
-# testf1 = df.Constant(41)
-# f1test = df.interpolate(testf1, V1)
-# testf2 = df.Constant(42)
-# f2test = df.interpolate(testf2, V2)
-#
-# f1vec = f1test.compute_vertex_values()
-# for i in range(len(f1vec)):
-# f1vec[i] = i
-#
-# # f1test.vector()[:] = f1vec
-# f1nodal = f1test.vector().get_local()
-# # for i in range(len(f1nodal)):
-# # f1nodal[i] = i
-#
-# print("nodal values of f1nodal\n", f1nodal)
-#
-# #
-# # print('dofs on interface \n', f1vec[dofs_on_interface])
-# # f1test.vector()[:] = f1vec
-# # print("f1test.vector()[:]\n", f1test.vector()[:])
-# # print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
-#
-# print("constant function on dom1\n", f1vec)
-#
-# for i in dofs_on_interface1:
-# print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
-# # f1test.vector() is enumerated by dof not by vertex
-# f1test.vector()[i] = i
-# # this enumerates according to the dof numbering
-# print("testfunction f1test with set interface values\n", f1test.vector()[:])
-# print("testfunction f1test with get_local\n", f1test.vector().get_local())
-# # f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
-# # .vector()[:] and .vector().get_local() order values according to the vertex numbering
-# print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
-# print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
-# i=0
-# for x in coordinates1:
-# print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
-# 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_vertices = interface[0]._vertex_indices(subdomain_boundary_marker1, 1, print_vertex_indices = True)
-#
-# interface[0].write_dofs(fem_function = f1test, #
-# interface_dofs = dofs_on_interface1,#
-# dof_to_vert_map = dom1_d2v,#
-# local_to_parent_vertex_map = dom1_parent_vertex_indices)
-#
-# interface[0].read_dofs(fem_function = f2test, #
-# interface_dofs = dofs_on_interface2,#
-# dof_to_vert_map = dom2_d2v,#
-# local_to_parent_vertex_map = dom2_parent_vertex_indices)
-#
-# i=0
-# for x in coordinates2:
-# print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
-# i += 1
-# 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_vertices = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
-# print("Parent indices", dom2_parent_vertex_indices)
-#
-#
-# u1 = df.TrialFunction(V1)
-# v1 = df.TestFunction(V1)
-#
-# u2 = df.TrialFunction(V2)
-# v2 = df.TestFunction(V2)
-#
-# def relative_permeability(s, subdomain_index):
-# if subdomain_index == 1:
-# # relative permeabilty on subdomain1
-# return s**2
-# else:
-# # relative permeabilty on subdomain2
-# return s**3
-#
-# def saturation(pressure, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-#
-# # exact solution
-# # cell = df.Cell("triangle", 2)
-# # x1 = df.SpatialCoordinate(mesh_subdomain[1])
-# # x2 = df.SpatialCoordinate(mesh_subdomain[2])
-#
-# p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
-# p1_e = df.interpolate(p1_exact, V1)
-# p2_e = df.interpolate(p2_exact, V2)
-# #
-# # initial conditions
-# p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
-# p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
-# 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)
-#
-# outerBC1 = df.DirichletBC(V1, p1_exact, subdomain_boundary_marker1, 0)
-# outerBC2 = df.DirichletBC(V2, p2_exact, subdomain_boundary_marker2, 0)
-#
-#
-# ### source terms
-# source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
-# source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
-# source1h = df.interpolate(source1, V1)
-# source2h = df.interpolate(source2, V2)
-#
-# L1 = 0.25
-# L2 = L1
-# timestep = 0.1
-# lambda1 = 4
-# lambda2 = lambda1
-# # initial iteartion
-# p1_i = p1_0
-# p2_i = p2_0
-# # initial g_i
-# g1_i = -lambda1*p1_i
-# g2_i = -lambda1*p2_i
+# print("\nvertices of subdomain2: \n", coordinates2)
+submesh2_data = mesh_subdomain[2].data()
+dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
+# print("'parent_vertex_indices',1 \n")
+#
+#
+# print("on domain: \n")
+# 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)
+
+
+dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_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)
+print("Dofs on Interface dom1", dofs_on_interface2)
+
+# markermesh = subdomain_boundary_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())
+
+dom1_d2v = df.dof_to_vertex_map(V1)
+dom1_v2d = df.vertex_to_dof_map(V1)
+print("dof to vertex map for V1:\n", dom1_d2v)
+print("vertex to dof map for V1:\n", dom1_v2d)
+dom2_d2v = df.dof_to_vertex_map(V2)
+dom2_v2d = df.vertex_to_dof_map(V2)
+print("dof to vertex map for V2:\n", dom2_d2v)
+print("vertex to dof map for V2:\n", dom2_v2d)
+
+testf1 = df.Constant(41)
+f1test = df.interpolate(testf1, V1)
+testf2 = df.Constant(42)
+f2test = df.interpolate(testf2, V2)
+
+f1vec = f1test.compute_vertex_values()
+for i in range(len(f1vec)):
+ f1vec[i] = i
+
+# f1test.vector()[:] = f1vec
+f1nodal = f1test.vector().get_local()
+# for i in range(len(f1nodal)):
+# f1nodal[i] = i
+
+print("nodal values of f1nodal\n", f1nodal)
+
#
-# 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
-# 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)
-# 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)
-# p1_i.vector()[:] = U1
-# interface[0].write_dofs(fem_function = p1_i, #
-# interface_dofs = dofs_on_interface1,#
-# dof_to_vert_map = dom1_d2v,#
-# local_to_parent_vertex_map = dom1_parent_vertex_indices)
-# #
-# # Save mesh to file
+# print('dofs on interface \n', f1vec[dofs_on_interface])
+# f1test.vector()[:] = f1vec
+# print("f1test.vector()[:]\n", f1test.vector()[:])
+# print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
+
+print("constant function on dom1\n", f1vec)
+
+for i in dofs_on_interface1:
+ print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
+ # f1test.vector() is enumerated by dof not by vertex
+ f1test.vector()[i] = i
+# this enumerates according to the dof numbering
+print("testfunction f1test with set interface values\n", f1test.vector()[:])
+print("testfunction f1test with get_local\n", f1test.vector().get_local())
+# f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
+# .vector()[:] and .vector().get_local() order values according to the vertex numbering
+print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
+print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
+i=0
+for x in coordinates1:
+ print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
+ 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)
+
+interface[0].write_dofs(fem_function = f1test, #
+ interface_dofs = dofs_on_interface1,#
+ dof_to_vert_map = dom1_d2v,#
+ local_to_parent_vertex_map = dom1_parent_vertex_indices)
+
+interface[0].read_dofs(fem_function = f2test, #
+ interface_dofs = dofs_on_interface2,#
+ dof_to_vert_map = dom2_d2v,#
+ local_to_parent_vertex_map = dom2_parent_vertex_indices)
+
+i=0
+for x in coordinates2:
+ print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
+ i += 1
+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)
+print("Parent indices", dom2_parent_vertex_indices)
+
+
+u1 = df.TrialFunction(V1)
+v1 = df.TestFunction(V1)
+
+u2 = df.TrialFunction(V2)
+v2 = df.TestFunction(V2)
+
+def relative_permeability(s, subdomain_index):
+ if subdomain_index == 1:
+ # relative permeabilty on subdomain1
+ return s**2
+ else:
+ # relative permeabilty on subdomain2
+ return s**3
+
+def saturation(pressure, subdomain_index):
+ # inverse capillary pressure-saturation-relationship
+ return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+
+
+# exact solution
+# cell = df.Cell("triangle", 2)
+# x1 = df.SpatialCoordinate(mesh_subdomain[1])
+# x2 = df.SpatialCoordinate(mesh_subdomain[2])
+
+p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
+p1_e = df.interpolate(p1_exact, V1)
+p2_e = df.interpolate(p2_exact, V2)
+#
+# initial conditions
+p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
+p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
+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)
+
+outerBC1 = df.DirichletBC(V1, p1_exact, subdomain_boundary_marker1, 0)
+outerBC2 = df.DirichletBC(V2, p2_exact, subdomain_boundary_marker2, 0)
+
+
+### source terms
+source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
+source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
+source1h = df.interpolate(source1, V1)
+source2h = df.interpolate(source2, V2)
+
+L1 = 0.25
+L2 = L1
+timestep = 0.1
+lambda1 = 4
+lambda2 = lambda1
+# initial iteartion
+p1_i = p1_0
+p2_i = p2_0
+# initial g_i
+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
+ 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)
+ 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)
+ p1_i.vector()[:] = U1
+ interface[0].write_dofs(fem_function = p1_i, #
+ interface_dofs = dofs_on_interface1,#
+ dof_to_vert_map = dom1_d2v,#
+ local_to_parent_vertex_map = dom1_parent_vertex_indices)
+#
+# Save mesh to file
df.File('./test_domain_layered_soil.xml.gz') << 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.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