Update

feat: 調整網格細化參數。
feat: 統一了 BJT 和 Opto 模型的程式風格網格細化設定。
2025-12-26 16:04:34 +08:00 · 2025-12-26 16:02:08 +08:00 · 2025-12-26 15:44:41 +08:00 · 2025-12-26 11:52:41 +08:00 · 2025-12-26 11:30:52 +08:00 · 2025-12-26 10:42:48 +08:00
43 changed files with 512 additions and 1426 deletions
@@ -1,446 +0,0 @@
-# BJT TCAD 模擬專案原理與流程解說
-
-本專案使用 **DEVSIM** 進行平面 NPN BJT 的元件級 TCAD 模擬。
-
->  **閱讀提示：** 本文件包含 Mermaid 圖表，建議使用以下方式閱讀：
-> ```bash
-> # 終端機閱讀
-> glow BJT_TCAD_GUIDE.md
-> 
-> # 瀏覽器閱讀 (支援 Mermaid 圖表)
-> pip install grip && grip BJT_TCAD_GUIDE.md
-> 
-> # VS Code 預覽 (需安裝 Mermaid 擴充套件)
-> # Ctrl+Shift+V (Windows/Linux) 或 Cmd+Shift+V (macOS)
-> ```
-
-> ***快速入門：** [README.md](README.md)
-
---
-
-## 一、物理基礎
-
-### 1.1 半導體基本方程
-
-DEVSIM 求解半導體器件的三個核心方程：
-
-```mermaid
-graph TD
-    Poisson[Poisson 方程] --> Psi[電位分佈 Ψ]
-    ElecCont[電子連續方程] --> n[電子濃度 n]
-    HoleCont[電洞連續方程] --> p[電洞濃度 p]
-    Psi --> Solve[自洽求解]
-    n --> Solve
-    p --> Solve
-```
-
-| 方程 | 物理意義 | 數學形式 |
-|------|----------|----------|
-| **Poisson** | 電荷分佈決定電場 | ∇²Ψ = -q(p - n + N_D - N_A)/ε |
-| **電子連續** | 電子流守恆 | ∂n/∂t = (1/q)∇·J_n + G - R |
-| **電洞連續** | 電洞流守恆 | ∂p/∂t = -(1/q)∇·J_p + G - R |
-
-### 1.2 載子傳輸模型
-
-專案使用 **Drift-Diffusion 模型**：
-
-```
-J_n = qμ_n·n·E + qD_n·∇n
-      ↑ 漂移       ↑ 擴散
-      (電場驅動)    (濃度梯度驅動)
-```
-
-| 參數 | 符號 | 物理意義 |
-|------|------|----------|
-| 電子遷移率 | μ_n | 電子在電場下的移動能力 |
-| 電洞遷移率 | μ_p | 電洞在電場下的移動能力 |
-| 擴散係數 | D_n, D_p | 愛因斯坦關係 D = μ·kT/q |
-
-### 1.3 PN 接面物理
-
-BJT 包含兩個 PN 接面：
-
-```
-      E-B 接面         B-C 接面
-         ↓                ↓
-    ┌────┬────┬────┬────┬────┐
-    │ N+ │    │ P  │    │ N  │  ← 正向偏壓時
-    │Emit│    │Base│    │Coll│    電子從 E→B→C 流動
-    └────┴────┴────┴────┴────┘
-           ↑         ↑
-        正偏        逆偏
-        (導通)      (收集)
-```
-
-**工作原理：**
-1. E-B 正偏 → 電子注入 Base
-2. Base 很薄 → 多數電子穿越到 Collector
-3. B-C 逆偏 → 電子被 Collector 收集
-4. 電流增益 β = I_C / I_B
-
---
-
-## 二、摻雜模型 (ERFC)
-
-### 2.1 ERFC 函數
-
-互補誤差函數 `erfc(x)` 模擬熱擴散摻雜分佈：
-
-```
-erfc(x) = 1 - erf(x) = (2/√π) ∫_x^∞ e^(-t²) dt
-
-         ↑
-    1.0  ├──●─────●
-         │     ╲
-    0.5  ├───────●─────
-         │         ╲
-    0.0  ├───────────●──●
-         └───┼───┼───┼───→ x
-            -2   0   2
-```
-
-### 2.2 摻雜公式
-
-```python
-# Emitter (N+) 摻雜
-emitter_doping
-* erfc((y - depth) / vdiff)           # 垂直方向：從表面衰減
-* erfc(-(x + 0.5*width - center) / hdiff)  # 左邊界
-* erfc((x - 0.5*width - center) / hdiff)   # 右邊界
-```
-
-| 參數 | 作用 |
-|------|------|
-| `depth` | 垂直 50% 濃度點 |
-| `vdiff` | 垂直擴散長度（控制深度過渡區寬度）|
-| `hdiff` | 水平擴散長度（控制側向過渡區寬度）|
-
---
-
-## 三、數值方法
-
-### 3.1 有限元素法 (FEM)
-
-將連續區域離散為三角形網格：
-
-```
-    ●───────●───────●
-   ╱ ╲     ╱ ╲     ╱ ╲
-  ╱   ╲   ╱   ╲   ╱   ╲
- ●─────●─●─────●─●─────●
-  ╲   ╱   ╲   ╱   ╲   ╱
-   ╲ ╱     ╲ ╱     ╲ ╱
-    ●───────●───────●
-```
-
-**優點：** 可自適應細化（PN 接面處網格更密）
-
-### 3.2 Newton-Raphson 迭代
-
-非線性方程組求解：
-
-```
-while |error| > tolerance:
-    J·Δx = -F(x)      # 求解線性系統
-    x = x + Δx         # 更新解
-```
-
-| 輸出 | 說明 |
-|------|------|
-| `RelError` | 相對誤差（收斂判據）|
-| `AbsError` | 絕對誤差 |
-| `Iteration` | 迭代次數 |
-
---
-
-## 四、模擬流程
-
-### 4.1 完整流程圖
-
-```mermaid
-graph TD
-    A[bjt.geo] -->|make mesh| B[bjt.msh]
-    B -->|make refine| C[bjt_refined.msh]
-    C -->|make init| D[bjt_dd_0.msh]
-    D -->|make gummel| E[gummel.csv]
-    D -->|make output| F[output_Vb0.70.csv]
-    D -->|make cb| G[cb_Vc0.50.csv]
-    D -->|make ac| H[ac_*.csv]
-```
-
-### 4.2 各步驟詳解
-
-#### Step 1: 網格生成 (`make mesh`)
-
-```bash
-gmsh -2 -format msh2 bjt.geo -o bjt.msh
-```
-
- **輸入:** `bjt.geo` (幾何定義)
- **輸出:** `bjt.msh` (初始網格)
- **原理:** Delaunay 三角化，~14000 節點
-
-#### Step 2: 自適應細化 (`make refine`)
-
-```bash
-python bjt_refine.py bjt.msh
-```
-
-**細化策略：**
-
-| 策略 | 模型 | 物理意義 |
-|------|------|----------|
-| E-field | `Emag` | 高電場區（PN 接面）細化 |
-| Contact | `SA` | 接觸面附近細化 |
-
- **輸出:** `bjt_refined.msh` (~38000 節點)
-
-#### Step 3: 初始化 (`make init`)
-
-```bash
-python bjt_dd.py
-```
-
-**求解順序：**
-1. **Potential Only** - 僅求解 Poisson 方程（建立初始電位）
-2. **Drift-Diffusion** - 耦合三方程求解（零偏壓平衡）
-
- **輸出:** `bjt_dd_0.msh` (DEVSIM 格式), `bjt_dd_0.tec` (ParaView 格式)
-
-#### Step 4: 電路模擬
-
-| 指令 | 腳本 | 掃描 | 輸出 |
-|------|------|------|------|
-| `make gummel` | bjt_circuit3.py | Vb: 0→0.8V | Ic, Ib vs Vb |
-| `make output` | bjt_circuit2.py | Vc: 0→2V (Vb=0.7V) | Ic vs Vc |
-| `make cb` | bjt_circuit4.py | Ve: 0→-1V | 共基極特性 |
-| `make ac` | bjt_circuit5.py | f: 1kHz→100GHz | AC 響應 |
-
---
-
-## 五、檔案結構與功能
-
-### 5.1 核心模組
-
-| 檔案 | 功能 | 說明 |
-|------|------|------|
-| `bjt_common.py` | **共用模組** | 摻雜參數、多執行緒設定、工具函數的單一真相來源 |
-| `bjt_device_setup.py` | 元件設定 | ERFC 摻雜分佈定義 |
-| `bjt_physics_model.py` | 物理模型 | DD 模型、接觸條件封裝 |
-
-### 5.2 模擬腳本
-
-| 檔案 | 功能 | 輸出 |
-|------|------|------|
-| `bjt_dd.py` | 零偏壓初始化 | `bjt_dd_0.msh`, `bjt_dd_0.tec` |
-| `bjt_refine.py` | 多策略網格細化 | `bjt_bg.pos` |
-| `bjt_circuit2.py` | 輸出特性 (Ic-Vc) | `output_*.csv` |
-| `bjt_circuit3.py` | Gummel 曲線 | `gummel.csv` |
-| `bjt_circuit4.py` | 共基極特性 | `cb_*.csv` |
-| `bjt_circuit5.py` | AC 小信號分析 | `ac_*.csv` |
-
-### 5.3 其他檔案
-
-| 檔案 | 功能 |
-|------|------|
-| `bjt.geo` | Gmsh 幾何定義 (45×14.5 μm²) |
-| `physics/` | 官方物理模組 (`new_physics.py`, `model_create.py`) |
-| `Makefile` | 建置自動化 |
-| `refinement_loop.sh` | 網格細化迴圈腳本 |
-| `sims.sh` | 批次模擬腳本 |
-
---
-
-## 5.4 多執行緒加速
-
-### 原理
-
-DEVSIM 支援平行運算，可利用多核心 CPU 加速模型計算：
-
-| 參數 | 說明 |
-|------|------|
-| `threads_available` | 可用的執行緒數量 |
-| `threads_task_size` | 最小任務大小（元素數 > 此值時才平行化）|
-
-### 使用方式
-
-本專案透過 `bjt_common.py` 自動啟用多執行緒：
-
-```bash
-# 預設使用全部 CPU 核心
-make all
-
-# 指定 4 個執行緒
-DEVSIM_THREADS=4 make all
-
-# 停用多執行緒 (除錯/效能比較)
-DEVSIM_THREADS=1 make all
-```
-
-### 動態調整
-
-在 Python 腳本中可隨時調整：
-
-```python
-import bjt_common
-
-# 檢查目前設定
-from devsim import get_parameter
-print(get_parameter(name="threads_available"))
-
-# 動態修改
-bjt_common.setup_threads(num_threads=4)
-bjt_common.setup_threads(num_threads=1)  # 關閉多執行緒
-```
-
---
-
-## 六、結果解讀
-
-### 6.1 摻雜分佈 (LogNetDoping)
-
-用 ParaView 檢視 `bjt_dd_0.tec`：
-
-| 顏色 | 值 | 區域 |
-|------|-----|------|
-| 紅 | +18~+19 | N+ Emitter |
-| 藍 | -16~-17 | P Base |
-| 橙 | +16~+18 | N Collector |
-
-### 6.2 Gummel 曲線
-
-```
-log(I)
-  ↑     ●
-        ●    ← Ic (指數增長)
-       ●
-      ●
-     ●   ●   ← Ib
-    ●   ●
-   ●   ●
-  ●───●───→ Vb
-  0   0.4  0.8
-```
-
-**電流增益 β = Ic/Ib ≈ 100~300**
-
---
-
-## 七、常見問題
-
-| 問題 | 原因 | 解決 |
-|------|------|------|
-| 收斂失敗 | 網格太粗/偏壓步進太大 | `make refine` / 減小步進 |
-| 電流過小 | 接觸電阻/摻雜不足 | 檢查濃度分佈 |
-| 結果不對稱 | 網格或幾何錯誤 | 用 Gmsh 檢視網格 |
-| 參數未定義錯誤 | circuit 檔案未同步 | 同步更新所有 DOPING_PARAMS |
-| 網格產生失敗 | .geo 語法錯誤 | `gmsh bjt.geo` 除錯 |
-
-### 自訂元件修改流程
-
-> ***重要：** 現行架構使用 `bjt_common.py` 作為摻雜參數的單一真相來源，修改參數只需編輯此檔案即可。
-
-| 步驟 | 檔案 | 必要性 | 說明 |
-|------|------|--------|------|
-| 1 | `bjt.geo` | ✓ 必須 | 定義幾何結構 |
-| 2 | `bjt_common.py` | ✓ 必須 | 修改 `DOPING_PARAMS` 字典 |
-| 3 | `bjt_device_setup.py` | ○ 進階 | ERFC 公式調整 (通常不需修改) |
-| 4 | `bjt_refine.py` | ○ 建議 | 調整細化策略 (若有特殊需求) |
-
-### 收斂性建議
-
-| 參數 | 建議值 | 說明 |
-|------|--------|------|
-| `hdiff`/`vdiff` | ≥ 1e-4 | 較大值 → 更平滑邊界 → 更穩定收斂 |
-| `emitter_doping` | ≤ 1e19 | 過高濃度梯度可能導致收斂困難 |
-
---
-
-## 八、與官方範例差異
-
-本專案 (`wisetop_bjt`) 與官方範例 (`devsim_bjt_example-main`) 的比較：
-
-### 8.1 元件結構
-
-| 項目 | 本專案 (Planar BJT) | 官方範例 (Vertical BJT) |
-|------|---------------------|-------------------------|
-| **尺寸** | 45 × 14.5 μm² | 27.5 × 5 μm² |
-| **Collector 位置** | 表面 (x=0-10μm) | 底部 (整個下方) |
-| **結構** | 平面結構 | 垂直結構 |
-| **Oxide** | 有 (隔離區) | 無 |
-
-```
-本專案 (Planar):              官方範例 (Vertical):
-┌──────────┬───────────┐       ┌───────────────────┐
-│Collector │ Oxide     │       │ Emitter │  Base   │
-├──────────┼───────────┤       ├─────────┴─────────┤
-│  N-Well  │   Base    │       │    Collector      │
-├──────────┴───────────┤       └───────────────────┘
-│      N-Sub           │
-└──────────────────────┘
-```
-
-### 8.2 摻雜模型
-
-| 參數 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **公式** | ERFC (已對齊) | ERFC |
-| **Emitter 濃度** | 1e19 | 1e19 |
-| **Base 濃度** | 1e17 | 1e17 |
-| **Collector** | N-Sub 1e16 + N-Well 1e18 | collector_doping 1e16 + sub_collector 1e19 |
-| **擴散長度** | hdiff/vdiff = 1e-4 | hdiff/vdiff = 1e-5 |
-
-> [!NOTE]
-> 本專案使用較大的擴散長度 (1e-4) 以確保在較大元件上的數值穩定性。
-
-### 8.3 網格細化
-
-| 策略 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **E-field** |  預設啟用 |  |
-| **Contact** |  預設啟用 |  |
-| **Potential gradient** | ○ 可選 |  |
-| **Doping gradient** | ○ 可選 |  |
-| **Emag 公式** | `(EField_x²+EField_y²)^0.5` | 同左 |
-| **策略組合** | `max(Enorm, SA)` | 同左 |
-
-### 8.4 專案結構
-
-| 項目 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **建置系統** | Makefile | Shell script (sims.sh) |
-| **共用模組** | `bjt_common.py` (單一真相來源) | 各檔重複定義 |
-| **摻雜設定** | 獨立模組 `bjt_device_setup.py` | 內嵌於 `netdoping.py` |
-| **物理模型** | 封裝於 `bjt_physics_model.py` | 直接使用 `bjt_common.py` |
-| **電路模擬** | 分檔 `bjt_circuit[2-5].py` | 單檔 `bjt.py` |
-
-### 8.5 物理模型
-
-| 模型 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **遷移率** | 官方 `new_physics` | 同左 |
-| **複合** | SRH + Auger | 同左 |
-| **帶隙窄化** | BGN 模型 | 同左 |
-| **接觸** | Ohmic (金屬) | 同左 |
-
-### 8.6 主要差異總結
-
-```mermaid
-graph LR
-    A[本專案] --> B[Planar BJT]
-    A --> C[Makefile 自動化]
-    A --> D[模組化設計]
-    A --> E[hdiff=1e-4]
-    
-    F[官方範例] --> G[Vertical BJT]
-    F --> H[Shell Script]
-    F --> I[單檔整合]
-    F --> J[hdiff=1e-5]
-```
-
-| 優勢 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **易用性** |  Makefile 一鍵執行 | 需手動執行多檔 |
-| **可維護** |  模組化清晰 | 邏輯集中 |
-| **學習曲線** | 結構複雜但說明完整 | 精簡易讀 |
-| **擴展性** |  易新增功能 | 需重構 |
@@ -21,7 +21,7 @@
 # -----------------------------------------------------------------------------
 # 工具設定
 # -----------------------------------------------------------------------------
