Repository: PMEAL/OpenPNM
Pull request #832 (2018-03-22 01:48:40)

Breakup base
This is a backend change that I've been meaning to make for a while.  The GenericPhysics and GenericGeometry objects need an extra level of inheritance to reduce redundant code, pertaining to setting locations, and the getitem method.


Patch:
diff --git a/example_script.py b/example_script.py
index d93cc1e8a..ca6933cf5 100644
--- a/example_script.py
+++ b/example_script.py
@@ -7,16 +7,32 @@ pn = op.network.Cubic(shape=[5, 5, 5], spacing=0.0001, name='pn11')
 pn.add_boundary_pores()
 proj = pn.project
 
-geom = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
+Ps = pn.pores('internal')
+Ts = pn.throats('internal')
+geom = op.geometry.StickAndBall(network=pn, pores=Ps, throats=Ts,
+                                name='intern')
+
+Ps = pn.pores('*boundary')
+Ts = pn.throats('*boundary')
+boun = op.geometry.StickAndBall(network=pn, pores=Ps, throats=Ts, name='bound')
+
+pn['pore.inlets'] = pn['pore.top_boundary'].copy()
+pn['pore.outlets'] = pn['pore.bottom_boundary'].copy()
 
 air = op.phases.Air(network=pn)
 water = op.phases.Water(network=pn)
 water['throat.viscosity'] = water['pore.viscosity'][0]
 
 mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
-phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
-phys_water.add_model(propname='throat.conductance',
-                     model=mod, viscosity='throat.viscosity')
+phys_water_1 = op.physics.GenericPhysics(network=pn, phase=water,
+                                         geometry=geom)
+phys_water_1.add_model(propname='throat.conductance',
+                       model=mod,
+                       viscosity='throat.viscosity')
+
+phys_water_2 = op.physics.GenericPhysics(network=pn, phase=water,
+                                         geometry=boun)
+phys_water_2.models = phys_water_1.models.copy()
 
 water['pore.A'] = 1e-10
 water['pore.k'] = 2
@@ -25,19 +41,16 @@ water.add_model(propname='pore.reaction',
                 prefactor='pore.A', exponent='pore.k',
                 quantity='pore.pressure', regen_mode='deferred')
 
-water.add_model(propname='throat.diffusive_conductance',
-                model=op.models.physics.diffusive_conductance.bulk_diffusion)
-
 s = {'conductance': 'throat.conductance',
      'quantity': 'pore.pressure'}
 alg1 = op.algorithms.GenericTransport(network=pn, phase=water, settings=s)
-alg1.set_dirchlet_BC(pores=pn.pores('top'), values=1)
-alg1.set_dirchlet_BC(pores=pn.pores('bottom'), values=0)
+alg1.set_dirchlet_BC(pores=pn.pores('inlets'), values=1)
+alg1.set_dirchlet_BC(pores=pn.pores('outlets'), values=0)
 alg1.run()
 
 alg2 = op.algorithms.ReactiveTransport(network=pn, phase=water, settings=s)
-alg2.set_dirchlet_BC(pores=pn.pores('top'), values=1)
-alg2.set_source_term(propname='pore.reaction', pores=pn.pores('bottom'))
+alg2.set_dirchlet_BC(pores=pn.pores('inlets'), values=1)
+alg2.set_source_term(propname='pore.reaction', pores=pn.pores('outlets'))
 alg2.run()
 water.update(alg2.results())
 
@@ -49,14 +62,14 @@ alg3.settings.update({'t_initial': 0,
                       't_step': 0.25,
                       'conductance': 'throat.conductance',
                       'quantity': 'pore.pressure'})
-alg3.set_dirchlet_BC(pores=pn.pores('top'), values=1)
+alg3.set_dirchlet_BC(pores=pn.pores('inlets'), values=1)
 alg3.set_IC(values=0)
 alg3.run()
 
 
 alg4 = op.algorithms.TransientReactiveTransport(network=pn, phase=water)
 alg4.settings.update(alg3.settings)  # Just copy settings from another alg
