LoginSignup
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@7of9

geometry (3D) > 正四角錐の内部を格子状のdipoleで埋め尽くす > volfil_tetrahedron_180519.py > v0.1 | 結果の可視化 (Tinkercad: Shape Generator) | v0.2: Numpy用ファイル出力

Last updated at Posted at 2018-05-19
動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04.4 LTS desktop amd64
TensorFlow v1.7.0
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6
gnustep-gui-runtime v0.24.0-3.1
PyMieScatt v1.7.0

check_intersection_3D_180512.py v0.1: geometry (3D) > rayとtriangleの交差判定 > check_intersection_3D_180512.py > v0.1 > C: intersect3D_RayTriangle()からの移植, v0.2: 戻り値の変更

前回: geometry (3D) > rayとtriangleの交差判定 > geometry_intersection_3D_180513.py (v0.3) > テスト用check_intersection_3D_180512.py (v0.5)

処理概要

  • 正四角錐の内部を格子状のdipoleで埋め尽くす

importコード

geometry (3D) > rayとtriangleの交差判定 > geometry_intersection_3D_180513.py (v0.3) > テスト用check_intersection_3D_180512.py (v0.5)
の
geometry_intersection_3D_180513.py (v0.3)

  • rayとtriangleの交差判定
  • rayがtriangleに到達しない場合は「非交差」と判定

アルゴリズム

  • 原点(0, 0, 0)から各dipoleの位置までのrayを定義する
  • rayが正四角錐を構成する三角形と交差するか判定する
    • 交差するdipole: 除外
    • 交差しないdipole: 残す
  • 全てのdipoleに対して、上記の判定を行う

残ったdipoleは正四角錐の内部にあると思われる。

code v0.1

volfil_tetrahedron_180519.py
import numpy as np
import geometry_intersection_3D_180513 as GI3
import sys

'''
v0.1 May, 19, 2018
  - add check_vertex_intersecting_ray()
  - import [geometry_intersection_3D_180513]
  - add get_dipole_positions_in_rectangle()
  - add get_regular_tetrahedron_vertices()
'''


def get_regular_tetrahedron_vertices(expand=1.0):
    vrt1 = np.array([expand, expand, expand])  # vertex
    vrt2 = np.array([-expand, -expand, expand])
    vrt3 = np.array([-expand, expand, -expand])
    vrt4 = np.array([expand, -expand, -expand])

    ret = []
    ret += [np.array([vrt1, vrt2, vrt3])]
    ret += [np.array([vrt1, vrt3, vrt4])]
    ret += [np.array([vrt2, vrt3, vrt4])]
    ret += [np.array([vrt1, vrt2, vrt4])]

    return np.array(ret)


def get_dipole_positions_in_rectangle(expand=5.0):
    mx = np.linspace(-expand, expand, 9)
    my = np.linspace(-expand, expand, 9)
    mz = np.linspace(-expand, expand, 9)
    ret = []
    for idx in range(len(mx)):
        for idy in range(len(my)):
            for idz in range(len(mz)):
                ret += [np.array([mx[idx], my[idy], mz[idz]])]
    return np.array(ret)


def check_vertex_intersecting_ray(start, end, vertices):
    ray = np.array([start, end])
    for vertex in vertices:
        res, retI = GI3.check_intersection(ray, vertex)
        isIntsct = GI3.check_reachingIntersection(ray, vertex)
        if isIntsct:
            return False
    return True


tethdr = get_regular_tetrahedron_vertices(expand=10.0)
dpl = get_dipole_positions_in_rectangle(expand=12.0)

startary = np.array([0, 0, 0])
for elem in dpl:
    res = check_vertex_intersecting_ray(startary, elem, tethdr)
    if res:
        print("[%f, %f, %f]," % (elem[0], elem[1], elem[2]))

実行

$ python3 volfil_tetrahedron_180519.py > test_dipole_180519.tinkercad

上記のコマンドにてTinkercadで使うための配列を作る(下記)。

$ head test_dipole_180519.tinkercad 
[-9.000000, -9.000000, 9.000000],
[-9.000000, -6.000000, 6.000000],
[-9.000000, -3.000000, 3.000000],
[-9.000000, 0.000000, 0.000000],
[-9.000000, 3.000000, -3.000000],
[-9.000000, 6.000000, -6.000000],
[-9.000000, 9.000000, -9.000000],
[-6.000000, -9.000000, 6.000000],
[-6.000000, -6.000000, 3.000000],
[-6.000000, -6.000000, 6.000000],

可視化

TinkercardのShape Generatorのスクリプトを作り可視化する。

  • 上記のtest_dipole_180519.tinkercad をmain.jsのpositionsの定義に貼る
  • 244行目をfor (var oi = 0; oi < 119; oi++) { // oi: outer indexとする (119を指定)
    • 119: 残ったdipoleの個数

関連: Tinkercad > Shape Generator > main.js v0.5 > BCCA (Ballistic cluster cluster aggregate), N=8 , N=128 | BPCA (N=128)

main.js
// Convenience Declarations For Dependencies.
// 'Core' Is Configured In Libraries Section.
// Some of these may not be used by this example.
var Conversions = Core.Conversions;
var Debug = Core.Debug;
var Path2D = Core.Path2D;
var Point2D = Core.Point2D;
var Point3D = Core.Point3D;
var Matrix2D = Core.Matrix2D;
var Matrix3D = Core.Matrix3D;
var Mesh3D = Core.Mesh3D;
var Plugin = Core.Plugin;
var Tess = Core.Tess;
var Sketch2D = Core.Sketch2D;
var Solid = Core.Solid;
var Vector2D = Core.Vector2D;
var Vector3D = Core.Vector3D;

// v0.6 Mar. 19, 2018
//    - positions[] has geometry for regular tetrahedron
//       + obtained using [volfil_tetrahedron_180519.py] v0.1
// v0.5 Mar. 08, 2018
//    - positions[] has geometry for BCCA(N=8)
//       + BCCA: Ballistic cluster cluster aggregate
// v0.4 Mar. 07, 2018
//   - add positions[] to locate each sphere
// v0.3 Mar. 06, 2018
//   - add parameter [dist]
// v0.2 Mar. 06, 2018
//   - mesh.triangle() uses [vrt1], [vrt2], [vrt3]
// v0.1 Mar. 06, 2018
//   - branched from [Shape Generator] > [Sphere]

// Template Code:
params = [
    { "id": "r", "displayName": "Radius", "type": "length", "rangeMin": 1.0, "rangeMax": 100.0, "default": 10.0 },
    { "id": "dist", "displayName": "Distance", "type": "length", "rangeMin": 1.0, "rangeMax": 500.0, "default": 20.0 }
];

function process(params) {
    var r = params["r"];
    var dist = params["dist"];
    var ndivs = Tess.circleDivisions(r);
    var tau = 0.8506508084;
    var one = 0.5257311121;
    var lod = 0;

    while(ndivs > 6){
        lod++;
        ndivs /= 2;
    }

    Debug.log(Tess.circleDivisions(r) + " " + lod);

    var v = [
        [ tau, one, 0.0 ],
        [-tau, one, 0.0 ],
        [-tau,-one, 0.0 ],
        [ tau,-one, 0.0 ],
        [ one, 0.0, tau ],
        [ one, 0.0,-tau ],
        [-one, 0.0,-tau ],
        [-one, 0.0, tau ],
        [ 0.0, tau, one ],
        [ 0.0,-tau, one ],
        [ 0.0,-tau,-one ],
        [ 0.0, tau,-one ]
    ];

    var tris = [
        [ v[4], v[8], v[7] ],
        [ v[4], v[7], v[9] ],
        [ v[5], v[6],v[11] ],
        [ v[5],v[10], v[6] ],
        [ v[0], v[4], v[3] ],
        [ v[0], v[3], v[5] ],
        [ v[2], v[7], v[1] ],
        [ v[2], v[1], v[6] ],
        [ v[8], v[0],v[11] ],
        [ v[8],v[11], v[1] ],
        [ v[9],v[10], v[3] ],
        [ v[9], v[2],v[10] ],
        [ v[8], v[4], v[0] ],
        [v[11], v[0], v[5] ],
        [ v[4], v[9], v[3] ],
        [ v[5], v[3],v[10] ],
        [ v[7], v[8], v[1] ],
        [ v[6], v[1],v[11] ], 
        [ v[7], v[2], v[9] ],
        [ v[6],v[10], v[2] ]
    ];

    // defined in radius unit
    //    for touching bisphere, set [2]
    var positions = [
      // from [volfil_tetrahedron_180519.py] v0.1
[-9.000000, -9.000000, 9.000000],
[-9.000000, -6.000000, 6.000000],
[-9.000000, -3.000000, 3.000000],
[-9.000000, 0.000000, 0.000000],
[-9.000000, 3.000000, -3.000000],
[-9.000000, 6.000000, -6.000000],
[-9.000000, 9.000000, -9.000000],
[-6.000000, -9.000000, 6.000000],
[-6.000000, -6.000000, 3.000000],
[-6.000000, -6.000000, 6.000000],
[-6.000000, -6.000000, 9.000000],
[-6.000000, -3.000000, 0.000000],
[-6.000000, -3.000000, 3.000000],
[-6.000000, -3.000000, 6.000000],
[-6.000000, 0.000000, -3.000000],
[-6.000000, 0.000000, 0.000000],
[-6.000000, 0.000000, 3.000000],
[-6.000000, 3.000000, -6.000000],
[-6.000000, 3.000000, -3.000000],
[-6.000000, 3.000000, 0.000000],
[-6.000000, 6.000000, -9.000000],
[-6.000000, 6.000000, -6.000000],
[-6.000000, 6.000000, -3.000000],
[-6.000000, 9.000000, -6.000000],
[-3.000000, -9.000000, 3.000000],
[-3.000000, -6.000000, 0.000000],
[-3.000000, -6.000000, 3.000000],
[-3.000000, -6.000000, 6.000000],
[-3.000000, -3.000000, -3.000000],
[-3.000000, -3.000000, 0.000000],
[-3.000000, -3.000000, 3.000000],
[-3.000000, -3.000000, 6.000000],
[-3.000000, -3.000000, 9.000000],
[-3.000000, 0.000000, -6.000000],
[-3.000000, 0.000000, -3.000000],
[-3.000000, 0.000000, 0.000000],
[-3.000000, 0.000000, 3.000000],
[-3.000000, 0.000000, 6.000000],
[-3.000000, 3.000000, -9.000000],
[-3.000000, 3.000000, -6.000000],
[-3.000000, 3.000000, -3.000000],
[-3.000000, 3.000000, 0.000000],
[-3.000000, 3.000000, 3.000000],
[-3.000000, 6.000000, -6.000000],
[-3.000000, 6.000000, -3.000000],
[-3.000000, 6.000000, 0.000000],
[-3.000000, 9.000000, -3.000000],
[0.000000, -9.000000, 0.000000],
[0.000000, -6.000000, -3.000000],
[0.000000, -6.000000, 0.000000],
[0.000000, -6.000000, 3.000000],
[0.000000, -3.000000, -6.000000],
[0.000000, -3.000000, -3.000000],
[0.000000, -3.000000, 0.000000],
[0.000000, -3.000000, 3.000000],
[0.000000, -3.000000, 6.000000],
[0.000000, 0.000000, -9.000000],
[0.000000, 0.000000, -6.000000],
[0.000000, 0.000000, -3.000000],
[0.000000, 0.000000, 0.000000],
[0.000000, 0.000000, 3.000000],
[0.000000, 0.000000, 6.000000],
[0.000000, 0.000000, 9.000000],
[0.000000, 3.000000, -6.000000],
[0.000000, 3.000000, -3.000000],
[0.000000, 3.000000, 0.000000],
[0.000000, 3.000000, 3.000000],
[0.000000, 3.000000, 6.000000],
[0.000000, 6.000000, -3.000000],
[0.000000, 6.000000, 0.000000],
[0.000000, 6.000000, 3.000000],
[0.000000, 9.000000, 0.000000],
[3.000000, -9.000000, -3.000000],
[3.000000, -6.000000, -6.000000],
[3.000000, -6.000000, -3.000000],
[3.000000, -6.000000, 0.000000],
[3.000000, -3.000000, -9.000000],
[3.000000, -3.000000, -6.000000],
[3.000000, -3.000000, -3.000000],
[3.000000, -3.000000, 0.000000],
[3.000000, -3.000000, 3.000000],
[3.000000, 0.000000, -6.000000],
[3.000000, 0.000000, -3.000000],
[3.000000, 0.000000, 0.000000],
[3.000000, 0.000000, 3.000000],
[3.000000, 0.000000, 6.000000],
[3.000000, 3.000000, -3.000000],
[3.000000, 3.000000, 0.000000],
[3.000000, 3.000000, 3.000000],
[3.000000, 3.000000, 6.000000],
[3.000000, 3.000000, 9.000000],
[3.000000, 6.000000, 0.000000],
[3.000000, 6.000000, 3.000000],
[3.000000, 6.000000, 6.000000],
[3.000000, 9.000000, 3.000000],
[6.000000, -9.000000, -6.000000],
[6.000000, -6.000000, -9.000000],
[6.000000, -6.000000, -6.000000],
[6.000000, -6.000000, -3.000000],
[6.000000, -3.000000, -6.000000],
[6.000000, -3.000000, -3.000000],
[6.000000, -3.000000, 0.000000],
[6.000000, 0.000000, -3.000000],
[6.000000, 0.000000, 0.000000],
[6.000000, 0.000000, 3.000000],
[6.000000, 3.000000, 0.000000],
[6.000000, 3.000000, 3.000000],
[6.000000, 3.000000, 6.000000],
[6.000000, 6.000000, 3.000000],
[6.000000, 6.000000, 6.000000],
[6.000000, 6.000000, 9.000000],
[6.000000, 9.000000, 6.000000],
[9.000000, -9.000000, -9.000000],
[9.000000, -6.000000, -6.000000],
[9.000000, -3.000000, -3.000000],
[9.000000, 0.000000, 0.000000],
[9.000000, 3.000000, 3.000000],
[9.000000, 6.000000, 6.000000],
[9.000000, 9.000000, 9.000000]
    ];

    function norm(a) {
        var l = Math.sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
        return [ a[0]/l, a[1]/l, a[2]/l ];
    }

    function sum(a, b) {
        return [ (a[0]+b[0]), (a[1]+b[1]), (a[2]+b[2]) ];
    }

    function xf(a, s) {
        return [ a[0]*s, a[1]*s, a[2]*s + s ];
    }

    for (var i = 0; i < lod; i++) {
        var ntris = [];
        for(var j = 0; j < tris.length; j++) {
            var ma = norm(sum(tris[j][0], tris[j][1]));
            var mb = norm(sum(tris[j][1], tris[j][2]));
            var mc = norm(sum(tris[j][2], tris[j][0]));

            ntris.push([tris[j][0], ma, mc]);
            ntris.push([tris[j][1], mb, ma]);
            ntris.push([tris[j][2], mc, mb]);
            ntris.push([ma, mb, mc]);
        }
        tris = ntris;
    }

    var mesh = new Mesh3D();
    for (var oi = 0; oi < 119; oi++) { // oi: outer index
        for (var ii = 0; ii < tris.length; ii++) {
            var vrt1 = tris[ii][0];
            var vrt2 = tris[ii][1];
            var vrt3 = tris[ii][2];
            var xf1 = xf(vrt1,r);
            var xf2 = xf(vrt2,r);
            var xf3 = xf(vrt3,r);
            for (var pi = 0; pi < 3; pi++) { // pi: position index
                xf1[pi] += positions[oi][pi] * dist / 2; // 2: for positions is defined in radisu unit
                xf2[pi] += positions[oi][pi] * dist / 2;
                xf3[pi] += positions[oi][pi] * dist / 2;
            }
            mesh.triangle(xf1, xf2, xf3);
        }
    }

    var solid = Solid.make(mesh);

    return solid;
}

作ったShapeを貼ったのが以下。

[3Dを表示]を選択すると、回転させることができる。

qiita.png

関連

上記に対する2D版の処理は下記となる。
geometry > star shape(星形)の内部をdipoleでfillする > findDipolesInisideStarShape_180429 v0.2 > shapely使用 > 1秒で完了 (MESH_RESOL=50)

code v0.2 > Numpy用ファイル出力

  • Numpyでgeonfromtxt()で読めるファイル出力
    • [dipole_180520.res]
volfil_tetrahedron_180519.py
import numpy as np
import geometry_intersection_3D_180513 as GI3
import sys
import os

'''
v0.2 May, 20, 2018
  - output as file for reading with Numpy genfromtxt()
    + add [OUTPUT]
    + remove [OUTPUT] before file save
v0.1 May, 19, 2018
  - add check_vertex_intersecting_ray()
  - import [geometry_intersection_3D_180513]
  - add get_dipole_positions_in_rectangle()
  - add get_regular_tetrahedron_vertices()
'''

OUTFILE = 'dipole_180520.res'  # for Numpy genfromtxt()


def get_regular_tetrahedron_vertices(expand=1.0):
    vrt1 = np.array([expand, expand, expand])  # vertex
    vrt2 = np.array([-expand, -expand, expand])
    vrt3 = np.array([-expand, expand, -expand])
    vrt4 = np.array([expand, -expand, -expand])

    ret = []
    ret += [np.array([vrt1, vrt2, vrt3])]
    ret += [np.array([vrt1, vrt3, vrt4])]
    ret += [np.array([vrt2, vrt3, vrt4])]
    ret += [np.array([vrt1, vrt2, vrt4])]

    return np.array(ret)


def get_dipole_positions_in_rectangle(expand=5.0):
    mx = np.linspace(-expand, expand, 9)
    my = np.linspace(-expand, expand, 9)
    mz = np.linspace(-expand, expand, 9)
    ret = []
    for idx in range(len(mx)):
        for idy in range(len(my)):
            for idz in range(len(mz)):
                ret += [np.array([mx[idx], my[idy], mz[idz]])]
    return np.array(ret)


def check_vertex_intersecting_ray(start, end, vertices):
    ray = np.array([start, end])
    for vertex in vertices:
        res, retI = GI3.check_intersection(ray, vertex)
        isIntsct = GI3.check_reachingIntersection(ray, vertex)
        if isIntsct:
            return False
    return True


tethdr = get_regular_tetrahedron_vertices(expand=10.0)
dpl = get_dipole_positions_in_rectangle(expand=12.0)

startary = np.array([0, 0, 0])

os.remove(OUTFILE)
with open(OUTFILE, 'ab') as wfp:
    for elem in dpl:
        res = check_vertex_intersecting_ray(startary, elem, tethdr)
        if res:
            print("[%f, %f, %f]," % (elem[0], elem[1], elem[2]))
            wrk = [elem]  # to avoid element-wise-LF for 1-D array
            np.savetxt(wfp, wrk, fmt='%.7f', delimiter=' ')
0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?