From 3a7fc0d5bcce81b69a2e1f900b4049aff3558948 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 8 Mar 2019 19:50:00 +0100
Subject: [PATCH] add init_boundary_markers, add support for rel perm as
 callable object

---
 domainPatch.py | 102 +++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 81 insertions(+), 21 deletions(-)

diff --git a/domainPatch.py b/domainPatch.py
index 41e8e6c..fe76f08 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -17,7 +17,9 @@ Assumed is the use of python 3.6
 import dolfin as df
 import mshr
 import numpy as np
+# import domainPatch as dp
 import typing as tp
+# Interface = tp.NewType('Interface', tp.any)
 #import domainPatch as domainPatch
 
 
@@ -48,11 +50,14 @@ class DomainPatch(df.SubDomain):
                              present (isRichards == False) viscosity[1] is
                              assumed to be the viscosity of the non-wetting phase.
 
-    has_interface           #type tp.List[int]:  list of indices (w.r.t. the global
-                                                 numbering of interfaces in
-                                                 LDDsimulation.interfaces) indicating
-                                                 which interfaces are part of the
-                                                 subdomain.
+    interfaces:             #type tp.List[dp.Interface]   List of Interface class
+                                                objects from the LDDsimulation
+                                                class.
+    has_interface           #type tp.List[int]: list of indices (w.r.t. the global
+                                                numbering of interfaces in
+                                                LDDsimulation.interfaces) indicating
+                                                which interfaces are part of the
+                                                subdomain.
 
     L                   #type tp.List[float]    L parameters for the L-scheme. In
                                                 case isRichards == False (i.e
@@ -69,25 +74,31 @@ class DomainPatch(df.SubDomain):
 
     ## Public Variables
 
-    self.isRichards          #type bool             set by input, see above
-    self.mesh                #type df.Mesh          set by input, see above
-    self.porosity            #type float            set by input, see above
-    self.viscosity           #type tp.List[float]   set by input, see above
-    self.has_interface       #type tp.List[int]     set by input, see above
-    self.L                   #type tp.List[float]   set by input, see above
-    self.lambda_param        #type tp.List[float]   set by input, see above
-    self.function_space      #type tp.Dict[str: df.Function]    function space
-                                                                used for wetting
-                                                                and nonwetting
-                                                                pressures. This is
-                                                                set by
-                                                                self._init_function_space()
+    isRichards              #type bool              set by input, see above
+    mesh                    #type df.Mesh           set by input, see above
+    porosity                #type float             set by input, see above
+    viscosity               #type tp.List[float]    set by input, see above
+    interfaces:             #type tp.List[dp.Interface] set by input, see above
+    has_interface           #type tp.List[int]      set by input, see above
+    L                       #type tp.List[float]    set by input, see above
+    lambda_param            #type tp.List[float]    set by input, see above
+    function_space          #type tp.Dict[str: df.Function]    function space
+                                                    used for wetting
+                                                    and nonwetting
+                                                    pressures. This is
+                                                    set by
+                                                    self._init_function_space()
+    parent_mesh_index       #type np.array          parent vertex map of the submesh.
+    boundary_marker         #type df.MeshFunction   marker function to mark all
+                                                    boundaries of the subdomain.
+                                                    The marker value marking
+                                                    interface i with global index
+                                                    has_interface[i] will be
+                                                    has_interface[i] aswell.
     ## Public Methods
 
 
     """
-    # default class variable values
-    tol: float = df.DOLFIN_EPS
 
     ### constructor
     def __init__(self, #
@@ -95,9 +106,12 @@ class DomainPatch(df.SubDomain):
                  mesh: df.Mesh,#
                  porosity: float,#
                  viscosity: tp.List[float],#
+                 # tp.List[Interface] doesn't work, because it hasen't been defined yet.
+                 interfaces: tp.List[object],#
                  has_interface: tp.List[int],#
                  L: tp.List[float],#
                  lambda_param: tp.List[float],#
+                 relative_permeability: tp.Callable[...,None]#
                  ):
         # because we declare our own __init__ method for the derived class BoundaryPart,
         # we overwrite the __init__-method of the parent class df.SubDomain. However,
@@ -109,11 +123,15 @@ class DomainPatch(df.SubDomain):
         self.mesh = mesh
         self.porosity = porosity
         self.viscosity = viscosity
+        self.interface = interfaces
         self.has_interface = has_interface
         self.L = L
         self.lambda_param = lambda_param
+        self.relative_permeability = relative_permeability
 
         self._init_function_space()
+        self._init_dof_and_vertex_maps()
+        self._init_boundary_markers()
 
     # END constructor
 
@@ -124,7 +142,8 @@ class DomainPatch(df.SubDomain):
         """ create function space for solution and trial functions
 
         Note that P1 FEM is hard coded here, as it is assumed in other methods
-        aswell.
+        aswell. This method sets the class variable
+            self.function_space
         """
         if self.isRichards:
             self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
@@ -133,6 +152,47 @@ class DomainPatch(df.SubDomain):
                 'wetting'    : df.FunctionSpace(self.mesh, 'P', 1),#
                 'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
                 }
+
+    def _init_dof_and_vertex_maps(self) -> None:
+        """ calculate dof maps and the vertex to parent map.
+
+        This method sets the class Variables
+            self.parent_mesh_index
+            self.dof2vertex
+            self.vertex2dof
+        """
+        mesh_data = self.mesh.data()
+        self.parent_mesh_index = mesh_data.array('parent_vertex_indices',0)
+        if self.isRichards:
+            self.dof2vertex = {#
+                'wetting' :df.dof_to_vertex_map(self.function_space['wetting'])#
+                }
+            self.vertex2dof = {#
+                'wetting' :df.vertex_to_dof_map(self.function_space['wetting'])#
+                }
+        else:
+            self.dof2vertex = {#
+                'wetting'   :df.dof_to_vertex_map(self.function_space['wetting']),#
+                'nonwetting':df.dof_to_vertex_map(self.function_space['nonwetting'])#
+                }
+            self.vertex2dof = {#
+                'wetting'   :df.vertex_to_dof_map(self.function_space['wetting']),#
+                'nonwetting':df.vertex_to_dof_map(self.function_space['nonwetting'])#
+                }
+
+    def _init_boundary_markers(self) -> None:
+        """ define boundary markers
+
+        This method sets the class Variables
+            self.boundary_marker
+
+        """
+        self.boundary_marker = df.MeshFunction('size_t', self.mesh, self.mesh.topology().dim()-1)
+        self.boundary_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)
+
     # END is_Richards
 
 # END OF CLASS DomainPatch
-- 
GitLab