長方形をつくっているコードを変形してT字形にしたい

長方形の箱をつくっているコードがあってそれを編集しT字形にしたいのですがやり方がわかりません。
元のコーデで長方形をつくっている部分はおそらく以下の部分で

python
1from numpy import concatenate, where, array
2
3(中略)           
4
5        indices = concatenate( (where( y <= -3.0 )[0],
6                                where( y >= 3.0 )[0],
7                                where( x >= 3.0 )[0],
8                                where( x <= 0 )[0],
9                                where( z <= 0 )[0]) )

ここを以下のように書き換えました

python
1        indices = concatenate( (where(y>=0.5 and x < 1.0)[0],
2                                where(y<=-0.5 and x < 1.0)[0],
3                                where( y <= -3.0 and x >= 1.0)[0],
4                                where( y >= 3.0 and x >= 1.0)[0],
5                                where( x >= 3.0 )[0],
6                                where( x <= 0 )[0],
7                                where( z <= 0 )[0]) )

しかしこのようにした場合ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()と表示されてしまいます。
また、

python
1        indices = concatenate( (where(y>=0.5 , x < 1.0)[0],
2                                where(y<=-0.5 , x < 1.0)[0],
3                                where( y <= -3.0 , x >= 1.0)[0],
4                                where( y >= 3.0 , x >= 1.0)[0],
5                                where( x >= 3.0 )[0],
6                                where( x <= 0 )[0],
7                                where( z <= 0 )[0]) )

とした場合はValueError: either both or neither of x and y should be given
となります。T字とするのは不可能なのでしょうか？
どのようにすればよいのでしょうか？
一応オープンソースであるため以下にコード全体も貼っておきます。
長方形の水槽に水を置いた場合の挙動をシミュレーションするモデルの設計となっています。