-PYTHON    ?= python                               # Python 解譯器
+PYTHON    ?= python3                              # Python 解譯器
 GMSH      ?= gmsh                                 # Gmsh 路徑
 NPROC     ?= $(shell nproc 2>/dev/null || echo 4) # CPU 核心數 (自動偵測)
 GMSH_OPTS  = -nt $(NPROC)                         # Gmsh 多執行緒選項
@@ -38,7 +38,7 @@ BG_POS       = bjt_bg.pos
 # -----------------------------------------------------------------------------
 # 細化設定
 # -----------------------------------------------------------------------------
-LOOPS ?= 3
+LOOPS ?= 1

 # -----------------------------------------------------------------------------
 # 預設參數
@@ -1,155 +0,0 @@
-# wisetop_bjt 專案分析報告
-
-> 生成時間：2025-12-12
-
---
-
-## 一、模組依賴關係
-
-```mermaid
-graph TD
-    subgraph "模擬腳本"
-        C2[bjt_circuit2.py]
-        C3[bjt_circuit3.py]
-        C4[bjt_circuit4.py]
-        C5[bjt_circuit5.py]
-        DD[bjt_dd.py]
-        RF[bjt_refine.py]
-    end
-    
-    subgraph "核心模組"
-        BC[bjt_common.py]
-        DS[bjt_device_setup.py]
-        PM[bjt_physics_model.py]
-    end
-    
-    subgraph "物理模組"
-        NP[physics/new_physics.py]
-        MC[physics/model_create.py]
-    end
-    
-    C2 --> BC
-    C3 --> BC
-    C4 --> BC
-    C5 --> BC
-    DD --> BC
-    DD --> DS
-    DD --> PM
-    
-    RF --> DS
-    
-    DS --> PM
-    PM --> NP
-    PM --> MC
-    
-    BC -.->|重複定義| DS
-```
-
---
-
-## 二、發現的問題
-
-### 🔴 問題 1：DOPING_PARAMS 重複定義（高優先）
-
-| 檔案 | 狀態 |
-|------|------|
-| `bjt_common.py` 第 57 行 | ✅ 單一真相來源（應保留）|
-| `bjt_device_setup.py` 第 10 行 | ❌ 重複定義（應移除）|
-
-**現況：** 兩個檔案都定義了 `DOPING_PARAMS`，目前內容相同，但這違反了「單一真相來源」原則。
-
-**風險：** 未來修改時可能只改其中一處，導致不一致性錯誤。
-
-**建議修正：**
-```python
-# bjt_device_setup.py 應改為
-from bjt_common import DOPING_PARAMS
-```
-
---
-
-### 🟡 問題 2：restore_doping_params() 未被使用
-
-`bjt_common.py` 定義了 `restore_doping_params()` 函數，但所有 circuit 腳本都直接使用：
-
-```python
-for k, v in DOPING_PARAMS.items():
-    set_parameter(device=device, region=region, name=k, value=v)
-```
-
-**建議：** 統一使用 `restore_doping_params()` 或移除此函數。
-
---
-
-### 🟢 正常：多執行緒自動啟用
-
-`bjt_common.py` 在 import 時自動呼叫 `setup_threads()`，所有依賴它的腳本都會受益。
-
-```
-[DEVSIM] Multi-threading: 8 threads, task_size=1024
-```
-
---
-
-## 三、與官方範例比較
-
-### 3.1 架構差異
-
-| 項目 | 本專案 (wisetop_bjt) | 官方範例 (devsim_bjt_example) |
-|------|----------------------|------------------------------|
-| **元件結構** | Planar BJT (水平) | Vertical BJT (垂直) |
-| **尺寸** | 45 × 14.5 μm² | 27.5 × 5 μm² |
-| **建置系統** | Makefile | Shell script |
-| **摻雜管理** | `DOPING_PARAMS` 字典 | `netdoping.py` 模組 |
-| **多執行緒** | ✅ 自動啟用 | ❌ 未設定 |
-
-### 3.2 摻雜參數比較
-
-| 參數 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **emitter_doping** | 1e19 | 1e19 ✅ |
-| **base_doping** | 1e17 | 1e17 ✅ |
-| **hdiff/vdiff** | 1e-4 (1μm) | 1e-5 (0.1μm) ⚠️ |
-| **Collector** | N-Sub + N-Well | collector_doping + sub_collector |
-
-> ⚠️ 本專案使用較大的擴散長度（1μm vs 0.1μm），這是為了在較大元件上確保數值穩定性。
-
-### 3.3 物理模型比較
-
-| 模型 | 本專案 | 官方範例 |
-|------|--------|----------|
-| **遷移率** | Arora + HF | ✅ 同左 |
-| **複合** | SRH | ✅ 同左 |
-| **BGN** | 帶隙窄化 | ✅ 同左 |
-| **接觸** | Ohmic | ✅ 同左 |
-
---
-
-## 四、模組功能摘要
-
-| 模組 | 功能 | 依賴 |
-|------|------|------|
-| `bjt_common.py` | 多執行緒、DOPING_PARAMS、工具函數 | devsim |
-| `bjt_device_setup.py` | 載入網格、ERFC 摻雜設定 | bjt_physics_model |
-| `bjt_physics_model.py` | DD 模型封裝、接觸條件 | physics/ |
-| `bjt_dd.py` | 零偏壓初始化 | bjt_common, bjt_device_setup, bjt_physics_model |
-| `bjt_refine.py` | 網格細化策略 | bjt_device_setup |
-| `bjt_circuit[2-5].py` | DC/AC 電路模擬 | bjt_common, physics/ |
-
---
-
-## 五、建議修正項目
-
-| 優先度 | 項目 | 說明 |
-|--------|------|------|
-| 🔴 高 | 移除 bjt_device_setup.py 中的 DOPING_PARAMS | 改為 `from bjt_common import DOPING_PARAMS` |
-| 🟡 中 | 統一使用 restore_doping_params() | 或移除此函數 |
-| 🟢 低 | 考慮減小 hdiff/vdiff | 如需更銳利的接面可改為 1e-5 |
-
---
-
-## 六、結論
-
-本專案架構設計合理，採用模組化設計將摻雜、物理模型、電路模擬分離。主要問題是 **DOPING_PARAMS 重複定義**，建議移除 `bjt_device_setup.py` 中的副本以維持單一真相來源。
-
-與官方範例的主要差異在於元件結構（Planar vs Vertical）和擴散長度設定，這些是設計選擇而非錯誤。
@@ -2,18 +2,6 @@

 DEVSIM 平面 NPN BJT 元件級模擬專案，支援完整的 DC/AC 電氣特性分析。

-> **閱讀提示：** 本文件使用 Markdown 格式，建議使用以下方式閱讀以獲得最佳體驗：
-> ```bash
-> # 終端機閱讀 (需安裝 glow)
-> glow README.md
-> 
-> # 瀏覽器閱讀 (需安裝 grip)
-> pip install grip && grip README.md
-> 
-> # VS Code 預覽
-> # 開啟檔案後按 Ctrl+Shift+V (或 Cmd+Shift+V)
-> ```
-
 ---

 ## 專案架構
@@ -32,7 +20,8 @@ wisetop_bjt/
 ├── bjt_circuit5.py          # AC 小信號分析
 ├── physics/                 # 官方物理模組
 │   ├── new_physics.py
-│   └── model_create.py
+│   ├── model_create.py
+│   └── ramp2.py
 ├── refinement_loop.sh       # 網格細化腳本
 ├── sims.sh                  # 批次模擬腳本
 └── Makefile                 # 建置自動化
@@ -40,7 +29,7 @@ wisetop_bjt/

 ---

-##  元件結構
+## 元件結構

 ```
 y (μm)  x=0      x=10     x=20     x=30     x=40    x=45
@@ -55,7 +44,7 @@ y (μm)  x=0      x=10     x=20     x=30     x=40    x=45

 ### 摻雜參數 (ERFC 模型)

-摻雜參數定義於 `bjt_common.py`，為專案的單一真相來源：
+摻雜參數定義於 `bjt_common.py`，為專案的**單一真相來源**：

 | 區域 | X 範圍 | Y 深度 | 濃度 (cm⁻³) | 類型 |
 |------|--------|--------|-------------|------|
@@ -64,8 +53,6 @@ y (μm)  x=0      x=10     x=20     x=30     x=40    x=45
 | **N-Well** | 0-10 μm | 0-5.5 μm | 1×10¹⁸ | N |
 | **N-Sub** | 全寬 | 全深 | 1×10¹⁶ | N- |

-> ***詳細說明：** [技術手冊](BJT_TCAD_GUIDE.md) — 物理原理、參數修改指南
-
 ---

 ## 快速開始
@@ -81,38 +68,24 @@ y (μm)  x=0      x=10     x=20     x=30     x=40    x=45
 ### 安裝步驟

 ```bash
-# 1. 下載專案
-git clone <repository-url>
-cd devsim
-
-# 2. 安裝系統套件
+# 1. 安裝系統套件
 sudo apt install python3 python3-venv gmsh    # Ubuntu/Debian
-# brew install python gmsh                     # macOS

-# 3. 建立並啟用虛擬環境
+# 2. 建立並啟用虛擬環境
 python3 -m venv denv
 source denv/bin/activate

-# 4. 安裝 Python 套件
+# 3. 安裝 Python 套件
 pip install devsim numpy
-
-# 5. 驗證安裝
-python -c "import devsim; print('DEVSIM OK')"
 ```

 ### 執行模擬

 ```bash
 cd devsim/wisetop_bjt
-make all    # 執行完整模擬流程 (已啟用多執行緒加速)
+make all    # 執行完整模擬流程
 ```

-> **效能提示：** 專案預設使用系統全部 CPU 核心進行平行運算，可透過環境變數調整：
-> ```bash
-> DEVSIM_THREADS=4 make all   # 指定使用 4 核心
-> DEVSIM_THREADS=1 make all   # 停用多執行緒 (除錯用)
-> ```
-
 ---

 ## Make 指令
@@ -122,7 +95,7 @@ make all    # 執行完整模擬流程 (已啟用多執行緒加速)
 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
 | `make all` | 執行完整模擬流程 | 所有 CSV 檔案 |
-| `make batch` | 批次參數掃描 | `data/` 目錄下的結果 |
+| `make quick` | 快速模擬 (無網格細化) | — |
 | `make clean` | 清除所有生成檔案 | — |

 ### 個別步驟
@@ -130,7 +103,7 @@ make all    # 執行完整模擬流程 (已啟用多執行緒加速)
 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
 | `make mesh` | 產生初始網格 | `bjt.msh` |
-| `make refine` | 網格細化 (3 次迭代) | `bjt_refined.msh` |
+| `make refine` | 網格細化 (2 次迭代) | `bjt_refined.msh` |
 | `make init` | 零偏壓初始化 | `bjt_dd_0.msh`, `bjt_dd_0.tec` |
 | `make gummel` | Gummel 曲線 | `gummel.csv` |
 | `make output` | 輸出特性 (Vb=0.7V) | `output_Vb0.70.csv` |
@@ -144,20 +117,41 @@ make all    # 執行完整模擬流程 (已啟用多執行緒加速)
 ```
 ┌─────────┐    ┌─────────┐    ┌─────────┐    ┌───────────────┐
 │  mesh   │───▶│ refine  │───▶│  init   │───▶│ gummel/output │
-│ (Gmsh)  │    │ (3迴圈) │    │ (DC=0V) │    │   /cb/ac      │
+│ (Gmsh)  │    │ (2迴圈) │    │ (DC=0V) │    │   /cb/ac      │
 └─────────┘    └─────────┘    └─────────┘    └───────────────┘
   bjt.msh    bjt_refined.msh  bjt_dd_0.msh      *.csv
 ```

 ---

+## 物理模型
+
+### 半導體基本方程
+
+| 方程 | 物理意義 |
+|------|----------|
+| **Poisson** | 電荷分佈決定電場：∇²Ψ = -q(p - n + N_D - N_A)/ε |
+| **電子連續** | 電子流守恆：∂n/∂t = (1/q)∇·J_n + G - R |
+| **電洞連續** | 電洞流守恆：∂p/∂t = -(1/q)∇·J_p + G - R |
+
+### 採用的模型
+
+| 模型 | 說明 |
+|------|------|
+| **遷移率** | Arora + High-Field |
+| **複合** | SRH + Auger |
+| **帶隙窄化** | BGN 模型 |
+| **接觸** | Ohmic (金屬) |
+
+---
+
 ## 結果檢視

 ```bash
 # 網格視覺化
 gmsh bjt_refined.msh

-# 摻雜/電位分佈 (需安裝 ParaView)
+# 摻雜/電位分佈
 paraview bjt_dd_0.tec
 # 載入後選擇 Variable: LogNetDoping 或 Potential
 ```
@@ -168,8 +162,6 @@ paraview bjt_dd_0.tec

 ### 多執行緒設定

-專案預設使用系統全部 CPU 核心進行平行運算（通過 `bjt_common.py` 自動啟用）：
-
 ```bash
 # 預設使用全部 CPU 核心
 make all
@@ -177,33 +169,10 @@ make all
 # 指定執行緒數量
 DEVSIM_THREADS=4 make all

-# 停用多執行緒 (除錯或效能比較用)
+# 停用多執行緒
 DEVSIM_THREADS=1 make all
 ```

-在 Python 腳本中也可動態調整：
-
-```python
-import bjt_common
-bjt_common.setup_threads(num_threads=4)  # 使用 4 核心
-bjt_common.setup_threads(num_threads=1)  # 關閉多執行緒
-```
-
-### 指定 Python 解譯器
-
-本專案使用**動態路徑計算**，無需修改硬編碼路徑：
-
-```bash
-# 方法 1：啟用虛擬環境 (推薦)
-source denv/bin/activate && make all
-
-# 方法 2：環境變數
-export PYTHON=/path/to/python && make all
-
-# 方法 3：命令列參數
-make PYTHON=/path/to/python all
-```
-
 ### 修改摻雜參數

 編輯 `bjt_common.py` 中的 `DOPING_PARAMS` 字典：
@@ -224,10 +193,8 @@ DOPING_PARAMS = {
 | 收斂失敗 | `make clean && make refine && make init` |
 | 模組找不到 | 確認已啟用虛擬環境：`source denv/bin/activate` |
 | `devsim` 未安裝 | `pip install devsim` |
-| Gmsh 找不到 | `sudo apt install gmsh` 或 `brew install gmsh` |
+| Gmsh 找不到 | `sudo apt install gmsh` |
 | 摻雜分佈異常 | 用 ParaView 檢視 `bjt_dd_0.tec` 的 `LogNetDoping` |
-| Shell 腳本權限不足 | `chmod +x *.sh` |
-| 模擬時間過長 | 預設已啟用多執行緒，可用 `DEVSIM_THREADS=N` 調整核心數 |

 ---

@@ -235,6 +202,7 @@ DOPING_PARAMS = {

 | 項目 | 說明 |
 |------|------|
-| **更新日期** | 2025-12-12 |
+| **更新日期** | 2025-12-24 |
+| **Python** | python3 |
+| **網格細化** | LOOPS=2 |
 | **架構特點** | 共用模組化設計 (`bjt_common.py`) |
-| **路徑相容性** | 動態計算，跨平台可攜 |
@@ -6,7 +6,8 @@

 Mesh.CharacteristicLengthExtendFromBoundary = 0;
 Mesh.Algorithm = 5;
-Mesh.CharacteristicLengthMax = 2.5e-5;
+Mesh.CharacteristicLengthMin = 2.0e-6;   // 2 µm (for refinement)
+Mesh.CharacteristicLengthMax = 2.5e-5;   // 25 µm (initial mesh)

 cl = 2.5e-5;
 sf = 1.0e-4;
@@ -1,10 +1,13 @@
+#!/usr/bin/env python3
 """
 bjt_circuit2.py - Output Characteristics (Ic-Vc sweep at fixed Vb)
 """
 import sys
 import os
 import csv
-from devsim import *
+from devsim import (
+    load_devices, set_parameter, get_parameter, solve
+)

 from bjt_common import DOPING_PARAMS, get_contact_currents
 from physics.new_physics import GetContactBiasName, SetSiliconParameters
@@ -1,10 +1,13 @@
+#!/usr/bin/env python3
 """
 bjt_circuit3.py - Gummel Plot (Ic, Ib vs Vb at fixed Vc)
 """
 import sys
 import os
 import csv
-from devsim import *
+from devsim import (
+    load_devices, set_parameter, get_parameter, solve
+)

 from bjt_common import DOPING_PARAMS, get_contact_currents
 from physics.new_physics import GetContactBiasName, SetSiliconParameters
@@ -1,10 +1,13 @@
+#!/usr/bin/env python3
 """
 bjt_circuit4.py - Common-Base Characteristics (sweep Ve at fixed Vc)
 """
 import sys
 import os
 import csv
-from devsim import *
+from devsim import (
+    load_devices, set_parameter, get_parameter, solve
+)

 from bjt_common import DOPING_PARAMS, get_contact_currents
 from physics.new_physics import GetContactBiasName, SetSiliconParameters
@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 """
 bjt_circuit5.py - AC Small-Signal Analysis (frequency sweep)
 """
@@ -5,7 +6,10 @@ import sys
 import os
 import csv
 import math
-from devsim import *
+from devsim import (
+    load_devices, set_parameter, solve,
+    circuit_element, circuit_alter, get_circuit_node_value
+)

 from bjt_common import DOPING_PARAMS
 from physics.new_physics import GetContactBiasName, SetSiliconParameters, CreateSiliconDriftDiffusionContact
@@ -1,33 +1,24 @@
 """
 bjt_common.py - Shared module for BJT simulation.
-
-Features:
-    - Multi-threading: auto-enabled via setup_threads()
-    - Doping params: DOPING_PARAMS dict (single source of truth)
-    - Utilities: get_contact_currents(), restore_doping_params()
-
-Environment:
-    DEVSIM_THREADS: Override thread count (default: all CPUs)
 """
 import sys
 import os
 from devsim import get_contact_current, set_parameter


-# --- Multi-threading ---
 def setup_threads(num_threads=None, task_size=1024):
-    """Configure DEVSIM threading. Uses all CPUs if num_threads is None."""
+    """Configure DEVSIM threading."""
    if num_threads is None:
        env_threads = os.environ.get('DEVSIM_THREADS')
        num_threads = int(env_threads) if env_threads else (os.cpu_count() or 1)
    set_parameter(name="threads_available", value=num_threads)
    set_parameter(name="threads_task_size", value=task_size)
-    print(f"[DEVSIM] Multi-threading: {num_threads} threads, task_size={task_size}")
+    print(f"[DEVSIM] Threads: {num_threads}, task_size={task_size}")

-setup_threads()  # Auto-enable on import
+setup_threads()


-# --- Path setup ---
+# --- Path Setup ---
 SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
 PROJECT_ROOT = os.path.dirname(SCRIPT_DIR)
 if PROJECT_ROOT not in sys.path:
@@ -37,7 +28,7 @@ sys.path.insert(0, os.path.join(os.getcwd(), 'physics'))


 def _apply_model_create_patch():
-    """Patch to avoid TypeError in EnsureEdgeFromNodeModelExists."""
+    """Patch EnsureEdgeFromNodeModelExists to avoid TypeError."""
    try:
        import python_packages.model_create as mc
        _orig = mc.EnsureEdgeFromNodeModelExists
@@ -53,49 +44,46 @@ def _apply_model_create_patch():
 _apply_model_create_patch()


-# Doping parameters (ERFC format) - single source of truth
+# --- Doping Parameters (ERFC, unit: cm) ---
 DOPING_PARAMS = {
    'H_silicon': 14.5e-4,
    # Emitter (N+)
-    'emitter_doping': 1e19,
-    'emitter_center': 25.0e-4,
-    'emitter_width': 10.0e-4,
-    'emitter_depth': 1.0e-4,
-    'emitter_hdiff': 1.0e-4,
-    'emitter_vdiff': 1.0e-4,
+    'emitter_doping': 1e19, 'emitter_center': 25.0e-4, 'emitter_width': 10.0e-4,
+    'emitter_depth': 1.0e-4, 'emitter_hdiff': 1.0e-4, 'emitter_vdiff': 1.0e-4,
    # Base (P)
-    'base_doping': 1e17,
-    'base_center': 27.5e-4,
-    'base_width': 35.0e-4,
-    'base_depth': 5.5e-4,
-    'base_hdiff': 1.0e-4,
-    'base_vdiff': 1.0e-4,
+    'base_doping': 1e17, 'base_center': 27.5e-4, 'base_width': 35.0e-4,
+    'base_depth': 5.5e-4, 'base_hdiff': 1.0e-4, 'base_vdiff': 1.0e-4,
    # N-Well
-    'nwell_doping': 1e18,
-    'nwell_center': 5.0e-4,
-    'nwell_width': 10.0e-4,
-    'nwell_depth': 5.5e-4,
-    'nwell_hdiff': 1.0e-4,
-    'nwell_vdiff': 1.0e-4,
+    'nwell_doping': 1e18, 'nwell_center': 5.0e-4, 'nwell_width': 10.0e-4,
+    'nwell_depth': 5.5e-4, 'nwell_hdiff': 1.0e-4, 'nwell_vdiff': 1.0e-4,
    # N-Sub
    'nsub_doping': 1e16,
 }


 def get_contact_currents(device):
-    """Get terminal currents (Ic, Ib, Ie) for BJT device."""
+    """Get BJT terminal currents (Ic, Ib, Ie)."""
    currents = {}
    for contact in ["collector", "base", "emitter"]:
-        i_n = get_contact_current(device=device, contact=contact,
-                                  equation="ElectronContinuityEquation")
-        i_p = get_contact_current(device=device, contact=contact,
-                                  equation="HoleContinuityEquation")
+        i_n = get_contact_current(device=device, contact=contact, equation="ElectronContinuityEquation")
+        i_p = get_contact_current(device=device, contact=contact, equation="HoleContinuityEquation")
        currents[contact] = i_n + i_p
    return currents


 def restore_doping_params(device, region="Silicon"):
    """Set DOPING_PARAMS to device region."""
-    from devsim import set_parameter
    for key, value in DOPING_PARAMS.items():
        set_parameter(device=device, region=region, name=key, value=value)
+
+
+# --- Refinement Parameters ---
+MESH_SIZE_MIN = 2.0e-6   # 2 um
+MESH_SIZE_MAX = 1.0e-4   # 100 um
+
+# E-field thresholds
+E_VERY_HIGH = 1.0e5   # 100 kV/cm -> 1x
+E_HIGH = 1.0e4        # 10 kV/cm  -> 2x
+E_MEDIUM = 1.0e3      # 1 kV/cm   -> 4x
+E_LOW = 1.0e2         # 100 V/cm  -> 8x
+E_VERYLOW = 1.0e1     # 10 V/cm   -> 16x
@@ -1,11 +1,12 @@
+#!/usr/bin/env python3
 """
-bjt_dd.py - BJT Zero Bias Initialization (equilibrium solve)
+bjt_dd.py - BJT Zero Bias Initialization
 """
 import sys
 import os
-from devsim import *
+from devsim import solve, node_solution, edge_from_node_model, set_node_values, write_devices

-import bjt_common  # Auto-enables multi-threading
+import bjt_common
 import bjt_device_setup
 import bjt_physics_model

@@ -15,15 +16,14 @@ def run_init():
    mesh_file = "bjt_refined.msh" if os.path.exists("bjt_refined.msh") else "bjt.msh"
    device = "bjt"
    
-    # Load mesh and setup doping
    si_regions, _ = bjt_device_setup.setup_device_and_doping(device, mesh_file)

-    # Step 1: Potential only solve
-    print("--- 1. Potential Only Solve ---")
+    # Potential only solve
+    print("--- 1. Potential Solve ---")
    solve(type="dc", absolute_error=1.0, relative_error=1e-9, maximum_iterations=50)

-    # Step 2: Drift-diffusion initialization
-    print("--- 2. Drift Diffusion Setup ---")
+    # Drift-diffusion setup
+    print("--- 2. DD Setup ---")
    for region in si_regions:
        node_solution(device=device, region=region, name="Electrons")
        node_solution(device=device, region=region, name="Holes")
@@ -36,16 +36,15 @@ def run_init():
    contact_map = {"collector": "Silicon", "emitter": "Silicon", "base": "Silicon"}
    bjt_physics_model.setup_contacts_drift_diffusion(device, contact_map, dd_opts)

-    # Step 3: Zero bias DC solve (relaxed for ERFC model)
-    print("--- 3. Drift Diffusion (Zero Bias) ---")
+    # Zero bias solve
+    print("--- 3. DD Solve ---")
    for i in range(5):
        try:
            solve(type="dc", absolute_error=1e10, relative_error=1.0, maximum_iterations=100)
-            print(f"  DD Iteration {i+1} converged")
+            print(f"  Iteration {i+1} OK")
        except:
-            print(f"  DD Iteration {i+1} failed, continuing...")
+            print(f"  Iteration {i+1} failed")

-    # Output results
    write_devices(file="bjt_dd_0.msh", type="devsim")
    write_devices(file="bjt_dd_0.tec", type="tecplot")
    print("Done: bjt_dd_0.msh")
@@ -1,26 +1,15 @@
+#!/usr/bin/env python3
 """
 bjt_device_setup.py - Device and ERFC Doping Setup for Planar BJT
 """
 import os
-from devsim import *
+from devsim import (
+    create_gmsh_mesh, add_gmsh_region, add_gmsh_contact,
+    finalize_mesh, create_device, get_device_list,
+    load_devices, set_parameter, node_model
+)
 import bjt_physics_model
-
-# --- Doping Parameters (ERFC, unit: cm) ---
-# Layout: Collector(0-10), Oxide(10-20), Emitter(20-30), Oxide(30-40), Base(40-45)
-DOPING_PARAMS = {
-    'H_silicon': 14.5e-4,
-    # Emitter (N+)
-    'emitter_doping': 1e19, 'emitter_center': 25.0e-4, 'emitter_width': 10.0e-4,
-    'emitter_depth': 1.0e-4, 'emitter_hdiff': 1.0e-4, 'emitter_vdiff': 1.0e-4,
-    # Base (P)
-    'base_doping': 1e17, 'base_center': 27.5e-4, 'base_width': 35.0e-4,
-    'base_depth': 5.5e-4, 'base_hdiff': 1.0e-4, 'base_vdiff': 1.0e-4,
-    # N-Well
-    'nwell_doping': 1e18, 'nwell_center': 5.0e-4, 'nwell_width': 10.0e-4,
-    'nwell_depth': 5.5e-4, 'nwell_hdiff': 1.0e-4, 'nwell_vdiff': 1.0e-4,
-    # N-Sub
-    'nsub_doping': 1e16,
-}
+import bjt_common


 def setup_device_and_doping(device, mesh_file):
@@ -47,42 +36,36 @@ def setup_device_and_doping(device, mesh_file):
        finalize_mesh(mesh=device)
        create_device(mesh=device, device=device)

-    for key, value in DOPING_PARAMS.items():
+    # Set doping parameters
+    for key, value in bjt_common.DOPING_PARAMS.items():
        set_parameter(device=device, region="Silicon", name=key, value=value)

-    # === ERFC Doping Profiles (Official Style) ===
-    
-    # Emitter (N+): ERFC profile
+    # ERFC doping profiles
    node_model(device=device, region="Silicon", name="nD_emit", equation='''
-        emitter_doping
-        * erfc((y - emitter_depth) / emitter_vdiff)
+        emitter_doping * erfc((y - emitter_depth) / emitter_vdiff)
        * erfc(-(x + 0.5 * emitter_width - emitter_center) / emitter_hdiff)
        * erfc((x - 0.5 * emitter_width - emitter_center) / emitter_hdiff)
    ''')
    
-    # Base (P): ERFC profile
    node_model(device=device, region="Silicon", name="nA_base", equation='''
-        base_doping
-        * erfc((y - base_depth) / base_vdiff)
+        base_doping * erfc((y - base_depth) / base_vdiff)
        * erfc(-(x + 0.5 * base_width - base_center) / base_hdiff)
        * erfc((x - 0.5 * base_width - base_center) / base_hdiff)
    ''')
    
-    # N-Well (under collector): ERFC profile
    node_model(device=device, region="Silicon", name="nD_nwell", equation='''
-        nwell_doping
-        * erfc((y - nwell_depth) / nwell_vdiff)
+        nwell_doping * erfc((y - nwell_depth) / nwell_vdiff)
        * erfc(-(x + 0.5 * nwell_width - nwell_center) / nwell_hdiff)
        * erfc((x - 0.5 * nwell_width - nwell_center) / nwell_hdiff)
    ''')
    
-    # Combine doping: N-Sub background + N-Well + Emitter
+    # Combine doping
    node_model(device=device, region="Silicon", name="Donors", equation="nsub_doping + nD_emit + nD_nwell")
    node_model(device=device, region="Silicon", name="Acceptors", equation="1e10 + nA_base")
    node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors")
    node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping/2)/log(10)")

-    # Initialize potential model
+    # Initialize potential
    bjt_physics_model.setup_silicon_potential_only(device, si_regions)
    contact_map = {"collector": "Silicon", "emitter": "Silicon", "base": "Silicon"}
    bjt_physics_model.setup_contacts_potential_only(device, contact_map)
@@ -1,11 +1,15 @@
+#!/usr/bin/env python3
 """
 bjt_physics_model.py - Physics Model Wrapper (potential/drift-diffusion)
 """
 import sys
 import os
-from devsim import *
+from devsim import (
+    set_parameter, get_parameter, node_solution,
+    edge_from_node_model, edge_model, equation
+)

-# --- Load Physics Module ---
+# Load physics module
 physics_path = os.path.join(os.path.dirname(__file__), 'physics')
 if physics_path not in sys.path:
    sys.path.insert(0, physics_path)
@@ -22,7 +26,7 @@ try:
            CreateEField, CreateDField, GetContactBiasName
        )
        from physics.model_create import CreateSolution
-        print("Physics: Official (new_physics)")
+        print("Physics: Official")
    else:
        raise ImportError("Fallback needed")
 except ImportError as e:
@@ -36,7 +40,7 @@ except ImportError as e:
    GetContactBiasName = simple_physics.GetContactBiasName
    CreateSiliconDriftDiffusionContact = getattr(simple_physics, 'CreateSiliconDriftDiffusionContact', None)

-# Oxide physics model (fallback definition)
+# Oxide physics (fallback)
 try:
    import python_packages.simple_physics as sp_oxide
    SetOxideParameters = sp_oxide.SetOxideParameters
@@ -56,7 +60,7 @@ except ImportError:


 def setup_silicon_potential_only(device, regions):
-    """Setup potential equation for silicon regions"""
+    """Setup potential equation for silicon regions."""
    for region in regions:
        SetSiliconParameters(device, region)
        if USE_OFFICIAL:
@@ -65,12 +69,7 @@ def setup_silicon_potential_only(device, regions):


 def setup_silicon_drift_diffusion(device, regions):
-    """
-    Setup drift-diffusion equations for silicon regions
-
-    Returns:
-        dict: Contains Jn, Jp current model names
-    """
+    """Setup drift-diffusion equations for silicon regions."""
    opts = {'Jn': 'Jn', 'Jp': 'Jp'}
    for region in regions:
        if USE_OFFICIAL:
@@ -85,14 +84,14 @@ def setup_silicon_drift_diffusion(device, regions):


 def setup_contacts_potential_only(device, contact_map):
-    """Setup potential boundary conditions for contacts"""
+    """Setup potential boundary conditions for contacts."""
    for contact, region in contact_map.items():
        set_parameter(device=device, name=GetContactBiasName(contact), value=0.0)
        CreateSiliconPotentialOnlyContact(device, region, contact)


 def setup_contacts_drift_diffusion(device, contact_map, opts=None):
-    """Setup drift-diffusion boundary conditions for contacts"""
+    """Setup drift-diffusion boundary conditions for contacts."""
    if opts is None:
        opts = {'Jn': 'Jn', 'Jp': 'Jp'}
    for contact, region in contact_map.items():
@@ -1,73 +1,67 @@
+#!/usr/bin/env python3
 """
-bjt_refine.py - Mesh Refinement (E-field and contact-based strategies)
+bjt_refine.py - Mesh Refinement based on E-field distribution
+
+Generates Gmsh background field (.pos) for adaptive mesh refinement.
 """
 import sys
-from devsim import *
+from devsim import (
+    solve, get_node_model_values, get_element_model_values,
+    element_from_edge_model, element_model, element_from_node_model,
+    edge_from_node_model, edge_model
+)
 import bjt_device_setup
-
-# --- Refinement Parameters ---
-MESH_SIZE_MIN = 2.0e-6  # 20nm (junction)
-MESH_SIZE_MAX = 1.0e-4  # 1um (bulk)
-
-# Strategy toggles
-ENABLE_EFIELD = True
-ENABLE_CONTACT = True
-ENABLE_POTENTIAL = False
-ENABLE_DOPING = False
-POTENTIAL_DIFF = 0.025
-DOPING_DIFF = 1.0
+import bjt_common


-# =============================================================================
-# Refinement Strategy Functions (Official Style)
-# =============================================================================
+# --- Refinement Strategy Functions ---

 def emag_refinement(device, region):
-    """Refine high E-field regions."""
+    """Create E-field magnitude refinement model."""
    element_from_edge_model(edge_model="EField", device=device, region=region)
-    element_model(device=device, region=region, name="Emag", 
+    element_model(device=device, region=region, name="Emag",
                  equation="(EField_x^2 + EField_y^2)^(0.5)")
-    element_model(device=device, region=region, name="Enorm", equation='''
-        ifelse(Emag > 1.0e5, 1.0,
-          ifelse(Emag > 1.0e4, 2.0,
-            ifelse(Emag > 1.0e3, 4.0,
-              ifelse(Emag > 1.0e2, 8.0,
-                if(Emag > 1.0e1, 16)))))
+    element_model(device=device, region=region, name="Enorm", equation=f'''
+        ifelse(Emag > {bjt_common.E_VERY_HIGH}, 1.0,
+          ifelse(Emag > {bjt_common.E_HIGH}, 2.0,
+            ifelse(Emag > {bjt_common.E_MEDIUM}, 4.0,
+              ifelse(Emag > {bjt_common.E_LOW}, 8.0,
+                if(Emag > {bjt_common.E_VERYLOW}, 16)))))
    ''')
    return "Enorm"


 def contact_refinement(device, region):
-    """Refine near-contact regions (4x finer)."""
+    """Create contact proximity refinement model (4x finer)."""
    element_from_node_model(node_model="SurfaceArea", device=device, region=region)
-    element_model(device=device, region=region, name="SA", 
+    element_model(device=device, region=region, name="SA",
                  equation="if((SurfaceArea@en0 + SurfaceArea@en1 + SurfaceArea@en2) > 0.0, 4.0)")
    return "SA"


 def potential_refinement(device, region, pdiff):
-    """Refine regions with large potential gradient."""
+    """Create potential gradient refinement model."""
    element_from_node_model(node_model="Potential", device=device, region=region)
-    element_model(device=device, region=region, name="potential_norm", equation='''
-        (abs(Potential@en0-Potential@en1) > %s) ||
-        (abs(Potential@en0-Potential@en2) > %s) ||
-        (abs(Potential@en1-Potential@en2) > %s)
-    ''' % (pdiff, pdiff, pdiff))
+    element_model(device=device, region=region, name="potential_norm", equation=f'''
+        (abs(Potential@en0-Potential@en1) > {pdiff}) ||
+        (abs(Potential@en0-Potential@en2) > {pdiff}) ||
+        (abs(Potential@en1-Potential@en2) > {pdiff})
+    ''')
    return "potential_norm"


 def doping_refinement(device, region, ldiff):
-    """Refine regions with large doping gradient."""
+    """Create doping gradient refinement model."""
    element_from_node_model(node_model="LogNetDoping", device=device, region=region)
-    element_model(device=device, region=region, name="lognetdoping_norm", equation='''
-        (abs(LogNetDoping@en0-LogNetDoping@en1) > %s) ||
-        (abs(LogNetDoping@en0-LogNetDoping@en2) > %s) ||
-        (abs(LogNetDoping@en1-LogNetDoping@en2) > %s)
-    ''' % (ldiff, ldiff, ldiff))
+    element_model(device=device, region=region, name="lognetdoping_norm", equation=f'''
+        (abs(LogNetDoping@en0-LogNetDoping@en1) > {ldiff}) ||
+        (abs(LogNetDoping@en0-LogNetDoping@en2) > {ldiff}) ||
+        (abs(LogNetDoping@en1-LogNetDoping@en2) > {ldiff})
+    ''')
    return "lognetdoping_norm"


-# --- Main Refinement Logic ---
+# --- Gmsh Output ---

 def write_gmsh_pos(filename, device, region, x, y, node_cl):
    """Write Gmsh background field .pos file."""
@@ -88,70 +82,57 @@ def write_gmsh_pos(filename, device, region, x, y, node_cl):
        print('};', file=fh)


+# --- Main ---
+
 def run_refinement():
-    """Execute single iteration of mesh refinement (Official Style)"""
+    """Execute mesh refinement based on E-field distribution."""
    mesh_file = sys.argv[1] if len(sys.argv) > 1 else "bjt.msh"
    device = "bjt"

+    # Setup device and solve for potential
    si_regions, _ = bjt_device_setup.setup_device_and_doping(device, mesh_file)
    region = si_regions[0]

-    # Solve potential for E-field
-    print("--- Solving Potential (for refinement) ---")
+    print("--- Solving Potential ---")
    try:
        solve(type="dc", absolute_error=1e10, relative_error=1e-1, maximum_iterations=30)
    except:
        print("Warning: Solve failed, using initial guess")

+    # Get node coordinates
    x = get_node_model_values(device=device, region=region, name="x")
    y = get_node_model_values(device=device, region=region, name="y")
    num_nodes = len(x)

-    # === Setup E-field for refinement ===
+    # Setup E-field model
    edge_from_node_model(device=device, region=region, node_model="Potential")
    edge_model(device=device, region=region, name="EField",
               equation="(Potential@n0 - Potential@n1) * EdgeInverseLength")

-    # === Apply refinement strategies (Official Style) ===
+    # Apply refinement strategies
    strategies = []
    
-    if ENABLE_EFIELD:
-        name = emag_refinement(device, region)
-        strategies.append(name)
-        print(f"  [✓] E-field refinement enabled ({name})")
+    strategy = emag_refinement(device, region)
+    strategies.append(strategy)
+    print(f"  [OK] E-field refinement ({strategy})")
    
-    if ENABLE_CONTACT:
-        name = contact_refinement(device, region)
-        strategies.append(name)
-        print(f"  [✓] Contact refinement enabled ({name})")
-    
-    if ENABLE_POTENTIAL:
-        name = potential_refinement(device, region, POTENTIAL_DIFF)
-        strategies.append(name)
-        print(f"  [✓] Potential gradient refinement enabled ({name})")
-    
-    if ENABLE_DOPING:
-        name = doping_refinement(device, region, DOPING_DIFF)
-        strategies.append(name)
-        print(f"  [✓] Doping gradient refinement enabled ({name})")
-    
-    # === Combine strategies using max() (Official Style) ===
-    if len(strategies) == 0:
-        print("Warning: No refinement strategy enabled!")
-        clen_eq = "1.0"
-    elif len(strategies) == 1:
+    strategy = contact_refinement(device, region)
+    strategies.append(strategy)
+    print(f"  [OK] Contact refinement ({strategy})")
+
+    # Combine strategies using max()
+    if len(strategies) == 1:
        clen_eq = strategies[0]
    else:
-        # Build nested max() expression: max(A, max(B, max(C, ...)))
        clen_eq = strategies[-1]
        for s in reversed(strategies[:-1]):
            clen_eq = f"max({s}, {clen_eq})"
    
-    print(f"  Combined strategy: {clen_eq}")
+    print(f"  Combined: {clen_eq}")
    element_model(device=device, region=region, name="clen", equation=clen_eq)
    cl = get_element_model_values(device=device, region=region, name='clen')
    
-    # Map element values to nodes (Official Style)
+    # Map element values to nodes
    element_from_node_model(node_model="node_index", device=device, region=region)
    en0 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en0')]
    en1 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en1')]
@@ -162,21 +143,21 @@ def run_refinement():
        v = cl[i] if i < len(cl) else 0
        ni0, ni1, ni2 = en0[i], en1[i], en2[i]
        if v > 0:
-            target = MESH_SIZE_MIN * v
-            # Use minimum size for high-field regions
-            if node_cl[ni0] == 0.0 or target < node_cl[ni0]:
-                node_cl[ni0] = target
-            if node_cl[ni1] == 0.0 or target < node_cl[ni1]:
-                node_cl[ni1] = target
-            if node_cl[ni2] == 0.0 or target < node_cl[ni2]:
-                node_cl[ni2] = target
+            target = bjt_common.MESH_SIZE_MIN * v
+            for ni in [ni0, ni1, ni2]:
+                if node_cl[ni] == 0.0 or target < node_cl[ni]:
+                    node_cl[ni] = target
    
-    # Fill remaining nodes with max size
+    # Fill remaining with max size
    for i in range(num_nodes):
        if node_cl[i] == 0.0:
-            node_cl[i] = MESH_SIZE_MAX
+            node_cl[i] = bjt_common.MESH_SIZE_MAX
    
    write_gmsh_pos("bjt_bg.pos", device, region, x, y, node_cl)
+    
+    # Stats
+    refined_count = sum(1 for c in node_cl if c < bjt_common.MESH_SIZE_MAX)
+    print(f"  Refined nodes: {refined_count} / {num_nodes}")
    print("Done - bjt_bg.pos written")


@@ -49,12 +49,12 @@ RESULT_3DF      = opto_3df_out
 # -----------------------------------------------------------------------------
 # 細化設定
 # -----------------------------------------------------------------------------
-LOOPS ?= 2
+LOOPS ?= 3

 # -----------------------------------------------------------------------------
 # PHONY 目標宣告
 # -----------------------------------------------------------------------------
-.PHONY: help 2d 2df 3d 3df clean
+.PHONY: help 2d 3d clean

 # =============================================================================
 # HELP - 顯示可用指令
@@ -64,13 +64,9 @@ help:
 	@echo "  DEVSIM Optocoupler Simulation Makefile"
 	@echo "=============================================================="
 	@echo ""
-	@echo "  完整模擬 (含網格細化):"
-	@echo "    make 2d             執行 2D 完整模擬"
-	@echo "    make 3d             執行 3D 完整模擬"
-	@echo ""
-	@echo "  快速模擬 (無網格細化):"
-	@echo "    make 2df            執行 2D 快速模擬"
-	@echo "    make 3df            執行 3D 快速模擬"
+	@echo "  模擬指令:"
+	@echo "    make 2d             執行 2D 模擬 (含網格細化)"
+	@echo "    make 3d             執行 3D 模擬 (含網格細化)"
 	@echo ""
 	@echo "  其他:"
 	@echo "    make clean          清除所有生成檔案"
@@ -196,23 +192,6 @@ sim2d: $(REFINED_MESH_2D)
 	@echo "    結果檔案:   $(RESULT_2D).msh, $(RESULT_2D).tec"
 	@echo "=========================================="

-# -----------------------------------------------------------------------------
-# 2df: 2D 快速模擬 (無網格細化)
-# 流程: mesh2d → 直接模擬
-# 指令:
-#   1. gmsh -nt N -2 -format msh2 opto_simplified.geo -o opto_simplified.msh
-#   2. python3 opto_simplified_run.py opto_simplified.msh
-# -----------------------------------------------------------------------------
-2df: mesh2d
-	@echo ">>> [2df] 執行 2D 快速模擬 (無細化)..."
-	@echo "    指令: $(PYTHON) opto_simplified_run.py $(MESH_2D) $(RESULT_2DF)"
-	$(PYTHON) opto_simplified_run.py $(MESH_2D) $(RESULT_2DF)
-	@echo ""
-	@echo "=========================================="
-	@echo ">>> [2DF] 快速模擬完成！"
-	@echo "    網格檔案: $(MESH_2D)"
-	@echo "    結果檔案: $(RESULT_2DF).msh, $(RESULT_2DF).tec"
-	@echo "=========================================="

 # =============================================================================
 # 3D 模擬流程
@@ -328,23 +307,6 @@ sim3d: $(REFINED_MESH_3D)
 	@echo "    結果檔案:   $(RESULT_3D).msh, $(RESULT_3D).tec"
 	@echo "=========================================="

-# -----------------------------------------------------------------------------
-# 3df: 3D 快速模擬 (無網格細化)
-# 流程: mesh3d → 直接模擬
-# 指令:
-#   1. gmsh -nt N -3 -format msh2 opto_simplified_3d.geo -o opto_simplified_3d.msh
-#   2. python3 opto_3d_run.py opto_simplified_3d.msh
-# -----------------------------------------------------------------------------
-3df: mesh3d
-	@echo ">>> [3df] 執行 3D 快速模擬 (無細化)..."
-	@echo "    指令: $(PYTHON) opto_3d_run.py $(MESH_3D) $(RESULT_3DF)"
-	$(PYTHON) opto_3d_run.py $(MESH_3D) $(RESULT_3DF)
-	@echo ""
-	@echo "=========================================="
-	@echo ">>> [3DF] 快速模擬完成！"
-	@echo "    網格檔案: $(MESH_3D)"
-	@echo "    結果檔案: $(RESULT_3DF).msh, $(RESULT_3DF).tec"
-	@echo "=========================================="

 # =============================================================================
 # CLEAN - 清除所有生成檔案
@@ -1,169 +0,0 @@
-# Optocoupler Simulation 專案結構
-
-## 模型概述
-
-本專案包含兩種模擬模型：
-
-| 特性 | 完整金屬模型 (Full) | 簡化模型 (Simplified) |
-|------|---------------------|----------------------|
-| 金屬處理 | 作為獨立區域 (Regions) | 作為挖空區域 (Cutouts) |
-| 邊界條件 | Equipotential constraints | Dirichlet BC on surfaces |
-| 幾何檔案 | `opto.geo` | `opto_simplified.geo` / `opto_simplified_3d.geo` |
-| 維度 | 2D | 2D / 3D |
-| 狀態 | 已棄用 | **目前使用** |
-
---
-
-## 簡化模型 (Simplified) - 目前使用
-
-### 執行方式
-
-```bash
-make 2d      # 2D 模擬（含網格細化）
-make 3d      # 3D 模擬
-```
-
-### 檔案依賴關係
-
-```
-┌─────────────────────────────────────────────────────────────────┐
-│                        Makefile                                  │
-│  (make 2d / make 3d)                                            │
-└────────────────┬─────────────────────────┬─────────────────────┘
-                 │                         │
-                 ▼                         ▼
-┌────────────────────────────┐  ┌────────────────────────────────┐
-│ 2D Pipeline                │  │ 3D Pipeline                    │
-│                            │  │                                │
-│ opto_simplified.geo        │  │ opto_simplified_3d.geo         │
-│         │                  │  │         │                      │
-│         ▼ (gmsh -2)        │  │         ▼ (gmsh -3)            │
-│ opto_simplified.msh        │  │ opto_simplified_3d.msh         │
-│         │                  │  │         │                      │
-│         ▼                  │  │         ▼                      │
-│ opto_simplified_run.py ────┼──┼─ opto_3d_run.py               │
-│         │                  │  │         │                      │
-│         ├── opto_common.py ◄──┼─────────┤                      │
-│         │                  │  │         │                      │
-│         └── opto_physics_model.py ◄─────┘                      │
-│         │                  │  │         │                      │
-│         ▼                  │  │         ▼                      │
-│ opto_simplified_result.*   │  │ opto_3d_result.*               │
-└────────────────────────────┘  └────────────────────────────────┘
-```
-
-### 2D 模擬流程
-
-```
-1. make 2d
-   ├── mesh2d: gmsh -2 opto_simplified.geo → opto_simplified.msh
-   ├── sim2d:  python opto_simplified_run.py
-   │   ├── 載入網格 (create_gmsh_mesh)
-   │   ├── 建立區域: region_dielectric, region_led, region_encap
-   │   ├── 建立接觸: contact_gnd (0V), contact_hv (100V)
-   │   │             contact_gnd_ext (0V), contact_hv_ext (100V)
-   │   ├── 設定介面連續性
-   │   └── 求解 Poisson 方程
-   └── refine2d: python opto_refine.py → 細化網格後重新模擬
-```
-
-### 3D 模擬流程
-
-```
-1. make 3d
-   ├── mesh3d: gmsh -3 opto_simplified_3d.geo → opto_simplified_3d.msh
-   ├── sim3d:  python opto_3d_run.py
-   │   ├── 載入 3D 網格
-   │   ├── 建立區域: region_dielectric, region_led, region_encap
-   │   ├── 建立接觸: contact_gnd, contact_hv, contact_gnd_ext, contact_hv_ext
-   │   ├── 建立介面 (含 Z 方向)
-   │   └── 求解 3D Poisson 方程
-   └── refine3d: python opto_refine_3d.py → 細化網格後重新模擬
-```
-
---
-
-## 完整金屬模型 (Full) - 已棄用
-
-> **注意**: 舊模型相關檔案已移至 `legacy/` 目錄
-
-### 檔案依賴關係
-
-```
-legacy/opto.geo.bak (需重新命名為 opto.geo 使用)
-    │
-    ▼ (gmsh -2)
-opto.msh
-    │
-    ▼
-legacy/opto_run.py
-    │
-    ├── legacy/opto_device_setup.py  ← 只被此模型使用
-    ├── legacy/opto_device.py
-    ├── opto_common.py
-    └── opto_physics_model.py
-    │
-    ▼
-opto_result.msh / .tec
-```
-
-### 如何恢復使用舊模型
-1. 複製 `legacy/opto.geo.bak` 到主目錄並重新命名為 `opto.geo`
-2. 複製 `legacy/` 內的 Python 檔案到主目錄
-
-### 差異說明
-
-| 項目 | 完整模型 | 簡化模型 |
-|------|----------|----------|
-| 金屬區域 | region_cu_base, region_conductor, region_led_metal, region_hv_pad | 不存在（作為挖空） |
-| 邊界條件 | setup_virtual_contact() 設定整個區域電位 | setup_contact() 設定表面邊界 |
-| 複雜度 | 高（需處理金屬區域的 Poisson 方程） | 低（金屬作為固定電位邊界） |
-| 效率 | 較低 | **較高** |
-
---
-
-## 共用模組
-
-### opto_common.py
- 多執行緒設定
- 物理常數 (eps_0, PERMITTIVITY)
- 介面定義 (INTERFACES)
- 工具函式 (generate_mesh, generate_mesh_3d)
-
-### opto_physics_model.py
- setup_poisson(): Poisson 方程設定
- setup_interface(): 介面連續性
- setup_contact(): Dirichlet 邊界條件
- setup_virtual_contact(): 整區域固定電位（完整模型用）
-
-### opto_refine.py
- 基於 E-field 的動態網格細化（2D / 3D）
-
---
-
-## 檔案清單
-
-### 核心檔案（需保留）
- `Makefile` - 建構腳本
- `opto_simplified.geo` - 2D 幾何
- `opto_simplified_3d.geo` - 3D 幾何
- `opto_simplified_run.py` - 2D 執行腳本
- `opto_3d_run.py` - 3D 執行腳本
- `opto_common.py` - 共用模組
- `opto_physics_model.py` - 物理模型
- `opto_refine.py` - 2D 網格細化
- `opto_refine_3d.py` - 3D 網格細化
- `PROJECT_STRUCTURE.md` - 本文件
-
-### legacy/ 目錄（舊模型檔案）
- `opto.geo.bak` - 完整金屬模型幾何
- `opto_run.py` - 完整模型執行腳本
- `opto_device_setup.py` - 完整模型設定
- `opto_full_run.py` - 完整模型執行腳本
- `opto_device.py` - 設備定義模組
-
-### 可重新生成的檔案
- `*.msh` - 網格檔案（執行 gmsh 可重建）
- `*_result.msh`, `*_result.tec` - 模擬結果
- `__pycache__/` - Python 快取
- `result/` - 舊結果資料
@@ -0,0 +1,171 @@
+# Optocoupler TCAD 模擬專案
+
+DEVSIM 光耦合器電場模擬專案，支援 2D/3D Poisson 方程求解與自適應網格細化。
+
+---
+
+## 專案架構
+
+```
+wisetop_opto/
+├── opto_simplified.geo       # 2D 幾何定義
+├── opto_simplified_3d.geo    # 3D 幾何定義
+├── opto_common.py            # 共用模組 (物理常數、介面定義)
+├── opto_physics_model.py     # 物理模型 (Poisson)
+├── opto_simplified_run.py    # 2D 執行腳本
+├── opto_3d_run.py            # 3D 執行腳本
+├── opto_refine.py            # 2D 網格細化
+├── opto_refine_3d.py         # 3D 網格細化
+├── opto_device.py            # 設備定義模組
+├── refinement_loop.sh        # 2D 網格細化腳本
+├── refinement_loop_3d.sh     # 3D 網格細化腳本
+├── legacy/                   # 舊模型檔案
+└── Makefile                  # 建置自動化
+```
+
+---
+
+## 模型概述
+
+| 特性 | 簡化模型 (目前使用) |
+|------|---------------------|
+| 金屬處理 | 作為挖空區域 (Cutouts) |
+| 邊界條件 | Dirichlet BC on surfaces |
+| 維度 | 2D / 3D |
+
+---
+
+## 快速開始
+
+### 系統需求
+
+| 軟體 | 版本 | 用途 |
+|------|------|------|
+| Python | ≥ 3.8 | 模擬核心 |
+| Gmsh | ≥ 4.0 | 網格生成 |
+| ParaView | 選用 | 結果視覺化 |
+
+### 執行模擬
+
+```bash
+cd devsim/wisetop_opto
+
+make 2d      # 2D 完整模擬 (含網格細化)
+make 3d      # 3D 完整模擬
+
+make 2df     # 2D 快速模擬 (無細化)
+make 3df     # 3D 快速模擬 (無細化)
+```
+
+---
+
+## Make 指令
+
+### 完整流程
+
+| 指令 | 說明 |
+|------|------|
+| `make 2d` | 2D 完整模擬 (mesh + refine + sim) |
+| `make 3d` | 3D 完整模擬 |
+| `make 2df` | 2D 快速模擬 (無網格細化) |
+| `make 3df` | 3D 快速模擬 |
+| `make clean` | 清除所有生成檔案 |
+
+### 個別步驟
+
+| 指令 | 說明 | 輸出檔案 |
+|------|------|----------|
+| `make mesh2d` | 生成 2D 初始網格 | `opto_2d.msh` |
+| `make refine2d` | 2D 網格細化迴圈 | `opto_2d_ref.msh` |
+| `make sim2d` | 執行 2D 模擬 | `opto_2d_out.msh` |
+| `make mesh3d` | 生成 3D 初始網格 | `opto_3d.msh` |
+| `make refine3d` | 3D 網格細化迴圈 | `opto_3d_ref.msh` |
+| `make sim3d` | 執行 3D 模擬 | `opto_3d_out.msh` |
+
+---
+
+## 模擬流程
+
+### 2D 模擬
+
+```
+1. make 2d
+   ├── mesh2d: gmsh -2 opto_simplified.geo → opto_2d.msh
+   ├── refine2d: python3 opto_refine.py → opto_2d_bg.pos
+   │   └── gmsh -bgm opto_2d_bg.pos → opto_2d_ref.msh
+   └── sim2d: python3 opto_simplified_run.py
+       ├── 載入網格 (create_gmsh_mesh)
+       ├── 建立區域: region_dielectric, region_led, region_encap
+       ├── 建立接觸: contact_gnd (0V), contact_hv (100V)
+       └── 求解 Poisson 方程
+```
+
+### 3D 模擬
+
+```
+1. make 3d
+   ├── mesh3d: gmsh -3 opto_simplified_3d.geo → opto_3d.msh
+   ├── refine3d: python3 opto_refine_3d.py → opto_3d_bg.pos
+   └── sim3d: python3 opto_3d_run.py → 求解 3D Poisson
+```
+
+---
+
+## 共用模組
+
+### opto_common.py
+- 多執行緒設定
+- 物理常數 (eps_0, PERMITTIVITY)
+- 介面定義 (INTERFACES)
+- 工具函式 (generate_mesh)
+
+### opto_physics_model.py
+- `setup_poisson()`: Poisson 方程設定
+- `setup_interface()`: 介面連續性
+- `setup_contact()`: Dirichlet 邊界條件
+
+---
+
+## 進階設定
+
+### 環境變數
+
+| 變數 | 說明 | 預設值 |
+|------|------|--------|
+| `LOOPS` | 細化迭代次數 | 2 |
+| `NPROC` | Gmsh 執行緒數 | 自動偵測 |
+
+```bash
+LOOPS=3 make 2d    # 3 次網格細化迭代
+NPROC=4 make 3d    # 指定 4 執行緒
+```
+
+---
+
+## 檔案清單
+
+### 核心檔案 (需保留)
+- `Makefile` - 建構腳本
+- `opto_simplified.geo` - 2D 幾何
+- `opto_simplified_3d.geo` - 3D 幾何
+- `opto_simplified_run.py` - 2D 執行腳本
+- `opto_3d_run.py` - 3D 執行腳本
+- `opto_common.py` - 共用模組
+- `opto_physics_model.py` - 物理模型
+- `opto_refine.py` - 2D 網格細化
+- `opto_refine_3d.py` - 3D 網格細化
+
+### 可重新生成的檔案
+- `*.msh` - 網格檔案
+- `*_out.msh`, `*_out.tec` - 模擬結果
+- `*.pos` - 背景場檔案
+
+---
+
+## 版本資訊
+
+| 項目 | 說明 |
+|------|------|
+| **更新日期** | 2025-12-24 |
+| **Python** | python3 |
+| **網格細化** | LOOPS=2 |
@@ -7,7 +7,6 @@ import subprocess
 from devsim import get_contact_charge, set_parameter


-# --- Multi-threading ---
 def setup_threads(num_threads=None, task_size=1024):
    """Configure DEVSIM threading."""
    if num_threads is None:
@@ -15,12 +14,12 @@ def setup_threads(num_threads=None, task_size=1024):
        num_threads = int(env_threads) if env_threads else (os.cpu_count() or 1)
    set_parameter(name="threads_available", value=num_threads)
    set_parameter(name="threads_task_size", value=task_size)
-    print(f"[DEVSIM] Multi-threading: {num_threads} threads, task_size={task_size}")
+    print(f"[DEVSIM] Threads: {num_threads}, task_size={task_size}")

-setup_threads()  # Auto-enable on import
+setup_threads()


-# --- Path setup ---
+# --- Path Setup ---
 SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
 PROJECT_ROOT = os.path.dirname(SCRIPT_DIR)
 if PROJECT_ROOT not in sys.path:
@@ -28,11 +27,11 @@ if PROJECT_ROOT not in sys.path:
 sys.path.append(os.path.join(PROJECT_ROOT, 'python_packages'))


-# --- Physical constants ---
-eps_0 = 8.854187817e-14  # Vacuum permittivity [F/cm]
+# --- Physical Constants ---
+eps_0 = 8.854187817e-14  # F/cm


-# Relative permittivity (eps_r)
+# --- Permittivity ---
 PERMITTIVITY = {
    'region_atmosphere': 1.0,
    'region_cu_base': 1.0,
@@ -45,8 +44,7 @@ PERMITTIVITY = {
 }


-# --- Interface definitions ---
-# Format: (interface_name, region0, region1)
+# --- Interfaces ---
 CORE_INTERFACES = [
    ("interface_di_led", "region_dielectric", "region_led"),
 ]
@@ -61,10 +59,21 @@ ENCAP_INTERFACES = [

 INTERFACES = CORE_INTERFACES + ENCAP_INTERFACES

+INTERFACES_3D = [
+    ("interface_di_led", "region_dielectric", "region_led"),
+    ("interface_encap_di", "region_encap", "region_dielectric"),
+    ("interface_encap_di_l", "region_encap", "region_dielectric"),
+    ("interface_encap_di_r", "region_encap", "region_dielectric"),
+    ("interface_encap_led_l", "region_encap", "region_led"),
+    ("interface_encap_led_r", "region_encap", "region_led"),
+    ("interface_encap_di_z", "region_encap", "region_dielectric"),
+    ("interface_encap_led_z", "region_encap", "region_led"),
+]

-# --- Utility functions ---
+
+# --- Utilities ---
 def get_contact_charges(device):
-    """Get contact charges for all contacts."""
+    """Get contact charges."""
    charges = {}
    for contact in ["contact_gnd", "contact_hv", "contact_gnd_ext", "contact_hv_ext"]:
        try:
@@ -78,44 +87,31 @@ def get_contact_charges(device):
 def generate_mesh(geo_file, mesh_file, dimension=2, force=False):
    """Generate mesh using Gmsh."""
    if not force and os.path.exists(mesh_file):
-        geo_mtime = os.path.getmtime(geo_file)
-        msh_mtime = os.path.getmtime(mesh_file)
-        if msh_mtime > geo_mtime:
-            print(f"Mesh {mesh_file} is up-to-date")
+        if os.path.getmtime(mesh_file) > os.path.getmtime(geo_file):
+            print(f"Mesh up-to-date: {mesh_file}")
            return mesh_file
    
-    dim_str = "3D " if dimension == 3 else ""
-    print(f"Generating {dim_str}mesh from {geo_file}...")
+    print(f"Generating {'3D ' if dimension == 3 else ''}mesh...")
    result = subprocess.run(
        ["gmsh", f"-{dimension}", geo_file, "-o", mesh_file, "-format", "msh2"],
        capture_output=True, text=True
    )
    if result.returncode != 0:
        print(f"Gmsh error: {result.stderr}")
-        raise RuntimeError(f"{dim_str}Mesh generation failed")
-    print(f"Generated {dim_str}mesh: {mesh_file}")
+        raise RuntimeError("Mesh generation failed")
+    print(f"Generated: {mesh_file}")
    return mesh_file


 # --- Refinement Parameters ---
-MESH_SIZE_MIN = 2.0e-4   # 200 um
-MESH_SIZE_MAX = 2.0e-3   # 2 mm
+MESH_SIZE_MIN_2D = 2.0e-4   # 2 um
+MESH_SIZE_MAX_2D = 5.0e-3   # 50 um
+MESH_SIZE_MIN_3D = 2.0e-4   # 2 um
+MESH_SIZE_MAX_3D = 5.0e-3   # 50 um

-E_VERY_HIGH = 1.0e4   # 10 kV/cm  → 1x mesh size
-E_HIGH = 5.0e3        # 5 kV/cm   → 2x mesh size
-E_MEDIUM = 1.0e3      # 1 kV/cm   → 4x mesh size
-E_LOW = 5.0e2         # 500 V/cm  → 8x mesh size
-E_VERYLOW = 1.0e2     # 100 V/cm  → 16x mesh size
-
-
-# --- 3D Interface Definitions ---
-INTERFACES_3D = [
-    ("interface_di_led", "region_dielectric", "region_led"),
-    ("interface_encap_di", "region_encap", "region_dielectric"),
-    ("interface_encap_di_l", "region_encap", "region_dielectric"),
-    ("interface_encap_di_r", "region_encap", "region_dielectric"),
-    ("interface_encap_led_l", "region_encap", "region_led"),
-    ("interface_encap_led_r", "region_encap", "region_led"),
-    ("interface_encap_di_z", "region_encap", "region_dielectric"),
-    ("interface_encap_led_z", "region_encap", "region_led"),
-]
+# E-field thresholds
+E_VERY_HIGH = 1.0e4   # 10 kV/cm  -> 1x
+E_HIGH = 5.0e3        # 5 kV/cm   -> 2x
+E_MEDIUM = 1.0e3      # 1 kV/cm   -> 4x
+E_LOW = 5.0e2         # 500 V/cm  -> 8x
+E_VERYLOW = 1.0e2     # 100 V/cm  -> 16x
@@ -3,7 +3,6 @@ opto_device.py - Unified Device Setup for Optocoupler Simulation
 """
 import os
 import sys
-
 sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

 from devsim import (
@@ -12,7 +11,6 @@ from devsim import (
    set_parameter, get_node_model_values, get_edge_model_values,
    solve, write_devices, element_from_edge_model
 )
-
 import opto_common
 import opto_physics_model

@@ -21,7 +19,6 @@ class OptoDevice:
    """Unified device setup for optocoupler simulation."""
    
    def __init__(self, name, dimension=2):
-        """Initialize device."""
        self.name = name
        self.dimension = dimension
        self.regions = []
@@ -29,9 +26,7 @@ class OptoDevice:
        
    def load_mesh(self, geo_file, mesh_file, force=False):
        """Generate and load mesh."""
-        # Generate mesh
        opto_common.generate_mesh(geo_file, mesh_file, self.dimension, force)
-        
        self.mesh_file = mesh_file
        print(f"--- Loading {self.dimension}D Mesh: {mesh_file} ---")
        create_gmsh_mesh(mesh=self.name, file=mesh_file)
@@ -40,18 +35,16 @@ class OptoDevice:
        """Add regions with material assignment."""
        for region, material in region_material_map.items():
            try:
-                add_gmsh_region(mesh=self.name, gmsh_name=region,
-                               region=region, material=material)
-                print(f"  Region added: {region}")
+                add_gmsh_region(mesh=self.name, gmsh_name=region, region=region, material=material)
+                print(f"  Region: {region}")
            except Exception as e:
-                print(f"  Warning: {region} failed: {e}")
+                print(f"  Warning: {region} - {e}")
                
    def setup_interfaces(self, interface_list):
        """Add interfaces between regions."""
        for iface_name, region0, region1 in interface_list:
            try:
-                add_gmsh_interface(mesh=self.name, gmsh_name=iface_name,
-                                  region0=region0, region1=region1, name=iface_name)
+                add_gmsh_interface(mesh=self.name, gmsh_name=iface_name, region0=region0, region1=region1, name=iface_name)
            except:
                pass
                
@@ -59,28 +52,23 @@ class OptoDevice:
        """Add contacts for boundary conditions."""
        for contact, region in contact_list:
            try:
-                add_gmsh_contact(mesh=self.name, gmsh_name=contact,
-                                region=region, name=contact, material="metal")
-                print(f"  Contact added: {contact} on {region}")
+                add_gmsh_contact(mesh=self.name, gmsh_name=contact, region=region, name=contact, material="metal")
+                print(f"  Contact: {contact}")
            except Exception as e:
-                print(f"  Warning: {contact} failed: {e}")
+                print(f"  Warning: {contact} - {e}")
                
    def finalize(self):
        """Finalize mesh and create device."""
        finalize_mesh(mesh=self.name)
        create_device(mesh=self.name, device=self.name)
        
-        # Get and store region list
        self.regions = list(get_region_list(device=self.name))
        print(f"Regions: {self.regions}")
        
-        # Set permittivity
        for region in self.regions:
            eps_r = opto_common.PERMITTIVITY.get(region, 1.0)
-            set_parameter(device=self.name, region=region,
-                         name="Permittivity", value=eps_r * opto_common.eps_0)
+            set_parameter(device=self.name, region=region, name="Permittivity", value=eps_r * opto_common.eps_0)
                         
-        # Setup Poisson physics
        opto_physics_model.setup_poisson(self.name, self.regions)
        
    def setup_interface_physics(self):
@@ -102,18 +90,18 @@ class OptoDevice:
                opto_physics_model.setup_contact(self.name, contact, voltage)
                print(f"  {contact}: {voltage} V")
            except Exception as e:
-                print(f"  Warning: {contact} failed: {e}")
+                print(f"  Warning: {contact} - {e}")
                
    def apply_equipotential(self, region_voltage_map):
-        """Apply equipotential constraints to entire regions."""
-        print("--- Metal Equipotential ---")
+        """Apply equipotential constraints."""
+        print("--- Equipotential ---")
        for region, (bias_name, voltage) in region_voltage_map.items():
            try:
                set_parameter(device=self.name, name=bias_name, value=voltage)
                opto_physics_model.setup_virtual_contact(self.name, region, bias_name)
-                print(f"  {region}: {voltage} V (equipotential)")
+                print(f"  {region}: {voltage} V")
            except Exception as e:
-                print(f"  Warning: {region} failed: {e}")
+                print(f"  Warning: {region} - {e}")
                
    def solve(self):
        """Solve Poisson equation."""
@@ -138,14 +126,14 @@ class OptoDevice:
                try:
                    e = get_edge_model_values(device=self.name, region=region, name="ElectricField")
                    e_max = max(abs(min(e)), abs(max(e)))
-                    print(f"  {region}: E_max = {e_max:.0f} V/cm = {e_max/100:.1f} kV/m")
+                    print(f"  {region}: E_max = {e_max:.0f} V/cm")
                except:
                    pass
                    
    def create_element_models(self):
-        """Create element models for visualization (3D)."""
+        """Create element models for 3D visualization."""
        if self.dimension == 3:
-            print("--- Creating 3D Element Models ---")
+            print("--- 3D Element Models ---")
            for region in self.regions:
                try:
                    element_from_edge_model(device=self.name, region=region, edge_model="ElectricField")
@@ -154,8 +142,6 @@ class OptoDevice:
                    
    def output(self, prefix):
        """Write output files."""
-        msh_file = f"{prefix}.msh"
-        tec_file = f"{prefix}.tec"
-        write_devices(file=msh_file, type="devsim")
-        write_devices(file=tec_file, type="tecplot")
-        print(f"Done: {msh_file}, {tec_file}")
+        write_devices(file=f"{prefix}.msh", type="devsim")
+        write_devices(file=f"{prefix}.tec", type="tecplot")
+        print(f"Done: {prefix}.msh, {prefix}.tec")
@@ -10,35 +10,25 @@ from devsim import (
 )


-# Metal regions (treated as cutouts in simplified model)
-METAL_REGIONS = ['region_cu_base', 'region_conductor', 'region_led_metal', 'region_hv_pad']
-
-
 def setup_poisson(device, regions):
    """Setup Poisson equation for dielectric regions."""
    for region in regions:
        node_solution(device=device, region=region, name="Potential")
        edge_from_node_model(device=device, region=region, node_model="Potential")
        
-        # Electric field: E = -dV/dx
+        # E = -dV/dx
        edge_model(device=device, region=region, name="ElectricField",
                   equation="(Potential@n0 - Potential@n1) * EdgeInverseLength")
-        edge_model(device=device, region=region, name="ElectricField:Potential@n0",
-                   equation="EdgeInverseLength")
-        edge_model(device=device, region=region, name="ElectricField:Potential@n1",
-                   equation="-EdgeInverseLength")
+        edge_model(device=device, region=region, name="ElectricField:Potential@n0", equation="EdgeInverseLength")
+        edge_model(device=device, region=region, name="ElectricField:Potential@n1", equation="-EdgeInverseLength")
        
-        # Displacement field: D = eps * E
-        edge_model(device=device, region=region, name="DField",
-                   equation="Permittivity * ElectricField")
-        edge_model(device=device, region=region, name="DField:Potential@n0",
-                   equation="Permittivity * ElectricField:Potential@n0")
-        edge_model(device=device, region=region, name="DField:Potential@n1",
-                   equation="Permittivity * ElectricField:Potential@n1")
+        # D = eps * E
+        edge_model(device=device, region=region, name="DField", equation="Permittivity * ElectricField")
+        edge_model(device=device, region=region, name="DField:Potential@n0", equation="Permittivity * ElectricField:Potential@n0")
+        edge_model(device=device, region=region, name="DField:Potential@n1", equation="Permittivity * ElectricField:Potential@n1")
        
-        # Poisson equation
-        equation(device=device, region=region, name="PotentialEquation",
-                 variable_name="Potential", edge_model="DField", variable_update="default")
+        equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+                 edge_model="DField", variable_update="default")
        
        try:
            element_from_edge_model(device=device, region=region, edge_model="ElectricField")
@@ -48,28 +38,20 @@ def setup_poisson(device, regions):

 def setup_interface(device, interface_name):
    """Setup potential continuity at interface."""
-    interface_model(device=device, interface=interface_name,
-                   name="continuousPotential", equation="Potential@r0 - Potential@r1")
-    interface_model(device=device, interface=interface_name,
-                   name="continuousPotential:Potential@r0", equation="1")
-    interface_model(device=device, interface=interface_name,
-                   name="continuousPotential:Potential@r1", equation="-1")
-    interface_equation(device=device, interface=interface_name,
-                      name="PotentialEquation", interface_model="continuousPotential",
-                      type="continuous")
+    interface_model(device=device, interface=interface_name, name="continuousPotential", equation="Potential@r0 - Potential@r1")
+    interface_model(device=device, interface=interface_name, name="continuousPotential:Potential@r0", equation="1")
+    interface_model(device=device, interface=interface_name, name="continuousPotential:Potential@r1", equation="-1")
+    interface_equation(device=device, interface=interface_name, name="PotentialEquation",
+                       interface_model="continuousPotential", type="continuous")


 def setup_contact(device, contact_name, bias_value=0.0):
-    """Setup Dirichlet boundary condition for contact."""
+    """Setup Dirichlet boundary condition."""
    bias_name = f"{contact_name}_bias"
    set_parameter(device=device, name=bias_name, value=bias_value)
-    
-    contact_node_model(device=device, contact=contact_name,
-                      name=f"{contact_name}_bc", equation=f"Potential - {bias_name}")
-    contact_node_model(device=device, contact=contact_name,
-                      name=f"{contact_name}_bc:Potential", equation="1")
-    contact_equation(device=device, contact=contact_name, name="PotentialEquation",
-                    node_model=f"{contact_name}_bc")
+    contact_node_model(device=device, contact=contact_name, name=f"{contact_name}_bc", equation=f"Potential - {bias_name}")
+    contact_node_model(device=device, contact=contact_name, name=f"{contact_name}_bc:Potential", equation="1")
+    contact_equation(device=device, contact=contact_name, name="PotentialEquation", node_model=f"{contact_name}_bc")


 def setup_contacts(device, contact_map):
@@ -77,21 +59,15 @@ def setup_contacts(device, contact_map):
    for contact, region in contact_map.items():
        bias_name = f"{contact}_bias"
        set_parameter(device=device, name=bias_name, value=0.0)
-        
-        contact_node_model(device=device, contact=contact,
-                          name=f"{contact}_bc", equation=f"Potential - {bias_name}")
-        contact_node_model(device=device, contact=contact,
-                          name=f"{contact}_bc:Potential", equation="1")
+        contact_node_model(device=device, contact=contact, name=f"{contact}_bc", equation=f"Potential - {bias_name}")
+        contact_node_model(device=device, contact=contact, name=f"{contact}_bc:Potential", equation="1")
        contact_equation(device=device, contact=contact, name="PotentialEquation",
-                        node_model=f"{contact}_bc", edge_charge_model="DField")
+                         node_model=f"{contact}_bc", edge_charge_model="DField")


 def setup_virtual_contact(device, region, bias_name):
    """Setup virtual contact for entire region."""
-    node_model(device=device, region=region,
-               name="virtual_contact_bc", equation=f"Potential - {bias_name}")
-    node_model(device=device, region=region,
-               name="virtual_contact_bc:Potential", equation="1")
-    equation(device=device, region=region, name="PotentialEquation",
-             variable_name="Potential", node_model="virtual_contact_bc",
-             variable_update="default")
+    node_model(device=device, region=region, name="virtual_contact_bc", equation=f"Potential - {bias_name}")
+    node_model(device=device, region=region, name="virtual_contact_bc:Potential", equation="1")
+    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+             node_model="virtual_contact_bc", variable_update="default")
@@ -1,35 +1,24 @@
 #!/usr/bin/env python3
 """
-opto_refine.py - Dynamic Mesh Refinement based on E-field
+opto_refine.py - 2D Mesh Refinement based on E-field
 """
-
 import sys
 import os
 from devsim import (
    create_gmsh_mesh, add_gmsh_region, add_gmsh_interface, add_gmsh_contact,
    finalize_mesh, create_device, get_region_list, get_interface_list,
-    set_parameter, solve,
-    get_node_model_values, get_element_model_values,
+    set_parameter, solve, get_node_model_values, get_element_model_values,
    element_from_edge_model, element_model, element_from_node_model
 )
-
-# Import shared modules
 import opto_common
 import opto_physics_model


-# --- Refinement Strategy Functions ---
-
 def emag_refinement(device, region):
-    """Refine based on E-field magnitude."""
-    # Create element E-field from edge model
+    """Create E-field magnitude refinement model."""
    element_from_edge_model(edge_model="ElectricField", device=device, region=region)
-    
-    # Calculate E-field magnitude
    element_model(device=device, region=region, name="Emag",
                  equation="(ElectricField_x^2 + ElectricField_y^2)^(0.5)")
-    
-    # Refinement factor based on E-field magnitude
    element_model(device=device, region=region, name="Enorm", equation=f'''
        ifelse(Emag > {opto_common.E_VERY_HIGH}, 1.0,
          ifelse(Emag > {opto_common.E_HIGH}, 2.0,
@@ -40,17 +29,14 @@ def emag_refinement(device, region):
    return "Enorm"


-# --- Main Refinement Logic ---
-
 def run_refinement():
-    """Execute mesh refinement based on E-field distribution."""
+    """Execute 2D mesh refinement."""
    mesh_file = sys.argv[1] if len(sys.argv) > 1 else "opto_simplified.msh"
    device = "opto_simplified"
    
-    print(f"=== OPTO Dynamic Mesh Refinement ===")
-    print(f"Input mesh: {mesh_file}")
+    print(f"=== 2D Mesh Refinement ===")
+    print(f"Input: {mesh_file}")
    
-    # --- Load mesh ---
    create_gmsh_mesh(mesh=device, file=mesh_file)
    
    region_material = {
@@ -65,39 +51,30 @@ def run_refinement():
        except:
            pass
    
-    # Add interfaces
    for iface_name, region0, region1 in opto_common.INTERFACES:
        try:
-            add_gmsh_interface(mesh=device, gmsh_name=iface_name,
-                              region0=region0, region1=region1, name=iface_name)
+            add_gmsh_interface(mesh=device, gmsh_name=iface_name, region0=region0, region1=region1, name=iface_name)
        except:
            pass
    
-    # Add contacts
-    add_gmsh_contact(mesh=device, gmsh_name="contact_gnd", region="region_dielectric",
-                    name="contact_gnd", material="metal")
-    add_gmsh_contact(mesh=device, gmsh_name="contact_hv", region="region_led",
-                    name="contact_hv", material="metal")
+    add_gmsh_contact(mesh=device, gmsh_name="contact_gnd", region="region_dielectric", name="contact_gnd", material="metal")
+    add_gmsh_contact(mesh=device, gmsh_name="contact_hv", region="region_led", name="contact_hv", material="metal")
    for contact in ["contact_hv_ext", "contact_gnd_ext"]:
        try:
-            add_gmsh_contact(mesh=device, gmsh_name=contact, region="region_encap",
-                            name=contact, material="metal")
+            add_gmsh_contact(mesh=device, gmsh_name=contact, region="region_encap", name=contact, material="metal")
        except:
            pass
    
    finalize_mesh(mesh=device)
    create_device(mesh=device, device=device)
    
-    # --- Set permittivity ---
    regions = [r for r in get_region_list(device=device) if r != "region_led_metal"]
-    print(f"Regions for refinement: {regions}")
+    print(f"Regions: {regions}")
    
    for region in regions:
        eps_r = opto_common.PERMITTIVITY.get(region, 1.0)
-        set_parameter(device=device, region=region, 
-                     name="Permittivity", value=eps_r * opto_common.eps_0)
+        set_parameter(device=device, region=region, name="Permittivity", value=eps_r * opto_common.eps_0)
    
-    # --- Setup physics ---
    opto_physics_model.setup_poisson(device, regions)
    
    for iface in get_interface_list(device=device):
@@ -106,7 +83,6 @@ def run_refinement():
        except:
            pass
    
-    # Setup contacts
    opto_physics_model.setup_contact(device, "contact_gnd", 0.0)
    opto_physics_model.setup_contact(device, "contact_hv", 100.0)
    for contact, voltage in [("contact_gnd_ext", 0.0), ("contact_hv_ext", 100.0)]:
@@ -115,50 +91,34 @@ def run_refinement():
        except:
            pass
    
-    # --- Solve Poisson ---
-    print("--- Solving Poisson (for refinement) ---")
+    print("--- Solving ---")
    try:
        solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=30)
        print("  Converged!")
    except Exception as e:
-        print(f"  Warning: Solve issue - {e}")
+        print(f"  Warning: {e}")
    
-    # --- Apply refinement to each region ---
    all_pos_data = []
    
    for region in regions:
-        print(f"\n--- Processing: {region} ---")
+        print(f"--- {region} ---")
        
        x = get_node_model_values(device=device, region=region, name="x")
        y = get_node_model_values(device=device, region=region, name="y")
        num_nodes = len(x)
        
        if num_nodes == 0:
-            print(f"  Skipped (no nodes)")
            continue
        
-        # E-field refinement
        try:
            strategy = emag_refinement(device, region)
-            print(f"  [OK] E-field refinement applied")
        except Exception as e:
-            print(f"  [FAIL] E-field refinement failed: {e}")
+            print(f"  E-field failed: {e}")
            continue
        
-        # Get refinement values
        element_model(device=device, region=region, name="clen", equation=strategy)
        cl = get_element_model_values(device=device, region=region, name='clen')
        
-        # E-field statistics
-        try:
-            emag = get_element_model_values(device=device, region=region, name='Emag')
-            e_max = max(emag) if emag else 0
-            e_avg = sum(emag) / len(emag) if emag else 0
-            print(f"  E_max = {e_max:.0f} V/cm, E_avg = {e_avg:.0f} V/cm")
-        except:
-            pass
-        
-        # Map element values to nodes
        element_from_node_model(node_model="node_index", device=device, region=region)
        en0 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en0')]
        en1 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en1')]
@@ -169,26 +129,21 @@ def run_refinement():
            v = cl[i] if i < len(cl) else 0
            ni0, ni1, ni2 = en0[i], en1[i], en2[i]
            if v > 0:
-                target = opto_common.MESH_SIZE_MIN * v
+                target = opto_common.MESH_SIZE_MIN_2D * v
                for ni in [ni0, ni1, ni2]:
                    if node_cl[ni] == 0.0 or target < node_cl[ni]:
                        node_cl[ni] = target
        
-        # Fill remaining with max size
        for i in range(num_nodes):
            if node_cl[i] == 0.0:
-                node_cl[i] = opto_common.MESH_SIZE_MAX
+                node_cl[i] = opto_common.MESH_SIZE_MAX_2D
        
-        # Collect for combined .pos file
        all_pos_data.append((region, x, y, node_cl, en0, en1, en2))
-        
-        # Stats
-        refined_count = sum(1 for c in node_cl if c < opto_common.MESH_SIZE_MAX)
-        print(f"  Refined nodes: {refined_count} / {num_nodes}")
+        refined_count = sum(1 for c in node_cl if c < opto_common.MESH_SIZE_MAX_2D)
+        print(f"  Refined: {refined_count}/{num_nodes}")
    
-    # --- Write combined .pos file ---
    pos_file = "opto_2d_bg.pos"
-    print(f"\n--- Writing: {pos_file} ---")
+    print(f"Writing: {pos_file}")
    
    with open(pos_file, 'w') as fh:
        print('View "background mesh" {', file=fh)
@@ -200,13 +155,8 @@ def run_refinement():
                    node_cl[ni0], node_cl[ni1], node_cl[ni2]), file=fh)
        print('};', file=fh)
    
-    print(f"\n=== Done ===")
-    print(f"Background field: {pos_file}")
-    print(f"To regenerate mesh:")
-    print(f"  gmsh -2 -format msh2 opto_simplified.geo -bgm {pos_file} -o opto_2d_ref.msh")
+    print(f"Done: {pos_file}")


-# --- Entry Point ---
-
 if __name__ == "__main__":
    run_refinement()
@@ -1,35 +1,24 @@
 #!/usr/bin/env python3
 """
-opto_refine_3d.py - Dynamic 3D Mesh Refinement based on E-field
+opto_refine_3d.py - 3D Mesh Refinement based on E-field
 """
-
 import sys
 import os
 from devsim import (
    create_gmsh_mesh, add_gmsh_region, add_gmsh_interface, add_gmsh_contact,
    finalize_mesh, create_device, get_region_list, get_interface_list,
-    set_parameter, solve,
-    get_node_model_values, get_element_model_values,
+    set_parameter, solve, get_node_model_values, get_element_model_values,
    element_from_edge_model, element_model, element_from_node_model
 )
-
-# Import shared modules
 import opto_common
 import opto_physics_model


-# --- Refinement Strategy Functions ---
-
 def emag_refinement_3d(device, region):
-    """Refine based on 3D E-field magnitude."""
-    # Create element E-field from edge model
+    """Create 3D E-field magnitude refinement model."""
    element_from_edge_model(edge_model="ElectricField", device=device, region=region)
-    
-    # Calculate 3D E-field magnitude (including Z component)
    element_model(device=device, region=region, name="Emag",
                  equation="(ElectricField_x^2 + ElectricField_y^2 + ElectricField_z^2)^(0.5)")
-    
-    # Refinement factor based on E-field magnitude
    element_model(device=device, region=region, name="Enorm", equation=f'''
        ifelse(Emag > {opto_common.E_VERY_HIGH}, 1.0,
          ifelse(Emag > {opto_common.E_HIGH}, 2.0,
@@ -40,38 +29,29 @@ def emag_refinement_3d(device, region):
    return "Enorm"


-# --- Gmsh Background Field Output (3D) ---
-
 def write_gmsh_pos_3d(filename, all_pos_data):
    """Write Gmsh 3D background field .pos file."""
    print(f"Writing: {filename}")
-    
    with open(filename, 'w') as fh:
        print('View "background mesh" {', file=fh)
        for region, x, y, z, node_cl, en0, en1, en2, en3 in all_pos_data:
            for i in range(len(en0)):
                ni0, ni1, ni2, ni3 = en0[i], en1[i], en2[i], en3[i]
-                # SS = Scalar Tetrahedron in Gmsh .pos format
                print("SS(%g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g) {%g, %g, %g, %g};" % (
-                    x[ni0], y[ni0], z[ni0],
-                    x[ni1], y[ni1], z[ni1],
-                    x[ni2], y[ni2], z[ni2],
-                    x[ni3], y[ni3], z[ni3],
+                    x[ni0], y[ni0], z[ni0], x[ni1], y[ni1], z[ni1],
+                    x[ni2], y[ni2], z[ni2], x[ni3], y[ni3], z[ni3],
                    node_cl[ni0], node_cl[ni1], node_cl[ni2], node_cl[ni3]), file=fh)
        print('};', file=fh)


-# --- Main Refinement Logic ---
-
 def run_refinement():
-    """Execute 3D mesh refinement based on E-field distribution."""
+    """Execute 3D mesh refinement."""
    mesh_file = sys.argv[1] if len(sys.argv) > 1 else "opto_simplified_3d.msh"
    device = "opto_3d"
    
-    print(f"=== OPTO 3D Dynamic Mesh Refinement ===")
-    print(f"Input mesh: {mesh_file}")
+    print(f"=== 3D Mesh Refinement ===")
+    print(f"Input: {mesh_file}")
    
-    # --- Load mesh ---
    create_gmsh_mesh(mesh=device, file=mesh_file)
    
    region_material = {
@@ -85,39 +65,30 @@ def run_refinement():
        except:
            pass
    
-    # Add 3D interfaces
    for iface_name, region0, region1 in opto_common.INTERFACES_3D:
        try:
-            add_gmsh_interface(mesh=device, gmsh_name=iface_name,
-                              region0=region0, region1=region1, name=iface_name)
+            add_gmsh_interface(mesh=device, gmsh_name=iface_name, region0=region0, region1=region1, name=iface_name)
        except:
            pass
    
-    # Add contacts
-    add_gmsh_contact(mesh=device, gmsh_name="contact_gnd", region="region_dielectric",
-                    name="contact_gnd", material="metal")
-    add_gmsh_contact(mesh=device, gmsh_name="contact_hv", region="region_led",
-                    name="contact_hv", material="metal")
+    add_gmsh_contact(mesh=device, gmsh_name="contact_gnd", region="region_dielectric", name="contact_gnd", material="metal")
+    add_gmsh_contact(mesh=device, gmsh_name="contact_hv", region="region_led", name="contact_hv", material="metal")
    for contact in ["contact_hv_ext", "contact_gnd_ext"]:
        try:
-            add_gmsh_contact(mesh=device, gmsh_name=contact, region="region_encap",
-                            name=contact, material="metal")
+            add_gmsh_contact(mesh=device, gmsh_name=contact, region="region_encap", name=contact, material="metal")
        except:
            pass
    
    finalize_mesh(mesh=device)
    create_device(mesh=device, device=device)
    
-    # --- Set permittivity ---
    regions = list(get_region_list(device=device))
-    print(f"Regions for refinement: {regions}")
+    print(f"Regions: {regions}")
    
    for region in regions:
        eps_r = opto_common.PERMITTIVITY.get(region, 1.0)
-        set_parameter(device=device, region=region, 
-                     name="Permittivity", value=eps_r * opto_common.eps_0)
+        set_parameter(device=device, region=region, name="Permittivity", value=eps_r * opto_common.eps_0)
    
-    # --- Setup physics ---
    opto_physics_model.setup_poisson(device, regions)
    
    for iface in get_interface_list(device=device):
@@ -126,7 +97,6 @@ def run_refinement():
        except:
            pass
    
-    # Setup contacts
    opto_physics_model.setup_contact(device, "contact_gnd", 0.0)
    opto_physics_model.setup_contact(device, "contact_hv", 100.0)
    for contact, voltage in [("contact_gnd_ext", 0.0), ("contact_hv_ext", 100.0)]:
@@ -135,19 +105,17 @@ def run_refinement():
        except:
            pass
    
-    # --- Solve Poisson ---
-    print("--- Solving 3D Poisson (for refinement) ---")
+    print("--- Solving ---")
    try:
        solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=30)
        print("  Converged!")
    except Exception as e:
-        print(f"  Warning: Solve issue - {e}")
+        print(f"  Warning: {e}")
    
-    # --- Apply refinement to each region ---
    all_pos_data = []
    
    for region in regions:
-        print(f"\n--- Processing: {region} ---")
+        print(f"--- {region} ---")
        
        x = get_node_model_values(device=device, region=region, name="x")
        y = get_node_model_values(device=device, region=region, name="y")
@@ -155,31 +123,17 @@ def run_refinement():
        num_nodes = len(x)
        
        if num_nodes == 0:
-            print(f"  Skipped (no nodes)")
            continue
        
-        # E-field refinement (3D)
        try:
            strategy = emag_refinement_3d(device, region)
-            print(f"  [OK] 3D E-field refinement applied")
        except Exception as e:
-            print(f"  [FAIL] E-field refinement failed: {e}")
+            print(f"  E-field failed: {e}")
            continue
        
-        # Get refinement values
        element_model(device=device, region=region, name="clen", equation=strategy)
        cl = get_element_model_values(device=device, region=region, name='clen')
        
-        # E-field statistics
-        try:
-            emag = get_element_model_values(device=device, region=region, name='Emag')
-            e_max = max(emag) if emag else 0
-            e_avg = sum(emag) / len(emag) if emag else 0
-            print(f"  E_max = {e_max:.0f} V/cm, E_avg = {e_avg:.0f} V/cm")
-        except:
-            pass
-        
-        # Map element values to nodes (3D tetrahedron: 4 nodes per element)
        element_from_node_model(node_model="node_index", device=device, region=region)
        en0 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en0')]
        en1 = [int(v) for v in get_element_model_values(device=device, region=region, name='node_index@en1')]
@@ -191,35 +145,23 @@ def run_refinement():
            v = cl[i] if i < len(cl) else 0
            ni0, ni1, ni2, ni3 = en0[i], en1[i], en2[i], en3[i]
            if v > 0:
-                target = opto_common.MESH_SIZE_MIN * v
+                target = opto_common.MESH_SIZE_MIN_3D * v
                for ni in [ni0, ni1, ni2, ni3]:
                    if node_cl[ni] == 0.0 or target < node_cl[ni]:
                        node_cl[ni] = target
        
-        # Fill remaining with max size
        for i in range(num_nodes):
            if node_cl[i] == 0.0:
-                node_cl[i] = opto_common.MESH_SIZE_MAX
+                node_cl[i] = opto_common.MESH_SIZE_MAX_3D
        
-        # Collect for combined .pos file
        all_pos_data.append((region, x, y, z, node_cl, en0, en1, en2, en3))
-        
-        # Stats
-        refined_count = sum(1 for c in node_cl if c < opto_common.MESH_SIZE_MAX)
-        print(f"  Refined nodes: {refined_count} / {num_nodes}")
+        refined_count = sum(1 for c in node_cl if c < opto_common.MESH_SIZE_MAX_3D)
+        print(f"  Refined: {refined_count}/{num_nodes}")
    
-    # --- Write combined 3D .pos file ---
    pos_file = "opto_3d_bg.pos"
-    print(f"\n--- Writing: {pos_file} ---")
    write_gmsh_pos_3d(pos_file, all_pos_data)
-    
-    print(f"\n=== Done ===")
-    print(f"Background field: {pos_file}")
-    print(f"To regenerate mesh:")
-    print(f"  gmsh -3 -format msh2 opto_simplified_3d.geo -bgm {pos_file} -o opto_3d_ref.msh")
+    print(f"Done: {pos_file}")


-# --- Entry Point ---
-
 if __name__ == "__main__":
    run_refinement()
@@ -62,12 +62,23 @@ BooleanDifference(10) = { Surface{1}; Delete; }{ Surface{2,3,4,7}; Delete; };
 v() = BooleanFragments{ Surface{10,5,6}; Delete; }{};

 /* === Physical Surfaces === */
-Physical Surface("region_dielectric") = {5};
-Physical Surface("region_led")        = {6};
-Physical Surface("region_encap")      = {7};
-
 d = 0.001;  // BoundingBox tolerance

+// Dielectric region
+di_surf() = Surface In BoundingBox{x_di_l-d, y_di_bot-d, -d, x_di_r+d, y_di_top+d, d};
+Physical Surface("region_dielectric") = {di_surf()};
+
+// LED region
+led_surf() = Surface In BoundingBox{x_led_l-d, y_led_bot-d, -d, x_led_r+d, y_led_top+d, d};
+Physical Surface("region_led") = {led_surf()};
+
+// Encapsulation (all surfaces in encap range, EXCLUDING nested regions)
+enc_surf_candidates() = Surface In BoundingBox{x_enc_l-d, y_enc_bot-d, -d, x_enc_r+d, y_enc_top+d, d};
+enc_surf() = enc_surf_candidates();
+enc_surf() -= di_surf();
+enc_surf() -= led_surf();
+Physical Surface("region_encap") = {enc_surf()};
+
 /* === Contacts === */
 contact_gnd() = Curve In BoundingBox{x_di_l-d, y_di_bot-d, -d, x_di_r+d, y_di_bot+d, d};
 Physical Curve("contact_gnd") = {contact_gnd()};
@@ -125,65 +136,9 @@ Physical Curve("interface_encap_di_r") = {ifc_enc_di_r()};
 Physical Curve("interface_encap_led_l") = {ifc_enc_led_l()};
 Physical Curve("interface_encap_led_r") = {ifc_enc_led_r()};

-/* === Mesh Refinement === */
-Field[1] = Attractor;
-Field[1].CurvesList = {ifc_di_led(), ifc_enc_di_top_l(), ifc_enc_di_top_r(),
-                       ifc_enc_di_l(), ifc_enc_di_r()};
+/* === Mesh Settings === */
+// Initial mesh uses lc_base, refinement can go down to lc_fine
+Mesh.CharacteristicLengthMin = lc_fine;
+Mesh.CharacteristicLengthMax = lc_base;

-Field[2] = Attractor;
-Field[2].CurvesList = {contact_gnd_ext_cu_t_l(), contact_gnd_ext_cu_t_r(),
-                       contact_gnd_ext_cu_r(), contact_gnd_ext_cu_l(), contact_gnd_ext_cu_b(),
-                       contact_gnd_ext_cond_l(), contact_gnd_ext_cond_r(),
-                       contact_gnd_ext_cond_t_l(), contact_gnd_ext_cond_t_r()};

-Field[3] = Attractor;
-Field[3].CurvesList = {contact_hv_ext_r(), contact_hv_ext_t(), contact_hv_ext_b()};
-
-Field[4] = Threshold;
-Field[4].InField = 1;
-Field[4].SizeMin = lc_fine;
-Field[4].SizeMax = lc_base;
-Field[4].DistMin = 0.001;
-Field[4].DistMax = 0.02;
-
-Field[5] = Threshold;
-Field[5].InField = 2;
-Field[5].SizeMin = lc_fine;
-Field[5].SizeMax = lc_base;
-Field[5].DistMin = 0.001;
-Field[5].DistMax = 0.02;
-
-Field[6] = Threshold;
-Field[6].InField = 3;
-Field[6].SizeMin = lc_fine;
-Field[6].SizeMax = lc_base;
-Field[6].DistMin = 0.001;
-Field[6].DistMax = 0.02;
-
-Field[7] = Box;
-Field[7].VIn = lc_fine * 2;
-Field[7].VOut = lc_base;
-Field[7].XMin = x_di_l;  Field[7].XMax = x_di_r;
-Field[7].YMin = y_di_bot; Field[7].YMax = y_di_top;
-
-Field[8] = Box;
-Field[8].VIn = lc_fine;
-Field[8].VOut = lc_base;
-Field[8].XMin = x_led_l; Field[8].XMax = x_led_r;
-Field[8].YMin = y_led_bot; Field[8].YMax = y_led_top;
-
-Field[9] = Box;
-Field[9].VIn = lc_fine;
-Field[9].VOut = lc_base;
-Field[9].XMin = x_led_l; Field[9].XMax = x_led_r;
-Field[9].YMin = y_metal_bot; Field[9].YMax = y_metal_top;
-
-Field[11] = Box;
-Field[11].VIn = lc_fine;
-Field[11].VOut = lc_base;
-Field[11].XMin = x_di_l; Field[11].XMax = x_di_r;
-Field[11].YMin = y_di_top - 0.005; Field[11].YMax = y_di_top + 0.005;
-
-Field[10] = Min;
-Field[10].FieldsList = {4, 5, 6, 7, 8, 9, 11};
-Background Field = 10;
@@ -18,7 +18,7 @@ Mesh.Algorithm3D = 10;
 sf = 0.1;                         // 1 mm = 0.1 cm
 um = 1.0e-4;                      // 1 µm = 1e-4 cm
 lc_fine = 20.0 * um;
-lc_base = 200.0 * um;
+lc_base = 50.0 * um;

 /* XY Dimensions (SAME as 2D version) */
 W_cu = 3.5*sf;    H_cu = 0.2*sf;      cu_ext = 1.0*sf;
@@ -213,21 +213,7 @@ ifc_enc_led_z_front() = Surface In BoundingBox{x_led_l-d, y_led_bot-d, z_led-d,
 ifc_enc_led_z_back() = Surface In BoundingBox{x_led_l-d, y_led_bot-d, -z_led-d, x_led_r+d, y_led_top+d, -z_led+d};
 Physical Surface("interface_encap_led_z") = {ifc_enc_led_z_front(), ifc_enc_led_z_back()};

-/* === Mesh Refinement === */
-Field[1] = Box;
-Field[1].VIn = lc_fine;
-Field[1].VOut = lc_base;
-Field[1].XMin = x_led_l - 0.01; Field[1].XMax = x_led_r + 0.01;
-Field[1].YMin = y_led_bot - 0.01; Field[1].YMax = y_metal_top + 0.01;
-Field[1].ZMin = -z_led; Field[1].ZMax = z_led;
-
-Field[2] = Box;
-Field[2].VIn = lc_fine * 2;
-Field[2].VOut = lc_base;
-Field[2].XMin = x_di_l;  Field[2].XMax = x_di_r;
-Field[2].YMin = y_di_bot; Field[2].YMax = y_di_top;
-Field[2].ZMin = -z_di; Field[2].ZMax = z_di;
-
-Field[3] = Min;
-Field[3].FieldsList = {1, 2};
-Background Field = 3;
+/* === Mesh Settings === */
+// Initial mesh uses lc_base, refinement can go down to lc_fine
+Mesh.CharacteristicLengthMin = lc_fine;
+Mesh.CharacteristicLengthMax = lc_base;
Author	SHA1	Message	Date
Dean Huang	cd4c78dcdc	Update	2025-12-26 16:04:34 +08:00
Dean Huang	37a2c32fde	feat: 調整網格細化參數。	2025-12-26 16:02:08 +08:00
Dean Huang	baf9ac1f4a	feat: 統一了 BJT 和 Opto 模型的程式風格網格細化設定。	2025-12-26 15:44:41 +08:00
Dean Huang	4ea453c75c	Update	2025-12-26 11:52:41 +08:00
Dean Huang	7b09cf5d56	Update	2025-12-26 11:30:52 +08:00
Dean Huang	c73fb0cb16	移除設定求解器的功能。	2025-12-26 10:42:48 +08:00
Dean Huang	610fa597ba	feat: 調整求解器預設值為 `superlu` 放寬網格尺寸以減少元素數量。	2025-12-26 10:35:18 +08:00
Dean Huang	d8ab2db6a5	feat: 新增 DEVSIM 線性求解器設定功能並預設啟用 MKL Pardiso	2025-12-26 10:24:49 +08:00
Dean Huang	90db272e25	feat: 區分 2D 與 3D 網格細化參數並更新相關引用	2025-12-26 09:58:39 +08:00
Dean Huang	fc5f69e4f7	更新 wisetop_opto 與 wisetop_bjt 模組的編譯設定。	2025-12-26 08:26:09 +08:00
Dean Huang	88b04500c4	docs: 更新README。	2025-12-24 17:29:43 +08:00
Dean Huang	af07cef7c8	feat: 調整 Python 腳本。	2025-12-24 16:32:29 +08:00
Dean Huang	2b770a77a0	Merge commit 'a008bb23ca76f0d8152c431ae606015ad34da2b8'	2025-12-24 16:22:04 +08:00