-alg4.set_dirchlet_BC(pores=pn.pores('top'), values=1)
+alg4.set_dirchlet_BC(pores=pn.pores('inlets'), values=1)
 alg4.set_IC(values=0)
 alg4.set_source_term(propname='pore.reaction', pores=pn.pores('bottom'))
 alg4.run()
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index fc0a6c666..74df8ca51 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -128,7 +128,7 @@ class GenericTransport(GenericAlgorithm):
         network = self.project.network
         phase = self.project.phases()[self.settings['phase']]
         g = phase[self.settings['conductance']]
-        am = network.create_adjacency_matrix(weights=-g, fmt='coo')
+        am = network.create_adjacency_matrix(weights=g, fmt='coo')
         A = spgr.laplacian(am)
         self.A = A
         return A
diff --git a/openpnm/core/Base.py b/openpnm/core/Base.py
index 1de3c5170..14503a816 100644
--- a/openpnm/core/Base.py
+++ b/openpnm/core/Base.py
@@ -12,8 +12,11 @@ class Base(dict):
     """
 
     def __new__(cls, *args, **kwargs):
-        cls.settings = SettingsDict()
         instance = super(Base, cls).__new__(cls, *args, **kwargs)
+        # The SettingsDict implements the __missing__ magic method, which
+        # returns None instead of KeyError.  This is useful for checking the
+        # value of a settings without first ensuring it exists.
+        instance.settings = SettingsDict()
         return instance
 
     def __init__(self, Np=0, Nt=0, name=None, project=None):
@@ -732,8 +735,12 @@ class Base(dict):
         sizes = [sp.size(a) for a in arrs]
         if all([item is None for item in arrs]):  # prop not found anywhere
             raise KeyError(prop)
-        if sp.any([i is None for i in arrs]):  # prop not found everywhere
-            logger.warning('\''+prop+'\' not found on at least one object')
+#        if sp.any([i is None for i in arrs]):  # prop not found everywhere
+#            logger.warning('\''+prop+'\' not found on at least one object')
+#        if sp.sum(sizes) < self._count(element):
+#            logger.warning('Not all '+element+'s are assigned to an object')
+#            N_missing = self._count(element) - sp.sum(sizes)
+#            arrs.append(sp.zeros(shape=(N_missing,), dtype=float)*sp.nan)
 
         # Check the general type of each array
         atype = []
@@ -760,10 +767,9 @@ class Base(dict):
                     temp_arr = sp.zeros((N, item.shape[1]), dtype=item.dtype)
                 temp_arr.fill(dummy_val[atype[0]])
 
-        # Convrert int arrays to float IF NaNs are expected
-        if (temp_arr.dtype.name.startswith('int') and
-            (sp.any([i is None for i in arrs]) or
-             sp.sum(sizes) != N)):
+        # Convert int arrays to float IF NaNs are expected
+        if temp_arr.dtype.name.startswith('int') and \
+           (sp.any([i is None for i in arrs]) or sp.sum(sizes) != N):
             temp_arr = temp_arr.astype(float)
             temp_arr.fill(sp.nan)
 
diff --git a/openpnm/core/Project.py b/openpnm/core/Project.py
index 76f141478..eb635886a 100644
--- a/openpnm/core/Project.py
+++ b/openpnm/core/Project.py
@@ -58,47 +58,76 @@ class Project(list):
         return obj
 
     def find_phase(self, obj):
-        mro = obj._mro()
-        if 'GenericPhase'in mro:
+        # If received phase, just return self
+        if obj._isa('phase'):
             return obj
-        if 'GenericAlgorithm' in mro:
+        # If phase happens to be in settings (i.e. algorithm), look it up
+        if 'phase' in obj.settings.keys():
             phase = self.phases()[obj.settings['phase']]
             return phase
-        for g in self.geometries().values():
-            for p in self.phases().values():
-                if obj.name == self.grid[g.name][p.name]:
-                    return p
+        # Otherwise find it using bottom-up approach (i.e. look in phase keys)
+        for phase in self.phases().values():
+            if ('pore.'+obj.name in phase) or ('throat.'+obj.name in phase):
+                return phase
+        # If all else fails, throw an exception
+        raise Exception('Cannot find a phase associated with '+obj.name)
 
     def find_geometry(self, physics):
-        for g in self.geometries().values():
-            for p in self.phases().values():
-                if physics.name == self.grid[g.name][p.name]:
-                    return g
+        # If geometry happens to be in settings, look it up directly
+        if 'geometry' in physics.settings.keys():
+            geom = self.geometries()[physics.settings['geometry']]
+            return geom
+        # Otherwise, use the bottom-up approach
+        for geo in self.geometries().values():
+            phys = self.find_physics(geometry=geo)
+            if phys is physics:
+                return geo
+        # If all else fails, throw an exception
+        raise Exception('Cannot find a geometry associated with '+physics.name)
 
     def find_physics(self, geometry=None, phase=None):
         if geometry and phase:
-            name = self.grid[geometry.name][phase.name]
-            phys = self[name]
+            physics = self.find_physics(geometry=geometry)
+            phases = list(self.phases().values())
+            phys = physics[phases.index(phase)]
+            return phys
         elif geometry:
-            row = self.grid.row(geometry.name)
-            phys = [self.physics().get(i, None) for i in row]
+            result = []
+            net = self.network
+            geoPs = net['pore.'+geometry.name]
+            geoTs = net['throat.'+geometry.name]
+            for phase in self.phases().values():
+                physics = self.find_physics(phase=phase)
+                temp = None
+                for phys in physics:
+                    Ps = phase.map_pores(phys['pore._id'])
+                    physPs = phase.tomask(pores=Ps)
+                    Ts = phase.map_throats(phys['throat._id'])
+                    physTs = phase.tomask(throats=Ts)
+                    if np.all(geoPs == physPs) and np.all(geoTs == physTs):
+                        temp = phys
+                result.append(temp)
+            return result
         elif phase:
-            col = self.grid.col(phase.name)
-            phys = [self.physics().get(i, None) for i in col]
+            names = set(self.physics().keys())
+            keys = set([item.split('.')[-1] for item in phase.keys()])
+            hits = names.intersection(keys)
+            phys = [self.physics().get(i, None) for i in hits]
+            return phys
         else:
-            raise Exception('Must specify at least one of geometry or phase')
-        if phys == ['']:
             phys = []
-        return phys
+            for geom in self.geometries().values():
+                phys.append(self.find_physics(geometry=geom))
+            return phys
 
     def _validate_name(self, name):
         names = [i.name for i in self]
         if name in names:
-            raise Exception('An object already exists named ' + name)
+            raise Exception('An object already exists named '+name)
         for item in self:
             for key in item.keys():
                 if key.split('.')[1] == name:
-                    raise Exception('A property/label is already named ' + name)
+                    raise Exception('A property/label is already named '+name)
 
     def _generate_name(self, obj):
         prefix = obj.settings['prefix']
@@ -191,40 +220,29 @@ class Project(list):
             self[obj_name][propname] = f[item]
         f.close()
 
-    def _get_net(self):
-        for item in self:
-            if 'GenericNetwork' in item._mro():
-                return item
-
-    network = property(fget=_get_net)
+    @property
+    def network(self):
+        net = list(self._get_objects_of_type('network').values())
+        if len(net) > 0:
+            net = net[0]
+        else:
+            net = None
+        return net
 
     def geometries(self):
-        _dict = {}
-        for item in self:
-            if 'GenericGeometry' in item._mro():
-                _dict.update({item.name: item})
-        return _dict
+        return self._get_objects_of_type('geometry')
 
     def phases(self):
-        _dict = {}
-        for item in self:
-            if 'GenericPhase' in item._mro():
-                _dict.update({item.name: item})
-        return _dict
+        return self._get_objects_of_type('phase')
 
     def physics(self):
-        _dict = {}
-        for item in self:
-            if 'GenericPhysics' in item._mro():
-                _dict.update({item.name: item})
-        return _dict
+        return self._get_objects_of_type('physics')
 
     def algorithms(self):
-        _dict = {}
-        for item in self:
-            if 'GenericAlgorithm' in item._mro():
-                _dict.update({item.name: item})
-        return _dict
+        return self._get_objects_of_type('algorithm')
+
+    def _get_objects_of_type(self, objtype):
+        return {item.name: item for item in self if item._isa(objtype)}
 
     def _set_comments(self, string):
         if hasattr(self, '_comments') is False:
@@ -240,20 +258,18 @@ class Project(list):
     comments = property(fget=_get_comments, fset=_set_comments)
 
     def _get_grid(self):
-        net = self.network
         grid = {}
         row = {phase: '' for phase in self.phases().keys()}
-        for geo in self.geometries().keys():
-            grid[geo] = row.copy()
+        for geo in self.geometries().values():
+            grid[geo.name] = row.copy()
             for phase in self.phases().values():
-                for phys in self.physics().keys():
-                    if phys in [n.split('.')[1] for n in phase.keys()]:
-                        geo_mask = net['pore.'+geo]
-                        phys_mask = phase['pore.'+phys]
-                        # TODO: This could be more or less strict
-                        if np.sum(geo_mask*phys_mask) > 0:
-                            grid[geo][phase.name] = phys
-        grid = ProjectGrid(self.name, grid)
+                phys = self.find_physics(phase=phase, geometry=geo)
+                if phys is None:
+                    phys = ''
+                else:
+                    phys = phys.name
+                grid[geo.name][phase.name] = phys
+        grid = ProjectGrid(self.network.name, grid)
         return grid
 
     grid = property(fget=_get_grid)
@@ -313,6 +329,10 @@ class Grid(dict):
 
 
 class ProjectGrid(Grid):
+    r"""
+    This is a subclass of Grid, which adds the ability to lookup by geometries
+    and phases, as more specific versions of rows and cols
+    """
 
     def geometries(self):
         return self.index()
diff --git a/openpnm/core/Subdomain.py b/openpnm/core/Subdomain.py
new file mode 100644
index 000000000..ac5875f45
--- /dev/null
+++ b/openpnm/core/Subdomain.py
@@ -0,0 +1,106 @@
+from openpnm.core import Base
+import numpy as np
+
+
+class Subdomain(Base):
+
+    def __getitem__(self, key):
+        element = key.split('.')[0]
+        # Find boss object (either phase or network)
+        if self._isa('physics'):
+            boss = self.project.find_phase(self)
+        else:
+            boss = self.project.network
+        # Deal with a few special key items
+        if key.split('.')[-1] == '_id':
+            inds = boss._get_indices(element=element, labels=self.name)
+            return boss[element+'._id'][inds]
+        # Convert self.name into 'all'
+        elif key.split('.')[-1] in [self.name]:
+            return self[element+'.all']
+        # Now get values if present, or regenerate them
+        vals = self.get(key)
+        # If still not found, check with boss object (interleave data)
+        if vals is None:
+            inds = boss._get_indices(element=element, labels=self.name)
+            vals = boss[key][inds]
+        return vals
+
+    def __setitem__(self, key, value):
+        if self.project:
+            # Find boss object (either phase or network)
+            if self._isa('phase'):
+                boss = self.project.find_phase(self)
+            else:
+                boss = self.project.network
+            if key in set(boss.keys()).difference({'pore.all', 'throat.all'}):
+                raise Exception(key + ' already exists on ' + boss.name)
+        super().__setitem__(key, value)
+
+    def add_locations(self, pores=[], throats=[]):
+        r"""
+        """
+        boss = self._get_boss()
+        pores = boss._parse_indices(pores)
+        throats = boss._parse_indices(throats)
+        if len(pores) > 0:
+            self._set_locations(element='pore', indices=pores, mode='add')
+        if len(throats) > 0:
+            self._set_locations(element='throat', indices=throats, mode='add')
+
+    def drop_locations(self, pores=[], throats=[]):
+        r"""
+        """
+        boss = self._get_boss()
+        pores = boss._parse_indices(pores)
+        throats = boss._parse_indices(throats)
+        if len(pores) > 0:
+            self._set_locations(element='pore', indices=pores, mode='drop')
+        if len(throats) > 0:
+            self._set_locations(element='throat', indices=throats, mode='drop')
+
+    def _set_locations(self, element, indices=[], mode='add'):
+        r"""
+        """
+        boss = self._get_boss()
+        element = self._parse_element(element=element, single=True)
+
+        # Make sure label array exists in boss
+        if (element+'.'+self.name) not in boss.keys():
+            boss[element+'.'+self.name] = False
+
+        # Check to ensure indices aren't already assigned
+        if mode == 'add':
+            if self._isa('geometry'):
+                objs = self.project.geometries().keys()
+            else:
+                objs = self.project.physics().keys()
+            for name in objs:
+                if element+'.'+name in boss.keys():
+                    if np.any(boss[element+'.'+name][indices]):
+                        raise Exception('Given indices are already assigned ' +
+                                        'to ' + name)
+
+        # Find mask of existing locations (network indexing)
+        mask = boss[element+'.'+self.name]
+        # Update mask with new locations (either add or remove)
+        if mode == 'add':
+            mask = mask + boss._tomask(indices=indices, element=element)
+        elif mode == 'drop':
+            mask = mask * (~boss._tomask(indices=indices, element=element))
+        # Change size of all arrays on self
+        for item in self.keys(element=element, mode='all'):
+            self.update({item: boss[item][mask]})
+        # Update label array in network
+        boss[element+'.'+self.name] = mask
+        # Remove label from boss if ALL locations are removed
+        if mode == 'drop':
+            if ~np.any(boss[element+'.'+self.name]):
+                del boss[element+'.'+self.name]
+
+    def _get_boss(self):
+        if self._isa('physics'):
+            boss = self.project.find_phase(self)
+        if self._isa('geometry'):
+            boss = self.project.network
+        return boss
diff --git a/openpnm/core/__init__.py b/openpnm/core/__init__.py
index f1813a599..90161df7c 100644
--- a/openpnm/core/__init__.py
+++ b/openpnm/core/__init__.py
@@ -3,3 +3,4 @@ from .Workspace import Workspace
 from .Project import Project
 from .ModelsMixin import ModelsMixin, ModelsDict
 from .Base import Base
+from .Subdomain import Subdomain
diff --git a/openpnm/geometry/GenericGeometry.py b/openpnm/geometry/GenericGeometry.py
index a1afb7605..86206a884 100644
--- a/openpnm/geometry/GenericGeometry.py
+++ b/openpnm/geometry/GenericGeometry.py
@@ -1,13 +1,11 @@
 import scipy as sp
-import numpy as np
 import matplotlib.pyplot as plt
-from openpnm.core import Base, ModelsMixin
-from openpnm.core import logging, Workspace
+from openpnm.core import Subdomain, ModelsMixin, Workspace, logging
 logger = logging.getLogger(__name__)
 ws = Workspace()
 
 
-class GenericGeometry(Base, ModelsMixin):
+class GenericGeometry(Subdomain, ModelsMixin):
     r"""
     GenericGeometry - Base class to construct a Geometry object
 
