diff --git a/HySoP/hysop/domain/subsets/control_box.py b/HySoP/hysop/domain/control_box.py
similarity index 98%
rename from HySoP/hysop/domain/subsets/control_box.py
rename to HySoP/hysop/domain/control_box.py
index 53e161aeb3bc60dca47729acac67dc6ad4148473..2c3bd2214030c074db8de62de079d1017820892b 100644
--- a/HySoP/hysop/domain/subsets/control_box.py
+++ b/HySoP/hysop/domain/control_box.py
@@ -1,5 +1,5 @@
"""Define a volume of control inside a domain (volume + all faces)"""
-from hysop.domain.subsets.boxes import SubBox
+from hysop.domain.subsets import SubBox
import numpy as np
import hysop.tools.numpywrappers as npw
@@ -11,7 +11,6 @@ class ControlBox(SubBox):
and set of indices belonging to surfaces of this domain (slices members).
Useful to define control volume to perform integration.
"""
-
def __init__(self, **kwds):
super(ControlBox, self).__init__(**kwds)
diff --git a/HySoP/hysop/domain/subsets/porous.py b/HySoP/hysop/domain/porous.py
similarity index 98%
rename from HySoP/hysop/domain/subsets/porous.py
rename to HySoP/hysop/domain/porous.py
index d8fe97bab1070041f863712a80e6ec28bf9838d0..f57fe564374ede00bb0e3e43653a92c60131262f 100644
--- a/HySoP/hysop/domain/subsets/porous.py
+++ b/HySoP/hysop/domain/porous.py
@@ -1,10 +1,8 @@
"""Pre-defined subsets as combination of basic geometries.
"""
-from hysop.domain.subsets.subset import Subset
-from hysop.domain.subsets.boxes import SubBox
-from hysop.domain.subsets.sphere import Sphere, HemiSphere
-from hysop.domain.subsets.cylinder import Cylinder, HemiCylinder
+from hysop.domain.subsets import Subset, Sphere, HemiSphere
+from hysop.domain.subsets import Cylinder, HemiCylinder, SubBox
import hysop.tools.numpywrappers as npw
from hysop.mpi.topology import Cartesian
import numpy as np
diff --git a/HySoP/hysop/domain/subsets/submesh.py b/HySoP/hysop/domain/submesh.py
similarity index 100%
rename from HySoP/hysop/domain/subsets/submesh.py
rename to HySoP/hysop/domain/submesh.py
diff --git a/HySoP/hysop/domain/subsets/subset.py b/HySoP/hysop/domain/subsets.py
similarity index 53%
rename from HySoP/hysop/domain/subsets/subset.py
rename to HySoP/hysop/domain/subsets.py
index bacc5eeea0dceffd0120a4cf85b5c618e1f35639..572620b8f2ae60aad420419403bfc66ce230a7ff 100644
--- a/HySoP/hysop/domain/subsets/subset.py
+++ b/HySoP/hysop/domain/subsets.py
@@ -1,4 +1,4 @@
-""" subset of a given domain.
+""" subsets of a given domain (sphere ...)
"""
from hysop.domain.domain import Domain
import numpy as np
@@ -6,6 +6,7 @@ from hysop.mpi.topology import Cartesian
from hysop.fields.discrete import DiscreteField
import hysop.tools.numpywrappers as npw
from hysop.fields.continuous import Field
+from hysop.domain.submesh import SubMesh
class Subset(object):
@@ -18,8 +19,6 @@ class Subset(object):
data[subset.slices] = ...
or
data[subset.tab] = ...
- or
- subset.apply(data) = ...
slices : a list of python slices defining the grid points inside the subset
ind : a tuple of arrays of indices
@@ -86,6 +85,8 @@ class Subset(object):
# where on_proc[topo] is true. Useful to reduce collective comm
# on the subset
self.subcomm = {}
+ # A box "around" the subset.
+ self._box = None
def chi(self, *args):
"""Indicator function of the subset
@@ -113,7 +114,16 @@ class Subset(object):
assert isinstance(topo, Cartesian)
if topo not in self.ind:
dim = self._parent.dimension
- self.ind[topo] = [np.where(self.chi(*topo.mesh.coords))]
+ if self._box is None:
+ self.ind[topo] = [np.where(self.chi(*topo.mesh.coords))]
+ else:
+ self._box.discretize(topo)
+ mesh = self._box.mesh[topo]
+ shift = mesh.local_start
+ self.ind[topo] = [np.where(self.chi(*mesh.coords))]
+ self.ind[topo][0] = tuple([self.ind[topo][0][d] + shift[d]
+ for d in xrange(dim)])
+
self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
for i in xrange(dim)])
!= 0).all()
@@ -128,13 +138,6 @@ class Subset(object):
subgroup = gtopo.Incl(rks)
self.subcomm[topo] = topo.comm.Create(subgroup)
- def apply(self, field):
- """
- return field on the current subset
- not yet implemented
- """
- raise Exception('not yet implemented')
-
@staticmethod
def intersection(slist, topo):
"""Compute the intersection of subsets
@@ -460,3 +463,343 @@ class Subset(object):
result = npw.real_sum(field[component][self.ind[topo][0]])
result *= dvol
return result
+
+
+class Sphere(Subset):
+ """
+ Spherical domain.
+ """
+
+ def __init__(self, origin, radius=1.0, **kwds):
+ """
+ Parameters
+ ----------
+ origin : list or array
+ position of the center
+ radius : double, optional
+ kwds : base class parameters
+
+ """
+ def dist(*args):
+ size = len(args)
+ return npw.asarray(np.sqrt(sum([(args[d] - self.origin[d]) ** 2
+ for d in xrange(size)]))
+ - self.radius)
+
+ super(Sphere, self).__init__(chi=dist, **kwds)
+ # Radius of the sphere
+ self.radius = radius
+ # Center position
+ self.origin = npw.asrealarray(origin).copy()
+ dim = self._parent.dimension
+ # a subbox used to reduced the size of coords
+ # during dist computation.
+ self._box = SubBox(parent=self._parent,
+ length=[2 * self.radius, ] * dim,
+ origin=self.origin - self.radius)
+
+ def __str__(self):
+ s = self.__class__.__name__ + ' of radius ' + str(self.radius)
+ s += ' and center position ' + str(self.origin)
+ return s
+
+
+class HemiSphere(Sphere):
+ """
+ HemiSpherical domain.
+ Area defined by the intersection of a sphere and a box.
+ The box is defined with :
+ - cutdir, normal direction to a plan
+ - cutpos, position of this plan along the 'cutdir" axis
+ - all points of the domain where x < xs.
+ """
+ def __init__(self, cutpos=None, cutdir=0, **kwds):
+ """
+ Parameters
+ ----------
+ cutpos : list or array of coordinates
+ position of the cutting plane
+ cutdir : real, optional
+ direction of the normal to the cutting plane.
+ Default = x-direction (0).
+ """
+ super(HemiSphere, self).__init__(**kwds)
+ # direction of normal to the cutting plane
+ self.cutdir = cutdir
+ if cutpos is None:
+ cutpos = self.origin[self.cutdir]
+
+ def left_box(x):
+ return x - cutpos
+ self.LeftBox = left_box
+
+ def chi(self, *args):
+ return np.logical_and(self._is_inside(*args) <= 0,
+ self.LeftBox(args[self.cutdir]) <= 0)
+
+
+class Cylinder(Subset):
+ """
+ Cylinder-like domain.
+ """
+
+ def __init__(self, origin, radius=1.0, axis=1, **kwds):
+ """
+ Parameters
+ ----------
+ origin : list or array
+ coordinates of the center
+ radius : double, optional
+ default = 1.
+ axis : int, optional
+ direction of the main axis of the cylinder, default=1 (y)
+
+ """
+
+ def dist(*args):
+ size = len(self._dirs)
+ return npw.asarray(np.sqrt(sum([(args[self._dirs[d]] -
+ self.origin[self._dirs[d]]) ** 2
+ for d in xrange(size)]))
+ - self.radius)
+
+ super(Cylinder, self).__init__(chi=dist, **kwds)
+ ## Radius of the cylinder
+ self.radius = radius
+ ## Main axis position
+ self.origin = npw.asrealarray(origin).copy()
+ ## direction of the main axis of the cylinder
+ self.axis = axis
+ dim = self._parent.dimension
+ dirs = np.arange(dim)
+ self._dirs = np.where(dirs != self.axis)[0]
+ length = [2 * self.radius] * dim
+ length[self.axis] = self._parent.length[self.axis]
+ origin = self.origin - self.radius
+ origin[self.axis] = self._parent.origin[self.axis]
+ # a subbox used to reduced the size of coords
+ # during dist computation.
+ self._box = SubBox(parent=self._parent,
+ length=length, origin=origin)
+
+ def chi(self, *args):
+ return np.logical_and(self._is_inside(*args) <= 0,
+ args[self.axis] ==
+ args[self.axis])
+
+ def __str__(self):
+ s = self.__class__.__name__ + ' of radius ' + str(self.radius)
+ s += ' and center position ' + str(self.origin)
+ return s
+
+
+class HemiCylinder(Cylinder):
+ """
+ Half cylinder domain.
+ """
+ def __init__(self, cutpos=None, cutdir=0, **kwds):
+ """A cylinder cut by a plane (normal to one axis dir).
+
+ Parameters
+ ----------
+ cutpos : list or array of coordinates
+ position of the cutting plane
+ cutdir : real, optional
+ direction of the normal to the cutting plane.
+ Default = x-direction (0).
+
+ """
+ super(HemiCylinder, self).__init__(**kwds)
+ # direction of normal to the cutting plane
+ self.cutdir = cutdir
+ if cutpos is None:
+ cutpos = self.origin[self.cutdir]
+
+ def left_box(x):
+ return x - cutpos
+ self.LeftBox = left_box
+
+ def chi(self, *args):
+ return (np.logical_and(
+ np.logical_and(self._is_inside(*args) <= 0,
+ self.LeftBox(args[self.cutdir]) <= 0),
+ args[self.axis] == args[self.axis]))
+
+
+class SubBox(Subset):
+ """
+ A rectangle (in 2 or 3D space), defined by the coordinates of
+ its lowest point, its lenghts and its normal.
+ """
+ def __init__(self, origin, length, normal=1, **kwds):
+ """
+ Parameters
+ ----------
+ origin : list or array of double
+ position of the lowest point of the box
+ length : list or array of double
+ lengthes of the sides of the box
+ normal : int = 1 or -1, optional
+ direction of the outward normal. Only makes
+ sense when the 'box' is a plane.
+ **kwds : extra args for parent class
+
+ """
+ super(SubBox, self).__init__(**kwds)
+ ## Dictionnary of hysop.domain.mesh.Submesh, keys=topo, values=mesh
+ self.mesh = {}
+ ## coordinates of the lowest point of this subset
+ self.origin = npw.asrealarray(origin).copy()
+ ## length of this subset
+ self.length = npw.asrealarray(length).copy()
+ ## coordinates of the 'highest' point of this subset
+ self.end = self.origin + self.length
+ ## directions of coordinate axes belonging to the subset
+ self.t_dir = np.where(self.length != 0)[0]
+ ## directions of coordinate axes normal to the subset
+ self.n_dir = np.where(self.length == 0)[0]
+ ## direction of the outward unit normal (+ or - 1)
+ ## Useful when the surface belong to a control box.
+ self.normal = normal
+ msg = 'Subset error, the origin is outside of the domain.'
+ assert ((self.origin - self._parent.origin) >= 0).all(), msg
+ msg = 'Subset error, the subset is too large for the domain.'
+ assert ((self._parent.end - self.end) >= 0).all(), msg
+ ## dict of coordinates of the lengthes of the subset (values)
+ ## for a given topo (keys). If the origin/length does not
+ ## fit exactly with the mesh, there may be a small shift of
+ ## these values.
+ self.real_length = {}
+ ## dict of coordinates of the origin of the subset (values)
+ ## for a given topo (keys). If the origin/length does not
+ ## fit exactly with the mesh, there may be a small shift of
+ ## these values.
+ self.real_orig = {}
+
+ self.chi = None
+
+ def discretize(self, topo):
+ """
+ Compute a local submesh on the subset, for a given topology
+
+ :param topo: `~hysop.mpi.topology.Cartesian`
+ """
+ assert isinstance(topo, Cartesian)
+ # Find the indices of starting/ending points in the global_end
+ # grid of the mesh of topo.
+ gstart = topo.mesh.global_indices(self.origin)
+ gend = topo.mesh.global_indices(self.origin + self.length)
+ # create a submesh
+ self.mesh[topo] = SubMesh(topo.mesh, gstart, gend)
+ self.real_length[topo] = self.mesh[topo].global_length
+ self.real_orig[topo] = self.mesh[topo].global_origin
+ self.on_proc[topo] = self.mesh[topo].on_proc
+ self.ind[topo] = [self.mesh[topo].iCompute]
+ self._reduce_topology(topo)
+ return self.ind[topo]
+
+ def integrate_field_on_proc(self, field, topo, component=0):
+ """Field integration
+
+ Parameters
+ ----------
+ field : `~hysop.field.continuous.Field`
+ a field to be integrated on the box
+ topo : `~hysop.mpi.topology.Cartesian`
+ set mesh/topology used for integration
+ component : int, optional
+ number of the field component to be integrated
+
+ Returns
+ --------
+ integral of a component of the input field over the current subset,
+ on the current process
+ (i.e. no mpi reduce over all processes) for the discretization
+ given by topo..
+ """
+ assert isinstance(field, Field)
+ assert isinstance(topo, Cartesian)
+ df = field.discretize(topo)
+ dvol = npw.prod(topo.mesh.space_step[self.t_dir])
+ ic = self.mesh[topo].ind4integ
+ result = npw.real_sum(df[component][ic])
+ result *= dvol
+ return result
+
+ def integrate_func(self, func, topo, nbc, root=None):
+ """Function integration
+
+ Parameters
+ ----------
+ func: python function of space coordinates
+ topo: `~hysop.mpi.topology.Cartesian`
+ set mesh/topology used for integration
+ nbc : int
+ number of components of the return value from func
+
+ Returns
+ -------
+ integral of a component of the input field
+ over the current subset, for the discretization given by
+ topo
+ """
+ res = npw.zeros(nbc)
+ gres = npw.zeros(nbc)
+ for i in xrange(res.size):
+ res[i] = self.integrate_func_on_proc(func, topo)
+ if root is None:
+ topo.comm.Allreduce(res, gres)
+ else:
+ topo.comm.Reduce(res, gres, root=root)
+ return gres
+
+ def integrate_func_on_proc(self, func, topo):
+ """Function local integration
+
+ Parameters
+ ----------
+ func: python function of space coordinates
+ topo: `~hysop.mpi.topology.Cartesian`
+ set mesh/topology used for integration
+
+ Returns
+ --------
+ integral of the function on the subset on the current process
+ (i.e. no mpi reduce over all processes)
+ """
+ assert hasattr(func, '__call__')
+ assert isinstance(topo, Cartesian)
+ vfunc = np.vectorize(func)
+ if self.mesh[topo].on_proc:
+ result = npw.real_sum(vfunc(*self.mesh[topo].coords4int))
+ else:
+ result = 0.
+ dvol = npw.prod(topo.mesh.space_step)
+ result *= dvol
+ return result
+
+ def integrate_dfield_on_proc(self, field, component=0):
+ """Discrete field local integration
+
+ Parameters
+ ----------
+ field : `~hysop.field.discrete.DiscreteField`
+ a field to be integrated on the box
+ component : int, optional
+ number of the field component to be integrated
+ integrate the field on the subset on the current process
+ (i.e. no mpi reduce over all processes) for the discretization
+ given by topo.
+
+ Returns
+ --------
+ integral of the field on the subset on the current process
+ (i.e. no mpi reduce over all processes)
+ """
+ assert isinstance(field, DiscreteField)
+ topo = field.topology
+ dvol = npw.prod(topo.mesh.space_step[self.t_dir])
+ ic = self.mesh[topo].ind4integ
+ result = npw.real_sum(field[component][ic])
+ result *= dvol
+ return result
diff --git a/HySoP/hysop/domain/subsets/__init__.py b/HySoP/hysop/domain/subsets/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/HySoP/hysop/domain/subsets/boxes.py b/HySoP/hysop/domain/subsets/boxes.py
deleted file mode 100644
index 41b665250c99e73e97b82b39127921c8fa38c09e..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/subsets/boxes.py
+++ /dev/null
@@ -1,187 +0,0 @@
-"""2D rectangle or 3D Rectangular Cuboid subsets """
-from hysop.domain.subsets.subset import Subset
-from hysop.domain.subsets.submesh import SubMesh
-import hysop.tools.numpywrappers as npw
-from hysop import Field
-from hysop.fields.discrete import DiscreteField
-import numpy as np
-from hysop.mpi.topology import Cartesian
-
-
-class SubBox(Subset):
- """
- A rectangle (in 2 or 3D space), defined by the coordinates of
- its lowest point, its lenghts and its normal.
- """
- def __init__(self, origin, length, normal=1, **kwds):
- """
- Parameters
- ----------
- origin : list or array of double
- position of the lowest point of the box
- length : list or array of double
- lengthes of the sides of the box
- normal : int = 1 or -1, optional
- direction of the outward normal. Only makes
- sense when the 'box' is a plane.
- **kwds : extra args for parent class
-
- """
- super(SubBox, self).__init__(**kwds)
- ## Dictionnary of hysop.domain.mesh.Submesh, keys=topo, values=mesh
- self.mesh = {}
- ## coordinates of the lowest point of this subset
- self.origin = npw.asrealarray(origin).copy()
- ## length of this subset
- self.length = npw.asrealarray(length).copy()
- ## coordinates of the 'highest' point of this subset
- self.end = self.origin + self.length
- ## directions of coordinate axes belonging to the subset
- self.t_dir = np.where(self.length != 0)[0]
- ## directions of coordinate axes normal to the subset
- self.n_dir = np.where(self.length == 0)[0]
- ## direction of the outward unit normal (+ or - 1)
- ## Useful when the surface belong to a control box.
- self.normal = normal
- msg = 'Subset error, the origin is outside of the domain.'
- assert ((self.origin - self._parent.origin) >= 0).all(), msg
- msg = 'Subset error, the subset is too large for the domain.'
- assert ((self._parent.end - self.end) >= 0).all(), msg
- ## dict of coordinates of the lengthes of the subset (values)
- ## for a given topo (keys). If the origin/length does not
- ## fit exactly with the mesh, there may be a small shift of
- ## these values.
- self.real_length = {}
- ## dict of coordinates of the origin of the subset (values)
- ## for a given topo (keys). If the origin/length does not
- ## fit exactly with the mesh, there may be a small shift of
- ## these values.
- self.real_orig = {}
-
- self.chi = None
-
- def discretize(self, topo):
- """
- Compute a local submesh on the subset, for a given topology
-
- :param topo: `~hysop.mpi.topology.Cartesian`
- """
- assert isinstance(topo, Cartesian)
- # Find the indices of starting/ending points in the global_end
- # grid of the mesh of topo.
- gstart = topo.mesh.global_indices(self.origin)
- gend = topo.mesh.global_indices(self.origin + self.length)
- # create a submesh
- self.mesh[topo] = SubMesh(topo.mesh, gstart, gend)
- self.real_length[topo] = self.mesh[topo].global_length
- self.real_orig[topo] = self.mesh[topo].global_origin
- self.on_proc[topo] = self.mesh[topo].on_proc
- self.ind[topo] = [self.mesh[topo].iCompute]
- self._reduce_topology(topo)
- return self.ind[topo]
-
- def integrate_field_on_proc(self, field, topo, component=0):
- """Field integration
-
- Parameters
- ----------
- field : `~hysop.field.continuous.Field`
- a field to be integrated on the box
- topo : `~hysop.mpi.topology.Cartesian`
- set mesh/topology used for integration
- component : int, optional
- number of the field component to be integrated
-
- Returns
- --------
- integral of a component of the input field over the current subset,
- on the current process
- (i.e. no mpi reduce over all processes) for the discretization
- given by topo..
- """
- assert isinstance(field, Field)
- assert isinstance(topo, Cartesian)
- df = field.discretize(topo)
- dvol = npw.prod(topo.mesh.space_step[self.t_dir])
- ic = self.mesh[topo].ind4integ
- result = npw.real_sum(df[component][ic])
- result *= dvol
- return result
-
- def integrate_func(self, func, topo, nbc, root=None):
- """Function integration
-
- Parameters
- ----------
- func: python function of space coordinates
- topo: `~hysop.mpi.topology.Cartesian`
- set mesh/topology used for integration
- nbc : int
- number of components of the return value from func
-
- Returns
- -------
- integral of a component of the input field
- over the current subset, for the discretization given by
- topo
- """
- res = npw.zeros(nbc)
- gres = npw.zeros(nbc)
- for i in xrange(res.size):
- res[i] = self.integrate_func_on_proc(func, topo)
- if root is None:
- topo.comm.Allreduce(res, gres)
- else:
- topo.comm.Reduce(res, gres, root=root)
- return gres
-
- def integrate_func_on_proc(self, func, topo):
- """Function local integration
-
- Parameters
- ----------
- func: python function of space coordinates
- topo: `~hysop.mpi.topology.Cartesian`
- set mesh/topology used for integration
-
- Returns
- --------
- integral of the function on the subset on the current process
- (i.e. no mpi reduce over all processes)
- """
- assert hasattr(func, '__call__')
- assert isinstance(topo, Cartesian)
- vfunc = np.vectorize(func)
- if self.mesh[topo].on_proc:
- result = npw.real_sum(vfunc(*self.mesh[topo].coords4int))
- else:
- result = 0.
- dvol = npw.prod(topo.mesh.space_step)
- result *= dvol
- return result
-
- def integrate_dfield_on_proc(self, field, component=0):
- """Discrete field local integration
-
- Parameters
- ----------
- field : `~hysop.field.discrete.DiscreteField`
- a field to be integrated on the box
- component : int, optional
- number of the field component to be integrated
- integrate the field on the subset on the current process
- (i.e. no mpi reduce over all processes) for the discretization
- given by topo.
-
- Returns
- --------
- integral of the field on the subset on the current process
- (i.e. no mpi reduce over all processes)
- """
- assert isinstance(field, DiscreteField)
- topo = field.topology
- dvol = npw.prod(topo.mesh.space_step[self.t_dir])
- ic = self.mesh[topo].ind4integ
- result = npw.real_sum(field[component][ic])
- result *= dvol
- return result
diff --git a/HySoP/hysop/domain/subsets/cylinder.py b/HySoP/hysop/domain/subsets/cylinder.py
deleted file mode 100644
index ff73c96cf9696696dd2e6e1eabf8d3225137db15..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/subsets/cylinder.py
+++ /dev/null
@@ -1,119 +0,0 @@
-"""Spherical or hemispherical subset. """
-from hysop.domain.subsets.subset import Subset
-import numpy as np
-import hysop.tools.numpywrappers as npw
-from hysop.mpi.topology import Cartesian
-from hysop.domain.subsets.boxes import SubBox
-
-
-class Cylinder(Subset):
- """
- Cylinder-like domain.
- """
-
- def __init__(self, origin, radius=1.0, axis=1, **kwds):
- """
- Parameters
- ----------
- origin : list or array
- coordinates of the center
- radius : double, optional
- default = 1.
- axis : int, optional
- direction of the main axis of the cylinder, default=1 (y)
-
- """
-
- def dist(*args):
- size = len(self._dirs)
- return npw.asarray(np.sqrt(sum([(args[self._dirs[d]] -
- self.origin[self._dirs[d]]) ** 2
- for d in xrange(size)]))
- - self.radius)
-
- super(Cylinder, self).__init__(chi=dist, **kwds)
- ## Radius of the cylinder
- self.radius = radius
- ## Main axis position
- self.origin = npw.asrealarray(origin).copy()
- ## direction of the main axis of the cylinder
- self.axis = axis
- dim = self._parent.dimension
- dirs = np.arange(dim)
- self._dirs = np.where(dirs != self.axis)[0]
- length = [2 * self.radius] * dim
- length[self.axis] = self._parent.length[self.axis]
- origin = self.origin - self.radius
- origin[self.axis] = self._parent.origin[self.axis]
- # a subbox used to reduced the size of coords
- # during dist computation.
- self._box = SubBox(parent=self._parent,
- length=length, origin=origin)
-
- def discretize(self, topo):
- assert isinstance(topo, Cartesian)
- if topo not in self.ind:
- self._box.discretize(topo)
- mesh = self._box.mesh[topo]
- shift = mesh.local_start
- dim = self._parent.dimension
- tmp = np.where(self.chi(*mesh.coords))
- self.ind[topo] = [tuple([tmp[d] + shift[d] for d in xrange(dim)])]
- self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
- for i in xrange(dim)])
- != 0).all()
- self._reduce_topology(topo)
- return self.ind[topo]
-
- def chi(self, *args):
- return np.logical_and(self._is_inside(*args) <= 0,
- args[self.axis] ==
- args[self.axis])
-
- def __str__(self):
- s = self.__class__.__name__ + ' of radius ' + str(self.radius)
- s += ' and center position ' + str(self.origin)
- return s
-
-
-class HemiCylinder(Cylinder):
- """
- Half cylinder domain.
- """
- def __init__(self, cutpos=None, cutdir=0, **kwds):
- """
- A cylinder cut by a plane (normal to one axis dir).
- @param cutpos : position of the cutting plane
- @param cutdir : direction of the normal to the cutting plane.
- Default = x-direction (0).
- """
- super(HemiCylinder, self).__init__(**kwds)
- ## direction of normal to the cutting plane
- self.cutdir = cutdir
- if cutpos is None:
- cutpos = self.origin[self.cutdir]
-
- def LeftBox(x):
- return x - cutpos
- self.LeftBox = LeftBox
-
- def chi(self, *args):
- return (np.logical_and(
- np.logical_and(self._is_inside(*args) <= 0,
- self.LeftBox(args[self.cutdir]) <= 0),
- args[self.axis] == args[self.axis]))
-
- def discretize(self, topo):
- if topo not in self.ind:
- self._box.discretize(topo)
- mesh = self._box.mesh[topo]
- shift = mesh.local_start
- dim = self._parent.dimension
- tmp = np.where(self.chi(*mesh.coords))
- self.ind[topo] = [tuple([tmp[d] + shift[d] for d in xrange(dim)])]
- self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
- for i in xrange(dim)])
- != 0).all()
- self._reduce_topology(topo)
-
- return self.ind[topo]
diff --git a/HySoP/hysop/domain/subsets/sphere.py b/HySoP/hysop/domain/subsets/sphere.py
deleted file mode 100644
index c187584d14773c44383c5210456b31346e4a85d1..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/subsets/sphere.py
+++ /dev/null
@@ -1,115 +0,0 @@
-"""Spherical or hemispherical subsets in 2 or 3D domains.
-
-"""
-from hysop.domain.subsets.subset import Subset
-import numpy as np
-import hysop.tools.numpywrappers as npw
-from hysop.mpi.topology import Cartesian
-from hysop.domain.subsets.boxes import SubBox
-
-
-class Sphere(Subset):
- """
- Spherical domain.
- """
-
- def __init__(self, origin, radius=1.0, **kwds):
- """
- Parameters
- ----------
- origin : list or array
- position of the center
- radius : double, optional
- kwds : base class parameters
-
- """
- def dist(*args):
- size = len(args)
- return npw.asarray(np.sqrt(sum([(args[d] - self.origin[d]) ** 2
- for d in xrange(size)]))
- - self.radius)
-
- super(Sphere, self).__init__(chi=dist, **kwds)
- # Radius of the sphere
- self.radius = radius
- # Center position
- self.origin = npw.asrealarray(origin).copy()
- dim = self._parent.dimension
- # a subbox used to reduced the size of coords
- # during dist computation.
- self._box = SubBox(parent=self._parent,
- length=[2 * self.radius, ] * dim,
- origin=self.origin - self.radius)
-
- def discretize(self, topo):
- assert isinstance(topo, Cartesian)
- if topo not in self.ind:
- self._box.discretize(topo)
- mesh = self._box.mesh[topo]
- shift = mesh.local_start
- dim = self._parent.dimension
- tmp = np.where(self.chi(*mesh.coords))
- self.ind[topo] = [tuple([tmp[d] + shift[d]
- for d in xrange(dim)])]
- self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
- for i in xrange(dim)])
- != 0).all()
- self._reduce_topology(topo)
-
- return self.ind[topo]
-
- def __str__(self):
- s = self.__class__.__name__ + ' of radius ' + str(self.radius)
- s += ' and center position ' + str(self.origin)
- return s
-
-
-class HemiSphere(Sphere):
- """
- HemiSpherical domain.
- Area defined by the intersection of a sphere and a box.
- The box is defined with :
- - cutdir, normal direction to a plan
- - cutpos, position of this plan along the 'cutdir" axis
- - all points of the domain where x < xs.
- """
- def __init__(self, cutpos=None, cutdir=0, **kwds):
- """
- Parameters
- ----------
- cutpos : list or array of coordinates
- position of the cutting plane
- cutdir : real, optional
- direction of the normal to the cutting plane.
- Default = x-direction (0).
- """
- super(HemiSphere, self).__init__(**kwds)
- # direction of normal to the cutting plane
- self.cutdir = cutdir
- if cutpos is None:
- cutpos = self.origin[self.cutdir]
-
- def left_box(x):
- return x - cutpos
- self.LeftBox = left_box
-
- def chi(self, *args):
- return np.logical_and(self._is_inside(*args) <= 0,
- self.LeftBox(args[self.cutdir]) <= 0)
-
- def discretize(self, topo):
- assert isinstance(topo, Cartesian)
- if topo not in self.ind:
- self._box.discretize(topo)
- mesh = self._box.mesh[topo]
- dim = self._parent.dimension
- shift = mesh.local_start
- self.ind[topo] = [np.where(self.chi(*mesh.coords))]
- self.ind[topo][0] = tuple([self.ind[topo][0][d] + shift[d]
- for d in xrange(dim)])
- self.on_proc[topo] = (np.asarray([len(self.ind[topo][0][i])
- for i in xrange(dim)])
- != 0).all()
- self._reduce_topology(topo)
-
- return self.ind[topo]
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
index 43577a5f96fcd31c858cbff477c78597e9b43357..73cc589ceefaa54b5ddbb5ded91c0806d277af62 100644
--- a/HySoP/hysop/domain/tests/test_subset.py
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -9,8 +9,6 @@ from hysop.mpi.topology import Cartesian
import hysop.tools.numpywrappers as npw
import numpy as np
import math
-from hysop.mpi import main_size
-import mpi4py
Nx = 128
@@ -53,10 +51,10 @@ def check_subset(discr, fileref, class_name):
assert hasattr(subs.chi, '__call__')
subs.discretize(topo)
assert subs.ind.keys()[0] == topo
- assert len(subs.ind[topo]) == dom.dimension
+ assert len(subs.ind[topo][0]) == dom.dimension
scal = Field(domain=topo.domain, is_vector=False)
sd = scal.discretize(topo)
- sd[0][subs.ind[topo]] = 2.
+ sd[0][subs.ind[topo][0]] = 2.
iCompute = topo.mesh.iCompute
return np.allclose(sd[0][iCompute], vd[0][iCompute])
@@ -123,11 +121,11 @@ def init_obst_list(discr, fileref):
return np.allclose(sd[0][iCompute], vd[0][iCompute])
-def test_union_3D():
+def test_union_3d():
assert init_obst_list(discr3D, 'multi_obst_3d')
-def test_union_2D():
+def test_union_2d():
assert init_obst_list(discr2D, 'multi_obst_2d')
@@ -148,11 +146,11 @@ def init_subtract(discr, fileref):
return np.allclose(sd[0][iCompute], vd[0][iCompute])
-def test_subtract_3D():
+def test_subtract_3d():
assert init_subtract(discr3D, 'multi_obst_subs_3d')
-def test_subtract_2D():
+def test_subtract_2d():
assert init_subtract(discr2D, 'multi_obst_subs_2d')
@@ -185,11 +183,11 @@ def init_subtract_lists(discr, fileref):
return np.allclose(sd[0][topo.mesh.iCompute], vd[0][topo.mesh.iCompute])
-def test_subtract_list_3D():
+def test_subtract_list_3d():
assert init_subtract_lists(discr3D, 'multi_obst_subslist_3d')
-def test_subtract_list_2D():
+def test_subtract_list_2d():
assert init_subtract_lists(discr2D, 'multi_obst_subslist_2d')
@@ -201,9 +199,9 @@ if __name__ == "__main__":
test_hemicylinder_3d()
test_sphere_2d()
test_hemisphere_2d()
- test_union_2D()
- test_union_3D()
- test_subtract_3D()
- test_subtract_2D()
- test_subtract_list_3D()
- test_subtract_list_2D()
+ test_union_2d()
+ test_union_3d()
+ test_subtract_3d()
+ test_subtract_2d()
+ test_subtract_list_3d()
+ test_subtract_list_2d()
diff --git a/HySoP/hysop/operator/tests/test_penalization.py b/HySoP/hysop/operator/tests/test_penalization.py
index b0cfb869253d6df3f9d913ec28c6d5caffaacd78..50e5915585f1ea65902153f7f0409721b3c403c5 100644
--- a/HySoP/hysop/operator/tests/test_penalization.py
+++ b/HySoP/hysop/operator/tests/test_penalization.py
@@ -1,7 +1,6 @@
# -*- coding: utf-8 -*-
-from hysop.domain.subsets.sphere import HemiSphere, Sphere
-from hysop.domain.subsets.cylinder import Cylinder
-from hysop.domain.subsets.porous import Porous
+from hysop.domain.subsets import HemiSphere, Sphere, Cylinder
+from hysop.domain.porous import Porous
from hysop.operator.penalization import Penalization
from hysop.operator.penalize_vorticity import PenalizeVorticity
from hysop.problem.simulation import Simulation
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index dc78957d1a26ac017168d7de0d389be088c71836..b2489f6c80d3d38cc255be212f218ec1d3b479e2 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -13,7 +13,6 @@ name = '@PYPACKAGE_NAME@'
# List of modules (directories) to be included
packages = ['hysop',
'hysop.domain',
- 'hysop.domain.subsets',
'hysop.fields',
'hysop.operator',
'hysop.operator.discrete',