From 0a8c01374c4bd18f1f685beece4923d6c3b898af Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Wed, 22 May 2013 13:14:39 +0000
Subject: [PATCH] Update obstacles + 'plates' in 2D domain and update
penalization example
---
Examples/testPenalization.py | 15 ++-
HySoP/hysop/domain/obstacle/cylinder2d.py | 118 +++++++++---------
HySoP/hysop/domain/obstacle/plates.py | 17 ++-
HySoP/hysop/domain/obstacle/sphere.py | 33 ++---
HySoP/hysop/operator/discrete/penalization.py | 2 +-
5 files changed, 99 insertions(+), 86 deletions(-)
diff --git a/Examples/testPenalization.py b/Examples/testPenalization.py
index 41a63ac12..7aebc8612 100644
--- a/Examples/testPenalization.py
+++ b/Examples/testPenalization.py
@@ -4,7 +4,7 @@ from parmepy.fields.continuous import Field
from parmepy.operator.analytic import Analytic
from parmepy.domain.obstacle.sphere import Sphere, HemiSphere
from parmepy.domain.obstacle.cylinder2d import Cylinder2D, SemiCylinder2D
-from parmepy.domain.obstacle.plates import ZPlates
+from parmepy.domain.obstacle.plates import Plates
pi = math.pi
from parmepy.operator.penalization import Penalization
from parmepy.operator.monitors.printer import Printer
@@ -49,24 +49,27 @@ sphere = Sphere(dom, position=[0., 0., 0.], radius=0.5, porousLayers=[0.13])
hsphere = HemiSphere(dom, position=[0., 0., 0.],
radius=0.5, porousLayers=[0.13])
-zplates = ZPlates(dom, epsilon=0.1)
+plates = Plates(dom, normal_dir=0, epsilon=0.1)
+plates2 = Plates(dom2, normal_dir=1, 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 = Penalization(scal, [hsphere, plates], coeff=[1e6, 10],
+ resolutions={scal: resol3D})
+penal2 = Penalization(scal2, [hcyl, plates2], coeff=[1e6, 10],
+ resolutions={scal2: resol2D})
penal.setUp()
-
+penal2.setUp()
printer = Printer(fields=[scal2], frequency=1)
printer.setUp()
print scal.norm()
penal.apply(dt=0.1)
-
+penal2.apply(dt=0.1)
printer.apply(ite=1)
print "print ..."
print scal.norm()
diff --git a/HySoP/hysop/domain/obstacle/cylinder2d.py b/HySoP/hysop/domain/obstacle/cylinder2d.py
index a9a769462..5d7f8b978 100644
--- a/HySoP/hysop/domain/obstacle/cylinder2d.py
+++ b/HySoP/hysop/domain/obstacle/cylinder2d.py
@@ -2,12 +2,12 @@
@file cylinder2d.py
Rigid disk.
"""
-from parmepy.domain.obstacle.obstacle import Obstacle
+from parmepy.domain.obstacle.sphere import Sphere, HemiSphere
import math
import numpy as np
-class Cylinder2D(Obstacle):
+class Cylinder2D(Sphere):
"""
Disk in a 2D domain.
"""
@@ -25,14 +25,9 @@ class Cylinder2D(Obstacle):
@param vd : velocity of the disk (considered as a rigid body),
default = 0.
"""
- Obstacle.__init__(self, domain, vd=vd)
+ Sphere.__init__(self, domain, position, radius, vd, porousLayers)
assert self.domain.dimension == 2
- ## Radius of the disk
- self.radius = radius
- ## Center position
- self.position = np.asarray(position)
-
def dist(x, y, R):
"""
"""
@@ -40,28 +35,27 @@ class Cylinder2D(Obstacle):
+ (y - self.position[1]) ** 2) - R
self.chi = [np.vectorize(dist)]
## List of thicknesses for porous layers
- 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] = []
- # 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 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"""
@@ -70,7 +64,7 @@ class Cylinder2D(Obstacle):
return s
-class SemiCylinder2D(Cylinder2D):
+class SemiCylinder2D(HemiSphere):
"""
Half disk in a 2D domain.
"""
@@ -84,38 +78,46 @@ class SemiCylinder2D(Cylinder2D):
@param vd : velocity of the disk (considered as a rigid body),
default = 0.
"""
- Cylinder2D.__init__(self, domain, position, radius, vd, porousLayers)
+ HemiSphere.__init__(self, domain, position, radius, vd, porousLayers)
+ assert self.domain.dimension == 2
+
+ def dist(x, y, R):
+ """
+ """
+ return math.sqrt((x - self.position[0]) ** 2
+ + (y - self.position[1]) ** 2) - R
+ self.chi = [np.vectorize(dist)]
- def HalfPlane(x, y):
+ def LeftBox(x, y):
return x - self.position[0]
- self.HalfPlane = np.vectorize(HalfPlane)
+ 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] = []
- # 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 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)
diff --git a/HySoP/hysop/domain/obstacle/plates.py b/HySoP/hysop/domain/obstacle/plates.py
index fe7f54747..79632643e 100644
--- a/HySoP/hysop/domain/obstacle/plates.py
+++ b/HySoP/hysop/domain/obstacle/plates.py
@@ -6,7 +6,7 @@ from parmepy.domain.obstacle.obstacle import Obstacle
import numpy as np
-class ZPlates(Obstacle):
+class Plates(Obstacle):
"""
Top and down (z-axis) plate-like sub-domains.
Defines two areas on the min and max z position of
@@ -16,7 +16,7 @@ class ZPlates(Obstacle):
\ or \ z_{min} - \epsilon \leq z \leq z_{min}\f$
"""
- def __init__(self, domain, epsilon=0.1, vd=0.0):
+ def __init__(self, domain, normal_dir, epsilon=0.1, vd=0.0):
"""
Description of a sphere in a domain.
@param domain : the physical domain that contains the sphere.
@@ -28,12 +28,16 @@ class ZPlates(Obstacle):
assert epsilon > 0.0, 'Plate thickness must be positive'
## Thickness/2
self.epsilon = epsilon
+ ## Direction of the normal to the plate (0:x, 1:y, 2:z))
+ self.normal = normal_dir
+ assert self.normal < domain.dimension and self.normal >= 0
- def Upper(x, y, z):
- return domain.max[-1] - self.epsilon - z
+ def Upper(*args):
+ return domain.max[self.normal] - self.epsilon - args[self.normal]
- def Lower(x, y, z):
- return z - self.epsilon - domain.origin[-1]
+ def Lower(*args):
+ return args[self.normal] - self.epsilon - \
+ domain.origin[self.normal]
self.Upper = np.vectorize(Upper)
self.Lower = np.vectorize(Lower)
@@ -56,3 +60,4 @@ class ZPlates(Obstacle):
s = 'Upper and lower (z-axis) plane areas of thickness '
s += str(2 * self.epsilon) + '.'
return s
+
diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py
index b02063b34..99fa9cb31 100644
--- a/HySoP/hysop/domain/obstacle/sphere.py
+++ b/HySoP/hysop/domain/obstacle/sphere.py
@@ -50,19 +50,20 @@ class Sphere(Obstacle):
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
+ 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)
+ currentRadius = currentRadius + thickness
return self.ind[topo]
@@ -107,22 +108,24 @@ class HemiSphere(Sphere):
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
+ 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))
+ currentRadius = currentRadius + thickness
+
return self.ind[topo]
def __str__(self):
diff --git a/HySoP/hysop/operator/discrete/penalization.py b/HySoP/hysop/operator/discrete/penalization.py
index a83a649cc..eeca52502 100644
--- a/HySoP/hysop/operator/discrete/penalization.py
+++ b/HySoP/hysop/operator/discrete/penalization.py
@@ -46,7 +46,7 @@ class Penalization_d(DiscreteOperator):
else:
lsize = min(len(self.cond), len(cond))
for c in xrange(lsize):
- self.cond[c] = cond[c] == self.cond[c]
+ np.logical_or(cond[c], self.cond[c], self.cond[c])
for c in xrange(lsize + 1, len(cond)):
self.cond[c] = cond[c]
assert self.factor.size == len(self.cond)
--
GitLab