動作環境

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を表示]を選択すると、回転させることができる。

関連

上記に対する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=' ')