@@ -51,89 +49,6 @@ class GenericGeometry(Base, ModelsMixin):
 
         self.add_locations(pores=pores, throats=throats)
 
-    def __getitem__(self, key):
-        # Find boss object (either phase or network)
-        element = key.split('.')[0]
-        if self._isa('phase'):
-            boss = self.project.find_phase(self)
-        else:
-            boss = self.project.network
-        # Deal with a few special key items
-        if key.split('.')[-1] == '_id':
-            inds = boss._get_indices(element=element, labels=self.name)
-            return boss[element+'._id'][inds]
-        # Convert self.name into 'all'
-        elif key.split('.')[-1] in [self.name]:
-            return self[element+'.all']
-        # Now get values if present, or regenerate them
-        vals = self.get(key)
-        if vals is None:
-            inds = boss._get_indices(element=element, labels=self.name)
-            vals = boss[key][inds]
-        return vals
-
-    def add_locations(self, pores=[], throats=[]):
-        r"""
-        """
-        if len(pores) > 0:
-            self._set_locations(element='pore', indices=pores, mode='add')
-        if len(throats) > 0:
-            self._set_locations(element='throat', indices=throats, mode='add')
-
-    def drop_locations(self, pores=[], throats=[]):
-        r"""
-        """
-        if len(pores) > 0:
-            self._set_locations(element='pore', indices=pores, mode='drop')
-        if len(throats) > 0:
-            self._set_locations(element='throat', indices=throats, mode='drop')
-
-    def _set_locations(self, element, indices=[], mode='add'):
-        r"""
-        """
-        element = self._parse_element(element=element, single=True)
-        # Use the network's _parse_indices, since indicies could be 'network'
-        # length boolean masks
-        network = self.project.network
-        indices = network._parse_indices(indices)
-
-        net = self.project.network
-        proj = self.project
-
-        if element+'.'+self.name not in net.keys():
-            net[element+'.'+self.name] = False
-
-        if mode == 'add':
-            # Ensure indices are not already assigned to another object
-            temp = sp.zeros(shape=[net._count(element=element), ], dtype=bool)
-            for item in proj.geometries().keys():
-                temp += net[element+'.'+item]
-            if sp.any(temp[indices]):
-                raise Exception('Some of the given '+element+' are already ' +
-                                'assigned to an existing object')
-            set_flag = True
-        elif mode == 'drop':
-            set_flag = False
-
-        # Change lables of all associated physics in their respective phases
-        for phase in proj.phases().values():
-            phys = proj.find_physics(geometry=self, phase=phase)
-            if phys:
-                if element+'.'+phys.name not in phase.keys():
-                    phase[element+'.'+phys.name] = False
-                phase[element+'.'+phys.name][indices] = set_flag
-                temp = sp.ones(shape=(sp.sum(phase[element+'.'+phys.name]),),
-                               dtype=bool)
-                phys.update({element+'.all': temp})
-
-        # Change labels in the network
-        net[element+'.'+self.name][indices] = set_flag
-        temp = sp.ones(shape=(sp.sum(net[element+'.'+self.name]),), dtype=bool)
-        self.update({element+'.all': temp})
-        inds = net._get_indices(element=element, labels=self.name)
-        for item in self.props(element=element):
-            self.update({item: net[item][inds]})
-
     def show_hist(self, props=['pore.diameter'], bins=20, fig=None, **kwargs):
         if fig is None:
             fig = plt.figure()
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index d8eb7a785..05b6e65a4 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -108,7 +108,8 @@ class Cubic(GenericNetwork):
 
         super().__init__(Np=points.shape[0], Nt=pairs.shape[0], name=name,
                          project=project)
-
+        self['pore.internal'] = True
+        self['throat.internal'] = True
         self._label_surfaces()
 
     def _label_surfaces(self):
@@ -121,11 +122,6 @@ class Cubic(GenericNetwork):
         for label in labels:
             if 'pore.'+label not in self.keys():
                 self['pore.'+label] = False
-        if 'pore.boundary' in self.keys():
-            internal = -self['pore.boundary']
-        else:
-            internal = self['pore.all']
-        self['pore.internal'] = internal
         self['pore.front'][x <= x.min()] = True
         self['pore.back'][x >= x.max()] = True
         self['pore.left'][y <= y.min()] = True
diff --git a/openpnm/network/GenericNetwork.py b/openpnm/network/GenericNetwork.py
index b3366a03d..d8637b957 100644
--- a/openpnm/network/GenericNetwork.py
+++ b/openpnm/network/GenericNetwork.py
@@ -22,8 +22,8 @@ class GenericNetwork(Base, ModelsMixin):
         self._im = {}
         self._am = {}
 
-    def __setitem__(self, prop, value):
-        if prop == 'throat.conns':
+    def __setitem__(self, key, value):
+        if key == 'throat.conns':
             if sp.shape(value)[1] != 2:
                 logger.error('Wrong size for throat conns!')
             else:
@@ -31,7 +31,12 @@ class GenericNetwork(Base, ModelsMixin):
                     logger.warning('Converting throat.conns to be upper \
                                     triangular')
                     value = sp.sort(value, axis=1)
-        super().__setitem__(prop, value)
+        if self.project:
+            for item in self.project.geometries().values():
+                exclude = {'pore.all', 'throat.all'}
+                if key in set(item.keys()).difference(exclude):
+                    raise Exception(key+' already exists on '+item.name)
+        super().__setitem__(key, value)
 
     def __getitem__(self, key):
         # Deal with special keys first
diff --git a/openpnm/phases/GenericPhase.py b/openpnm/phases/GenericPhase.py
index 6b200b1e8..dc351f85c 100644
--- a/openpnm/phases/GenericPhase.py
+++ b/openpnm/phases/GenericPhase.py
@@ -44,6 +44,14 @@ class GenericPhase(Base, ModelsMixin):
         self['pore.temperature'] = 298.0
         self['pore.pressure'] = 101325.0
 
+    def __setitem__(self, key, value):
+        if self.project:
+            for item in self.project.find_physics(phase=self):
+                exclude = {'pore.all', 'throat.all'}
+                if key in set(item.keys()).difference(exclude):
+                    raise Exception(key+' already exists on '+item.name)
+        super().__setitem__(key, value)
+
     def __getitem__(self, key):
         element = key.split('.')[0]
         # Deal with special keys first
