長方形（四角柱）オブジェをつくるコードに手を加えて円柱をつくりたい

元のコードは以下の通りです。

python
1class Geometry(object):
2    def __init__(
3        obstacle_center_x=2.5, obstacle_center_y=0,
4        obstacle_length=0.16, obstacle_height=0.161, obstacle_width=0.4,
5        with_obstacle=True, dx=0.02, hdx=1.2, rho0=1000.0):
6
7        # save the geometry details
8        
9        self.obstacle_center_x = obstacle_center_x
10        self.obstacle_center_y = obstacle_center_y
11
12        self.obstacle_width=obstacle_width
13        self.obstacle_length=obstacle_length
14        self.obstacle_height=obstacle_height
15
16
17　　　　 self.hdx = hdx
18        self.rho0 = rho0
19        self.with_obstacle = with_obstacle
20
21
22    def create_particles(self, **kwargs):
23
24        obstacle_height = self.obstacle_height
25        obstacle_length = self.obstacle_length
26        obstacle_width = self.obstacle_width
27
28        obstacle_center_x = self.obstacle_center_x
29        obstacle_center_y = self.obstacle_center_y
30
31        nboundary_layers = self.nboundary_layers
32        dx = self.dx
33
34        # get the domain limits
35        ghostlims = nboundary_layers * dx
36
37        # create all particles
38        eps = 0.1 * dx
39        xx, yy, zz = numpy.mgrid[xmin:xmax+eps:dx,
40                                 ymin:ymax+eps:dx,
41                                 zmin:zmax+eps:dx]
42
43        x = xx.ravel(); y = yy.ravel(); z = zz.ravel()
44
45        # create a dummy particle array from which we'll sort
46        pa = get_particle_array_wcsph(name='block', x=x, y=y, z=z)
47
48        # get the individual arrays
49        findices = []
50        oindices = []
51
52        obw2 = 0.5 * obstacle_width
53        obl2 = 0.5 * obstacle_length
54        obh = obstacle_height
55        ocx = obstacle_center_x
56        ocy = obstacle_center_y
57
58        for i in range(x.size):
59            xi = x[i]; yi = y[i]; zi = z[i]
60
61            if ( (ocx-obl2 <= xi <= ocx+obl2) and \
62                     (ocy-obw2 <= yi <= ocy+obw2) and \
63                     (0 < zi <= obh) ):
64
65                oindices.append(i)
66
67        # extract the individual arrays
68        if self.with_obstacle:
69            oa = LongArray(len(oindices)); oa.set_data(numpy.array(oindices))
70            obstacle = pa.extract_particles(oa)
71            obstacle.set_name('obstacle')
72   
73
74
75        # create the particles
76        if self.with_obstacle:
77            particles = [fluid, boundary, obstacle]
78        else:
79            particles = [fluid, boundary]
80
81        # set up particle properties
82        h0 = self.hdx * dx
83
84        volume = dx**3
85        m0 = self.rho0 * volume
86
87        for pa in particles:
88            pa.m[:] = m0
89            pa.h[:] = h0
90
91            pa.rho[:] = self.rho0
92
93        nf = fluid.num_real_particles
94        nb = boundary.num_real_particles
95

本来はもっと長いコードなのですが、オブジェ関係の部分のみ抜粋しています。（よって不可解な部分があるかもしれません）
オブジェクト指向で書かれているためこのコード内でオブジェの作成が完結していることはなんとなく推測していますが、どの部分で長方形（四角柱）を設計しているのかわかりません。
あまりにも未熟な質問であることは重々理解していますが設計箇所だけでもご教授ねがえないでしょうか？
ちなみにこれはSPH法とよばれる方法で、指定した空間にことにより任意の圧力をもつ粒子をしきつめることにより流体・個体などのシミュレーションをするというものです。

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

回答1件

ベストアンサー

        self.obstacle_center_x = obstacle_center_x
        self.obstacle_center_y = obstacle_center_y

        self.obstacle_width=obstacle_width
        self.obstacle_length=obstacle_length
        self.obstacle_height=obstacle_height

ここで四角柱の大きさと位置を定義して

        for i in range(x.size):
            xi = x[i]; yi = y[i]; zi = z[i]

            if ( (ocx-obl2 <= xi <= ocx+obl2) and \
                     (ocy-obw2 <= yi <= ocy+obw2) and \
                     (0 < zi <= obh) ):

                oindices.append(i)

ここで領域内のパーティクルだけ受け入れているように見えるので

前者のwidth length の代わりに半径を定義して

後者の側で、xi,yi,ziが定義された円柱の中に収まっているかどうか
という判定を行えば円柱内のパーティクルだけ受け入れられるようになりそうですね。

点が円の中にあるかどうかは

(ocx - xi)**2 + (ocy - yi)**2 <= 半径**2

で判定できますよね。

投稿2018/11/05 03:52