python
1import numpy
2from numpy import concatenate, where, array
3from pysph.base.utils import get_particle_array_wcsph, get_particle_array_iisph
4from cyarray.api import LongArray
5class DamBreak3DGeometry(object):
6    def __init__(
7        self, container_height=1.0, container_width=3.0, container_length=5,
8        fluid_column_height=0.8, fluid_column_width=3.0, fluid_column_length=1,
9        obstacle_center_x=2.5, obstacle_center_y=0,
10        obstacle_length=0.16, obstacle_height=0.161, obstacle_width=0.4,
11        nboundary_layers=5, with_obstacle=True, dx=0.02, hdx=1.2, rho0=1000.0):
12
13        # save the geometry details
14        self.container_width = container_width
15        self.container_length = container_length
16        self.container_height = container_height
17
18        self.fluid_column_length=fluid_column_length
19        self.fluid_column_width=fluid_column_width
20        self.fluid_column_height=fluid_column_height
21
22        self.obstacle_center_x = obstacle_center_x
23        self.obstacle_center_y = obstacle_center_y
24
25        self.obstacle_width=obstacle_width
26        self.obstacle_length=obstacle_length
27        self.obstacle_height=obstacle_height
28
29        self.nboundary_layers=nboundary_layers
30        self.dx=dx
31
32        self.hdx = hdx
33        self.rho0 = rho0
34        self.with_obstacle = with_obstacle
35
36    def get_max_speed(self, g=9.81):
37        return numpy.sqrt( 2 * g * self.fluid_column_height )
38
39    def create_particles(self, **kwargs):
40        fluid_column_height=self.fluid_column_height
41        fluid_column_width=self.fluid_column_width
42        fluid_column_length=self.fluid_column_length
43
44        container_height = self.container_height
45        container_length = self.container_length
46        container_width = self.container_width
47
48        obstacle_height = self.obstacle_height
49        obstacle_length = self.obstacle_length
50        obstacle_width = self.obstacle_width
51
52        obstacle_center_x = self.obstacle_center_x
53        obstacle_center_y = self.obstacle_center_y
54
55        nboundary_layers = self.nboundary_layers
56        dx = self.dx
57
58        # get the domain limits
59        ghostlims = nboundary_layers * dx
60
61        xmin, xmax = 0.0 -ghostlims, container_length + ghostlims
62        zmin, zmax = 0.0 - ghostlims, container_height + ghostlims
63
64        cw2 = 0.5 * container_width
65        ymin, ymax = -cw2 - ghostlims, cw2 + ghostlims
66
67        # create all particles
68        eps = 0.1 * dx
69        xx, yy, zz = numpy.mgrid[xmin:xmax+eps:dx,
70                                 ymin:ymax+eps:dx,
71                                 zmin:zmax+eps:dx]
72
73        x = xx.ravel(); y = yy.ravel(); z = zz.ravel()
74
75        # create a dummy particle array from which we'll sort
76        pa = get_particle_array_wcsph(name='block', x=x, y=y, z=z)
77
78        # get the individual arrays
79        indices = []
80        findices = []
81        oindices = []
82
83        obw2 = 0.5 * obstacle_width
84        obl2 = 0.5 * obstacle_length
85        obh = obstacle_height
86        ocx = obstacle_center_x
87        ocy = obstacle_center_y
88
89        for i in range(x.size):
90            xi = x[i]; yi = y[i]; zi = z[i]
91
92            # fluid
93            if ( (0 < xi <= fluid_column_length) and \
94                     (-cw2 < yi < cw2) and \
95                     (0 < zi <= fluid_column_height) ):
96
97                findices.append(i)
98
99            # obstacle
100            if ( (ocx-obl2 <= xi <= ocx+obl2) and \
101                     (ocy-obw2 <= yi <= ocy+obw2) and \
102                     (0 < zi <= obh) ):
103
104                oindices.append(i)
105
106        # extract the individual arrays
107        fa = LongArray(len(findices)); fa.set_data(numpy.array(findices))
108        fluid = pa.extract_particles(fa)
109        fluid.set_name('fluid')
110
111        if self.with_obstacle:
112            oa = LongArray(len(oindices)); oa.set_data(numpy.array(oindices))
113            obstacle = pa.extract_particles(oa)
114            obstacle.set_name('obstacle')
115
116        #上に乗せた部分はここです。cw2などは編集しました。
117        indices = concatenate( (where(y>=0.5 and x < 1.0)[0],
118                                where(y<=-0.5 and x < 1.0)[0],
119                                where( y <= -cw2 and x >= 1.0)[0],
120                                where( y >= cw2 and x >= 1.0)[0],
121                                where( x >= container_length )[0],
122                                where( x <= 0 )[0],
123                                where( z <= 0 )[0]) )
124
125        # remove duplicates
126        indices = array(list(set(indices)))
127
128        wa = LongArray(indices.size); wa.set_data(indices)
129        boundary = pa.extract_particles(wa)
130        boundary.set_name('boundary')
131
132        # create the particles
133        if self.with_obstacle:
134            particles = [fluid, boundary, obstacle]
135        else:
136            particles = [fluid, boundary]
137
138        # set up particle properties
139        h0 = self.hdx * dx
140
141        volume = dx**3
142        m0 = self.rho0 * volume
143
144        for pa in particles:
145            pa.m[:] = m0
146            pa.h[:] = h0
147
148            pa.rho[:] = self.rho0
149
150        nf = fluid.num_real_particles
151        nb = boundary.num_real_particles
152
153        if self.with_obstacle:
154            no = obstacle.num_real_particles
155            print("3D dam break with %d fluid, %d boundary, %d obstacle particles"%(nf, nb, no))
156        else:
157            print("3D dam break with %d fluid, %d boundary particles"%(nf, nb))
158
159
160        # load balancing props for the arrays
161        #fluid.set_lb_props(['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'h', 'm', 'gid',
162        #                    'x0', 'y0', 'z0', 'u0', 'v0', 'w0', 'rho0'])
163        fluid.set_lb_props( list(fluid.properties.keys()) )
164
165        #boundary.set_lb_props(['x', 'y', 'z', 'rho', 'h', 'm', 'gid', 'rho0'])
166        #obstacle.set_lb_props(['x', 'y', 'z', 'rho', 'h', 'm', 'gid', 'rho0'])
167        boundary.set_lb_props( list(boundary.properties.keys()) )
168
169        # boundary and obstacle particles can do with a reduced list of properties
170        # to be saved to disk since they are fixed
171        boundary.set_output_arrays( ['x', 'y', 'z', 'rho', 'm', 'h', 'p', 'tag', 'pid', 'gid'] )
172
173        if self.with_obstacle:
174            obstacle.set_lb_props( list(obstacle.properties.keys()) )
175            obstacle.set_output_arrays( ['x', 'y', 'z', 'rho', 'm', 'h', 'p', 'tag', 'pid', 'gid'] )
176
177        return particles

行動規範の内容に同意します

回答1件

ベストアンサー

処理詳細は確認していません。numpyの条件演算についてのみの回答です。
y>=0.5かつx < 1.0である条件で抽出したい場合は
where(y>=0.5 and x < 1.0)→where((y>=0.5) & (x < 1.0))と記述してください。
以下、検証コードです。

Python
1import numpy
2from numpy import concatenate, where, array
3
4a = numpy.array([10,11,12,13,14])
5b = numpy.array([15,16,17,18,19])
6
7#print( where( (a>10) and (b<=18) )) # ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
8print( where( (a>10)  &  (b<=18) )) # (array([1, 2, 3], dtype=int64),)