diff --git a/HySoP/hysop/domain/obstacle/__init__.py b/HySoP/hysop/domain/obstacle/__init__.py
deleted file mode 100644
index be5911bbcc91dc597c5c18e7479042e779dfd672..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/__init__.py
+++ /dev/null
@@ -1,57 +0,0 @@
-## @package parmepy.domain.obstacle
-# Obstacles description (geometry).
-#
-#
-# An 'obstacle' is the description of a sub-domain
-# of the main physical domain with one or more user-defined
-# functions of the space coordinates (at least).
-#
-# What for? \n
-# Mainly to provide some 'index sets' to penalization operator.\n
-# For a given obstacle, and a given discrete field, we must be able to call:
-# \code
-# cond = obstacle.discretize(topo)
-# field[cond] = 1.9
-# \endcode
-# This example means that the field (discretized on topology topo)
-# will be set to 1.9 everywhere inside the obstacle.
-#
-#
-# Obviouslvy the index sets will depend on the discretization
-# of the domain (the underlying topology indeed).
-# So each obstacle handles a dictionnary of boolean arrays. The keys
-# of the dictionnaries are the topologies and the values some boolean arrays
-# which values are true on and inside the object and false outside.
-# \code
-# obstacle.ind[topo] = someBooleanArray
-# \endcode
-# Each component of the dictionnary is created using the method 'discretize':
-# \code
-# # Create the boolean array that represents the obstacle for topology topo:
-# someBooleanArray = obstacle.discretize(topo)
-# # So for a field already discretized on topo, we may call
-# field[someBooleanArray]
-# \endcode
-#
-# A more complete example : initialize a scalar field with one inside a sphere
-# and zero everywhere else.
-#
-# \code
-# Lx = Ly = Lz = 2
-# dom = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[-1., -1., -1.])
-# # Definition of a sphere in the middle of the domain
-# sphere = Sphere(dom, position=[0., 0., 0.], radius=0.5)
-# # A topology
-# topo = Cartesian(dom, 3, [33, 33, 33])
-# # A scalar field on the domain :
-# scal = pp.Field(domain=dom, name='Scalar')
-# # Discretization on topo:
-# scal_discr = scal.discretize(topo)
-# # scal set to 1. inside the obstacle:
-# condition = sphere.discretize(topo)
-# scal.discreteFields[topo][condition] = 1.
-# # equivalent to
-# scal_discr[condition] = 1.
-# # to set a value everywhere except in the sphere
-# scal_discr[not condition] = 8.
-# \endcode
diff --git a/HySoP/hysop/domain/obstacle/controlBox.py b/HySoP/hysop/domain/obstacle/controlBox.py
deleted file mode 100644
index 9c982da5668fc5fbc563c0ca490a847ea00b5c64..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/controlBox.py
+++ /dev/null
@@ -1,300 +0,0 @@
-"""
-@file controlBox.py
-Define a sub-domain with a box-liked shape.
-"""
-from parmepy.domain.obstacle.obstacle import Obstacle
-from parmepy.domain.obstacle.planes import SubSpace, SubPlane
-from parmepy.mpi.mesh import SubMesh
-import numpy as np
-import parmepy.tools.numpywrappers as npw
-
-
-class ControlBox(Obstacle):
- """
- Build a sub-domain, box-shaped
- ==> define set of indices inside this domain (ind member)
- and set of indices belonging to surfaces of this domain (slices members).
- Useful to define control volume to perform integration.
- """
-
- def __init__(self, origin, lengths, **kwds):
- """
- Build the volume of control
- @param origin : coordinates of the lowest point in the sub-domain
- @param lengths : lengths of box sides.
- """
- super(ControlBox, self).__init__(**kwds)
-
- ## Lowest point of the box
- self.origin = npw.asrealarray(origin)
- ## Box's sides dimension
- self.lengths = npw.asrealarray(lengths)
- ## Dictionnary of local meshes, such that
- ## mesh[topo] is the restriction of topo.mesh
- ## to the current control box.
- self.mesh = {}
- self.upper = None
- self.lower = None
- self.upperS = None
- self.lowerS = None
- self.slices = {}
- self.indReduced = {}
- self._boxCreated = False
- ## Check if the defined box contains points
- ## for a given topology. self.isEmpty[topo] = False
- ## if some grid points are inside the box on
- ## the current processor for topo discretization.
- self.isEmpty = {}
- ## Dict of local coordinates for a given topology
- self.coords = {}
- ## Global resolution of the obstacle (plane, sub space ...)
- ## At the time, it's computed only for subspaces, by calling
- ## globalResolution method.
- self.gRes = None
- ## Global index in the original topology of the 'lowest' point
- ## of the obstacle. Only for subspaces.
- self.gstart = None
-
- def _createVolumeAndSides(self, spaceStep):
- """
- @param[in] array of size self._dim, space step size in each direction
- This value will be used to compute a tolerance and detect
- points inside the box.
- """
- # Build Half-spaces indices list in all directions
- normalUp = np.identity(self._dim)
- normalDown = np.identity(self._dim) * -1
- pointsUp = npw.zeros((self._dim, self._dim))
- # Control box will be used for integration, so we remove the
- # last point in the grid.
- boxlengths = self.lengths - spaceStep
- tol = spaceStep * 0.5
-
- for i in xrange(self._dim):
- pointsUp[:, i] = self.origin
- pointsUp.flat[::self._dim + 1] += self.lengths
- # -- Control volume : union of two halfspaces --
- if self.upper is None:
- self.upper = [SubSpace(domain=self.domain, normal=normalUp[:, i],
- point=pointsUp[:, i],
- lengths=boxlengths,
- epsilon=tol[i])
- for i in xrange(self._dim)]
- if self.lower is None:
- self.lower = [SubSpace(domain=self.domain, normal=normalDown[:, i],
- point=self.origin, lengths=boxlengths,
- epsilon=tol[i])
- for i in xrange(self._dim)]
-
- # Create objects to describe the sides of the box
- if self.upperS is None:
- self.upperS = [SubPlane(domain=self.domain, normal=normalUp[:, i],
- point=pointsUp[:, i],
- lengths=boxlengths,
- epsilon=tol[i])
- for i in xrange(self._dim)]
-
- if self.lowerS is None:
- self.lowerS = [SubPlane(domain=self.domain,
- normal=normalDown[:, i],
- point=self.origin, lengths=boxlengths,
- epsilon=tol[i])
- for i in xrange(self._dim)]
- self._boxCreated = True
-
- def discretize(self, topo):
- """
- Discretize the box volume and its surfaces.
- @param topo : the topology that described the discretization.
- """
- # Check if already done. If so, this function has no effect.
- if topo not in self.ind.keys():
- spaceStep = topo.mesh.space_step
- # -- Control volume : union of two halfspaces --
- if not self._boxCreated:
- self._createVolumeAndSides(spaceStep)
-
- # Discretize all volume and surfaces of
- # the box for topo
- for i in xrange(self._dim):
- self.lower[i].discretize(topo)
- self.upper[i].discretize(topo)
- self.lowerS[i].discretize(topo)
- self.upperS[i].discretize(topo)
-
- # 1 -- Compute list of indices inside the box,
- # for topo --> ind[topo]
- self.ind[topo] = []
-
- self.ind[topo].append(np.logical_and(self.upper[0].ind[topo][0],
- self.lower[0].ind[topo][0]))
- for i in xrange(1, self._dim):
- cond = np.logical_and(self.upper[i].ind[topo][0],
- self.lower[i].ind[topo][0])
- self.ind[topo][0] = np.logical_and(self.ind[topo][0], cond)
-
- ind = np.where(self.ind[topo][0])
-
- # 2 -- Convert ind[topo] (array of bool) to slices
- # which may be more convenient for computations
- # --> slices[topo]
- # + mesh[topo], a parmepy.mpi.SubMesh, useful
- # to get local coordinates and so on
- if ind[0].size == 0:
- self.slices[topo] = [slice(0, 0) for i in xrange(self._dim)]
- self.mesh[topo] = None
- self.isEmpty[topo] = True
- else:
- self.isEmpty[topo] = False
- ic = topo.mesh.iCompute
- lstart = [ind[i].min() if ind[i].size > 0 else None
- for i in xrange(self._dim)]
- lstart = npw.asintarray([max(lstart[i], ic[i].start)
- for i in xrange(self._dim)])
- end = [ind[i].max() for i in xrange(self._dim)]
- end = npw.asintarray([min(end[i], ic[i].stop - 1)
- for i in xrange(self._dim)])
- # slice(start,end) --> end not included, so +1
- end += 1
- resol = end - lstart + 2 * topo.mesh.discretization.ghosts
- gstart = lstart + topo.mesh.start()
- gstart -= topo.mesh.discretization.ghosts
- self.mesh[topo] = SubMesh(self.domain,
- topo.mesh.discretization,
- gstart, resol)
- self.slices[topo] = [slice(lstart[i], end[i])
- for i in xrange(self._dim)]
- coords = []
- for i in xrange(self._dim):
- cc = topo.mesh.coords[i].flat[self.slices[topo][i]]
- coords.append(cc)
- coords = tuple(coords)
- self.coords[topo] = np.ix_(*coords)
-
- # --> self.ind[topo][0] components are True
- # for points inside the volume
- # --> self.slices[topo] represent the same thing
- # but using slices of numpy arrays.
- # Usage (vd being a numpy array discretized
- # on the whole domain, cb a control box):
- # # Set values to all points inside the control box
- # vd[cb.ind[topo][0]] = ...
- # # Get a sub-array of vd representing the control box
- # # and use it
- # result[...] = vd[cb.slices] + ...
- # The important difference between slices and ind is:
- # 1 - vd[ind] returns a 1D array whatever vd shape is.
- # 2 - vd[slices] return an array of the same dim as vd,
- # with shape given by slices.
-
- return self.ind[topo]
-
- def sub(self, obstacle, topo):
- """
- Remove all points corresponding to the input obstacle from
- the current control box
- """
- if topo not in self.ind.keys():
- obstacle.discretize(topo)
- self.discretize(topo)
- self.indReduced[topo] = []
- # Warning : obstacle may have several layers
- cond = obstacle.ind[topo][0]
- for i in xrange(1, len(obstacle.ind[topo])):
- cond = npw.asboolarray(np.logical_or(cond,
- obstacle.ind[topo][i]))
- cond = np.logical_not(cond)
- self.indReduced[topo].append(np.logical_and(self.ind[topo][0],
- cond))
- return self.indReduced[topo][-1]
-
- def integrate_on_proc(self, field, topo, useSlice=True, component=0):
- """
- integrate field on the box
- """
- if useSlice:
- cond = self.slices[topo]
- else:
- iC = topo.mesh.iCompute
- cond = self.ind[topo][0][iC]
- dvol = npw.prod(topo.mesh.space_step)
- result = npw.real_sum(field.discretize(topo)[component][cond])
- result *= dvol
- return result
-
- def integrate(self, field, topo, useSlice=True,
- component=0, root=0, mpiall=True):
- res = self.integrate_on_proc(field, topo, useSlice, component)
- if mpiall:
- return topo.comm.allreduce(res)
- else:
- return topo.comm.reduce(res, root=root)
-
- def integrateOnSurface(self, field, topo, normalDir=0, up=True,
- useSlice=True, component=0, root=0, mpiall=True):
- """
- integrate field on top (if up is True) or down surface
- normal to a direction
- """
- res = self.integrateOnSurf_proc(field, topo, normalDir, up, useSlice,
- component)
- if mpiall:
- return topo.comm.allreduce(res)
- else:
- return topo.comm.reduce(res, root=root)
-
- def integrateOnSurf_proc(self, field, topo, normalDir=0,
- up=True, useSlice=True, component=0):
- """
- integrate field on top and down surfaces normal to a direction
- """
- if up:
- surf = self.upperS[normalDir]
- else:
- surf = self.lowerS[normalDir]
- if useSlice:
- cond = surf.slices[topo]
- else:
- iC = topo.mesh.iCompute
- cond = surf.ind[topo][0][iC]
- dirs = np.logical_not(np.arange(self._dim) == normalDir)
- dS = npw.prod(topo.mesh.space_step[dirs])
- result = npw.real_sum(field.discretize(topo)[component][cond])
- result *= dS
- return result
-
- def globalResolution(self, parent_topo):
- """
- Compute 'global resolution' of the subplane
- """
- # We create a false topology, with only one proc
- # to get the global resolution for the plane.
- # This could also be done with local computation
- # sum but that would need a lot of communications.
- if parent_topo.rank == 0:
- color = 0
- else:
- color = 1
- subcomm = parent_topo.comm.Split(color)
- dimension = self.domain.dimension
- tmp = None
- if parent_topo.rank == 0:
- resolution = parent_topo.globalMeshResolution
- ghosts = parent_topo.mesh.discretization.ghosts
- topo = self.domain.getOrCreateTopology(3, resolution,
- ghosts=ghosts, comm=subcomm)
- self.discretize(topo)
- sl = self.slices[topo]
- self.gRes = [sl[i].stop - sl[i].start for i in xrange(dimension)]
- self.gstart = [sl[i].start for i in xrange(dimension)]
- # if the topology has been created just to
- # get the global resolution, we can remove it
- if topo.isNew:
- self.domain.remove(topo)
- self.slices.pop(topo)
- self.ind.pop(topo)
- tmp = self.gRes + self.gstart
- tmp = parent_topo.comm.bcast(tmp)
- self.gRes = tmp[:dimension]
- self.gstart = tmp[dimension:]
- return self.gRes
diff --git a/HySoP/hysop/domain/obstacle/disk.py b/HySoP/hysop/domain/obstacle/disk.py
deleted file mode 100644
index f57cbcc9284fac1aeab3a38190cc1b8f8490694e..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/disk.py
+++ /dev/null
@@ -1,64 +0,0 @@
-"""
-@file disk.py
-Rigid disk (2D)
-"""
-from parmepy.domain.obstacle.sphere import Sphere, HemiSphere
-import numpy as np
-import parmepy.tools.numpywrappers as npw
-
-
-class Disk(Sphere):
- """
- Disk in a 2D domain.
- """
-
- def __init__(self, **kwds):
- """
- Description of a disk in a domain.
- @param domain : the physical domain that contains the sphere.
- @param position : position of the center
- @param radius : sphere radius, default = 1
- @param porousLayers : a list of thicknesses
- for successive porous layers
- radius is the inside sphere radius and thicknesses are given from
- inside layer to outside one.
- @param vd : velocity of the disk (considered as a rigid body),
- default = 0.
- """
- super(Disk, self).__init__(**kwds)
- assert self.domain.dimension == 2
-
- def dist(x, y, R):
- return npw.asarray(np.sqrt((x - self.position[0]) ** 2
- + (y - self.position[1]) ** 2) - R)
- self.chi = [dist]
-
-
-class HalfDisk(HemiSphere):
- """
- Half disk in a 2D domain.
- """
- def __init__(self, **kwds):
- """
- Constructor for the semi-disk.
- @param domain : the physical domain that contains the sphere.
- @param position : position of the center
- @param radius : sphere radius, default = 1
- (if box ...)
- @param vd : velocity of the disk (considered as a rigid body),
- default = 0.
- """
- super(HalfDisk, self).__init__(**kwds)
- assert self.domain.dimension == 2
-
- def dist(x, y, R):
- """
- """
- return npw.asarray(np.sqrt((x - self.position[0]) ** 2
- + (y - self.position[1]) ** 2) - R)
- self.chi = [dist]
-
- def LeftBox(x, y):
- return x - self.position[0]
-
- self.LeftBox = LeftBox
diff --git a/HySoP/hysop/domain/obstacle/obstacle.py b/HySoP/hysop/domain/obstacle/obstacle.py
deleted file mode 100644
index c052ed7aa218ec544d66958a060d55611592a888..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/obstacle.py
+++ /dev/null
@@ -1,96 +0,0 @@
-"""@file obstacle.py
-
-General interface to define a new geometry
-inside a domain (sphere, control box ...)
-"""
-import numpy as np
-
-
-class Obstacle(object):
- """Ddescription of a physical obstacle.
- An obstacle is the geometrical description of
- a physical sub-domain.
- """
- def __init__(self, domain, formula=None, vd=0.0):
- """ Constructor
- @param domain : the domain that contains this obstacle.
- @param formula : a list of functions that describe
- the geometry of the obstacle.
- @param vd : velocity of the obstacle (considered as a rigid body),
- default = 0.
- """
- ## Domain.
- self.domain = domain
- from parmepy.domain.box import Box
- assert isinstance(domain, Box),\
- 'Obstacle only implemented for box-like domains'
- ## Obstacle dimension.
- self._dim = domain.dimension
- ## A function that describe the geometry of the obstacle.
- ## see parmepy.domain.obstacle.
- self.chi = []
- if formula is not None:
- if isinstance(formula, list):
- for func in formula:
- self.chi.append = np.vectorize(func)
- else:
- self.chi = [np.vectorize(formula)]
- ## A dictionnary of lists of indices ...
- ## ind[topo][i] represents the set of points of the domain
- ## discretized with topo that are in the area defined with chi[i].
- self.ind = {}
- ## Velocity of the center of mass of the obstacle
- ## (considered as a rigid body)
- self.vd = vd
- ## Check if some grid points are present inside the current object
- ## for the current mpi proc. If not, isEmpty[topo] = True.
- self.isEmpty = {}
-
-
- def discretize(self, topo):
- """
- For a given topology, computes the list of points in the domain
- that belongs to the obstacle.
- Add a new index into self.indices.
- @param topo : topology specification for discretization
- @return an array of bool such that array[i,j,k] = True if
- point(i,j,k) is in the obstacle.
-
- Note FP : there are two ways to 'save' which points are
- in the obstacle : either we set a test function and
- fill a boolean numpy array (case A) or we compute domain.dimension
- lists of indice (B).
- - case A : indices[topo] is a numpy array of the same
- size as topo.mesh.
- +++ : very fast to compute, very fast to apply
- --- : needs more memory (size of bool * topo.mesh.size)
- - case B : indices[topo] is a tupple of lists, like (ix, iy, iz).
- A point inside the obstacle is thus given
- by the indices ix[i], iy[i], iz[i]
- +++ : needs less memory
- --- : very slow to compute/apply compared to case A
- Default choice = case A
- \todo : provide a way for user to choose between case A and B.
- Note FP 2: maybe that would be better to save indices in
- operator (penalization) and not in obstacle, to save memory?
- """
- # first check if we have already compute indices for
- # this topology
- if topo not in self.ind[0].keys():
- # for each indicator function
- self.ind[topo] = []
- for func in self.chi:
- # current function position
- i = self.chi.index(func)
- # apply indicator function on local mesh for the required topo
- self.ind[topo].append(self.chi[i](*topo.mesh.coords) <= 0)
-
- return self.ind[topo]
-
- def _isempty(self, topo):
- ilist = np.where(self.ind[topo])
- if ilist[0].size == 0:
- self.isEmpty[topo] = True
- else:
- self.isEmpty[topo] = False
-
diff --git a/HySoP/hysop/domain/obstacle/planes.py b/HySoP/hysop/domain/obstacle/planes.py
deleted file mode 100644
index e5f2edbdb29bec969222b413441d27be727708fa..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/planes.py
+++ /dev/null
@@ -1,359 +0,0 @@
-"""
-@file planes.py
-Plate-like sub-domains at boundaries, normal
-to a given direction.
-"""
-from parmepy.domain.obstacle.obstacle import Obstacle
-import numpy as np
-import parmepy.tools.numpywrappers as npw
-
-
-class HalfSpace(Obstacle):
- """
- Divide domain into two sub-spaces, on each side of a plane
- defined by its normal and a point.
- Indices of this.ind describe the half-space below the plane,
- 'normal' being the outward normal of the plane.
- """
-
- def __init__(self, normal, point, epsilon=1e-2, **kwds):
- """
- Half space define by the points of the domain on one side
- of a plane.
- @param domain : the physical domain that contains the plane
- @param normal : outward normal
- @param point : coordinates of a point of the plane.
- @param epsilon : tolerance
- """
- super(HalfSpace, self).__init__(**kwds)
- assert epsilon > 0.0, 'Tolerance value must be positive'
- ## Tolerance used to considered that points at the boundary are
- ## in the subspace. Good choice may be grid space_step / 2.
- self.epsilon = epsilon
- ## Direction of the normal to the plate (0:x, 1:y, 2:z))
- ## normal is the 'outer' normal of the 'in' subspace.
- self.normal = npw.asintarray(normal)
- self.point = point
- self.origin = npw.asrealarray(point)
-
- def Outside(*coords):
- return sum([(coords[i] - self.point[i]) * self.normal[i]
- for i in xrange(self.domain.dimension)])
-
- ## Test function for half-space.
- ## Positive value if outside subdomain else negative
- self.chi = [Outside]
- self.slices = {}
- ## Global resolution of the obstacle (plane, sub space ...)
- ## At the time, it's computed only for subspaces, by calling
- ## globalResolution method.
- self.gRes = None
- ## Global index in the original topology of the 'lowest' point
- ## of the obstacle. Only for subspaces.
- self.gstart = None
-
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
-
- if topo not in self.ind.keys():
- self.ind[topo] = []
- # apply indicator function on topo local mesh
- cond = npw.asarray(self.chi[0](*topo.mesh.coords) <= self.epsilon)
- self.ind[topo].append(cond)
- self._isempty(topo)
- return self.ind[topo]
-
- def __str__(self):
- s = 'Plane normal to vector' + str(self.normal)
- s += ' going through point ' + str(self.point)
- return s
-
-
-class Plane(HalfSpace):
- """
- A plane in the domain, defined by its normal and a point.
- Indices of plane.ind describe the points belonging to the plane.
- """
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
-
- if topo not in self.ind.keys():
- self.ind[topo] = []
- # apply indicator function on topo local mesh
- cond = npw.abs(self.chi[0](*topo.mesh.coords)) < self.epsilon
- self.ind[topo].append(cond)
- self._isempty(topo)
-
- # assert that the plane is a real surface, i.e.
- # only one value for coords[normalDir].
- # The expr is a bit tricky but it works ...
- ndir = np.where(self.normal != 0)[0][0]
- assert assertSubPlane(ndir, self.ind[topo][0], *topo.mesh.coords),\
- 'Your plane is not a surface but a volume.\
- Please reduce epsilon value.'
-
- return self.ind[topo]
-
-
-class SubSpace(HalfSpace):
- """
- Define a rectangular space in a plane normal to one
- coord. axis and the subspace below this surface.
- 'Below' = direction opposite to the outward normal of the plane
- (input param)
- """
- def __init__(self, lengths, **kwds):
- """
- @param domain : the physical domain that contains the space
- @param normal : outward normal
- @param point : coordinates of a point of the plane.
- @param lengths : lengths of the subplane
- @param epsilon : tolerance
- @param vd : velocity of the obstacle (considered as a rigid body),
- default = 0.
- """
- super(SubSpace, self).__init__(**kwds)
-
- def dist(cdir, val, *coords):
- return coords[cdir] - val
-
- self.dist = dist
- self.max = self.origin + npw.asrealarray(lengths)
- self.lengths = npw.asrealarray(lengths)
- ndir = np.where(self.normal != 0)[0][0]
- if self.normal[ndir] > 0:
- self.max[ndir] = self.origin[ndir]
- elif self.normal[ndir] < 0:
- self.max[ndir] = self.domain.max[ndir]
- # Only implemented for planes orthogonal to coord. axes
- assert len(self.normal[self.normal == 0]) == self.domain.dimension - 1
- self.coords = {}
-
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
- condMax = [0] * self.domain.dimension
- condMin = [0] * self.domain.dimension
- if topo not in self.ind.keys():
- self.ind[topo] = []
- # apply indicator function on topo local mesh
- coords = topo.mesh.coords
- cond = npw.asboolarray(self.chi[0](*coords) < self.epsilon)
- indices = np.where(self.normal == 0)[0]
-
- for i in indices:
- condMax[i] = self.dist(i, self.max[i], *coords) < self.epsilon
- condMin[i] = self.dist(i, self.origin[i], *coords) > - self.epsilon
- condMin[i] = np.logical_and(condMax[i], condMin[i])
- cond = npw.asarray(np.logical_and(cond, condMin[i]))
-
- self.ind[topo].append(cond)
- self._isempty(topo)
- return self.ind[topo]
-
-
-class SubPlane(SubSpace):
- """
- Define a rectangular surf in a plane normal to one
- coord. axis.
- """
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
- dim = self.domain.dimension
- condMax = [0] * dim
- condMin = [0] * dim
- if topo not in self.ind.keys():
- self.ind[topo] = []
- # apply indicator function on topo local mesh
- coords = topo.mesh.coords
- cond = npw.abs(self.chi[0](*coords)) < self.epsilon
- indices = np.where(self.normal == 0)[0]
- for i in indices:
- condMax[i] = self.dist(i, self.max[i], *coords) < self.epsilon
- condMin[i] = self.dist(i, self.origin[i], *coords) > -self.epsilon
- condMin[i] = np.logical_and(condMax[i], condMin[i])
- cond = npw.asboolarray(np.logical_and(cond, condMin[i]))
-
- self.ind[topo].append(cond)
- ilist = np.where(cond)
- if ilist[0].size == 0:
- self.slices[topo] = [slice(0, 0) for i in xrange(dim)]
- self.isEmpty[topo] = True
- else:
- self.isEmpty[topo] = False
- start = [ilist[i].min() for i in xrange(dim)]
- # Ghost points must not be included into surf. points
- ic = topo.mesh.iCompute
- start = [max(start[i], ic[i].start) for i in xrange(dim)]
- end = [ilist[i].max() for i in xrange(dim)]
- end = npw.asintarray([min(end[i], ic[i].stop - 1)
- for i in xrange(dim)])
- end += 1
- ndir = np.where(self.normal != 0)[0][0]
- end[ndir] = start[ndir] + 1
- self.slices[topo] = [slice(start[i], end[i])
- for i in xrange(dim)]
- assert assertSubPlane(ndir, cond, *topo.mesh.coords),\
- 'Your plane is not a surface but a volume.\
- Please reduce epsilon value.'
- subcoords = []
- # !! Warning : slices will be used for integration,
- # so the last point in each dir is not included.
- # Same thing for coords.
- for i in xrange(dim):
- subcoords.append(coords[i].flat[self.slices[topo][i]])
- subcoords = tuple(subcoords)
- self.coords[topo] = np.ix_(*subcoords)
- return self.ind[topo]
-
- def globalResolution(self, parent_topo):
- """
- Compute 'global resolution' of the subplane
- """
- # We create a false topology, with only one proc
- # to get the global resolution for the plane.
- # This could also be done with local computation
- # sum but that would need a lot of communications.
- if parent_topo.rank == 0:
- color = 0
- else:
- color = 1
- subcomm = parent_topo.comm.Split(color)
- dimension = self.domain.dimension
- tmp = None
- if parent_topo.rank == 0:
- resolution = parent_topo.globalMeshResolution
- ghosts = parent_topo.ghosts
- topo = self.domain.getOrCreateTopology(3, resolution,
- ghosts=ghosts, comm=subcomm)
- self.discretize(topo)
- sl = self.slices[topo]
- self.gRes = [sl[i].stop - sl[i].start for i in xrange(dimension)]
- self.gstart = [sl[i].start for i in xrange(dimension)]
- # if the topology has been created just to
- # get the global resolution, we can remove it
- if topo.isNew:
- self.domain.remove(topo)
- self.slices.pop(topo)
- self.ind.pop(topo)
- tmp = self.gRes + self.gstart
- tmp = parent_topo.comm.bcast(tmp)
- self.gRes = tmp[:dimension]
- self.gstart = tmp[dimension:]
- return self.gRes
-
-
-class PlaneBoundaries(Obstacle):
- """
- Defines top and down (meaning for min and max value in
- a given direction) planes at boundaries.
- All points in the spaces above the top plane and below the down plane
- will be included in the PlaneBoundaries list of indices.
- Thickness of the top/down areas is given as an input param.
- Example for z dir:
- \f$ \{x,y,z\} \ for \ z_{max} - \epsilon \leq z \leq z_{max} + \epsilon
- \ or \ z_{min} - \epsilon \leq z \leq z_{min}\f$
- """
-
- def __init__(self, normal_dir, thickness=0.1, **kwds):
- """
- Description of a sphere in a domain.
- @param domain : the physical domain that contains the sphere.
- @param thickness : thickness of boundary areas
- @param vd : velocity of obstacle (considered as a rigid body),
- default = 0.
- """
- super(PlaneBoundaries, self).__init__(**kwds)
- assert thickness > 0.0, 'Plate thickness must be positive'
- ## Thickness/2
- self.thickness = thickness
- ## Direction of the normal to the plate (0:x, 1:y, 2:z))
- normalUp = np.zeros((self.domain.dimension))
- normalUp[normal_dir] = -1
- pointUp = npw.zeros((self.domain.dimension))
- pointUp[normal_dir] = self.domain.max[normal_dir] - thickness
- self.upper = HalfSpace(domain=self.domain, normal=normalUp, point=pointUp,
- epsilon=1e-3)
- normalDown = np.zeros((self.domain.dimension))
- normalDown[normal_dir] = 1
- pointDown = npw.zeros((self.domain.dimension))
- pointDown[normal_dir] = self.domain.origin[normal_dir] + thickness
- self.lower = HalfSpace(domain=self.domain, normal=normalDown,
- point=pointDown, epsilon=1e-3)
-
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
-
- self.lower.discretize(topo)
- self.upper.discretize(topo)
- if topo not in self.ind.keys():
- # Warning FP : ind[topo] must be a list to be coherent
- # with sphere definition, where porous layers are allowed.
- # todo if required : add porous layers for planes.
- self.ind[topo] = []
- self.ind[topo].append(np.logical_or(self.upper.ind[topo][0],
- self.lower.ind[topo][0]))
- self._isempty(topo)
- return self.ind[topo]
-
-
-def assertSubPlane(ndir, ind, *coords):
- dim = len(coords)
- if dim == 2:
- return assertline(ndir, ind, *coords)
- elif dim == 3:
- return assertsurface(ndir, ind, *coords)
-
-
-def assertsurface(nd, ind, *coords):
-
- dim = len(coords)
- shape = np.zeros(dim, dtype=np.int32)
- shape[:] = [coords[i].shape[i] for i in xrange(dim)]
- cshape = coords[nd].shape
- if nd == 0:
- return max([a.max() - a.min()
- for a in [coords[nd][ind[:, i, j]]
- for i in xrange(shape[1])
- for j in xrange(shape[2])
- if coords[nd][ind[:, i, j]].size
- > 0]] + [0]) == 0.
- elif nd == 1:
- return max([a.max() - a.min()
- for a in [coords[nd][ind[i, :, j].reshape(cshape)]
- for i in xrange(shape[0])
- for j in xrange(shape[2])
- if coords[nd][ind[i, :, j].reshape(cshape)].size
- > 0]] + [0]) == 0.
-
- else:
- return max([a.max() - a.min()
- for a in [coords[nd][ind[i, j, :].reshape(cshape)]
- for i in xrange(shape[0])
- for j in xrange(shape[1])
- if coords[nd][ind[i, j, :].reshape(cshape)].size
- > 0]] + [0]) == 0.
-
-
-def assertline(nd, ind, *coords):
-
- dim = len(coords)
- shape = np.zeros(dim, dtype=np.int32)
- shape[:] = [coords[i].shape[i] for i in xrange(dim)]
- cshape = coords[nd].shape
- if nd == 0:
- return max([a.max() - a.min()
- for a in [coords[nd][ind[:, i]]
- for i in xrange(shape[1])
- if coords[nd][ind[:, i]].size
- > 0]] + [0]) == 0.
- elif nd == 1:
- return max([a.max() - a.min()
- for a in [coords[nd][ind[i, :].reshape(cshape)]
- for i in xrange(shape[0])
- if coords[nd][ind[i, :].reshape(cshape)].size
- > 0]] + [0]) == 0.
diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py
deleted file mode 100644
index 1bb0d62c7f6ad0b9d8cb19803f1629a31e9f8702..0000000000000000000000000000000000000000
--- a/HySoP/hysop/domain/obstacle/sphere.py
+++ /dev/null
@@ -1,132 +0,0 @@
-"""
-@file sphere.py
-Spherical or hemispherical sub-domain.
-"""
-from parmepy.domain.obstacle.obstacle import Obstacle
-import numpy as np
-import parmepy.tools.numpywrappers as npw
-
-
-class Sphere(Obstacle):
- """
- Spherical domain.
- """
-
- def __init__(self, position, radius=1.0, porousLayers=None, **kwds):
- """
- Description of a sphere in a domain.
- @param domain : the physical domain that contains the sphere.
- @param position : position of the center
- @param radius : sphere radius, default = 1
- @param porousLayers : a list of thicknesses
- for successive porous layers
- radius is the inside sphere radius and thicknesses are given from
- inside layer to outside one.
- @param vd : velocity of the sphere (considered as a rigid body),
- default = 0.
- """
- super(Sphere, self).__init__(**kwds)
-
- ## Radius of the sphere
- self.radius = radius
- ## Center position
- self.position = np.asarray(position)
-
- def dist(x, y, z, R):
- """
- """
- return npw.asarray(np.sqrt((x - self.position[0]) ** 2
- + (y - self.position[1]) ** 2
- + (z - self.position[2]) ** 2) - R)
-
- self.chi = [dist]
- ## List of thicknesses for porous layers
- if porousLayers is None:
- porousLayers = []
- self.layers = porousLayers
-
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
-
- if topo not in self.ind.keys():
- currentRadius = self.radius
- self.ind[topo] = []
- # First, internal sphere
- args = (currentRadius,)
- self.ind[topo].append(self.chi[0](*(topo.mesh.coords + args)) <= 0)
- # Then each layer from inside to outside
- # for each indicator function
- for thickness in self.layers:
- # apply indicator function on topo local mesh
- args = (currentRadius,)
- condA = self.chi[0](*(topo.mesh.coords + args)) > 0
- args = (currentRadius + thickness,)
- condB = self.chi[0](*(topo.mesh.coords + args)) <= 0
- self.ind[topo].append(np.logical_and(condA, condB))
- # update current radius
- currentRadius = currentRadius + thickness
- self._isempty(topo)
- return self.ind[topo]
-
- def __str__(self):
- """ToString method"""
- s = self.__class__.__name__ + ' of radius ' + str(self.radius)
- s += ' and center position ' + str(self.position)
- return s
-
-
-class HemiSphere(Sphere):
- """
- HemiSpherical domain.
- Area defined by the intersection of a sphere and the volume where
- x < xs for xs == x position of the center of the sphere.
- """
- def __init__(self, **kwds):
- """
- Description of a sphere in a domain.
- @param domain : the physical domain that contains the sphere.
- @param position : position of the center
- @param radius : sphere radius, default = 1
- @param porousLayers : a list of thicknesses
- for successive porous layers
- radius is the inside sphere radius and thicknesses are given from
- inside layer to outside one.
- @param vd : velocity of the sphere (considered as a rigid body),
- default = 0.
- """
- super(HemiSphere, self).__init__(**kwds)
-
- def LeftBox(x, y, z):
- return x - self.position[0]
- self.LeftBox = LeftBox
-
- def discretize(self, topo):
- # first check if we have already compute indices for
- # this topology
- if topo not in self.ind.keys():
- currentRadius = self.radius
- self.ind[topo] = []
- # check if we are in the left half-box
- cond0 = self.LeftBox(*(topo.mesh.coords)) <= 0
- # First, internal sphere
- args = (currentRadius,)
- condA = self.chi[0](*(topo.mesh.coords + args)) <= 0
- self.ind[topo].append(np.logical_and(condA, cond0))
- # Then each layer from inside to outside
- # for each indicator function
- for thickness in self.layers:
- # apply indicator function on topo local mesh
- args = (currentRadius,)
- condA = self.chi[0](*(topo.mesh.coords + args)) > 0
- args = (currentRadius + thickness,)
- condB = self.chi[0](*(topo.mesh.coords + args)) <= 0
- np.logical_and(condA, condB, condA)
- np.logical_and(condA, cond0, condA)
- condA = npw.asarray(condA)
- self.ind[topo].append(condA)
- # update current radius
- currentRadius = currentRadius + thickness
- self._isempty(topo)
-
- return self.ind[topo]
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 86e72c608ea5a4dedca1115d60d7f4a820b306b2..148dae2686d67df2f265964223174ca3c6724c89 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -15,7 +15,6 @@ name = '@PYPACKAGE_NAME@'
packages = ['parmepy',
'parmepy.domain',
'parmepy.domain.subsets',
- 'parmepy.domain.obstacle',
'parmepy.fields',
'parmepy.operator',
'parmepy.operator.discrete',