diff --git a/openpnm/physics/GenericPhysics.py b/openpnm/physics/GenericPhysics.py
index dba9492d0..c9e48c3e5 100644
--- a/openpnm/physics/GenericPhysics.py
+++ b/openpnm/physics/GenericPhysics.py
@@ -1,10 +1,10 @@
-from openpnm.core import Base, Workspace, ModelsMixin, logging
+from openpnm.core import Subdomain, ModelsMixin, Workspace, logging
 from numpy import ones
 logger = logging.getLogger(__name__)
 ws = Workspace()
 
 
-class GenericPhysics(Base, ModelsMixin):
+class GenericPhysics(Subdomain, ModelsMixin):
     r"""
     Generic class to generate Physics objects
 
@@ -43,41 +43,10 @@ class GenericPhysics(Base, ModelsMixin):
 
         if phase is not None:
             self.settings['phase'] = phase.name
+            phase['pore.'+self.name] = False
+            phase['throat.'+self.name] = False
         if geometry is not None:
             self.settings['geometry'] = geometry.name
-            self.set_locations(geometry)
-
-    def set_locations(self, geometry):
-        phase = self.project.phases()[self.settings['phase']]
-        network = self.project.network
-
-        # Adjust array sizes
-        self['pore.all'] = ones((geometry.Np, ), dtype=bool)
-        self['throat.all'] = ones((geometry.Nt, ), dtype=bool)
-
-        # Initialize a label array in the associated phase
-        phase['pore.'+self.name] = False
-        phase['pore.'+self.name][network.pores(geometry.name)] = True
-        phase['throat.'+self.name] = False
-        phase['throat.'+self.name][network.throats(geometry.name)] = True
-
-    def __getitem__(self, key):
-        # Find boss object (either phase or network)
-        element = key.split('.')[0]
-        if self._isa('phase'):
-            boss = self.project.find_phase(self)
-        else:
-            boss = self.project.network
-        # Deal with a few special key items
-        if key.split('.')[-1] == '_id':
-            inds = boss._get_indices(element=element, labels=self.name)
-            return boss[element+'._id'][inds]
-        # Convert self.name into 'all'
-        elif key.split('.')[-1] in [self.name]:
-            return self[element+'.all']
-        # Now get values if present, or regenerate them
-        vals = self.get(key)
-        if vals is None:
-            inds = boss._get_indices(element=element, labels=self.name)
-            vals = boss[key][inds]
-        return vals
+            Ps = network.pores(geometry.name)
+            Ts = network.throats(geometry.name)
+            self.add_locations(pores=Ps, throats=Ts)


Test patch:
diff --git a/tests/unit/core/BaseTest.py b/tests/unit/core/BaseTest.py
index 511c2f953..61fb2797c 100644
--- a/tests/unit/core/BaseTest.py
+++ b/tests/unit/core/BaseTest.py
@@ -327,7 +327,7 @@ class BaseTest:
 
     def test_labels_on_throats(self):
         a = self.net.labels(element='throats')
-        b = ['throat.all', 'throat.'+self.geo.name]
+        b = ['throat.all', 'throat.internal', 'throat.'+self.geo.name]
         assert sorted(a) == sorted(b)
 
     def test_labels_on_foo(self):
@@ -343,7 +343,7 @@ class BaseTest:
 
     def test_labels_on_all_throats(self):
         a = self.net.labels(throats=self.net.Ts)
-        b = ['throat.all', 'throat.'+self.geo.name]
+        b = ['throat.all', 'throat.internal', 'throat.'+self.geo.name]
         assert sorted(a) == sorted(b)
 
     def test_labels_on_one_pore(self):
diff --git a/tests/unit/core/SubdomainTest.py b/tests/unit/core/SubdomainTest.py
new file mode 100644
index 000000000..0ef71137c
--- /dev/null
+++ b/tests/unit/core/SubdomainTest.py
@@ -0,0 +1,84 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class SubdomainTest:
+
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.settings['local_data'] = True
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['pore.diameter'] = sp.rand(self.net.Np)
+        self.geo.add_model(propname='pore.volume',
+                           model=op.models.geometry.pore_volume.sphere)
+        self.geo['throat.diameter'] = sp.rand(self.net.Nt)
+        self.geo.add_model(propname='throat.area',
+                           model=op.models.geometry.throat_area.cylinder)
+        self.geo.regenerate_models()
+        self.phase1 = op.phases.GenericPhase(network=self.net)
+        self.phase2 = op.phases.GenericPhase(network=self.net)
+        self.phys1 = op.physics.GenericPhysics(network=self.net,
+                                               geometry=self.geo,
+                                               phase=self.phase1)
+        self.phys1['pore.blah'] = 1.0
+        self.phys2 = op.physics.GenericPhysics(network=self.net,
+                                               geometry=self.geo,
+                                               phase=self.phase2)
+        self.phys2['pore.blah'] = 2.0
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+    def test_drop_locations_from_geom_successively_with_single_geometry(self):
+        assert self.geo.Np == 27
+        assert self.geo.Nt == 54
+        self.geo.drop_locations(pores=[0, 1, 2], throats=[0, 1, 2])
+        assert self.geo.Np == 24
+        assert self.geo.Nt == 51
+        self.geo.drop_locations(pores=[3, 4], throats=[3, 4])
+        assert self.geo.Np == 22
+        assert self.geo.Nt == 49
+        self.geo.add_locations(pores=[0, 1, 2, 3, 4], throats=[0, 1, 2, 3, 4])
+        assert self.geo.Np == 27
+        assert self.geo.Nt == 54
+
+    def test_drop_locations_from_physics_successively_with_two_physics(self):
+        assert self.phys1.Np == 27
+        assert self.phys1.Nt == 54
+        self.phys1.drop_locations(pores=[0, 1], throats=[0, 1])
+        assert self.phys1.Np == 25
+        assert self.phys1.Nt == 52
+        self.phys1.drop_locations(pores=[3, 4], throats=[3, 4])
+        assert self.phys1.Np == 23
+        assert self.phys1.Nt == 50
+        self.phys1.add_locations(pores=[0, 1, 3, 4], throats=[0, 1, 3, 4])
+        assert self.phys1.Np == 27
+        assert self.phys1.Nt == 54
+
+    def test_drop_locations_completely(self):
+        assert self.phys1.Np == 27
+        assert self.phys1.Nt == 54
+        assert 'pore.'+self.phys1.name in self.phase1.keys()
+        assert 'throat.'+self.phys1.name in self.phase1.keys()
+        self.phys1.drop_locations(pores=self.net.Ps)
+        assert 'pore.'+self.phys1.name not in self.phase1.keys()
+        assert 'throat.'+self.phys1.name in self.phase1.keys()
+        self.phys1.drop_locations(throats=self.net.Ts)
+        assert 'throat.'+self.phys1.name not in self.phase1.keys()
+        self.phys1.add_locations(pores=self.net.Ps, throats=self.net.Ts)
+
+
+if __name__ == '__main__':
+
+    t = SubdomainTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/io/PandasTest.py b/tests/unit/io/PandasTest.py
index 883d391fa..3081ed760 100644
--- a/tests/unit/io/PandasTest.py
+++ b/tests/unit/io/PandasTest.py
@@ -63,12 +63,12 @@ class PandasTest:
         df = Pandas.to_dataframe(network=self.net, phases=[self.phase_1],
                                  join=False)
         assert len(df.pore.keys()) == 21
-        assert len(df.throat.keys()) == 11
+        assert len(df.throat.keys()) == 12
 
     def test_to_dataframe_joined(self):
         df = Pandas.to_dataframe(network=self.net, phases=[self.phase_1],
                                  join=True)
-        assert len(df.keys()) == 32
+        assert len(df.keys()) == 33
 
     def test_save(self):
         with pytest.raises(NotImplementedError):