[yt-svn] commit/yt-3.0: MatthewTurk: Removing duplicated Hexahedral Mesh code. Speed up hexahedral_connectivity.
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Nov 26 06:32:18 PST 2013
1 new commit in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/commits/3ecddcf90d39/
Changeset: 3ecddcf90d39
Branch: yt-3.0
User: MatthewTurk
Date: 2013-11-25 22:03:33
Summary: Removing duplicated Hexahedral Mesh code. Speed up hexahedral_connectivity.
Affected #: 2 files
diff -r 6e029bcb0dbf8680278737524cce19e3365cbea1 -r 3ecddcf90d39dcb55ed7bb16cbc28076b326f536 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -861,145 +861,6 @@
return spf
-def hexahedral_connectivity(xgrid, ygrid, zgrid):
- nx = len(xgrid)
- ny = len(ygrid)
- nz = len(zgrid)
- coords = np.fromiter(chain.from_iterable(product(xgrid, ygrid, zgrid)),
- dtype=np.float64, count = nx*ny*nz*3)
- coords.shape = (nx*ny*nz, 3)
- cis = np.fromiter(chain.from_iterable(product([0,1], [0,1], [0,1])),
- dtype=np.int64, count = 8*3)
- cis.shape = (8, 3)
- cycle = np.fromiter(chain.from_iterable(product(*map(range, (nx-1, ny-1, nz-1)))),
- dtype=np.int64, count = (nx-1)*(ny-1)*(nz-1)*3)
- cycle.shape = ((nx-1)*(ny-1)*(nz-1), 3)
- off = cis + cycle[:, np.newaxis]
- connectivity = ((off[:,:,0] * ny) + off[:,:,1]) * nz + off[:,:,2]
- return coords, connectivity
-
-class StreamHexahedralMesh(SemiStructuredMesh):
- _connectivity_length = 8
- _index_offset = 0
-
-class StreamHexahedralHierarchy(UnstructuredGeometryHandler):
-
- def __init__(self, pf, data_style = None):
- self.stream_handler = pf.stream_handler
- super(StreamHexahedralHierarchy, self).__init__(pf, data_style)
-
- def _initialize_mesh(self):
- coords = self.stream_handler.fields.pop('coordinates')
- connec = self.stream_handler.fields.pop('connectivity')
- self.meshes = [StreamHexahedralMesh(0,
- self.hierarchy_filename, connec, coords, self)]
-
- def _setup_data_io(self):
- if self.stream_handler.io is not None:
- self.io = self.stream_handler.io
- else:
- self.io = io_registry[self.data_style](self.pf)
-
- def _detect_fields(self):
- self.field_list = list(set(self.stream_handler.get_fields()))
-
-class StreamHexahedralStaticOutput(StreamStaticOutput):
- _hierarchy_class = StreamHexahedralHierarchy
- _fieldinfo_fallback = StreamFieldInfo
- _fieldinfo_known = KnownStreamFields
- _data_style = "stream_hexahedral"
-
-def load_hexahedral_mesh(data, connectivity, coordinates,
- sim_unit_to_cm, bbox=None,
- sim_time=0.0, periodicity=(True, True, True)):
- r"""Load a hexahedral mesh of data into yt as a
- :class:`~yt.frontends.stream.data_structures.StreamHandler`.
-
- This should allow a semistructured grid of data to be loaded directly into
- yt and analyzed as would any others. This comes with several caveats:
- * Units will be incorrect unless the data has already been converted to
- cgs.
- * Some functions may behave oddly, and parallelism will be
- disappointing or non-existent in most cases.
- * Particles may be difficult to integrate.
-
- Particle fields are detected as one-dimensional fields. The number of particles
- is set by the "number_of_particles" key in data.
-
- Parameters
- ----------
- data : dict
- This is a dict of numpy arrays, where the keys are the field names.
- There must only be one.
- connectivity : array_like
- This should be of size (N,8) where N is the number of zones.
- coordinates : array_like
- This should be of size (M,3) where M is the number of vertices
- indicated in the connectivity matrix.
- sim_unit_to_cm : float
- Conversion factor from simulation units to centimeters
- bbox : array_like (xdim:zdim, LE:RE), optional
- Size of computational domain in units sim_unit_to_cm
- sim_time : float, optional
- The simulation time in seconds
- periodicity : tuple of booleans
- Determines whether the data will be treated as periodic along
- each axis
-
- """
-
- domain_dimensions = np.ones(3, "int32") * 2
- nprocs = 1
- if bbox is None:
- bbox = np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]], 'float64')
- domain_left_edge = np.array(bbox[:, 0], 'float64')
- domain_right_edge = np.array(bbox[:, 1], 'float64')
- grid_levels = np.zeros(nprocs, dtype='int32').reshape((nprocs,1))
-
- sfh = StreamDictFieldHandler()
-
- particle_types = set_particle_types(data)
-
- sfh.update({'connectivity': connectivity,
- 'coordinates': coordinates,
- 0: data})
- grid_left_edges = domain_left_edge
- grid_right_edges = domain_right_edge
- grid_dimensions = domain_dimensions.reshape(nprocs,3).astype("int32")
-
- # I'm not sure we need any of this.
- handler = StreamHandler(
- grid_left_edges,
- grid_right_edges,
- grid_dimensions,
- grid_levels,
- -np.ones(nprocs, dtype='int64'),
- np.zeros(nprocs, dtype='int64').reshape(nprocs,1), # Temporary
- np.zeros(nprocs).reshape((nprocs,1)),
- sfh,
- particle_types=particle_types,
- periodicity=periodicity
- )
-
- handler.name = "HexahedralMeshData"
- handler.domain_left_edge = domain_left_edge
- handler.domain_right_edge = domain_right_edge
- handler.refine_by = 2
- handler.dimensionality = 3
- handler.domain_dimensions = domain_dimensions
- handler.simulation_time = sim_time
- handler.cosmology_simulation = 0
-
- spf = StreamHexahedralStaticOutput(handler)
- spf.units["cm"] = sim_unit_to_cm
- spf.units['1'] = 1.0
- spf.units["unitary"] = 1.0
- box_in_mpc = sim_unit_to_cm / mpc_conversion['cm']
- for unit in mpc_conversion.keys():
- spf.units[unit] = mpc_conversion[unit] * box_in_mpc
-
- return spf
-
class StreamOctreeSubset(OctreeSubset):
domain_id = 1
_domain_offset = 1
@@ -1190,20 +1051,22 @@
return spf
+_cis = np.fromiter(chain.from_iterable(product([0,1], [0,1], [0,1])),
+ dtype=np.int64, count = 8*3)
+_cis.shape = (8, 3)
+
def hexahedral_connectivity(xgrid, ygrid, zgrid):
nx = len(xgrid)
ny = len(ygrid)
nz = len(zgrid)
- coords = np.fromiter(chain.from_iterable(product(xgrid, ygrid, zgrid)),
- dtype=np.float64, count = nx*ny*nz*3)
- coords.shape = (nx*ny*nz, 3)
- cis = np.fromiter(chain.from_iterable(product([0,1], [0,1], [0,1])),
- dtype=np.int64, count = 8*3)
- cis.shape = (8, 3)
- cycle = np.fromiter(chain.from_iterable(product(*map(range, (nx-1, ny-1, nz-1)))),
- dtype=np.int64, count = (nx-1)*(ny-1)*(nz-1)*3)
+ coords = np.zeros((nx, ny, nz, 3), dtype="float64", order="C")
+ coords[:,:,:,0] = xgrid[:,None,None]
+ coords[:,:,:,1] = ygrid[None,:,None]
+ coords[:,:,:,2] = zgrid[None,None,:]
+ coords.shape = (nx * ny * nz, 3)
+ cycle = np.rollaxis(np.indices((nx-1,ny-1,nz-1)), 0, 4)
cycle.shape = ((nx-1)*(ny-1)*(nz-1), 3)
- off = cis + cycle[:, np.newaxis]
+ off = _cis + cycle[:, np.newaxis]
connectivity = ((off[:,:,0] * ny) + off[:,:,1]) * nz + off[:,:,2]
return coords, connectivity
diff -r 6e029bcb0dbf8680278737524cce19e3365cbea1 -r 3ecddcf90d39dcb55ed7bb16cbc28076b326f536 yt/utilities/lib/mesh_utilities.pyx
--- a/yt/utilities/lib/mesh_utilities.pyx
+++ b/yt/utilities/lib/mesh_utilities.pyx
@@ -76,7 +76,7 @@
def smallest_fwidth(np.ndarray[np.float64_t, ndim=2] coords,
np.ndarray[np.int64_t, ndim=2] indices,
int offset = 0):
- cdef np.float64_t fwidth = 1e60
+ cdef np.float64_t fwidth = 1e60, pos
cdef int nc = indices.shape[0]
cdef int nv = indices.shape[1]
if nv != 8:
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list