質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

89.72%

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

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 489

aine_

score 20

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

class Geometry(object):
    def __init__(
        obstacle_center_x=2.5, obstacle_center_y=0,
        obstacle_length=0.16, obstacle_height=0.161, obstacle_width=0.4,
        with_obstacle=True, dx=0.02, hdx=1.2, rho0=1000.0):

        # save the geometry details

        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


     self.hdx = hdx
        self.rho0 = rho0
        self.with_obstacle = with_obstacle


    def create_particles(self, **kwargs):

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

        obstacle_center_x = self.obstacle_center_x
        obstacle_center_y = self.obstacle_center_y

        nboundary_layers = self.nboundary_layers
        dx = self.dx

        # get the domain limits
        ghostlims = nboundary_layers * dx

        # create all particles
        eps = 0.1 * dx
        xx, yy, zz = numpy.mgrid[xmin:xmax+eps:dx,
                                 ymin:ymax+eps:dx,
                                 zmin:zmax+eps:dx]

        x = xx.ravel(); y = yy.ravel(); z = zz.ravel()

        # create a dummy particle array from which we'll sort
        pa = get_particle_array_wcsph(name='block', x=x, y=y, z=z)

        # get the individual arrays
        findices = []
        oindices = []

        obw2 = 0.5 * obstacle_width
        obl2 = 0.5 * obstacle_length
        obh = obstacle_height
        ocx = obstacle_center_x
        ocy = obstacle_center_y

        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)

        # extract the individual arrays
        if self.with_obstacle:
            oa = LongArray(len(oindices)); oa.set_data(numpy.array(oindices))
            obstacle = pa.extract_particles(oa)
            obstacle.set_name('obstacle')



        # create the particles
        if self.with_obstacle:
            particles = [fluid, boundary, obstacle]
        else:
            particles = [fluid, boundary]

        # set up particle properties
        h0 = self.hdx * dx

        volume = dx**3
        m0 = self.rho0 * volume

        for pa in particles:
            pa.m[:] = m0
            pa.h[:] = h0

            pa.rho[:] = self.rho0

        nf = fluid.num_real_particles
        nb = boundary.num_real_particles


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

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

0

        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 15:45

    ありがとうございます。if文で領域を指定しているわけですね。ありがとうございます。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 89.72%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

同じタグがついた質問を見る