(REMOVED) support no longer contains property values

this is stored on the interpolator now and the can
be evaluated by providing an array to the support

Author: lachlangrose

LoopStructural/interpolators/discrete_interpolator.py

@@ -31,12 +31,9 @@ def __init__(self, support):
         GeologicalInterpolator.__init__(self)
         self.B = []
         self.support = support
-        self.region_function = None
-        self.region = np.arange(0, support.n_nodes)
-        self.region_map = np.zeros(support.n_nodes).astype(int)
+        self.region_function = lambda xyz : np.ones(xyz.shape[0],dtype=int)
         # self.region_map[self.region] = np.array(range(0,
         # len(self.region_map[self.region])))
-        self.nx = len(self.support.nodes[self.region])
         self.shape = 'rectangular'
         if self.shape == 'square':
             self.B = np.zeros(self.nx)
@@ -53,7 +50,21 @@ def __init__(self, support):
         self.constraints = {}
         self.interpolation_weights= {}
         logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.nx))
-        
+
+    @property
+    def nx(self):
+        return len(self.support.nodes[self.region])
+
+    @property
+    def region(self):
+        return self.region_function(self.support.nodes)
+
+    @property
+    def region_map(self):
+        region_map = np.zeros(self.support.n_nodes).astype(int)
+        region_map[self.region] = np.array(
+            range(0, len(region_map[self.region])))
+        return region_map
     def set_property_name(self, propertyname):
         """
         Set the property name attribute, this is usually used to
@@ -85,11 +96,6 @@ def set_region(self, region=None):
         # evaluate the region function on the support to determine
         # which nodes are inside update region map and degrees of freedom
         self.region_function = region
-        self.region = region(self.support.nodes)
-        self.region_map = np.zeros(self.support.n_nodes).astype(int)
-        self.region_map[self.region] = np.array(
-            range(0, len(self.region_map[self.region])))
-        self.nx = len(self.support.nodes[self.region])
         logger.info("Interpolation now uses region and has {} degrees of freedom".format(self.nx))
 
     def set_interpolation_weights(self, weights):
@@ -511,7 +517,7 @@ def _solve(self, solver='cg', **kwargs):
             logger.warning("Using external solver")
             self.c[self.region] = kwargs['external'](A, B)[:self.nx]
         # check solution is not nan
-        self.support.properties[self.propertyname] = self.c
+        # self.support.properties[self.propertyname] = self.c
         if np.all(self.c == np.nan):
             logger.warning("Solver not run, no scalar field")
         # if solution is all 0, probably didn't work
@@ -544,7 +550,7 @@ def evaluate_value(self, evaluation_points):
 
         if evaluation_points[~mask, :].shape[0] > 0:
             evaluated[~mask] = self.support.evaluate_value(
-                evaluation_points[~mask], self.propertyname)
+                evaluation_points[~mask], self.c)
         return evaluated
 
     def evaluate_gradient(self, evaluation_points):
@@ -561,5 +567,5 @@ def evaluate_gradient(self, evaluation_points):
         """
         if evaluation_points.shape[0] > 0:
             return self.support.evaluate_gradient(evaluation_points,
-                                                  self.propertyname)
+                                                  self.c)
         return np.zeros((0, 3))

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -33,7 +33,6 @@ def __init__(self, mesh):
         DiscreteInterpolator.__init__(self, mesh)
         # whether to assemble a rectangular matrix or a square matrix
         self.interpolator_type = 'PLI'
