apt-get update
apt-get install -y ca-certificates
update-ca-certificates

apt-get update
apt-get install -y software-properties-common

apt-get update
apt-get install -y –no-install-recommends apt-transport-https ca-certificates gnupg
update-ca-certificates

add-apt-repository ppa:fenics-packages/fenics
apt-get update

apt-get install -y fenics

python3 -c “import dolfin; print(dolfin.version)”

import math, random
from gengeo import Vector3

def random_euler_angles():
alpha = random.uniform(0, 2math.pi) # 回転Z beta = random.uniform(0, 2math.pi) # 回転Y
gamma = random.uniform(0, 2*math.pi) # 回転X
return alpha, beta, gamma

def rotation_matrix_from_euler(alpha, beta, gamma):
# 回転Z(alpha)
Rz = [
[math.cos(alpha), -math.sin(alpha), 0],
[math.sin(alpha), math.cos(alpha), 0],
[0, 0, 1]
]
# 回転Y(beta)
Ry = [
[math.cos(beta), 0, math.sin(beta)],
[0, 1, 0],
[-math.sin(beta), 0, math.cos(beta)]
]
# 回転X(gamma)
Rx = [
[1, 0, 0],
[0, math.cos(gamma), -math.sin(gamma)],
[0, math.sin(gamma), math.cos(gamma)]
]

# 合成回転行列 R = R_x * R_y * R_z
# 順番に行列積を行う関数を用意
def mat_mul(A, B):
    return [
        [
            A[i][0]*B[0][j] + A[i][1]*B[1][j] + A[i][2]*B[2][j] 
            for j in range(3)
        ] for i in range(3)
    ]
Rzy = mat_mul(Rz, Ry)
R = mat_mul(Rzy, Rx)  # 順番を変更して好みのオイラー順序にもできる
return R

def apply_rotation(v, R):
# vはVector3、Rは3×3リスト
x = v.X(); y = v.Y(); z = v.Z()
rx = R[0][0]x + R[0][1]y + R[0][2]z ry = R[1][0]x + R[1][1]y + R[1][2]z
rz = R[2][0]x + R[2][1]y + R[2][2]*z
return Vector3(rx, ry, rz)

クラスターオフセット例

r = 3e-6
cluster_offsets = [
Vector3(0,0,0),
Vector3(2r,0,0), Vector3(r, math.sqrt(3)r,0)
]

ランダム回転適用例

alpha, beta, gamma = random_euler_angles()
R = rotation_matrix_from_euler(alpha, beta, gamma)

rotated_offsets = [apply_rotation(off, R) for off in cluster_offsets]

これでrotated_offsetsはランダムな配向になったクラスター座標を表す

from gengeo import *
import math, random

クラスターオフセット（初期状態で原点近傍に定義）

r = 3e-6
cluster_offsets = [
Vector3(0,0,0),
Vector3(2r,0,0), Vector3(r, math.sqrt(3)r,0)
]

def random_rotation_transform():
# ランダムなEuler角度を生成（0～360度）
alpha = random.uniform(0, 2math.pi) beta = random.uniform(0, 2math.pi)
gamma = random.uniform(0, 2*math.pi)

# Transformオブジェクトを作成
t = Transform()
# 順序はZ-Y-X回りなど任意
t.rotateZ(alpha)
t.rotateY(beta)
t.rotateX(gamma)
return t

def is_overlapping(new_spheres, existing_spheres, margin=0.0):
# オーバーラップチェック
for ns in new_spheres:
nc = ns.Centre()
nr = ns.Radius()
if (nc.X()-nr < 0 or nc.X()+nr > xdim or
nc.Y()-nr < 0 or nc.Y()+nr > ydim or
nc.Z()-nr < 0 or nc.Z()+nr > zdim):
return True
for es in existing_spheres:
ec = es.Centre()
er = es.Radius()
dist = (nc-ec).norm()
if dist < nr+er+margin:
return True
return False

パラメータ設定(例)

xdim = 50e-6
ydim = 50e-6
zdim = 50e-6
minPoint = Vector3(0,0,0)
maxPoint = Vector3(xdim, ydim, zdim)
mntable = MNTable3D(minPoint,maxPoint,2.5*r,1)
box = BoxWithPlanes3D(minPoint,maxPoint)

… 境界面の定義と追加は省略

num_clusters = 10
existing_spheres = []
random_seed = 12345
random.seed(random_seed)

for i in range(num_clusters):
max_tries = 1000
inserted = False
for attempt in range(max_tries):
cx = random.uniform(0, xdim)
cy = random.uniform(0, ydim)
cz = random.uniform(0, zdim)

    # クラスターオフセットにランダム回転適用
    t = random_rotation_transform()
    rotated_offsets = [t.transform(off) for off in cluster_offsets]

    cluster_spheres = []
    for off in rotated_offsets:
        center = Vector3(cx+off.X(), cy+off.Y(), cz+off.Z())
        s = Sphere(center, r)
        s.setTag(1)
        cluster_spheres.append(s)

    if not is_overlapping(cluster_spheres, existing_spheres):
        for s in cluster_spheres:
            mntable.insert(s)
            existing_spheres.append(s)
        inserted = True
        break
if not inserted:
    print(f"Failed to insert cluster {i}.")

ボンド生成やファイル出力など

mntable.generateBonds(0, 1.0e-5, 1)
mntable.write(“clustered_box_random_orient.geo”,1)
mntable.write(“clustered_box_random_orient.vtu”,2)
print(“Inserted clusters with random orientations.”)

from gengeo import *
import math, random

パラメータ設定

xdim = 50e-6
ydim = 50e-6
zdim = 50e-6
minPoint = Vector3(0,0,0)
maxPoint = Vector3(xdim, ydim, zdim)
mntable = MNTable3D(minPoint, maxPoint, 2.5*(3e-6), 1) # 検索範囲は適宜変更
box = BoxWithPlanes3D(minPoint,maxPoint)

… BoxWithPlanes3D に各平面をaddPlane

クラスター定義（3粒子）

3粒子を正三角形に配置する例

r = 3e-6 # 半径

原点近傍に3粒子クラスターを定義

cluster_offsets = [
Vector3(0,0,0),
Vector3(2r,0,0), Vector3(r, math.sqrt(3)r,0)
]

def is_overlapping(new_spheres, existing_spheres, margin=0.0):
# 簡易オーバーラップチェック
for ns in new_spheres:
nc = ns.Centre()
nr = ns.Radius()
# 境界チェック
if (nc.X()-nr < 0 or nc.X()+nr > xdim or
nc.Y()-nr < 0 or nc.Y()+nr > ydim or
nc.Z()-nr < 0 or nc.Z()+nr > zdim):
return True
# 既存粒子との重なりチェック
for es in existing_spheres:
ec = es.Centre()
er = es.Radius()
dist = (nc-ec).norm()
if dist < nr+er+margin:
return True
return False

挿入処理

num_clusters = 10 # 10クラスター挿入例
existing_spheres = [] # mntableから取得可能だが、簡易的にリストで管理
random_seed = 12345
random.seed(random_seed)

for i in range(num_clusters):
max_tries = 1000
inserted = False
for attempt in range(max_tries):
# クラスター中心をランダムに選ぶ
cx = random.uniform(0, xdim)
cy = random.uniform(0, ydim)
cz = random.uniform(0, zdim)
cluster_spheres = []
for off in cluster_offsets:
center = Vector3(cx+off.X(), cy+off.Y(), cz+off.Z())
s = Sphere(center, r)
s.setTag(1) # クラスター粒子のタグを1で揃えるなど
cluster_spheres.append(s)

    if not is_overlapping(cluster_spheres, existing_spheres):
        # 挿入可能
        for s in cluster_spheres:
            mntable.insert(s)
            existing_spheres.append(s)
        inserted = True
        break
if not inserted:
    print(f"Failed to insert cluster {i} after {max_tries} attempts.")

最後にボンド生成など可能

mntable.generateBonds(0, 1.0e-5, 1)

mntable.write(“clustered_box.geo”,1)
mntable.write(“clustered_box.vtu”,2)

print(“Insertion of 3-particle clusters complete.”)

# 粒子挿入ループ maxLoops = 10 previous_count = 0 for loop_i in range(maxLoops): packer.generatePacking(box, mntable, 0, 1) spheres = mntable.getSphereListFromGroup(0) current_count = len(spheres) if current_count <= previous_count: # 粒子数が増えなくなったら終了 break else: # ボンド生成 (距離3e-6mで近接粒子間にボンド) # 3e-6 mは粒子半径と同等スケールで、近接粒子にボンドがつきやすい mntable.generateBonds(0, 3e-6, 0) previous_count = current_count print(f”Loop {loop_i}: particle count = {current_count}”)

!/usr/bin/env python

— geometry setup script for block with smooth sides —

from gengeo import *

– input parameters —

block dimensions (50µm角)

xdim = 50e-6
ydim = 50e-6
zdim = 50e-6

particle size range

minRadius = 0.2
maxRadius = 1.0

———————

corner points

minPoint = Vector3(0.0,0.0,0.0)
maxPoint = Vector3(xdim,ydim,zdim)

neighbour table

mntable = MNTable3D(minPoint,maxPoint,2.5*maxRadius,1)

block volume

box = BoxWithPlanes3D(minPoint,maxPoint)

boundary planes

bottomPlane = Plane(minPoint,Vector3(0.0,1.0,0.0))
leftPlane = Plane(minPoint,Vector3(1.0,0.0,0.0))
frontPlane = Plane(minPoint,Vector3(0.0,0.0,1.0))
topPlane = Plane(maxPoint,Vector3(0.0,-1.0,0.0))
rightPlane = Plane(maxPoint,Vector3(-1.0,0.0,0.0))
backPlane = Plane(maxPoint,Vector3(0.0,0.0,-1.0))

add them to the box

box.addPlane(bottomPlane)
box.addPlane(leftPlane)
box.addPlane(frontPlane)
box.addPlane(topPlane)
box.addPlane(rightPlane)
box.addPlane(backPlane)

— setup packer —

iteration parameters

insertFails = 1000
maxIter = 1000
tol = 1.0e-6

packer

packer = InsertGenerator3D(minRadius,maxRadius,insertFails,maxIter,tol,False)

pack particles into volume

packer.generatePacking(box,mntable,0,1)

create bonds between neighbouring particles:

mntable.generateBonds(0,1.0e-5,0)

calculate and print the porosity:

volume = xdimydimzdim
porosity = (volume – mntable.getSumVolume(groupID=0))/volume
print “Porosity: “, porosity

write a geometry file

mntable.write(“smooth_box.geo”, 1)
mntable.write(“smooth_box.vtu”, 2)

— geometry setup script for 50µm box —

from gengeo import *

パラメータ設定(50µm角)

xdim = 50e-6
ydim = 50e-6
zdim = 50e-6

粒子半径範囲(例：0.2～0.5 µm)

minRadius = 0.2e-6
maxRadius = 0.5e-6

minPoint = Vector3(0.0,0.0,0.0)
maxPoint = Vector3(xdim, ydim, zdim)
box = BoxWithPlanes3D(minPoint, maxPoint)

MNTable3D設定

mntable = MNTable3D(minPoint, maxPoint, 2.5*maxRadius, 1)

InsertGenerator3Dパラメータ

insertFails = 1000
maxIter = 500
tol = 1.0e-6
packer = InsertGenerator3D(minRadius, maxRadius, insertFails, maxIter, tol)

粒子挿入のループ処理

maxLoops = 5 # 最大5回挿入再試行
previous_count = 0

for loop_i in range(maxLoops):
packer.generatePacking(box, mntable)
spheres = mntable.getSphereListFromGroup(0)
current_count = len(spheres)

if current_count <= previous_count:
    print(f"No increase in particle count after loop {loop_i}, stopping loop.")
    break

# 粒子数が増えたらボンドを再生成して内部状態を安定化
mntable.generateBonds(0, 1.0e-5, 0)
print(f"After loop {loop_i}, particle count: {current_count}")
previous_count = current_count

ループ終了後にもう一度ボンド生成で状態安定化（必要に応じて）

mntable.generateBonds(0, 1.0e-5, 0)

ファイル出力 (VTU形式)

mntable.write(“box.geo”, 1)
mntable.write(“box.vtu”, 2)

final_count = len(mntable.getSphereListFromGroup(0))
print(“Geometry generation complete.”)
print(f”Final particle count in 50µm box: {final_count}”)
print(“Generated files: box.geo (ASCII), box.vtu (VTU). Try loading box.vtu in ParaView.”)

!/usr/bin/env python

— geometry setup script for simple box —

from gengeo import * # これによりVector3が利用可能になります

パラメータ

xdim = 50e-6
ydim = 50e-6
zdim = 50e-6

minRadius = 0.2e-6
maxRadius = 0.5e-6

Vector3を使用

minPoint = Vector3(0.0, 0.0, 0.0)
maxPoint = Vector3(xdim, ydim, zdim)

box = BoxWithPlanes3D(minPoint, maxPoint)

mntable = MNTable3D(minPoint, maxPoint, 2.5*maxRadius, 1)

insertFails = 10000
maxIter = 500
tol = 1.0e-6
packer = InsertGenerator3D(minRadius, maxRadius, insertFails, maxIter, tol)

maxLoops = 10
previous_count = 0

for loop_i in range(maxLoops):
packer.generatePacking(box, mntable)
sphere_list = mntable.getSphereListFromGroup(0)

current_count = len(sphere_list)
if current_count <= previous_count:
print(f”No increase in particle count after loop {loop_i}. Stopping.”)
break
else:
print(f”After loop {loop_i}, particle count: {current_count}”)
previous_count = current_count

mntable.generateBonds(0, 1.0e-5, 0)

mntable.write(“box.geo”, 1)
mntable.write(“box.vtu”, 2)

sphere_list = mntable.getSphereListFromGroup(0)

particle_count = len(sphere_list)

print(f”Complete. Final particle count: {particle_count}”)

// —————————-
// Groups and Regions
// —————————-
Group {
Region[1] = Region[Physical Surface(“MaterialRegion1”)]; // Upper circle
Region[2] = Region[Physical Surface(“MaterialRegion2”)]; // Lower circle
Region[3] = Region[Physical Surface(“RectangleRegion”)]; // Rectangle (air)
Boundary[1] = Region[Physical Line(“Bottomedge”)]; // Bottom edge
Boundary[3] = Region[Physical Line(“Topedge”)]; // Top edge
}

// —————————-
// Material Properties
// —————————-
Function {
MuSoft[] = 1e-2; // Soft magnetic material permeability
MuAir[] = 4 * Pi * 1e-7; // Air permeability
PhiBoundary1[] = 0; // Scalar potential at Bottomedge
PhiBoundary3[] = 1; // Scalar potential at Topedge
}

// —————————-
// Function Spaces
// —————————-
FunctionSpace {
{ Name H1 ; // Define the H1 space
Type Form0 ; // Scalar field
BasisFunction {
{ Name sn ; NameOfCoef nodal ; Function BF_Node ; }
}
Constraint {
{ NameOfConstraint BoundaryConditions ; }
}
}
}

// —————————-
// Formulation
// —————————-
Formulation {
{ Name MagneticScalarFormulation ;
Type FemEquation ;
Quantity {
{ Name MagneticScalarPotential ; Type Local ; NameOfSpace H1 ; }
}
Equation {
// Weak form of the scalar magnetic potential equation
Galerkin { [ 1/MuSoft[] * Dof{MagneticScalarPotential}, {MagneticScalarPotential} ]; In Region[1, 2]; }
Galerkin { [ 1/MuAir[] * Dof{MagneticScalarPotential}, {MagneticScalarPotential} ]; In Region[3]; }
}
}
}

// —————————-
// Boundary Conditions
// —————————-
Constraint {
{ Name BoundaryConditions ;
Type Assign ;
Case {
{ Region Boundary[1]; Value PhiBoundary1[]; } // Bottom edge
{ Region Boundary[3]; Value PhiBoundary3[]; } // Top edge
}
}
}

// —————————-
// Resolution
// —————————-
Resolution {
{ Name MagneticScalarResolution ;
System {
{ NameOfFormulation MagneticScalarFormulation ; }
}
Operation {
Generate[] ;
Solve[] ;
SaveSolution[MagneticScalarPotential] ;
}
}
}

// —————————-
// Post-Processing
// —————————-
PostProcessing {
{ Name MagneticScalarField ;
Quantity {
{ Name MagneticScalarPotential ; Value { Local { [ {MagneticScalarPotential} ]; In Region[1, 2, 3]; }}}
}
}
}
—————————- —————————- —————————-
Esys

Docker のインストール　省略

Docker コンテナ作成

docker run -it -d –name esysparticle-container ubuntu

必要な依存関係のインストール

apt update
apt upgrade

apt install -y build-essential cmake git python3 python3-pip wget
apt install libtool automake build-essential bzr cmake mpi-default-bin mpi-default-dev libboost-all-dev python-is-python3

ソースコードの取得

dionの親切に乗りたかったがNG

https://launchpad.net/esys-particle/trunk/2.3/+download

mkdir -p ESyS-Particle
cd ESyS-Particle
bzr branch lp:esys-particle esys-particle-trunk
cd esys-particle-trunk

直DL

wget https://launchpad.net/esys-particle/+download#:~:text=esys%2Dparticle%2D3.0%2Dalpha.tar.gz　–no-check-certificate

アーカイブを展開

tar -xvzf esys-particle-2.3.tar.gz
cd esys-particle-2.3

MPI-3.0 以降では削除された関数が含まれているため対応

export OMPI_MCA_mpi_compatibility=true

apt remove –purge mpi-default-bin mpi-default-dev
apt install openmpi-bin=1.10.7-1build1 openmpi-common=1.10.7-1build1 libopenmpi-dev=1.10.7-1build1

ビルド

mkdir build
cd build

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

cmake ..

make -j$(nproc)

make install

環境変数

export PATH=/usr/local/bin:$PATH
export LIBRARY_PATH=/usr/local/lib:$LIBRARY_PATH
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export PYTHONPATH=/usr/local/lib/python3/dist-packages:$PYTHONPATH

永続化

echo ‘export PATH=/usr/local/bin:$PATH’ >> ~/.bashrc
echo ‘export LIBRARY_PATH=/usr/local/lib:$LIBRARY_PATH’ >> ~/.bashrc
echo ‘export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH’ >> ~/.bashrc
echo ‘export PYTHONPATH=/usr/local/lib/python3/dist-packages:$PYTHONPATH’ >> ~/.bashrc
source ~/.bashrc

sudo apt update
sudo apt install –reinstall ca-certificates
sudo update-ca-certificates

ビルドとインストール

MPI-3.0 以降では削除された関数が含まれているため対応

export OMPI_MCA_mpi_compatibility=true

apt remove –purge mpi-default-bin mpi-default-dev
apt install openmpi-bin=1.10.7-1build1 openmpi-common=1.10.7-1build1 libopenmpi-dev=1.10.7-1build1

mkdir -p ~/build/esys-particle
cd ~/build/esys-particle
cmake -S ~/src/ESyS-Particle/esys-particle-trunk
make
make install

環境変数の設定 

 export PATH=/usr/local/bin/:$PATH
export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
export PYTHONPATH=/usr/local/lib/python3.8/dist-packages:/usr/local/lib/python3/dist-packages:$PYTHONPATH

mpirun -np 2 esysparticle ../Doc/Examples/two_particle.py

警告無視

export OMPI_ALLOW_RUN_AS_ROOT=1
export OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1

一般ユーザ追加

adduser mpiuser
su – mpiuser
mpirun -np 2 ./your_program

インストールの確認

mkdir ~/scratch
cd ~/scratch
mpirun -np 2 esysparticle ~/src/ESyS-Particle/esys-particle-trunk/Doc/Examples/two_particle.py

https://launchpad.net/esys-particle/+download#:~:text=ESyS%2DParticle%2D2.3.5.tar.gz

wget https://launchpad.net/esys-particle/+download#:~:text=ESyS%2DParticle%2D2.3.5.tar.gz

cd /usr/BUILD/sources
wget https://boostorg.jfrog.io/artifactory/main/release/1.52.0/source/boost_1_52_0.tar.bz2
tar –bzip2 -xf boost_1_52_0.tar.bz2
cd boost_1_52_0

./bootstrap.sh –prefix=/usr/BUILD –with-libraries=filesystem,python,regex,system
./b2
./b2 install

cd /usr/BUILD/sources
wget https://github.com/Kitware/CMake/releases/download/v3.22.0/cmake-3.22.0.tar.gz
tar xzf cmake-3.22.0.tar.gz
cd cmake-3.22.0

export CXXFLAGS=”-std=c++11″ export CFLAGS=”-std=c++11″
./configure –prefix=/usr/BUILD CXXFLAGS=”-std=c++11″
make
make install

/usr/BUILD/bin/cmake –version

sudo apt-get install libncurses5-dev

cd /usr/BUILD/sources
wget https://www.vtk.org/files/release/6.1/VTK-6.1.0.tar.gz
tar xzf VTK-6.1.0.tar.gz
mkdir /usr/BUILD/sources/VTK-build
cd /usr/BUILD/sources/VTK-build

cmake ../VTK-6.1.0 \
-DCMAKE_BUILD_TYPE=Release \
-DCMAKE_INSTALL_PREFIX=/usr/BUILD \
-DBUILD_SHARED_LIBS=ON \
-DVTK_WRAP_PYTHON=ON \
-DModule_vtkPython=ON \
-DModule_vtkPythonInterpreter=ON \
-DPYTHON_EXECUTABLE=/usr/bin/python2.7 \
-DPYTHON_INCLUDE_DIR=/usr/include/python2.7 \
-DPYTHON_LIBRARY=/usr/lib/x86_64-linux-gnu/libpython2.7.so \
-DVTK_INSTALL_PYTHON_MODULE_DIR=/usr/BUILD/lib/python2.7/site-packages

sudo apt-get install -y libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev

ls /usr/include/python2.7/patchlevel.h

sudo apt-get install -y python2.7-dev

ls /usr/local/include/python2.6/patchlevel.h

apt-get install libxt-dev

make -j$(nproc)
make install

sudo apt-get install -y software-properties-common

sudo apt-get update
sudo apt-get install -y ca-certificates
sudo update-ca-certificates

cd /home/ubuntu/build/sources
wget https://cmake.org/files/v2.8/cmake-2.8.12.2.tar.gz

tar xfz cmake-2.8.12.2.tar.gz
cd cmake-2.8.12.2ß

./configure –prefix=/home/ubuntu/build
make
make install

sudo apt-get install -y gcc-7 g++-7

sudo apt-get install -y build-essential libgmp-dev libmpfr-dev libmpc-dev texinfo

wget https://ftp.gnu.org/gnu/gcc/gcc-7.5.0/gcc-7.5.0.tar.gz tar xzf gcc-7.5.0.tar.gz cd gcc-7.5.0 ./contrib/download_prerequisites

../configure –prefix=/usr/local/gcc-7.5 –enable-languages=c,c++ –disable-multilib make -j$(nproc) sudo make install

export PATH=/usr/local/gcc-7.5/bin:$PATH

gcc –version

../configure –prefix=/home/ubuntu/build/gcc-7.5.0 \
–enable-languages=c,c++ \
–disable-multilib \
–disable-libsanitizer

export PATH=/home/ubuntu/build/cmake-3.10.2/bin:$PATH

cmake ../sources/VTK-6.2.0 \
-DCMAKE_BUILD_TYPE=Release \
-DCMAKE_INSTALL_PREFIX=/home/ubuntu/build/VTK-6.2 \
-DBUILD_SHARED_LIBS=ON \
-DVTK_USE_X=ON \
-DVTK_Group_StandAlone=ON \
-DVTK_Group_Rendering=ON \
-DVTK_USE_GL2PS=ON \
-DCMAKE_CXX_FLAGS=”-std=c++11 -Wno-register”

wget https://launchpad.net/esys-particle/trunk/3.0-alpha/+download/gengeo-163.tar.gz
tar -xzf gengeo-163.tar.gz
cd gengeo

bzr branch lp:esys-particle/gengeo

sudo apt-get update
sudo apt-get install -y autoconf automake libtool

from gengeo import RandomBoxPacker, Vec3, VtkSphereWriter

Define the simulation box dimensions (50 µm x 50 µm x 50 µm)

box_min = Vec3(0, 0, 0) # Minimum corner of the box
box_max = Vec3(50, 50, 50) # Maximum corner of the box

Define particle properties

particle_count = 50 # Number of particles to generate
particle_radius_min = 0.5 # Minimum particle radius (µm)
particle_radius_max = 1.0 # Maximum particle radius (µm)

Initialize the random particle packer

packer = RandomBoxPacker(
radiusMin=particle_radius_min, # Minimum particle radius
radiusMax=particle_radius_max, # Maximum particle radius
region=box_max – box_min, # Box dimensions
tolerance=0.01, # Tolerance for packing
numSpheres=particle_count # Number of particles to generate
)

Generate the particles

packer.generate()

Retrieve the generated particle data

particles = packer.getParticles()

Create a VTK writer to output the particles

output_file = “particles.vtk”
writer = VtkSphereWriter(output_file)

Add each particle to the VTK writer

for particle in particles:
writer.addSphere(
position=particle.position, # Position of the particle
radius=particle.radius # Radius of the particle
)

Write the particle data to the VTK file

writer.write()

print(f”Particle data has been saved to {output_file}.”)

apt-get update && apt-get install -y \
build-essential \
cmake \
libboost-all-dev \
python3 \
python3-dev \
python3-pip

import gmsh
import numpy as np

Gmsh 初期化

gmsh.initialize()

モデルのロード

gmsh.open(“simple_electromagnetic.geo”)

メッシュ生成 (2D)

gmsh.model.mesh.generate(2)

磁界の簡易計算 (単純なフィールド例)

注: 実際の電磁気解析は追加の数値解法ライブラリが必要

nodes = gmsh.model.mesh.getNodes()
x, y, z = nodes[1].reshape(-1, 3).T

簡単な磁界の計算（例: 磁場が z 軸方向に一様である場合）

Bx = np.zeros_like(x)
By = np.zeros_like(y)
Bz = np.ones_like(z) # 一様な z 軸方向磁場

メッシュノードごとに磁場を ParaView 用に保存

vtk_file = “simple_magnetic_field.vtk”
with open(vtk_file, “w”) as f:
# VTK ヘッダ
f.write(“# vtk DataFile Version 3.0\n”)
f.write(“Magnetic field data\n”)
f.write(“ASCII\n”)
f.write(“DATASET UNSTRUCTURED_GRID\n”)
f.write(f”POINTS {len(x)} float\n”)
for xi, yi, zi in zip(x, y, z):
f.write(f”{xi} {yi} {zi}\n”)
f.write(“\n”)
f.write(f”POINT_DATA {len(x)}\n”)
f.write(“VECTORS MagneticField float\n”)
for bxi, byi, bzi in zip(Bx, By, Bz):
f.write(f”{bxi} {byi} {bzi}\n”)

print(f”Magnetic field data saved to {vtk_file}”)

Gmsh の終了

gmsh.finalize()

// Gmsh スクリプト: 単純な 2D 平面で電磁気解析用メッシュ生成

// モデルの設定
SetFactory(“OpenCASCADE”);

// 四角形導体の作成 (幅 1×1 の領域)
Rectangle(1) = {0, 0, 0, 1, 1};

// 周囲の空間 (解析領域) を作成
Rectangle(2) = {-2, -2, 0, 5,å 5};

// メッシュ生成の設定
Physical Surface(“Conductor”) = {1};
Physical Surface(“Air”) = {2};

// メッシュを作成
Mesh 2;

// 結果を ParaView 用にエクスポート (VTK)
Save “simple_electromagnetic.vtk”;å

// モデルの設定
SetFactory(“OpenCASCADE”);

// 導体 (1×1 四角形)
Rectangle(1) = {0, 0, 0, 1, 1};

// 周囲の空気領域 (解析領域: 5×5)
Rectangle(2) = {-2, -2, 0, 5, 5};

// 導体を空気領域から引き算 (導体が空気の一部を占める構造にする)
BooleanDifference{ Surface{2}; Delete; }{ Surface{1}; Delete; }

// 物理領域を定義
Physical Surface(“Conductor”) = {1}; // 導体領域
Physical Surface(“Air”) = {2}; // 空気領域

// メッシュを生成
Mesh 2;

// ParaView 用に出力
Save “two_regions.vtk”;

SetFactory(“OpenCASCADE”);

// モデル領域を作成
Rectangle(1) = {0, 0, 0, 1, 1};
Rectangle(2) = {-2, -2, 0, 5, 5};

// サーフェスの差を作成
out[] = BooleanDifference{ Surface{2}; Delete; }{ Surface{1}; Delete; };

// 結果のタグを表示して確認
Printf(“Resulting surface tags: %g and %g”, out[0], out[1]);

// 物理領域を再定義
Physical Surface(“Conductor”) = {out[1]}; // 導体領域
Physical Surface(“Air”) = {out[0]}; // 空気領域

// メッシュ生成
Mesh 2;

// VTK 出力
Save “two_regions.vtk”;

// モデルの設定
SetFactory(“Built-in”);

// ポイントを定義
// 空気領域の四隅
Point(1) = {-2, -2, 0, 1.0};
Point(2) = {3, -2, 0, 1.0};
Point(3) = {3, 3, 0, 1.0};
Point(4) = {-2, 3, 0, 1.0};

// 導体領域の四隅
Point(5) = {0, 0, 0, 1.0};
Point(6) = {1, 0, 0, 1.0};
Point(7) = {1, 1, 0, 1.0};
Point(8) = {0, 1, 0, 1.0};

// ラインを定義
// 空気領域
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

// 導体領域
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};

// ラインループを作成
Line Loop(1) = {1, 2, 3, 4}; // 空気領域
Line Loop(2) = {5, 6, 7, 8}; // 導体領域

// サーフェスを作成
Plane Surface(1) = {1}; // 空気領域
Plane Surface(2) = {2}; // 導体領域

// 物理領域を定義
Physical Surface(“Air”) = {1};
Physical Surface(“Conductor”) = {2};

// メッシュを生成
Mesh 2;

// VTK 出力
Save “line_based_structure.vtk”;

import gmsh
import numpy as np
import math

初期化

gmsh.initialize()

モデルを読み込む (Gmsh スクリプトで生成された .geo ファイル)

gmsh.open(“line_based_structure.geo”)

メッシュを生成 (2D)

gmsh.model.mesh.generate(2)

ノード情報を取得

node_tags, node_coords, _ = gmsh.model.mesh.getNodes()

座標を取得

node_coords = np.array(node_coords).reshape(-1, 3)
x_coords = node_coords[:, 0]
y_coords = node_coords[:, 1]

磁場計算: 簡易モデル (中心に単位電流を仮定)

mu_0 = 4 * math.pi * 1e-7 # 真空の透磁率
I = 1.0 # 電流 (1A と仮定)
Bx = np.zeros_like(x_coords)
By = np.zeros_like(y_coords)

for i, (x, y) in enumerate(zip(x_coords, y_coords)):
r = math.sqrt(x2 + y2)
if r > 0: # 自身の位置での無限大回避
B = mu_0 * I / (2 * math.pi * r) # 磁場の大きさ (円形電流の場合)
Bx[i] = -B * y / r # x成分
By[i] = B * x / r # y成分

結果を VTK ファイルとして保存 (ParaView 用)

output_filename = “magnetic_field.vtk”
with open(output_filename, “w”) as vtk_file:
vtk_file.write(“# vtk DataFile Version 3.0\n”)
vtk_file.write(“Magnetic field example\n”)
vtk_file.write(“ASCII\n”)
vtk_file.write(“DATASET UNSTRUCTURED_GRID\n”)
vtk_file.write(f”POINTS {len(x_coords)} float\n”)
for x, y in zip(x_coords, y_coords):
vtk_file.write(f”{x} {y} 0\n”)

vtk_file.write(f"\nPOINT_DATA {len(x_coords)}\n")
vtk_file.write("VECTORS MagneticField float\n")
for bx, by in zip(Bx, By):
    vtk_file.write(f"{bx} {by} 0\n")

print(f”Magnetic field data saved to {output_filename}”)

終了

gmsh.finalize()

import gmsh
import numpy as np

Gmsh 初期化

gmsh.initialize()
gmsh.open(“line_based_structure.geo”)

メッシュを生成

gmsh.model.mesh.generate(2)

電気伝導率を定義

conductivity = {} # 領域ごとの電気伝導率を格納

領域の物理名を取得して電気伝導率を割り当て

entities = gmsh.model.getEntities(dim=2) # 2次元エンティティ (サーフェス)
for entity in entities:
name = gmsh.model.getPhysicalName(2, entity[1])
if name == “Air”:
conductivity[entity[1]] = 1e-12 # 空気の仮想的な電気伝導率
elif name == “Conductor”:
conductivity[entity[1]] = 5.96e7 # 銅の電気伝導率

設定された電気伝導率を確認

print(“Conductivity settings for regions:”)
for tag, value in conductivity.items():
print(f”Region {tag}: Conductivity = {value} S/m”)

ノード情報を取得

node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
node_coords = np.array(node_coords).reshape(-1, 3)

ノードごとの電流密度と電場を計算

node_current_density = []
node_electric_field = []
for coord in node_coords:
# 仮想的に電場を線形分布と仮定
# 下: 0V, 上: 10V
if coord[1] == 3:
E = 10 / (3 – (-2)) # 電場の大きさ
elif coord[1] == -2:
E = 10 / (3 – (-2))
else:
E = 10 / (3 – (-2))

# 電流密度 = 電気伝導率 * 電場
J = 1e-12 * E  # 空気層での仮定値（適宜変更可）
node_electric_field.append(E)
node_current_density.append(J)

VTK ファイルに電流密度と電場を保存

output_filename = “conductivity_simulation.vtk”
with open(output_filename, “w”) as vtk_file:
vtk_file.write(“# vtk DataFile Version 3.0\n”)
vtk_file.write(“Conductivity simulation\n”)
vtk_file.write(“ASCII\n”)
vtk_file.write(“DATASET UNSTRUCTURED_GRID\n”)
# ノード座標を書き出し
vtk_file.write(f”POINTS {len(node_coords)} float\n”)
for coord in node_coords:
vtk_file.write(f”{coord[0]} {coord[1]} {coord[2]}\n”)
vtk_file.write(“\n”)

# ノードごとの電場データを書き出し
vtk_file.write(f"POINT_DATA {len(node_coords)}\n")
vtk_file.write("SCALARS ElectricField float 1\n")
vtk_file.write("LOOKUP_TABLE default\n")
for E in node_electric_field:
    vtk_file.write(f"{E}\n")

# ノードごとの電流密度データを書き出し
vtk_file.write("SCALARS CurrentDensity float 1\n")
vtk_file.write("LOOKUP_TABLE default\n")
for J in node_current_density:
    vtk_file.write(f"{J}\n")

print(f”Simulation completed and saved to {output_filename}”)

Gmsh を終了

gmsh.finalize()

import gmsh
import numpy as np

Gmsh 初期化

gmsh.initialize()
gmsh.open(“line_based_structure.geo”)

メッシュを生成

gmsh.model.mesh.generate(2)

電位の計算のための準備

下部のライン（AirBottom）を 0V、上部のライン（AirTop）を 10V と設定

voltage = {} # 境界ごとの電圧を格納

境界ラインのエンティティを取得

entities = gmsh.model.getEntities(dim=1) # 1次元エンティティ (ライン)
for entity in entities:
name = gmsh.model.getPhysicalName(1, entity[1])
if name == “AirBottom”: # 空気層の下界面
voltage[entity[1]] = 0 # 0V
elif name == “AirTop”: # 空気層の上界面
voltage[entity[1]] = 10 # 10V

設定された電圧を確認

print(“Voltage settings for boundaries:”)
for tag, value in voltage.items():
print(f”Boundary {tag}: Voltage = {value} V”)

ノード情報を取得

node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
node_coords = np.array(node_coords).reshape(-1, 3)

電位ポテンシャルの計算

node_potentials = []
for coord in node_coords:
# 下: 0V、上: 10V
if coord[1] == 3: # y=3 は AirTop
node_potentials.append(10)
elif coord[1] == -2: # y=-2 は AirBottom
node_potentials.append(0)
else:
# 線形補間で中間の電位を計算
potential_at_node = 10 * (coord[1] – (-2)) / (3 – (-2))
node_potentials.append(potential_at_node)

VTK ファイルに電位ポテンシャルを保存

output_filename = “conductivity_simulation_with_potential.vtk”
with open(output_filename, “w”) as vtk_file:
vtk_file.write(“# vtk DataFile Version 3.0\n”)
vtk_file.write(“Conductivity simulation with potential\n”)
vtk_file.write(“ASCII\n”)
vtk_file.write(“DATASET UNSTRUCTURED_GRID\n”)
# ノード座標を書き出し
vtk_file.write(f”POINTS {len(node_coords)} float\n”)
for coord in node_coords:
vtk_file.write(f”{coord[0]} {coord[1]} {coord[2]}\n”)
vtk_file.write(“\n”)

# ノードごとの電位ポテンシャルデータを書き出し
vtk_file.write(f"POINT_DATA {len(node_coords)}\n")
vtk_file.write("SCALARS Potential float 1\n")
vtk_file.write("LOOKUP_TABLE default\n")
for potential in node_potentials:
    vtk_file.write(f"{potential}\n")

print(f”Voltage potential simulation completed and saved to {output_filename}”)

Gmsh を終了

gmsh.finalize()