diff --git a/domainPatch.py b/domainPatch.py
index 3bb5286e61d32b44e8bba94774ef4958da87addd..b67873cbda2dad5c04390f129510b51a71555170 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -163,6 +163,11 @@ class DomainPatch(df.SubDomain):
# which is given by self.govering_problem().
#
self.gli_assembled = None
+ # Dictionary of FEM function of self.function_space holding the interface
+ # dofs of the previous iteration of the neighbouring pressure. This is used
+ # to assemble the gli terms in the method self.calc_gli_term().
+ self.neighbouring_interface_pressure = None
+
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_markers()
@@ -176,7 +181,7 @@ class DomainPatch(df.SubDomain):
""" save the interface values of self.pressure to all neighbouring interfaces"""
if self.isRichards:
for ind in self.has_interface:
- interface[ind].write_pressure_dofs(from_function = self.pressure['wetting'], #
+ interface[ind].read_pressure_dofs(from_function = self.pressure['wetting'], #
interface_dofs = self._dof_indices_of_interface[ind]['wetting'],#
dof_to_vert_map = self.dof2vertex['wetting'],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -185,7 +190,7 @@ class DomainPatch(df.SubDomain):
else:
for ind in self.has_interface:
for phase in ['wetting', 'nonwetting']:
- interface[ind].write_pressure_dofs(from_function = self.pressure[phase], #
+ interface[ind].read_pressure_dofs(from_function = self.pressure[phase], #
interface_dofs = self._dof_indices_of_interface[ind][phase],#
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -332,7 +337,7 @@ class DomainPatch(df.SubDomain):
neighbour_iter_num = interface.current_iteration[neighbour]
dt = self.timestep_size
ds = self.ds(ind)
- # needed for the write_gli_dofs() functions
+ # needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind],#
Lambda = self.lambda_param
@@ -365,10 +370,10 @@ class DomainPatch(df.SubDomain):
flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
self.gli_assembled = {
- 'wetting': df.assemble(dt*gl0*phi_w*ds)
+ 'wetting': df.assemble(dt*gl0*phi_w*ds).get_local()
}
# write glo to the current iteration dictionary of the interface
- interface.write_gli_dofs(from_array = self.gli_assembled['wetting'], #
+ interface.read_gli_dofs(from_array = self.gli_assembled['wetting'], #
interface_dofs = interface_dofs['wetting'],#
dof_to_vert_map = self.dof2vertex['wetting'],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -391,9 +396,9 @@ class DomainPatch(df.SubDomain):
fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
self.gli_assembled = {
- 'wetting': df.assemble(dt*gwl0*phi_w*ds),
+ 'wetting': df.assemble(dt*gwl0*phi_w*ds).get_local(),
# assemble the nonwetting phase gl0 term as zero functional
- 'nonwetting': df.assemble(df.Constant(0)*ds)
+ 'nonwetting': df.assemble(df.Constant(0)*ds).get_local()
}
else:
# here we need to assemble a gli term for the nonwetting phase
@@ -406,12 +411,12 @@ class DomainPatch(df.SubDomain):
gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
self.gli_assembled = {
- 'wetting': df.assemble(dt*gwl0*phi_w*ds),
- 'nonwetting': df.assemble(dt*gnwl0*phi_nw*ds)
+ 'wetting': df.assemble(dt*gwl0*phi_w*ds).get_local(),
+ 'nonwetting': df.assemble(dt*gnwl0*phi_nw*ds).get_local()
}
# write glo to the current iteration dictionary of the interface
# wetting
- interface.write_gli_dofs(from_array = self.gli_assembled['wetting'], #
+ interface.read_gli_dofs(from_array = self.gli_assembled['wetting'], #
interface_dofs = interface_dofs['wetting'],#
dof_to_vert_map = self.dof2vertex['wetting'],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -420,7 +425,7 @@ class DomainPatch(df.SubDomain):
previous_iter = False
)
# nonwetting
- interface.write_gli_dofs(from_array = self.gli_assembled['nonwetting'], #
+ interface.read_gli_dofs(from_array = self.gli_assembled['nonwetting'], #
interface_dofs = interface_dofs['nonwetting'],#
dof_to_vert_map = self.dof2vertex['nonwetting'],#
local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -430,37 +435,60 @@ class DomainPatch(df.SubDomain):
)
interface.current_iteration[subdomain] += 1
- # else: # iteration > 0
- # if neighbour_iter_num == iteration:
- # if self.isRichards:
- # interface.gli_term_prev[subdomain].update(#
- # {'wetting': interface.gli_term[neighbour]['wetting']}
- # )
- #
- # pw_prev = pw_prev_iter['wetting']
- # #calc gl0 from the flux term
- # flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
- # gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
- # interface.gli_term[subdomain].update({'wetting': gl0})
- # else:
- # mu_w = mu['wetting']
- # mu_nw = mu['nonwetting']
- # pw_prev = pw_prev_iter['wetting']
- # pnw_prev = pw_prev_iter['nonwetting']
- # kw = k['wetting']
- # knw = k['nonwetting']
- # #calc gl0 from the flux term
- # fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
- # fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
- # gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
- # gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
- # interface.gli_term[subdomain].update({'wetting': gwl0})
- # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
- #
- # else:
- # print('uiae')
- #
- # interface.current_iteration[subdomain] += 1
+ else: # iteration > 0
+ if neighbour_iter_num == iteration:
+ if self.isRichards:
+ # first, save the gli_terms of the neighbour to gli_prev_term
+ # of the current subdomain.
+ interface.gli_term_prev[subdomain].update(#
+ {'wetting': interface.gli_term[neighbour]['wetting']}
+ )
+ # then read gli_prev term from interface and write it to
+ # the subdomain gli_assembled array
+ interface.write_gli_dofs(to_array = self.gli_assembled['wetting'], #
+ interface_dofs = interface_dofs['wetting'],#
+ dof_to_vert_map = self.dof2vertex['wetting'],#
+ local_to_parent_vertex_map = self.parent_mesh_index,#
+ phase = 'wetting',#
+ subdomain_ind = subdomain,#
+ previous_iter = True
+ )
+
+ # read the neighbouring pressure values and write them to
+ # the placeholder fucntion self.neighbouring_interface_pressure
+ interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure['wetting'], #
+ interface_dofs = interface_dofs['wetting'],#
+ dof_to_vert_map = self.dof2vertex['wetting'],#
+ local_to_parent_vertex_map = self.parent_mesh_index,#
+ phase = 'wetting',#
+ subdomain_ind = subdomain,#
+ previous_iter = True)
+
+
+ pw_prev = pw_prev_iter['wetting']
+ #calc gl0 from the flux term
+ flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
+ gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
+ interface.gli_term[subdomain].update({'wetting': gl0})
+ else:
+ mu_w = mu['wetting']
+ mu_nw = mu['nonwetting']
+ pw_prev = pw_prev_iter['wetting']
+ pnw_prev = pw_prev_iter['nonwetting']
+ kw = k['wetting']
+ knw = k['nonwetting']
+ #calc gl0 from the flux term
+ fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
+ fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
+ gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+ gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
+ interface.gli_term[subdomain].update({'wetting': gwl0})
+ interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+
+ else:
+ print('uiae')
+
+ interface.current_iteration[subdomain] += 1
#### PRIVATE METHODS
@@ -477,7 +505,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_assembled = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
+ self.neighbouring_interface_pressure = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
else:
self.function_space = {#
'wetting' : df.FunctionSpace(self.mesh, 'P', 1),#
@@ -491,10 +519,10 @@ class DomainPatch(df.SubDomain):
'wetting' : df.TestFunction(self.function_space['wetting']),
'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
}
- # self.gli_assembled = {#
- # 'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
- # 'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
- # }
+ self.neighbouring_interface_pressure = {#
+ '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.
@@ -914,7 +942,7 @@ class Interface(BoundaryPart):
self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
# end init_interface_values
- def write_pressure_dofs(self, from_function: df.Function, #
+ def read_pressure_dofs(self, from_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
@@ -922,10 +950,10 @@ class Interface(BoundaryPart):
subdomain_ind: int,#
previous_iter: bool = False,#
) -> None:
- """ write dofs of from_function with indices interface_dofs to self.pressure_values
+ """ read dofs of from_function with indices interface_dofs and write to self.pressure_values
- Write the degrees of freedom of the function from_function with indices
- interface_dofs, i.e. the dofs of the interface, to one of the dictionaries
+ Read the degrees of freedom of the function from_function with indices
+ interface_dofs, i.e. the dofs of the interface, and write them to 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)
@@ -963,7 +991,7 @@ class Interface(BoundaryPart):
print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
- def write_gli_dofs(self, from_array: np.array, #
+ def read_gli_dofs(self, from_array: np.array, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
@@ -1012,7 +1040,7 @@ class Interface(BoundaryPart):
print("have written previous gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
- def read_pressure_dofs(self, to_function: df.Function, #
+ def write_pressure_dofs(self, to_function: df.Function, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
@@ -1070,7 +1098,7 @@ class Interface(BoundaryPart):
print("wrote these dofs to function \n", to_function.vector()[:])
- def read_gli_dofs(self, to_gli_array: np.array, #
+ def write_gli_dofs(self, to_array: np.array, #
interface_dofs: np.array,#
dof_to_vert_map: np.array,#
local_to_parent_vertex_map: np.array,#
@@ -1079,17 +1107,17 @@ class Interface(BoundaryPart):
previous_iter: bool = False,#
) -> None:
""" read dofs off self.gli_term and overwrite the dofs stored in the array
- to_gli_array with indices interface_dofs
+ to_array with indices interface_dofs
Read the degrees of freedom of one of the dictionaries of the interface
self.gli_term[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
- and overwrite the values (holding dofs) of the array to_gli_array with indices
+ and overwrite the values (holding dofs) of the array to_array with indices
interface_dofs, i.e. update the interface dofs of the gli terms.
.
# Parameters
- to_gli_array: #type: np.array, array to overwrite the
+ to_array: #type: np.array, array to overwrite the
dofs of.
interface_dofs: #type: np.array, index array of interface
dofs of from_function
@@ -1112,18 +1140,18 @@ class Interface(BoundaryPart):
d2v = dof_to_vert_map
if not previous_iter:
for dof_index in interface_dofs:
- # to_gli_array is orderd by dof and not by vertex numbering of the mesh.
+ # to_array 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_gli_array[dof_index] = self.gli_term[subdomain_ind][phase][parent[d2v[dof_index]]]
+ to_array[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 array \n", to_gli_array[:])
+ print("wrote these dofs to array \n", to_array[:])
else:
for dof_index in interface_dofs:
- # to_gli_array is orderd by dof and not by vertex numbering of the mesh.
+ # to_array 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_gli_array[dof_index] = self.gli_term_prev[subdomain_ind][phase][parent[d2v[dof_index]]]
+ to_array[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 array \n", to_gli_array[:])
+ print("wrote these dofs to array \n", to_array[:])
#### Private Methods ####