-        self.nx = len(self.support.nodes[self.region])
         self.support = mesh
 
         self.interpolation_weights = {'cgw': 0.1, 'cpw': 1., 'npw': 1.,

LoopStructural/interpolators/structured_grid.py

@@ -6,18 +6,22 @@
 
 import numpy as np
 
+from .base_structured_3d_support import BaseStructuredSupport
+
+
 from LoopStructural.utils import getLogger
 logger = getLogger(__name__)
 
 
-class StructuredGrid:
+class StructuredGrid(BaseStructuredSupport):
     """
 
     """
     def __init__(self,
                  origin=np.zeros(3),
                  nsteps=np.array([10, 10, 10]),
                  step_vector=np.ones(3),
+                 name=None
                  ):
         """
 
@@ -27,57 +31,15 @@ def __init__(self,
         nsteps - 3d list or numpy array of ints
         step_vector - 3d list or numpy array of int
         """
-
-        self.nsteps = np.array(nsteps)
-        self.step_vector = np.array(step_vector)
-        self.origin = np.array(origin)
-        self.maximum = origin+self.nsteps*self.step_vector
-
-        # self.nsteps+=1
-        self.n_nodes = self.nsteps[0] * self.nsteps[1] * self.nsteps[2]
-        # self.nsteps-=1
-        self.dim = 3
-        self.nsteps_cells = self.nsteps - 1
-        self.n_cell_x = self.nsteps[0] - 1
-        self.n_cell_y = self.nsteps[1] - 1
-        self.n_cell_z = self.nsteps[2] - 1
-        self.properties = {}
-        self.n_elements = self.n_cell_x * self.n_cell_y * self.n_cell_z
-
-        # calculate the node positions using numpy (this should probably not
-        # be stored as it defeats
-        # the purpose of a structured grid
-
-        # self.barycentre = self.cell_centres(np.arange(self.n_elements))
-
+        BaseStructuredSupport.__init__(self,origin,nsteps,step_vector)
         self.regions = {}
         self.regions['everywhere'] = np.ones(self.n_nodes).astype(bool)
-
-
-    @property
-    def nodes(self):
-        max = self.origin + self.nsteps_cells * self.step_vector
-        x = np.linspace(self.origin[0], max[0], self.nsteps[0])
-        y = np.linspace(self.origin[1], max[1], self.nsteps[1])
-        z = np.linspace(self.origin[2], max[2], self.nsteps[2])
-        xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
-        return np.array([xx.flatten(order='F'), yy.flatten(order='F'),
-                               zz.flatten(order='F')]).T
+        self.name = name
 
     def barycentre(self):
         return self.cell_centres(np.arange(self.n_elements))
 
-    # @property
-    # def barycentre(self):
-    #     return self.cell_centres(np.arange(self.n_elements))
 
-    def print_geometry(self):
-        print('Origin: %f %f %f' % (
-            self.origin[0], self.origin[1], self.origin[2]))
-        print('Cell size: %f %f %f' % (
-            self.step_vector[0], self.step_vector[1], self.step_vector[2]))
-        max = self.origin + self.nsteps_cells * self.step_vector
-        print('Max extent: %f %f %f' % (max[0], max[1], max[2]))
 
     def update_property(self, propertyname, values):
         """[summary]
@@ -373,7 +335,7 @@ def position_to_cell_corners(self, pos):
         globalidx[~inside] = -1
         return globalidx, inside
 
-    def evaluate_value(self, evaluation_points, property_name):
+    def evaluate_value(self, evaluation_points, property_array):
         """
         Evaluate the value of of the property at the locations.
         Trilinear interpolation dot corner values
@@ -387,25 +349,33 @@ def evaluate_value(self, evaluation_points, property_name):
         -------
 
         """
+        if property_array.shape[0] != self.n_nodes:
+            logger.error("Property array does not match grid")
+            raise BaseException
         idc, inside = self.position_to_cell_corners(evaluation_points)
         v = np.zeros(idc.shape)
         v[:, :] = np.nan
 
         v[inside, :] = self.position_to_dof_coefs(
             evaluation_points[inside, :]).T
-        v[inside, :] *= self.properties[property_name][idc[inside, :]]
+        
+        v[inside, :] *= property_array[idc[inside, :]]
+        
         return np.sum(v, axis=1)
 
-    def evaluate_gradient(self, evaluation_points, property_name):
+    def evaluate_gradient(self, evaluation_points, property_array):
+        if property_array.shape[0] != self.n_nodes:
+            logger.error("Property array does not match grid")
+            raise BaseException
         idc, inside = self.position_to_cell_corners(evaluation_points)
         T = np.zeros((idc.shape[0], 3, 8))
         T[inside, :, :] = self.calcul_T(evaluation_points[inside, :])
         # indices = np.array([self.position_to_cell_index(evaluation_points)])
         # idc = self.global_indicies(indices.swapaxes(0,1))
         # print(idc)
-        T[inside, 0, :] *= self.properties[property_name][idc[inside, :]]
-        T[inside, 1, :] *= self.properties[property_name][idc[inside, :]]
-        T[inside, 2, :] *= self.properties[property_name][idc[inside, :]]
+        T[inside, 0, :] *= property_array[idc[inside, :]]
+        T[inside, 1, :] *= property_array[idc[inside, :]]
+        T[inside, 2, :] *= property_array[idc[inside, :]]
         return np.array(
             [np.sum(T[:, 0, :], axis=1), np.sum(T[:, 1, :], axis=1) ,
              np.sum(T[:, 2, :], axis=1) ]).T

LoopStructural/interpolators/structured_tetra.py

@@ -85,11 +85,7 @@ def barycentre(self, elements = None):
                                  axis=1) / 4.
         return barycentre
 
-    def update_property(self, name, value):
-
-        self.properties[name] = value
-
-    def evaluate_value(self, pos, prop):
+    def evaluate_value(self, pos, property_array):
         """
         Evaluate value of interpolant
 
@@ -107,10 +103,10 @@ def evaluate_value(self, pos, prop):
         values = np.zeros(pos.shape[0])
         values[:] = np.nan
         vertices, c, tetras, inside = self.get_tetra_for_location(pos)
-        values[inside] = np.sum(c[inside,:]*self.properties[prop][tetras[inside,:]],axis=1)
+        values[inside] = np.sum(c[inside,:]*property_array[tetras[inside,:]],axis=1)
         return values
 
-    def evaluate_gradient(self, pos, prop):
+    def evaluate_gradient(self, pos, property_array):
         """
         Evaluate the gradient of an interpolant at the locations
 
@@ -131,7 +127,7 @@ def evaluate_gradient(self, pos, prop):
         vertices, element_gradients, tetras, inside = self.get_tetra_gradient_for_location(pos)
         vertex_vals = self.properties[prop][tetras]
         #grads = np.zeros(tetras.shape)
-        values[inside,:] = (element_gradients[inside,:,:]*self.properties[prop][tetras[inside,None,:]]).sum(2)
+        values[inside,:] = (element_gradients[inside,:,:]*property_array[tetras[inside,None,:]]).sum(2)
         length = np.sum(values[inside,:],axis=1)
         # values[inside,:] /= length[:,None]
         return values

LoopStructural/modelling/features/geological_feature_builder.py

@@ -14,7 +14,7 @@
     weight_name, gradient_vec_names, tangent_vec_names, interface_name
 from LoopStructural.modelling.features import GeologicalFeature
 from LoopStructural.utils.helper import get_data_bounding_box_map as get_data_bounding_box
-
+from LoopStructural.utils import get_data_axis_aligned_bounding_box
 
 class GeologicalFeatureInterpolator:
     """[summary]
@@ -310,13 +310,29 @@ def build(self, fold=None, fold_weights={}, data_region=None, **kwargs):
 
         """
 
-        if data_region is not None:
-            xyz = self.get_data_locations()
-            bb, region = get_data_bounding_box(xyz, data_region)
-            self.interpolator.set_region(region=region)
+
         if not self.data_added:
             self.add_data_to_interpolator(**kwargs)
-
+        if data_region is not None:
+            xyz = self.interpolator.get_data_locations()
+            bb, region = get_data_bounding_box(xyz, data_region)    
+            # if self.model.reuse_supports == False:
+            if np.any(np.min(bb[0,:])< self.interpolator.support.origin):
+                neworigin = np.min([self.interpolator.support.origin,bb[0,:]],axis=0)
+                logger.info("Changing origin of support for {} from {} {} {} to {} {} {}"\
+                                .format(self.name,self.interpolator.support.origin[0],\
+                                    self.interpolator.support.origin[1],self.interpolator.support.origin[2],\
+                                                neworigin[0],neworigin[1],neworigin[2]))
+                self.interpolator.support.origin = neworigin
+            if np.any(np.max(bb[0,:])< self.interpolator.support.maximum):
+                newmax = np.max([self.interpolator.support.maximum,bb[0,:]],axis=0)
+                logger.info("Changing origin of support for {} from {} {} {} to {} {} {}"\
+                        .format(self.name,self.interpolator.support.maximum[0],
+                        self.interpolator.support.maximum[1],self.interpolator.support.maximum[2],\
+                                                newmax[0],newmax[1],newmax[2]))
+
+                self.interpolator.support.maximum = newmax
+            self.interpolator.set_region(region=region)
         # moving this to init because it needs to be done before constraints
         # are added?
         if fold is not None: