diff --git a/Examples/testPenalization.py b/Examples/testPenalization.py index 8623e7a175e00397a6a8e81d058e51fb86fa37ee..41a63ac1291c4ab9237545b75e29d1416ea391cc 100644 --- a/Examples/testPenalization.py +++ b/Examples/testPenalization.py @@ -1,9 +1,10 @@ -import numpy as np import parmepy as pp import math from parmepy.fields.continuous import Field from parmepy.operator.analytic import Analytic -from parmepy.domain.obstacle.sphere import Sphere +from parmepy.domain.obstacle.sphere import Sphere, HemiSphere +from parmepy.domain.obstacle.cylinder2d import Cylinder2D, SemiCylinder2D +from parmepy.domain.obstacle.plates import ZPlates pi = math.pi from parmepy.operator.penalization import Penalization from parmepy.operator.monitors.printer import Printer @@ -12,18 +13,26 @@ from parmepy.operator.monitors.printer import Printer def vitesse(x, y, z): return x, y, z -nb = 65 + +def vitesse2(x, y): + return x, y + +nb = 33 Lx = Ly = Lz = 2 dom = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[-1., -1., -1.]) +dom2 = pp.Box(dimension=2, length=[Lx, Ly], origin=[-1., -1.]) resol3D = [nb, nb, nb] +resol2D = [nb, nb] # Fields scal = Field(domain=dom, name='Scalar') op = Analytic(scal, formula=vitesse, resolutions={scal: resol3D}) +scal2 = Field(domain=dom2, name='Scalar') +op2 = Analytic(scal2, formula=vitesse2, resolutions={scal2: resol2D}) op.setUp() - +op2.setUp() topo = dom.topologies[0] coords = topo.mesh.coords @@ -31,20 +40,33 @@ ff = scal.discreteFields[topo].data ff[...] = 128 +topo2 = dom2.topologies.values()[0] + +ff2 = scal2.discreteFields[topo2].data +ff2[...] = 128 sphere = Sphere(dom, position=[0., 0., 0.], radius=0.5, porousLayers=[0.13]) -ind = sphere.discretize(topo) -penal = Penalization(scal, sphere, coeff=[1e6, 10], - resolutions={scal: resol3D}) +hsphere = HemiSphere(dom, position=[0., 0., 0.], + radius=0.5, porousLayers=[0.13]) + +zplates = ZPlates(dom, epsilon=0.1) + +cyl = Cylinder2D(dom2, position=[0., 0.], radius=0.5, porousLayers=[0.13]) +hcyl = SemiCylinder2D(dom2, position=[0., 0.], radius=0.5, porousLayers=[0.13]) + +ind = hsphere.discretize(topo) + +penal = Penalization(scal2, hcyl, coeff=[1e6, 10.], + resolutions={scal: resol2D}) penal.setUp() -printer = Printer(fields=[scal], frequency=1) +printer = Printer(fields=[scal2], frequency=1) printer.setUp() print scal.norm() -penal.apply(0.1) +penal.apply(dt=0.1) -printer.apply(1) +printer.apply(ite=1) print "print ..." print scal.norm() diff --git a/HySoP/hysop/domain/obstacle/__init__.py b/HySoP/hysop/domain/obstacle/__init__.py index 69cfe680cd147965d7e314b3a40f1e8d1223388d..be5911bbcc91dc597c5c18e7479042e779dfd672 100644 --- a/HySoP/hysop/domain/obstacle/__init__.py +++ b/HySoP/hysop/domain/obstacle/__init__.py @@ -1,2 +1,57 @@ ## @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/cylinder2d.py b/HySoP/hysop/domain/obstacle/cylinder2d.py index a3b261afab691c3f0b7459067b5e461bb7dc380a..a9a76946230ee59eb5050d1b9a5eb84490955c66 100644 --- a/HySoP/hysop/domain/obstacle/cylinder2d.py +++ b/HySoP/hysop/domain/obstacle/cylinder2d.py @@ -7,34 +7,61 @@ import math import numpy as np -class Cylinder2d(Obstacle): +class Cylinder2D(Obstacle): """ - 2D disk. + Disk in a 2D domain. """ - def __init__(self, domain, position, radius=1.0, vd=0.0): + def __init__(self, domain, position, radius=1.0, vd=0.0, porousLayers=[]): """ 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 - (if box ...) + @param porousLayers : a list of thicknesses + for successive porous layers + radius is the outside sphere radius and thicknesses are given from + outside layer to inside one. @param vd : velocity of the disk (considered as a rigid body), default = 0. """ - ## Radius of the sphere + Obstacle.__init__(self, domain, vd=vd) + assert self.domain.dimension == 2 + + ## Radius of the disk self.radius = radius ## Center position self.position = np.asarray(position) - def dist(x, y): + def dist(x, y, R): """ """ return math.sqrt((x - self.position[0]) ** 2 - + (y - self.position[1]) ** 2) - self.radius + + (y - self.position[1]) ** 2) - R + self.chi = [np.vectorize(dist)] + ## List of thicknesses for porous layers + self.layers = porousLayers - Obstacle.__init__(self, domain, formula=dist, vd=vd) - assert self.domain.dimension == 2 + 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] = [] + # 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 + # and finally the 'internal' sphere + args = (currentRadius,) + self.ind[topo].append(self.chi[0](*(topo.mesh.coords + args)) <= 0) + return self.ind[topo] def __str__(self): """ToString method""" @@ -43,6 +70,59 @@ class Cylinder2d(Obstacle): return s +class SemiCylinder2D(Cylinder2D): + """ + Half disk in a 2D domain. + """ + def __init__(self, domain, position, radius=1.0, vd=0.0, porousLayers=[]): + """ + 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. + """ + Cylinder2D.__init__(self, domain, position, radius, vd, porousLayers) + + def HalfPlane(x, y): + return x - self.position[0] + + self.HalfPlane = np.vectorize(HalfPlane) + + 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 cond : we must be in the left half plane + cond0 = self.HalfPlane(*(topo.mesh.coords)) <= 0 + # 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) + self.ind[topo].append(condA) + # update current radius + currentRadius = currentRadius - thickness + # and finally the 'internal' sphere + args = (currentRadius,) + condA = self.chi[0](*(topo.mesh.coords + args)) <= 0 + self.ind[topo].append(np.logical_and(condA, cond0)) + return self.ind[topo] + + def __str__(self): + s = '2D semi-cylinder of radius ' + str(self.radius) + s += ' and center position ' + str(self.position) + return s + if __name__ == "__main__": print "This module defines the following classe:" - print "Sphere: ", Cylinder2d.__doc__ + print "Disk: ", Cylinder2D.__doc__ + print "Half-disk: ", SemiCylinder2D.__doc__ diff --git a/HySoP/hysop/domain/obstacle/hemisphere.py b/HySoP/hysop/domain/obstacle/hemisphere.py index 1c6e4aa3111c104c266d37732ee6a71e8013b70e..7bed176c22ba8883663a53a9ba431496447b2b4e 100644 --- a/HySoP/hysop/domain/obstacle/hemisphere.py +++ b/HySoP/hysop/domain/obstacle/hemisphere.py @@ -1,11 +1,13 @@ """ +@file sphere.py +Spherical sub-domain Hemisphere obstacle description """ from parmepy.constants import np -from parmepy.domain.obstacle.obstacle_3D import Obstacle3D +from parmepy.domain.obstacle.sphere import Sphere -class Hemisphere(Obstacle3D): +class Hemisphere(Sphere): """ Discrete obstacles representation. """ diff --git a/HySoP/hysop/domain/obstacle/obstacle.py b/HySoP/hysop/domain/obstacle/obstacle.py index 36077c5c4f7770dec9509fbc4e9e50370f717c2a..e764524d5bbda56b0abe4cafa6b97ae47c7f077a 100644 --- a/HySoP/hysop/domain/obstacle/obstacle.py +++ b/HySoP/hysop/domain/obstacle/obstacle.py @@ -21,6 +21,9 @@ class Obstacle(object): """ ## Domain. self.domain = domain + from parmepy.domain.box import Box + assert isinstance(domain, Box),\ + 'Obstacle only implemented for box-like domains' ## Obstacle dimension. self.dimension = domain.dimension ## A function that describe the geometry of the obstacle. diff --git a/HySoP/hysop/domain/obstacle/plates.py b/HySoP/hysop/domain/obstacle/plates.py new file mode 100644 index 0000000000000000000000000000000000000000..fe7f5474716ec745c4bb7ba8eaa2109f89d54163 --- /dev/null +++ b/HySoP/hysop/domain/obstacle/plates.py @@ -0,0 +1,58 @@ +""" +@file sphere.py +Plate-like sub-domains. +""" +from parmepy.domain.obstacle.obstacle import Obstacle +import numpy as np + + +class ZPlates(Obstacle): + """ + Top and down (z-axis) plate-like sub-domains. + Defines two areas on the min and max z position of + the main domain. + zplates: + \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, domain, epsilon=0.1, vd=0.0): + """ + Description of a sphere in a domain. + @param domain : the physical domain that contains the sphere. + @param epsilon : thickness/2 of the plates + @param vd : velocity of obstacle (considered as a rigid body), + default = 0. + """ + Obstacle.__init__(self, domain, vd=vd) + assert epsilon > 0.0, 'Plate thickness must be positive' + ## Thickness/2 + self.epsilon = epsilon + + def Upper(x, y, z): + return domain.max[-1] - self.epsilon - z + + def Lower(x, y, z): + return z - self.epsilon - domain.origin[-1] + + self.Upper = np.vectorize(Upper) + self.Lower = np.vectorize(Lower) + + 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 + condUpper = self.Upper(*topo.mesh.coords) <= 0 + condLower = self.Lower(*topo.mesh.coords) <= 0 + self.ind[topo].append(np.logical_or(condUpper, condLower)) + + return self.ind[topo] + + def __str__(self): + """ToString method""" + s = 'Upper and lower (z-axis) plane areas of thickness ' + s += str(2 * self.epsilon) + '.' + return s diff --git a/HySoP/hysop/domain/obstacle/semi_cylinder.py b/HySoP/hysop/domain/obstacle/semi_cylinder.py deleted file mode 100644 index c5bca741410696591b40898859874beba41b8af7..0000000000000000000000000000000000000000 --- a/HySoP/hysop/domain/obstacle/semi_cylinder.py +++ /dev/null @@ -1,110 +0,0 @@ -""" -Semi-cylinder obstacle description -""" -from parmepy.constants import np, PARMES_INTEGER -from parmepy.domain.obstacle.obstacle_2D import Obstacle2D - - -class SemiCylinder(Obstacle2D): - """ - Discrete obstacles representation. - """ - - def __init__(self, obst, radius=0.2): - """ - Creates the semi-cylinder obsctacle and y-boundaries and - returns 3 arrays corresponding to the characteristic - functions of three different porous media. - - @param obstacle2D : Two dimensional obstacle description. - @param radius : Semi-cylinder radius - """ - ## 2D parent obstacle - self.obst = obst - ## Radius of the semi-cylinder - self.radius = radius - ## Characteristic function (array) for y-boundaries - self.chiBoundary = None - ## Characteristic function (array) for the solid - self.chiSolid = None - ## Characteristic function (array) for the porous area - self.chiPorous = None - print 'obstacle =', self.obst.obstacleName - - def setUp(self, topology): - """ - Compute the characteristic functions associated - to y-boundaries and semi-cylinder - """ - ## Temporary arrays - chiBoundary_i = [] - chiBoundary_j = [] - chiSolid_i = [] - chiSolid_j = [] - chiPorous_i = [] - chiPorous_j = [] - - step = topology.mesh.space_step - ghosts = topology.ghosts - local_start = topology.mesh.local_start - local_end = topology.mesh.local_end - coord_start = topology.mesh.origin + (ghosts * step) - coord_end = topology.mesh.end - (ghosts * step) - layerMin = coord_start[1] + self.obst.zlayer - layerMax = coord_end[1] - self.obst.zlayer - if not (self.obst.porousLayerThickn <= self.radius): - raise ValueError("Error, porous layer thickness" + - "is higher than semi-cylinder radius.") - radiusMinuslayer = self.radius- self.obst.porousLayerThickn - - print 'step, ghosts, local_start, local_end, zlayer', \ - step, ghosts, local_start, local_end, self.obst.zlayer - print 'start, end, layerMin, layerMax, radiusMinuslayer', \ - coord_start, coord_end, layerMin, layerMax, radiusMinuslayer - - for j in xrange (local_start[1], local_end[1] + 1): - cy = coord_start[1] + j * step[1] - for i in xrange (local_start[0], local_end[0] + 1): - if (cy >= layerMax or cy <= layerMin): - # we are in the y-layer boundary: - chiBoundary_i.append(i) - chiBoundary_j.append(j) - else : - cx = coord_start[0] + i * step[0] - dist = np.sqrt((cx - self.obst.center[0]) ** 2 + - (cy - self.obst.center[1]) ** 2) - if (radiusMinuslayer < dist - and dist <= self.radius + 1E-12 - and cx <= self.obst.center[0] - and self.obst.porousLayerThickn != 0.): - # we are in the porous region of the semi-cylinder: - chiPorous_i.append(i) - chiPorous_j.append(j) - if (dist <= radiusMinuslayer + 1E-12 - and cx <= self.obst.center[0]): - # we are in the solid region of the semi-cylinder: - chiSolid_i.append(i) - chiSolid_j.append(j) - - ## Characteristic function of penalized boundaries - chiBoundary_i = np.asarray(chiBoundary_i, dtype=PARMES_INTEGER) - chiBoundary_j = np.asarray(chiBoundary_j, dtype=PARMES_INTEGER) - self.chiBoundary = tuple([chiBoundary_i, chiBoundary_j]) - - ## Characteristic function of solid areas - chiSolid_i = np.asarray(chiSolid_i, dtype=PARMES_INTEGER) - chiSolid_j = np.asarray(chiSolid_j, dtype=PARMES_INTEGER) - self.chiSolid = tuple([chiSolid_i, chiSolid_j]) - - ## Characteristic function of porous areas - chiPorous_i = np.asarray(chiPorous_i, dtype=PARMES_INTEGER) - chiPorous_j = np.asarray(chiPorous_j, dtype=PARMES_INTEGER) - self.chiPorous = tuple([chiPorous_i, chiPorous_j]) - - def __str__(self): - """ToString method""" - return "SemiCylinder" - -if __name__ == "__main__" : - print "This module defines the following classe:" - print "SemiCylinder: ", SemiCylinder.__doc__ diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py index 90e7fa186295ad96a58df2f4d116d975aff6a06e..b02063b3455b51d8a1c44a8498eb662bf141549c 100644 --- a/HySoP/hysop/domain/obstacle/sphere.py +++ b/HySoP/hysop/domain/obstacle/sphere.py @@ -1,6 +1,6 @@ """ @file sphere.py -Spherical sub-domain +Spherical or hemispherical sub-domain. """ from parmepy.domain.obstacle.obstacle import Obstacle import math @@ -40,7 +40,6 @@ class Sphere(Obstacle): + (z - self.position[2]) ** 2) - R self.chi = [np.vectorize(dist)] - self.ind = {} ## List of thicknesses for porous layers self.layers = porousLayers @@ -58,7 +57,7 @@ class Sphere(Obstacle): 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(condA == condB) + self.ind[topo].append(np.logical_and(condA, condB)) # update current radius currentRadius = currentRadius - thickness # and finally the 'internal' sphere @@ -74,6 +73,64 @@ class Sphere(Obstacle): 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, domain, position, radius=1.0, + vd=0.0, porousLayers=[]): + """ + 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 outside sphere radius and thicknesses are given from + outside layer to inside one. + @param vd : velocity of the sphere (considered as a rigid body), + default = 0. + """ + Sphere.__init__(self, domain, position, radius, vd, porousLayers) + + def LeftBox(x, y, z): + return x - self.position[0] + self.LeftBox = np.vectorize(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 + # 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) + self.ind[topo].append(condA) + # update current radius + currentRadius = currentRadius - thickness + # and finally the 'internal' sphere + args = (currentRadius,) + condA = self.chi[0](*(topo.mesh.coords + args)) <= 0 + self.ind[topo].append(np.logical_and(condA, cond0)) + return self.ind[topo] + + def __str__(self): + s = 'hemisphere of radius ' + str(self.radius) + s += ' and center position ' + str(self.position) + return s + if __name__ == "__main__": - print "This module defines the following classe:" + print "This module defines the following classes:" print "Sphere: ", Sphere.__doc__ + print "HemiSphere: ", HemiSphere.__doc__ diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py index 9d0a682ad3abdde05420c6035789f92035cb8e50..dd2ca52c2f7173b15ccadf018ed46bd8ae2647e8 100644 --- a/HySoP/hysop/operator/continuous.py +++ b/HySoP/hysop/operator/continuous.py @@ -85,7 +85,7 @@ class Operator(object): self.discreteOperator.finalize() @debug - def apply(self, t, dt, ite): + def apply(self, t=None, dt=None, ite=None): """ apply the operator on its variables. @param t : Current simulation time. diff --git a/HySoP/hysop/operator/discrete/penalization.py b/HySoP/hysop/operator/discrete/penalization.py index c87ce5a11aba46f9a36bb94d7d1db9f4253ad7e0..a83a649ccef8f74ce07ca2d84649939fe3fa7ede 100644 --- a/HySoP/hysop/operator/discrete/penalization.py +++ b/HySoP/hysop/operator/discrete/penalization.py @@ -56,7 +56,7 @@ class Penalization_d(DiscreteOperator): pass @debug - def apply(self, t, dt, ite): + def apply(self, t=None, dt=0.0, ite=None): self.compute_time = MPI.Wtime() diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py index 37e207ac70dab19642830233b932f37397303252..524a908ce9dc3246c9dc6b601f1385902003b1aa 100644 --- a/HySoP/hysop/operator/monitors/printer.py +++ b/HySoP/hysop/operator/monitors/printer.py @@ -6,6 +6,7 @@ Classes for handling ouputs. from parmepy.constants import np, S_DIR, PARMES_REAL from parmepy.operator.monitors.monitoring import Monitoring import evtk.hl as evtk +from parmepy.mpi import main_rank import time @@ -27,7 +28,7 @@ class Printer(Monitoring): """ Monitoring.__init__(self, fields, frequency) ## Filename prefix - self.outputPrefix = outputPrefix + self.outputPrefix = outputPrefix + str(main_rank) ## Method to collect data in case of distributed data self.get_data_method = None if self.freq != 0: @@ -39,7 +40,7 @@ class Printer(Monitoring): self.compute_time = 0. self.call_number = 0 - def apply(self, t, dt, ite): + def apply(self, t=None, dt=None, ite=None): self.step(ite) def _build_vtk_dict(self):