From 0c527b8aaa194ec82ccb52749eb653564effb369 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 2 Apr 2019 13:37:59 +0200
Subject: [PATCH] reorganise repo
---
.../LDDsimulation.py | 0
LDDsimulation/__init__.py | 0
.../domainPatch.py | 2 +
helpers.py => LDDsimulation/helpers.py | 0
RR-2-patch-test-case/RR-2-patch-test.py | 530 ++++++++++++++++++
RR-2-patch-test.py | 478 ----------------
Test_misc/taaaeeeeeeschd.py | 45 ++
.../layered_soil.py | 0
8 files changed, 577 insertions(+), 478 deletions(-)
rename LDDsimulation.py => LDDsimulation/LDDsimulation.py (100%)
create mode 100644 LDDsimulation/__init__.py
rename domainPatch.py => LDDsimulation/domainPatch.py (99%)
rename helpers.py => LDDsimulation/helpers.py (100%)
create mode 100755 RR-2-patch-test-case/RR-2-patch-test.py
delete mode 100755 RR-2-patch-test.py
create mode 100755 Test_misc/taaaeeeeeeschd.py
rename layered_soil.py => layered-soil-case/layered_soil.py (100%)
diff --git a/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
similarity index 100%
rename from LDDsimulation.py
rename to LDDsimulation/LDDsimulation.py
diff --git a/LDDsimulation/__init__.py b/LDDsimulation/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/domainPatch.py b/LDDsimulation/domainPatch.py
similarity index 99%
rename from domainPatch.py
rename to LDDsimulation/domainPatch.py
index 1f69ce3..47d72f1 100644
--- a/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -142,6 +142,8 @@ class DomainPatch(df.SubDomain):
# 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()
+ # List of objects of clas outerBoundary initialised by self._init_markers()
+ # can be NONE if no outer boundary is present.
self.outer_boundary = []
# measures are set by self._init_measures()
self.iteration_number = 0
diff --git a/helpers.py b/LDDsimulation/helpers.py
similarity index 100%
rename from helpers.py
rename to LDDsimulation/helpers.py
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
new file mode 100755
index 0000000..12938f2
--- /dev/null
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -0,0 +1,530 @@
+#!/usr/bin/python3
+import dolfin as df
+import mshr
+import numpy as np
+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 ####
+# global simulation domain domain
+sub_domain0_vertices = [df.Point(0.0,0.0), #
+ df.Point(1.0,0.0),#
+ df.Point(1.0,1.0),#
+ df.Point(0.0,1.0)]
+# interface between subdomain1 and subdomain2
+interface12_vertices = [df.Point(0.0, 0.5),
+ df.Point(1.0, 0.5) ]
+# subdomain1.
+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
+# 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_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_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]]
+isRichards = {
+ 1: True, #
+ 2: True
+ }
+
+Tmax = 2
+dt1 = 0.1
+dt2 = dt1
+# the timestep is taken as a dictionary to be able to run different timesteps on
+# different subdomains.
+timestep_size = {
+ 1: dt1,
+ 2: dt2
+}
+
+viscosity = {#
+# subdom_num : viscosity
+ 1 : {'wetting' :1}, #
+ 2 : {'wetting' :1}
+}
+
+porosity = {#
+# subdom_num : porosity
+ 1 : 1,#
+ 2 : 1
+}
+
+L = {#
+# subdom_num : subdomain L for L-scheme
+ 1 : {'wetting' :0.25},#
+ 2 : {'wetting' :0.25}
+}
+
+lambda_param = {#
+# subdom_num : lambda parameter for the L-scheme
+ 1 : {'wetting' :4},#
+ 2 : {'wetting' :4}
+}
+
+## relative permeabilty functions on subdomain 1
+def rel_perm1(s):
+ # relative permeabilty on subdomain1
+ return s**2
+
+_rel_perm1 = ft.partial(rel_perm1)
+
+subdomain1_rel_perm = {
+ 'wetting': _rel_perm1,#
+ 'nonwetting': None
+}
+## relative permeabilty functions on subdomain 2
+def rel_perm2(s):
+ # relative permeabilty on subdomain2
+ return s**3
+_rel_perm2 = ft.partial(rel_perm2)
+
+subdomain2_rel_perm = {
+ 'wetting': _rel_perm2,#
+ 'nonwetting': None
+}
+
+## dictionary of relative permeabilties on all domains.
+relative_permeability = {#
+ 1: subdomain1_rel_perm,
+ 2: subdomain2_rel_perm
+}
+
+def saturation(pressure, subdomain_index):
+ # inverse capillary pressure-saturation-relationship
+ return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+
+sat_pressure_relationship = {#
+ 1: ft.partial(saturation, subdomain_index = 1),#
+ 2: ft.partial(saturation, subdomain_index = 2)
+}
+
+source_expression = {
+ 1: {'wetting': '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]) )'},
+ 2: {'wetting': '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))'}
+}
+
+initial_condition = {
+ 1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
+ 2: {'wetting': '-x[1]*x[1]'}
+}
+
+exact_solution = {
+ 1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
+ 2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+}
+
+dirichletBC = exact_solution
+
+# def saturation(pressure, subdomain_index):
+# # inverse capillary pressure-saturation-relationship
+# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+#
+# sa
+
+mesh_resolution = 3
+
+# initialise LDD simulation class
+simulation = ldd.LDDsimulation()
+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 = mesh_resolution,#
+ viscosity = viscosity,#
+ porosity = porosity,#
+ L = L,#
+ lambda_param = lambda_param,#
+ relative_permeability = relative_permeability,#
+ saturation = sat_pressure_relationship,#
+ timestep_size = timestep_size,#
+ sources = source_expression,#
+ initial_conditions = initial_condition,#
+ outer_dirichletBC = dirichletBC,#
+ )
+
+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
+
+# 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_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_def_points, adjacent_subdomains)
+
+interface = simulation.interface
+interface_marker = simulation.interface_marker
+
+subdoms = simulation.subdomain
+
+for iterations in range(1):
+ #
+
+ a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*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(1)
+ rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+ extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+ splitrhs = df.assemble(rhssplit1)
+ print("splitrhs: \n", splitrhs.get_local())
+ extrat_assembled = df.assemble(extratrem)
+ print("extratrem: \n", extrat_assembled.get_local())
+ added = splitrhs + extrat_assembled
+ print("both added together:\n", added.get_local())
+
+ 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("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].read_pressure_dofs(from_function = p1_i, #
+ interface_dofs = dofs_on_interface1,#
+ dof_to_vert_map = dom1_d2v,#
+ local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+ phase = 'wetting',#
+ subdomain_ind = 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_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
+
+# 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
+#
+# 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(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])
+# dx2 = df.Measure('dx', domain = mesh_subdomain[2])
+#
+# # boundary measures
+# 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)
+#
+# 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()
+# # 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(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, 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, interface_marker2, 1)
+# # print("Dofs on Interface dom1", dofs_on_interface2)
+#
+# # 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())
+#
+# 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(interface_marker1, 1, print_coordinates = True)
+# dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
+#
+# interface[0].read_pressure_dofs(from_function = f1test, #
+# interface_dofs = dofs_on_interface1,#
+# dof_to_vert_map = dom1_d2v,#
+# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+# phase = 'wetting',
+# subdomain_ind = 1)
+#
+# interface[0].write_pressure_dofs(to_function = f2test, #
+# interface_dofs = dofs_on_interface2,#
+# dof_to_vert_map = dom2_d2v,#
+# local_to_parent_vertex_map = dom2_parent_vertex_indices,#
+# phase = 'wetting',
+# subdomain_ind = 2)
+#
+# 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(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)
+#
+#
+# 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, p1_exact, interface_marker1, 1)
+# p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
+#
+# outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
+# outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
+#
+#
+# ### 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)
+# source1 = df.Expression(source_expression[1]['wetting'], 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
+# n1 = df.FacetNormal(mesh_subdomain[1])
+# flux = -df.grad(p1_0)
+# gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
+# print('gl0 term assembliert', df.assemble(gl0).get_local())
+# g1_i = -lambda1*p1_i
+# g2_i = -lambda1*p2_i
+
+# for iterations in range(1):
+# #
+#
+# a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*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(1)
+# rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+# extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+# splitrhs = df.assemble(rhssplit1)
+# print("splitrhs: \n", splitrhs.get_local())
+# extrat_assembled = df.assemble(extratrem)
+# print("extratrem: \n", extrat_assembled.get_local())
+# added = splitrhs + extrat_assembled
+# print("both added together:\n", added.get_local())
+#
+# 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("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].read_pressure_dofs(from_function = p1_i, #
+# interface_dofs = dofs_on_interface1,#
+# dof_to_vert_map = dom1_d2v,#
+# local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+# phase = 'wetting',#
+# subdomain_ind = 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_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/RR-2-patch-test.py b/RR-2-patch-test.py
deleted file mode 100755
index df0fb80..0000000
--- a/RR-2-patch-test.py
+++ /dev/null
@@ -1,478 +0,0 @@
-#!/usr/bin/python3
-import dolfin as df
-import mshr
-import numpy as np
-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 ####
-# global simulation domain domain
-sub_domain0_vertices = [df.Point(0.0,0.0), #
- df.Point(1.0,0.0),#
- df.Point(1.0,1.0),#
- df.Point(0.0,1.0)]
-# interface between subdomain1 and subdomain2
-interface12_vertices = [df.Point(0.0, 0.5),
- df.Point(1.0, 0.5) ]
-# subdomain1.
-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
-# 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_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_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]]
-isRichards = {
- 1: True, #
- 2: True
- }
-
-Tmax = 2
-dt1 = 0.1
-dt2 = dt1
-# the timestep is taken as a dictionary to be able to run different timesteps on
-# different subdomains.
-timestep_size = {
- 1: dt1,
- 2: dt2
-}
-
-viscosity = {#
-# subdom_num : viscosity
- 1 : {'wetting' :1}, #
- 2 : {'wetting' :1}
-}
-
-porosity = {#
-# subdom_num : porosity
- 1 : 1,#
- 2 : 1
-}
-
-L = {#
-# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.25},#
- 2 : {'wetting' :0.25}
-}
-
-lambda_param = {#
-# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :4},#
- 2 : {'wetting' :4}
-}
-
-## relative permeabilty functions on subdomain 1
-def rel_perm1(s):
- # relative permeabilty on subdomain1
- return s**2
-
-_rel_perm1 = ft.partial(rel_perm1)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1,#
- 'nonwetting': None
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2(s):
- # relative permeabilty on subdomain2
- return s**3
-_rel_perm2 = ft.partial(rel_perm2)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2,#
- 'nonwetting': None
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
-}
-
-def saturation(pressure, subdomain_index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-
-sat_pressure_relationship = {#
- 1: ft.partial(saturation, subdomain_index = 1),#
- 2: ft.partial(saturation, subdomain_index = 2)
-}
-
-source_expression = {
- 1: {'wetting': '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]) )'},
- 2: {'wetting': '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))'}
-}
-
-initial_condition = {
- 1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
- 2: {'wetting': '-x[1]*x[1]'}
-}
-
-exact_solution = {
- 1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
- 2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
-}
-
-dirichletBC = exact_solution
-
-# def saturation(pressure, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-mesh_resolution = 3
-
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation()
-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 = mesh_resolution,#
- viscosity = viscosity,#
- porosity = porosity,#
- L = L,#
- lambda_param = lambda_param,#
- relative_permeability = relative_permeability,#
- saturation = sat_pressure_relationship,#
- timestep_size = timestep_size,#
- sources = source_expression,#
- initial_conditions = initial_condition,#
- outer_dirichletBC = dirichletBC,#
- )
-
-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
-
-# 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_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_def_points, adjacent_subdomains)
-
-interface = simulation.interface
-interface_marker = simulation.interface_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
-
-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(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])
-dx2 = df.Measure('dx', domain = mesh_subdomain[2])
-
-# boundary measures
-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)
-
-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()
-# 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(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, 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, interface_marker2, 1)
-# print("Dofs on Interface dom1", dofs_on_interface2)
-
-# 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())
-
-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(interface_marker1, 1, print_coordinates = True)
-dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
-
-interface[0].write_pressure_dofs(from_function = f1test, #
- interface_dofs = dofs_on_interface1,#
- dof_to_vert_map = dom1_d2v,#
- local_to_parent_vertex_map = dom1_parent_vertex_indices,#
- phase = 'wetting',
- subdomain_ind = 1)
-
-interface[0].read_pressure_dofs(to_function = f2test, #
- interface_dofs = dofs_on_interface2,#
- dof_to_vert_map = dom2_d2v,#
- local_to_parent_vertex_map = dom2_parent_vertex_indices,#
- phase = 'wetting',
- subdomain_ind = 2)
-
-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(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)
-
-
-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, p1_exact, interface_marker1, 1)
-p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
-
-outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
-outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
-
-
-### 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)
-source1 = df.Expression(source_expression[1]['wetting'], 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
-n1 = df.FacetNormal(mesh_subdomain[1])
-flux = -df.grad(p1_0)
-gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
-print('gl0 term assembliert', df.assemble(gl0).get_local())
-g1_i = -lambda1*p1_i
-g2_i = -lambda1*p2_i
-
-for iterations in range(1):
- # + +
- a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*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(1)
- rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
- extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
- splitrhs = df.assemble(rhssplit1)
- print("splitrhs: \n", splitrhs.get_local())
- extrat_assembled = df.assemble(extratrem)
- print("extratrem: \n", extrat_assembled.get_local())
- added = splitrhs + extrat_assembled
- print("both added together:\n", added.get_local())
-
- 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("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_pressure_dofs(from_function = p1_i, #
- interface_dofs = dofs_on_interface1,#
- dof_to_vert_map = dom1_d2v,#
- local_to_parent_vertex_map = dom1_parent_vertex_indices,#
- phase = 'wetting',#
- subdomain_ind = 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_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/Test_misc/taaaeeeeeeschd.py b/Test_misc/taaaeeeeeeschd.py
new file mode 100755
index 0000000..3d07eb3
--- /dev/null
+++ b/Test_misc/taaaeeeeeeschd.py
@@ -0,0 +1,45 @@
+#!/usr/bin/python3
+
+import dolfin as df
+
+mesh = df.UnitIntervalMesh(2)
+
+dx = df.dx
+ds = df.ds
+
+V = df.FunctionSpace(mesh, 'CG', 1)
+
+v = df.TestFunction(V)
+u = df.Function(V)
+
+# set k constant to 42
+# k = df.Constant(42)
+k = df.Expression('2*x[0]', domain = mesh, degree=1)
+k = df.interpolate(k, V)
+
+n = df.FacetNormal(mesh)
+print(k.vector().get_local())
+print(df.assemble(k*v*dx).get_local())
+print(df.assemble(df.dot(df.grad(k),n)*v*ds).get_local())
+
+
+# solve k*u'' = 0
+a = u.dx(0)*k*v.dx(0)*dx + u.dx(0)*k.dx(0)*v*dx
+
+print('Form: \n', df.assemble(a).get_local())
+bcL = df.DirichletBC(V,df.Constant(10.),"on_boundary && x[0] < 0.5")
+bcR = df.DirichletBC(V,df.Constant(1.),"on_boundary && x[0] > 0.5")
+bcs= [bcL,bcR]
+
+df.solve(a == 0,u,bcs=bcs)
+
+df.File('./taaeschd_u1.pvd') << u
+
+# change value of k at one boundary point
+k.vector()[0] = 1.e6
+print(k.vector().get_local())
+print(df.assemble(k*dx))
+print(df.assemble(k*ds))
+
+df.solve(a == 0,u,bcs=bcs)
+df.File('./taaeschd_u2.pvd') << u
diff --git a/layered_soil.py b/layered-soil-case/layered_soil.py
similarity index 100%
rename from layered_soil.py
rename to layered-soil-case/layered_soil.py
--
GitLab