feat: Stabilizing High-Voltage BJT Simulation

feat: Remove .msh file generation for Paraview
feat: Add PROJECT GUIDE.md.
2026-03-27 16:13:23 +08:00 · 2026-03-19 16:30:19 +08:00 · 2026-03-16 14:48:02 +08:00 · 2026-03-16 09:51:13 +08:00 · 2026-03-13 16:01:52 +08:00 · 2026-03-13 15:13:51 +08:00
66 changed files with 1535 additions and 1483 deletions
@@ -0,0 +1,280 @@
+# DEVSIM TCAD 專案 — 完整安裝與執行流程
+
+## 專案概覽
+
+**DEVSIM** 是一套開源的 TCAD（Technology Computer-Aided Design）元件模擬器，採用有限體積法（Finite Volume Method）求解半導體方程式，透過 Python 腳本驅動模擬。
+
+本專案路徑為 `/home/wsx52114/devsim`，目錄結構如下：
+
+```
+devsim/                           # 專案根目錄
+├── denv/                         # Python 虛擬環境（已建立）
+├── devsim/                       # DEVSIM 原始碼 + 自訂模擬
+│   ├── wisetop_bjt/              # ★ 自訂 BJT 模擬專案（主要工作區）
+│   ├── wisetop_opto/             # 光電元件模擬
+│   ├── examples/                 # 官方範例（diode, capacitance 等）
+│   ├── testing/                  # 官方測試腳本
+│   ├── src/                      # C++ 原始碼
+│   ├── scripts/                  # 建置腳本
+│   └── INSTALL.md, BUILD.md ...  # 官方文件
+└── devsim_bjt_example-main/      # 官方 BJT 範例（論文用）
+```
+
+---
+
+## 一、環境安裝
+
+### 1.1 系統需求
+
+| 軟體 | 最低版本 | 用途 | 目前已安裝版本 |
+|------|---------|------|---------------|
+| **Python** | ≥ 3.9 | 模擬核心引擎 | 3.12.3 |
+| **Gmsh** | ≥ 4.0 | 有限元素網格生成 | 4.15.0 |
+| **ParaView** | 選用 | `.tec` 結果視覺化 | — |
+
+### 1.2 安裝步驟（從零開始）
+
+```bash
+# ── Step 1：安裝系統套件 ──
+sudo apt update
+sudo apt install python3 python3-venv gmsh
+
+# ── Step 2：Clone 專案 ──
+git clone https://gitlab.com/wisetop/devsim-tcad/devsim.git
+cd devsim
+
+# ── Step 3：建立 Python 虛擬環境 ──
+python3 -m venv denv
+
+# ── Step 4：啟用虛擬環境 ──
+source denv/bin/activate
+
+# ── Step 5：安裝 Python 套件 ──
+pip install devsim numpy
+# 選用：安裝 MKL（提升線性代數運算效能）
+pip install mkl
+```
+
+> [!IMPORTANT]
+> **每次開啟新終端都必須先啟用虛擬環境**：`source denv/bin/activate`
+
+### 1.3 目前已安裝的套件
+
+本專案虛擬環境 `denv/` 中已安裝：
+
+| 套件 | 版本 | 說明 |
+|------|------|------|
+| `devsim` | 2.10.0 | TCAD 模擬核心 |
+| `gmsh` | 4.15.0 | Gmsh Python API |
+| `numpy` | 2.4.0 | 數值計算 |
+| `mkl` | 2025.3.0 | Intel Math Kernel Library |
+| `intel_openmp` | 2025.3.1 | 多執行緒支援 |
+
+### 1.4 驗證安裝
+
+```bash
+source denv/bin/activate
+
+# 測試 devsim
+python3 -c "import devsim; print('devsim OK')"
+
+# 測試 gmsh
+gmsh --version
+
+# 跑官方範例
+cd devsim/testing
+python3 cap2.py
+```
+
+---
+
+## 二、模擬執行流程
+
+### 2.1 主要工作區：`wisetop_bjt/`
+
+該目錄包含自訂的平面 NPN BJT 模擬，已封裝為 **Makefile** 自動化流程。
+
+```bash
+# 進入工作目錄
+cd /home/wsx52114/devsim/devsim/wisetop_bjt
+
+# 啟用虛擬環境（若尚未啟用）
+source /home/wsx52114/devsim/denv/bin/activate
+```
+
+### 2.2 模擬流程圖
+
+```
+┌─────────┐    ┌─────────┐    ┌─────────┐    ┌──────────────────┐
+│  mesh   │───▶│ refine  │───▶│  init   │───▶│ gummel / output  │
+│ (Gmsh)  │    │ (迴圈)  │    │ (DC=0V) │    │   / cb / ac      │
+└─────────┘    └─────────┘    └─────────┘    └──────────────────┘
+  bjt.msh    bjt_refined.msh  bjt_dd_0.msh       *.csv
+```
+
+### 2.3 Make 指令一覽
+
+#### 完整流程指令
+
+| 指令 | 說明 | 包含步驟 |
+|------|------|---------|
+| `make all` | **完整模擬**（推薦） | refine → init → gummel → output → cb → ac |
+| `make quick` | 快速模擬（無網格細化） | mesh → init → gummel |
+| `make clean` | 清除所有生成檔案 | 刪除 `*.msh *.pos *.tec *.csv *.log` |
+
+#### 個別步驟指令
+
+| 指令 | 說明 | 輸入 | 輸出 |
+|------|------|------|------|
+| `make mesh` | 生成初始網格 | `bjt.geo` | `bjt.msh` |
+| `make refine` | 網格細化迴圈 | `bjt.geo` | `bjt_refined.msh` |
+| `make init` | 零偏壓 DD 初始化 | `bjt_refined.msh` | `bjt_dd_0.msh`, `bjt_dd_0.tec` |
+| `make gummel` | Gummel 曲線 (Ic,Ib vs Vb) | `bjt_dd_0.msh` | `gummel.csv` |
+| `make output` | 輸出特性 (Ic vs Vc) | `bjt_dd_0.msh` | `output_Vb0.70.csv` |
+| `make cb` | 共基極特性 (Ic vs Ve) | `bjt_dd_0.msh` | `cb_Vc0.50.csv` |
+| `make ac` | AC 小信號分析 | `bjt_dd_0.msh` | `ac_Vc0.5_Ve-0.7.csv` |
+| `make batch` | 批次多參數掃描 | `bjt_dd_0.msh` | `data/*.log` |
+
+#### 環境變數調整
+
+```bash
+# 設定網格細化迭代次數（預設 1）
+make all LOOPS=3
+
+# 指定 Gmsh 執行緒數
+make all NPROC=8
+
+# 修改偏壓參數
+make output VB_DEFAULT=0.8
+make cb VC_DEFAULT=1.0
+```
+
+---
+
+## 三、各步驟詳細說明
+
+### 3.1 網格生成 (`make mesh`)
+
+```bash
+gmsh -nt <NPROC> -2 -format msh2 bjt.geo -o bjt.msh
+```
+
+- 讀取 `bjt.geo` 幾何定義檔
+- 生成 2D 網格，格式為 MSH2
+- `-nt` 啟用多執行緒
+
+### 3.2 網格細化 (`make refine`)
+
+迭代流程（重複 LOOPS 次）：
+
+1. 用 `bjt_refine.py` 分析電場，計算各區域需要的網格密度
+2. 輸出背景場檔 `bjt_bg.pos`
+3. Gmsh 依據背景場重新生成網格
+4. 最終結果複製為 `bjt_refined.msh`
+
+### 3.3 零偏壓初始化 (`make init`)
+
+```bash
+python3 bjt_dd.py
+```
+
+- 載入細化後的網格
+- 設定 ERFC 摻雜模型
+- 求解 Poisson 方程（Potential Only → Drift-Diffusion）
+- 零偏壓 DC 求解
+- 輸出 `bjt_dd_0.msh`（可續算狀態）和 `bjt_dd_0.tec`（視覺化用）
+
+### 3.4 電氣特性分析
+
+| 分析 | 腳本 | 說明 |
+|------|------|------|
+| **Gummel** | `bjt_circuit3.py` | 固定 Vc=2V，掃描 Vb 0→0.8V（20mV 步進） |
+| **輸出特性** | `bjt_circuit2.py <Vb>` | 固定 Vb，掃描 Vc 0→2V（0.1V 步進） |
+| **共基極** | `bjt_circuit4.py <Vc>` | 固定 Vc，掃描 Ve 0→-1V（-50mV 步進） |
+| **AC 分析** | `bjt_circuit5.py <Vc> <fmin> <fmax> <ppd>` | AC 小信號，頻率 1kHz→100GHz |
+
+---
+
+## 四、結果檢視
+
+```bash
+# 網格視覺化
+gmsh bjt_refined.msh
+
+# 摻雜/電位分佈（需 ParaView）
+paraview bjt_dd_0.tec
+# 載入後選擇 Variable: LogNetDoping 或 Potential
+
+# CSV 數據
+# 可用 Python/matplotlib 繪圖，或用試算表軟體開啟
+```
+
+---
+
+## 五、官方範例
+
+### 5.1 內建範例
+
+位於 [devsim/examples/](file:///home/wsx52114/devsim/devsim/examples)：
+- `diode/` — 二極體模擬
+- `capacitance/` — 電容模擬
+- `mobility/` — 遷移率模型
+- `bioapp1/` — 生物應用
+- `plotting/` — 繪圖範例
+
+### 5.2 BJT 論文範例
+
+位於 [devsim_bjt_example-main/](file:///home/wsx52114/devsim/devsim_bjt_example-main)，對應論文 *Semiconductor Device Simulation Using DEVSIM*。
+
+---
+
+## 六、從原始碼建置（進階）
+
+若需修改 DEVSIM 核心，可從原始碼建置：
+
+```bash
+# Clone
+git clone https://github.com/devsim/devsim
+cd devsim
+git submodule init
+git submodule update
+
+# Linux 建置
+bash scripts/build_manylinux_2_28.sh <version>
+
+# 或使用 Docker
+bash scripts/build_docker_manylinux_2_28.sh <version>
+
+# 建置後的 .whl 檔在 dist/ 目錄
+pip install dist/devsim-<version>*.whl
+```
+
+---
+
+## 七、疑難排解
+
+| 問題 | 解決方案 |
+|------|---------|
+| `ModuleNotFoundError: devsim` | 確認虛擬環境已啟用：`source denv/bin/activate` |
+| `devsim` 未安裝 | `pip install devsim` |
+| `gmsh: command not found` | `sudo apt install gmsh` 或 `pip install gmsh` |
+| 模擬收斂失敗 | `make clean && make refine && make init` 重新開始 |
+| 摻雜分佈異常 | 用 ParaView 檢視 `bjt_dd_0.tec` 的 `LogNetDoping` |
+| 記憶體不足 | 減少網格細化次數：`make all LOOPS=1` |
+
+---
+
+## 八、快速開始指令摘要
+
+```bash
+# ── 完整流程（一行搞定）──
+cd /home/wsx52114/devsim/devsim/wisetop_bjt
+source /home/wsx52114/devsim/denv/bin/activate
+make all
+
+# ── 快速測試 ──
+make quick
+
+# ── 清除重來 ──
+make clean && make all
+```
@@ -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 一鍵執行 | 需手動執行多檔 |
-| **可維護** |  模組化清晰 | 邏輯集中 |
-| **學習曲線** | 結構複雜但說明完整 | 精簡易讀 |
-| **擴展性** |  易新增功能 | 需重構 |
@@ -1,27 +1,27 @@
 # =============================================================================
 # Makefile for DEVSIM BJT Simulation
 # =============================================================================
-# 本 Makefile 詳細列出每個模擬步驟所需的指令
+# -----------------------------------------------------------------------------
+# 快速操作:
+#   make all        - 執行完整模擬：網格細化 → 零偏壓初始化 → 各項特性分析
+#   make quick      - 快速模擬驗證：跳過網格細化，直接分析
 #
-# 使用方式:
-#   make all        - 完整模擬流程 (refine → init → 所有分析)
-#   make quick      - 快速模擬 (無網格細化)
-#
-# 細部步驟:
-#   make mesh       - 生成初始網格
-#   make refine     - 執行網格細化迴圈 (LOOPS 次)
-#   make init       - 零偏壓初始化 (drift-diffusion)
-#   make gummel     - Gummel 曲線掃描
-#   make output     - 輸出特性掃描
-#   make cb         - 共基極特性掃描
-#   make ac         - AC 小信號分析
-#   make batch      - 批次參數掃描 (多組參數)
+# 個別分析項目:
+#   make mesh       - 生成初始網格 (.msh)
+#   make refine     - 依據電場強度執行網格細化迴圈 (次數由 LOOPS 參數控制)
+#   make init       - 零偏壓初始化 (Drift-Diffusion 模型)
+#   make gummel     - 繪製 Gummel 曲線 (Ic, Ib 對 Vb)
+#   make output     - 輸出特性曲線 (Ic 對 Vc)
+#   make cb         - 共基極特性曲線 (Ic 對 Ve)
+#   make ac         - 小信號交流響應分析 (頻率掃描)
+#   make static     - 靜態偏壓分析（輸出場分佈供 ParaView 檢視）
+#   make batch      - 批次執行所有模擬條件
 # =============================================================================

 # -----------------------------------------------------------------------------
 # 工具設定
 # -----------------------------------------------------------------------------
-PYTHON    ?= python                               # Python 解譯器
+PYTHON    ?= python3                              # Python 解譯器
 GMSH      ?= gmsh                                 # Gmsh 路徑
 NPROC     ?= $(shell nproc 2>/dev/null || echo 4) # CPU 核心數 (自動偵測)
 GMSH_OPTS  = -nt $(NPROC)                         # Gmsh 多執行緒選項
@@ -38,19 +38,25 @@ BG_POS       = bjt_bg.pos
 # -----------------------------------------------------------------------------
 # 細化設定
 # -----------------------------------------------------------------------------
-LOOPS ?= 3
+LOOPS ?= 1

 # -----------------------------------------------------------------------------
 # 預設參數
 # -----------------------------------------------------------------------------
 VB_DEFAULT ?= 0.7
 VC_DEFAULT ?= 0.5
+VE_DEFAULT ?= 0.0
 VC_GUMMEL  ?= 2.0

+# static 命令預設環境變數
+VB ?= 0.7
+VC ?= 30.0
+VE ?= 0.0
+
 # -----------------------------------------------------------------------------
 # PHONY 目標宣告
 # -----------------------------------------------------------------------------
-.PHONY: help all quick mesh refine init gummel output cb ac batch clean
+.PHONY: help all quick mesh refine init gummel output cb ac static batch clean

 # =============================================================================
 # HELP - 顯示可用指令
@@ -74,16 +80,17 @@ help:
 	@echo "    make output           輸出特性 (Ic vs Vc)"
 	@echo "    make cb               共基極特性 (Ic vs Ve)"
 	@echo "    make ac               AC 小信號分析"
+	@echo "    make static           靜態偏壓分析 (場分佈輸出)"
 	@echo "    make batch            批次參數掃描"
 	@echo ""
 	@echo "  其他:"
 	@echo "    make clean            清除所有生成檔案"
 	@echo ""
 	@echo "  環境變數:"
-	@echo "    LOOPS=N               細化迭代次數 (預設: 3)"
-	@echo "    NPROC=N               設定 Gmsh 執行緒數 (預設: 自動偵測)"
-	@echo "    VB_DEFAULT=x          output 的 Vb 值 (預設: 0.7)"
-	@echo "    VC_DEFAULT=x          cb/ac 的 Vc 值 (預設: 0.5)"
+	@echo "    LOOPS=N               網格細化迭代次數 (預設: 1，數值越大網格越密)"
+	@echo "    NPROC=N               設定 Gmsh 處理器核心數 (預設: 自動偵測系統規格)"
+	@echo "    VB_DEFAULT=x          設定 output 分析的基極電壓 Vb (預設: 0.7V)"
+	@echo "    VC_DEFAULT=x          設定 cb/ac 分析的集極電壓 Vc (預設: 0.5V)"
 	@echo ""
 	@echo "  目前設定: NPROC=$(NPROC), LOOPS=$(LOOPS)"
 	@echo ""
@@ -93,8 +100,8 @@ help:
 # =============================================================================

 # all: 完整模擬流程
-# 流程: refine → init → gummel → output → cb → ac
-all: refine init gummel output cb ac
+# 流程: refine → init → gummel → output → cb → ac → static
+all: refine init gummel output cb ac static
 	@echo "=========================================="
 	@echo ">>> [ALL] 完整模擬流程完成！"
 	@echo "=========================================="
@@ -171,7 +178,7 @@ refine: mesh
 # 指令: python bjt_dd.py
 # 輸入: bjt_refined.msh (或 bjt.msh)
 # 輸出: bjt_dd_0.msh, bjt_dd_0.tec
-# 說明:
+# 說明: 
 #   1. 載入網格並設定摻雜
 #   2. 求解 Potential Only
 #   3. 設定 Drift-Diffusion 模型
@@ -245,10 +252,10 @@ output:
 	@echo "    參數:"
 	@echo "      - Vb = $(VB_DEFAULT)V (固定)"
 	@echo "      - Vc: 0 → 2V (0.1V 步進)"
-	@echo "    輸出: output_Vb$(VB_DEFAULT).csv"
+	@echo "    輸出: output_Vb$(VB_DEFAULT)0.csv"
 	$(PYTHON) bjt_circuit2.py $(VB_DEFAULT)
 	@echo ""
-	@echo ">>> [output] 完成: output_Vb$(VB_DEFAULT).csv"
+	@echo ">>> [output] 完成: output_Vb$(VB_DEFAULT)0.csv"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -272,10 +279,10 @@ cb:
 	@echo "    參數:"
 	@echo "      - Vc = $(VC_DEFAULT)V (固定)"
 	@echo "      - Ve: 0 → -1V (-50mV 步進)"
-	@echo "    輸出: cb_Vc$(VC_DEFAULT).csv"
+	@echo "    輸出: cb_Vc$(VC_DEFAULT)0.csv"
 	$(PYTHON) bjt_circuit4.py $(VC_DEFAULT)
 	@echo ""
-	@echo ">>> [cb] 完成: cb_Vc$(VC_DEFAULT).csv"
+	@echo ">>> [cb] 完成: cb_Vc$(VC_DEFAULT)0.csv"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -307,6 +314,29 @@ ac:
 	@echo ">>> [ac] 完成"
 	@echo ""

+# -----------------------------------------------------------------------------
+# static: 靜態偏壓分析 (輸出場分佈 .tec)
+# 指令: python bjt_circuit6.py <Vb> <Vc> [Ve]
+# 輸入: bjt_dd_0.msh
+# 輸出: static_Vb*_Vc*_Ve*.tec
+# 說明:
+#   - 給定特定偏壓，DC 求解後輸出完整場分佈
+#   - 用 ParaView 檢視 Potential, Electrons, Holes, EField 等
+# -----------------------------------------------------------------------------
+static:
+	@if [ ! -f $(INIT_MSH) ]; then \
+		echo ">>> [static] 未找到初始狀態，先執行 init..."; \
+		$(MAKE) init; \
+	fi
+	@echo ">>> [static] 執行靜態偏壓分析..."
+	@echo "    指令: $(PYTHON) bjt_circuit6.py $(VB) $(VC) $(VE)"
+	@echo "    輸入: $(INIT_MSH)"
+	@echo "    參數: Vb=$(VB)V, Vc=$(VC)V, Ve=$(VE)V"
+	$(PYTHON) bjt_circuit6.py $(VB) $(VC) $(VE)
+	@echo ""
+	@echo ">>> [static] 完成: 用 paraview static_*.tec 檢視結果"
+	@echo ""
+
 # -----------------------------------------------------------------------------
 # batch: 批次參數掃描
 # 指令: bash sims.sh
@@ -341,9 +371,14 @@ batch:
 # CLEAN - 清除所有生成檔案
 # =============================================================================
 clean:
-	@echo ">>> [clean] 清除所有生成檔案..."
+	@echo ">>> [clean] 清除所有生成檔案與暫存..."
 	@echo "    rm -f *.msh *.pos *.tec *.csv *.log"
 	@echo "    rm -rf data"
+	@echo "    rm -rf __pycache__ physics/__pycache__"
+	@echo "    rm -f *.pyc physics/*.pyc"
 	rm -f *.msh *.pos *.tec *.csv *.log
 	rm -rf data
+	rm -rf __pycache__ physics/__pycache__
+	find . -name '*.pyc' -delete 2>/dev/null || true
+	find . -name '__pycache__' -type d -exec rm -rf {} + 2>/dev/null || true
 	@echo ">>> [clean] 完成"
@@ -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 指令
@@ -121,21 +94,37 @@ make all    # 執行完整模擬流程 (已啟用多執行緒加速)

 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
-| `make all` | 執行完整模擬流程 | 所有 CSV 檔案 |
-| `make batch` | 批次參數掃描 | `data/` 目錄下的結果 |
-| `make clean` | 清除所有生成檔案 | — |
+| `make all` | 執行完整模擬流程 | 所有 CSV 與 TEC 檔案 |
+| `make quick` | 快速模擬 (無網格細化) | — |
+| `make clean` | 清除所有生成檔案與暫存檔 (`__pycache__`, `*.pyc`) | — |

 ### 個別步驟

 | 指令 | 說明 | 輸出檔案 |
 |------|------|----------|
 | `make mesh` | 產生初始網格 | `bjt.msh` |
-| `make refine` | 網格細化 (3 次迭代) | `bjt_refined.msh` |
+| `make refine` | 網格細化 (LOOPS 次迭代) | `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` |
 | `make cb` | 共基極掃描 (Vc=0.5V) | `cb_Vc0.50.csv` |
 | `make ac` | AC 小信號分析 | `ac_Vc0.5_Ve-0.7.csv` |
+| `make static` | 靜態偏壓分析 (場分佈輸出) | `static_Vb*_Vc*_Ve*.tec`, `*.msh` |
+| `make batch` | 批次參數掃描 | `data/*.log` |
+
+### 靜態偏壓分析 (`bjt_circuit6.py`)
+
+給定 Vb、Vc、Ve，DC 求解後輸出場分佈 `.tec` 檔，用 ParaView 觀察電位、電場、載子濃度等（檔案內輸出與註解為全英文）。
+
+```bash
+make static VB=0.7 VC=2.0              # Vb=0.7V, Vc=2.0V, Ve=0V
+make static VB=0.7 VC=2.0 VE=-0.5      # 三端都指定
+make static                            # 使用預設值 (全部 0V)
+```
+
+輸出 `static_Vb*_Vc*_Ve*.tec`，用 `paraview static_*.tec` 開啟。
+
+可觀察變數：`Potential`、`Electrons`/`Holes`、`LogElectrons`/`LogHoles`、`NetDoping`、`LogNetDoping`、`EField`

 ---

@@ -144,20 +133,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 +178,6 @@ paraview bjt_dd_0.tec

 ### 多執行緒設定

-專案預設使用系統全部 CPU 核心進行平行運算（通過 `bjt_common.py` 自動啟用）：
-
 ```bash
 # 預設使用全部 CPU 核心
 make all
@@ -177,33 +185,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 +209,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 +218,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
@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+"""
+bjt_circuit6.py - Static DC Bias Analysis (export .tec for ParaView)
+"""
+import sys
+import os
+import math
+from devsim import (
+    load_devices, set_parameter, solve, write_devices,
+    node_model, edge_from_node_model, edge_model
+)
+
+from bjt_common import DOPING_PARAMS, get_contact_currents, restore_doping_params
+from physics.new_physics_c6 import GetContactBiasName, SetSiliconParameters, CreateNodeModel, CreateNodeModelDerivative
+
+
+def _ramp_contact(device, contact, target, step_v=0.1, max_retries=3):
+    """Ramp a single contact from 0 to target voltage in fixed step_v increments."""
+    if target == 0.0:
+        return
+    n_steps = max(1, math.ceil(abs(target) / step_v))
+    print(f"  Ramping {contact}: 0 -> {target}V ({n_steps} steps, ~{abs(target)/n_steps:.3f}V/step)")
+    for i in range(1, n_steps + 1):
+        v = target * i / n_steps
+        set_parameter(device=device, name=GetContactBiasName(contact), value=v)
+        try:
+            solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+        except Exception:
+            # Retry with finer sub-steps
+            ok = target * (i - 1) / n_steps
+            success = False
+            for retry in range(1, max_retries + 1):
+                n_sub = 2 ** retry
+                try:
+                    for j in range(1, n_sub + 1):
+                        sv = ok + (v - ok) * j / n_sub
+                        set_parameter(device=device, name=GetContactBiasName(contact), value=sv)
+                        solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+                    success = True
+                    break
+                except Exception:
+                    set_parameter(device=device, name=GetContactBiasName(contact), value=ok)
+                    try:
+                        solve(type="dc", absolute_error=1e6, relative_error=1.0, maximum_iterations=100)
+                    except Exception:
+                        pass
+            if not success:
+                raise RuntimeError(f"Failed to ramp {contact} to {v:.3f}V")
+    print(f"  {contact} -> {target}V OK")
+
+
+def run_static(Vb=0.7, Vc=2.0, Ve=0.0):
+    """Apply static DC bias and export field distributions."""
+    device = "bjt"
+    region = "Silicon"
+    init_file = "bjt_dd_0.msh"
+
+    if not os.path.exists(init_file):
+        sys.exit("Error: Run 'make init' first")
+
+    print(f"Loading: {init_file}")
+    load_devices(file=init_file)
+
+    # Override IntrinsicElectrons to prevent overflow at high voltage
+    # Clamp exp argument to ±80 (corresponds to 10^35 upper limit)
+    ie_exp = "NIE*exp(ifelse(Potential/V_t > 80, 80, ifelse(Potential/V_t < -80, -80, Potential/V_t)))"
+    CreateNodeModel(device, region, "IntrinsicElectrons", ie_exp)
+    CreateNodeModelDerivative(device, region, "IntrinsicElectrons", ie_exp, 'Potential')
+    CreateNodeModel(device, region, "IntrinsicHoles", "NIE^2/IntrinsicElectrons")
+    CreateNodeModelDerivative(device, region, "IntrinsicHoles", "NIE^2/IntrinsicElectrons", 'Potential')
+
+    restore_doping_params(device, region)
+    SetSiliconParameters(device, region)
+
+    # Initialize contact bias
+    for contact in ["collector", "emitter", "base"]:
+        set_parameter(device=device, name=GetContactBiasName(contact), value=0.0)
+
+    # Ramp bias sequentially: Vc -> Vb -> Ve
+    print(f"Ramping bias: Vc={Vc}V, Vb={Vb}V, Ve={Ve}V")
+    try:
+        _ramp_contact(device, "collector", Vc)
+        _ramp_contact(device, "base", Vb)
+        _ramp_contact(device, "emitter", Ve)
+    except RuntimeError as e:
+        print(f"FAILED: {e}")
+        sys.exit(1)
+    print("DC Solve OK")
+
+    # Setup additional models for visualization
+    node_model(device=device, region=region, name="BuiltinPotential", equation="V_t * asinh(NetDoping / (2 * NIE))")
+    node_model(device=device, region=region, name="AppliedPotential", equation="Potential - BuiltinPotential")
+    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")
+    node_model(device=device, region=region, name="LogElectrons", equation="log(Electrons)/log(10)")
+    node_model(device=device, region=region, name="LogHoles", equation="log(Holes)/log(10)")
+    node_model(device=device, region=region, name="NetCarrier", equation="Electrons - Holes")
+
+    try:
+        node_model(device=device, region=region, name="LogNetDoping", equation="asinh(NetDoping/2)/log(10)")
+    except Exception:
+        pass
+
+    # Print terminal currents
+    currents = get_contact_currents(device)
+    print("\n[Terminal Currents]")
+    print(f"Ic = {currents['collector']:+.6e} A/cm")
+    print(f"Ib = {currents['base']:+.6e} A/cm")
+    print(f"Ie = {currents['emitter']:+.6e} A/cm")
+    print(f"KCL: Ic+Ib+Ie = {sum(currents.values()):+.4e} A/cm")
+
+    if abs(currents['collector']) > 1e-30:
+        beta = abs(currents['collector'] / currents['base']) if abs(currents['base']) > 1e-30 else float('inf')
+        print(f"β ≈ {beta:.1f}")
+
+    # Output files
+    tag = f"static_Vb{Vb:.2f}_Vc{Vc:.2f}_Ve{Ve:.2f}"
+    tec_file = f"{tag}.tec"
+
+    write_devices(file=tec_file, type="tecplot")
+
+    print(f"\nDone: {tec_file}")
+
+
+if __name__ == "__main__":
+    if len(sys.argv) < 3:
+        print("Usage: python3 bjt_circuit6.py <Vb> <Vc> [Ve]")
+        sys.exit(0)
+
+    Vb = float(sys.argv[1])
+    Vc = float(sys.argv[2])
+    Ve = float(sys.argv[3]) if len(sys.argv) > 3 else 0.0
+    run_static(Vb=Vb, Vc=Vc, Ve=Ve)
+
@@ -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")


@@ -0,0 +1,484 @@
+# Copyright 2013 Devsim LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from devsim import *
+from .model_create import *
+
+def SetUniversalParameters(device, region):
+    universal = {
+        'q' :          1.6e-19,          #, 'coul'),
+      'k' :          1.3806503e-23,    #, 'J/K'),
+      'Permittivity_0' :  8.85e-14     #, 'F/cm^2')
+    }
+    for k, v in universal.items():
+        set_parameter(device=device, region=region, name=k, value=v)
+
+
+
+def SetSiliconParameters(device, region):
+    '''
+      Sets Silicon device parameters on the specified region.
+    '''
+
+    SetUniversalParameters(device, region)
+
+##D. B. M. Klaassen, J. W. Slotboom, and H. C. de Graaff, "Unified apparent bandgap narrowing in n- and p-type Silicon," Solid-State Electronics, vol. 35, no. 2, pp. 125-29, 1992.
+    par = {
+        'Permittivity'     : 11.1*get_parameter(device=device, region=region, name='Permittivity_0'),
+      'NC300'       : 2.8e19,  # '1/cm^3'
+      'NV300'       : 3.1e19, # '1/cm^3'
+      'EG300'       : 1.12,    # 'eV'
+      'EGALPH'      : 2.73e-4, # 'eV/K'
+      'EGBETA'      : 0      , # 'K'
+      'Affinity'    : 4.05      , # 'K'
+      # Canali model
+      'BETAN0'      : 2.57e-2, # '1'
+      'BETANE'      : 0.66,    # '1'
+      'BETAP0'      : 0.46,    # '1'
+      'BETAPE'      : 0.17,    # '1'
+      'VSATN0'      : 1.43e9,
+      'VSATNE'      : -0.87,
+      'VSATP0'      : 1.62e8,
+      'VSATPE'      : -0.52,
+      # Arora model
+      'MUMN'  : 88,
+      'MUMEN' : -0.57,
+      'MU0N'  : 7.4e8,
+      'MU0EN' : -2.33,
+      'NREFN' : 1.26e17,
+      'NREFNE' : 2.4,
+      'ALPHA0N' : 0.88,
+      'ALPHAEN' : -0.146,
+      'MUMP'  : 54.3,
+      'MUMEP' : -0.57,
+      'MU0P'  : 1.36e8,
+      'MU0EP' : -2.23,
+      'NREFP' : 2.35e17,
+      'NREFPE' : 2.4,
+      'ALPHA0P' : 0.88,
+      'ALPHAEP' : -0.146,
+      # SRH
+      "taun" : 1e-5,
+      "taup" : 1e-5,
+      "n1" : 1e10,
+      "p1" : 1e10,
+      # TEMP
+      "T" : 300
+    }
+
+    for k, v in par.items():
+        set_parameter(device=device, region=region, name=k, value=v)
+
+def CreateQuasiFermiLevels(device, region, electron_model, hole_model, variables):
+    '''
+    Creates the models for the quasi-Fermi levels.  Assuming Boltzmann statistics.
+    '''
+    eq = (
+        ('EFN', 'EC + V_t * log(%s/NC)' % electron_model, ('Potential', 'Electrons')),
+      ('EFP', 'EV - V_t * log(%s/NV)' % hole_model, ('Potential', 'Holes')),
+    )
+    for (model, equation, variable_list) in eq:
+        #print "MODEL: " + model + " equation " + equation
+        CreateNodeModel(device, region, model, equation)
+        vset = set(variable_list)
+        for v in variables:
+            if v in vset:
+                CreateNodeModelDerivative(device, region, model, equation, v)
+
+def CreateDensityOfStates(device, region, variables):
+    '''
+      Set up models for density of states.
+      Neglects Bandgap narrowing.
+    '''
+    eq = (
+        ('NC', 'NC300 * (T/300)^1.5', ('T',)),
+      ('NV', 'NV300 * (T/300)^1.5', ('T',)),
+      ('NTOT', 'Donors + Acceptors', ()),
+      # Band Gap Narrowing
+      ('DEG', '0', ()),
+      #('DEG', 'V0.BGN * (log(NTOT/N0.BGN) + ((log(NTOT/N0.BGN)^2 + CON.BGN)^(0.5)))', ()),
+      ('EG', 'EG300 + EGALPH*((300^2)/(300+EGBETA) - (T^2)/(T+EGBETA)) - DEG', ('T')),
+      ('NIE', '((NC * NV)^0.5) * exp(-EG/(2*V_t))*exp(DEG)', ('T')),
+      ('EC', '-Potential - Affinity - DEG/2', ('Potential',)),
+      ('EV', 'EC - EG + DEG/2', ('Potential', 'T')),
+      ('EI', '0.5 * (EC + EV + V_t*log(NC/NV))', ('Potential', 'T')),
+    )
+
+    for (model, equation, variable_list) in eq:
+        #print "MODEL: " + model + " equation " + equation
+        CreateNodeModel(device, region, model, equation)
+        vset = set(variable_list)
+        for v in variables:
+            if v in vset:
+                CreateNodeModelDerivative(device, region, model, equation, v)
+
+
+def GetContactBiasName(contact):
+    return "{0}_bias".format(contact)
+
+def GetContactNodeModelName(contact):
+    return "{0}nodemodel".format(contact)
+
+
+def CreateVT(device, region, variables):
+    '''
+      Calculates the thermal voltage, based on the temperature.
+      V_t : node model
+      V_t_edge : edge model from arithmetic mean
+    '''
+    CreateNodeModel(device, region, 'V_t', "k*T/q")
+    CreateArithmeticMean(device, region, 'V_t', 'V_t_edge')
+    if 'T' in variables:
+        CreateArithmeticMeanDerivative(device, region, 'V_t', 'V_t_edge', 'T')
+
+
+def CreateEField(device, region):
+    '''
+      Creates the EField and DField.
+    '''
+    edge_average_model(device=device, region=region, node_model="Potential",
+                       edge_model="EField", average_type="negative_gradient")
+    edge_average_model(device=device, region=region, node_model="Potential",
+                       edge_model="EField", average_type="negative_gradient", derivative="Potential")
+
+def CreateDField(device, region):
+    CreateEdgeModel(device, region, "DField", "Permittivity * EField")
+    CreateEdgeModel(device, region, "DField:Potential@n0", "Permittivity * EField:Potential@n0")
+    CreateEdgeModel(device, region, "DField:Potential@n1", "Permittivity * EField:Potential@n1")
+
+def CreateSiliconPotentialOnly(device, region):
+    '''
+      Creates the physical models for a Silicon region for equilibrium simulation.
+    '''
+
+    variables = ("Potential",)
+    CreateVT(device, region, variables)
+    CreateDensityOfStates(device, region, variables)
+
+    SetSiliconParameters(device, region)
+
+    # require NetDoping
+    for i in (
+        ("IntrinsicElectrons",       "NIE*exp(ifelse(Potential/V_t > 80, 80, ifelse(Potential/V_t < -80, -80, Potential/V_t)))"),
+         ("IntrinsicHoles",           "NIE^2/IntrinsicElectrons"),
+         ("IntrinsicCharge",          "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"),
+         ("PotentialIntrinsicCharge", "-q * IntrinsicCharge")
+    ):
+        n = i[0]
+        e = i[1]
+        CreateNodeModel(device, region, n, e)
+        CreateNodeModelDerivative(device, region, n, e, 'Potential')
+
+    CreateQuasiFermiLevels(device, region, 'IntrinsicElectrons', 'IntrinsicHoles', variables)
+
+    CreateEField(device, region)
+    CreateDField(device, region)
+
+    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+             node_model="PotentialIntrinsicCharge", edge_model="DField", variable_update="log_damp")
+
+def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
+    '''
+      Creates the potential equation at the contact
+      if is_circuit is true, than use node given by GetContactBiasName
+    '''
+    if not InNodeModelList(device, region, "contactcharge_node"):
+        CreateNodeModel(device, region, "contactcharge_node", "q*IntrinsicCharge")
+
+    celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
+    -V_t*log({1}/NIE), \
+    V_t*log({2}/NIE))".format(GetContactBiasName(contact), celec_model, chole_model)
+
+    contact_model_name = GetContactNodeModelName(contact)
+    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
+    if is_circuit:
+        CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")
+
+    if is_circuit:
+        contact_equation(device=device, contact=contact, name="PotentialEquation",
+                         node_model=contact_model_name, edge_model="",
+                         node_charge_model="contactcharge_node", edge_charge_model="DField",
+                         node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
+    else:
+        contact_equation(device=device, contact=contact, name="PotentialEquation",
+                         node_model=contact_model_name, edge_model="",
+                         node_charge_model="contactcharge_node", edge_charge_model="DField",
+                         node_current_model="", edge_current_model="")
+
+
+def CreateSRH(device, region, variables):
+    '''
+      Shockley Read hall recombination model in terms of generation.
+    '''
+    USRH="(Electrons*Holes - NIE^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
+    Gn = "-q * USRH"
+    Gp = "+q * USRH"
+    CreateNodeModel(device, region, "USRH", USRH)
+    CreateNodeModel(device, region, "ElectronGeneration", Gn)
+    CreateNodeModel(device, region, "HoleGeneration", Gp)
+    for i in ("Electrons", "Holes", "T"):
+        if i in variables:
+            CreateNodeModelDerivative(device, region, "USRH", USRH, i)
+            CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
+            CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)
+
+def CreateECE(device, region, Jn):
+    '''
+      Electron Continuity Equation using specified equation for Jn
+    '''
+    NCharge = "q * Electrons"
+    CreateNodeModel(device, region, "NCharge", NCharge)
+    CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")
+
+    equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
+             time_node_model = "NCharge",
+             edge_model=Jn, variable_update="positive", node_model="ElectronGeneration")
+
+def CreateHCE(device, region, Jp):
+    '''
+      Hole Continuity Equation using specified equation for Jp
+    '''
+    PCharge = "-q * Holes"
+    CreateNodeModel(device, region, "PCharge", PCharge)
+    CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")
+
+    equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
+             time_node_model = "PCharge",
+             edge_model=Jp, variable_update="positive", node_model="HoleGeneration")
+
+def CreatePE(device, region):
+    '''
+      Create Poisson Equation assuming the Electrons and Holes as solution variables
+    '''
+    pne = "-q*kahan3(Holes, -Electrons, NetDoping)"
+    CreateNodeModel(device, region,           "PotentialNodeCharge", pne)
+    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
+    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")
+
+    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
+             node_model="PotentialNodeCharge", edge_model="DField",
+             time_node_model="", variable_update="log_damp")
+
+
+def CreateSiliconDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p", Jn='Jn', Jp='Jp'):
+    '''
+      Instantiate all equations for drift diffusion simulation
+    '''
+    CreateDensityOfStates(device, region, ("Potential",))
+    CreateQuasiFermiLevels(device, region, "Electrons", "Holes", ("Electrons", "Holes", "Potential"))
+    CreatePE(device, region)
+    CreateSRH(device, region, ("Electrons", "Holes", "Potential"))
+    CreateECE(device, region, Jn)
+    CreateHCE(device, region, Jp)
+
+
+def CreateSiliconDriftDiffusionContact(device, region, contact, Jn, Jp, is_circuit=False): 
+    '''
+      Restrict electrons and holes to their equilibrium values
+      Integrates current into circuit
+    '''
+    CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit)
+
+    celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * NIE^2)^(0.5)))"
+    contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, NIE^2/{1})".format(celec_model, chole_model)
+    contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +NIE^2/{0})".format(celec_model, chole_model)
+    contact_electrons_name = "{0}nodeelectrons".format(contact)
+    contact_holes_name = "{0}nodeholes".format(contact)
+
+    CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")
+
+    CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
+    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")
+
+    if is_circuit:
+        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
+                         node_model=contact_electrons_name,
+                         edge_current_model=Jn, circuit_node=GetContactBiasName(contact))
+
+        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
+                         node_model=contact_holes_name,
+                         edge_current_model=Jp, circuit_node=GetContactBiasName(contact))
+
+    else:
+        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
+                         node_model=contact_electrons_name,
+                         edge_current_model=Jn)
+
+        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
+                         node_model=contact_holes_name,
+                         edge_current_model=Jp)
+
+
+def CreateBernoulliString (Potential="Potential", scaling_variable="V_t", sign=-1):
+    '''
+    Creates the Bernoulli function for Scharfetter Gummel
+    sign -1 for potential
+    sign +1 for energy
+    scaling variable should be V_t
+    Potential should be scaled by V_t in V
+    Ec, Ev should scaled by V_t in eV
+
+    returns the Bernoulli expression and its argument
+    Caller should understand that B(-x) = B(x) + x
+    '''
+
+    tdict = {
+        "Potential" : Potential,
+      "V_t"       : scaling_variable
+    }
+    #### test for requisite models here
+    if sign == -1:
+        vdiff="(%(Potential)s@n0 - %(Potential)s@n1)/%(V_t)s" % tdict
+    elif sign == 1:
+        vdiff="(%(Potential)s@n1 - %(Potential)s@n0)/%(V_t)s" % tdict
+    else:
+        raise NameError("Invalid Sign %s" % sign)
+
+    Bern01 = "B(%s)" % vdiff
+    return (Bern01, vdiff)
+
+
+def CreateElectronCurrent(device, region, mu_n, Potential="Potential", sign=-1, ElectronCurrent="ElectronCurrent", V_t="V_t_edge"):
+    '''
+    Electron current
+    mu_n = mobility name
+    Potential is the driving potential
+    '''
+    EnsureEdgeFromNodeModelExists(device, region, "Potential")
+    EnsureEdgeFromNodeModelExists(device, region, "Electrons")
+    EnsureEdgeFromNodeModelExists(device, region, "Holes")
+    if Potential == "Potential":
+        (Bern01, vdiff) = CreateBernoulliString(scaling_variable=V_t, Potential=Potential, sign=sign)
+    else:
+        raise NameError("Implement proper call")
+
+    tdict = {
+        'Bern01' : Bern01,
+      'vdiff'  : vdiff,
+      'mu_n'   : mu_n,
+      'V_t'    : V_t
+    }
+
+    Jn = "q*%(mu_n)s*EdgeInverseLength*%(V_t)s*kahan3(Electrons@n1*%(Bern01)s,  Electrons@n1*%(vdiff)s,  -Electrons@n0*%(Bern01)s)" % tdict
+
+    CreateEdgeModel(device, region, ElectronCurrent, Jn)
+    for i in ("Electrons", "Potential", "Holes"):
+        CreateEdgeModelDerivatives(device, region, ElectronCurrent, Jn, i)
+
+def CreateHoleCurrent(device, region, mu_p, Potential="Potential", sign=-1, HoleCurrent="HoleCurrent", V_t="V_t_edge"):
+    '''
+    Hole current
+    '''
+    EnsureEdgeFromNodeModelExists(device, region, "Potential")
+    EnsureEdgeFromNodeModelExists(device, region, "Electrons")
+    EnsureEdgeFromNodeModelExists(device, region, "Holes")
+    # Make sure the bernoulli functions exist
+    if Potential == "Potential":
+        (Bern01, vdiff) = CreateBernoulliString(scaling_variable=V_t, Potential=Potential, sign=sign)
+    else:
+        raise NameError("Implement proper call for " + Potential)
+
+    tdict = {
+        'Bern01' : Bern01,
+      'vdiff'  : vdiff,
+      'mu_p'   : mu_p,
+      'V_t'    : V_t
+    }
+
+    Jp ="-q*%(mu_p)s*EdgeInverseLength*%(V_t)s*kahan3(Holes@n1*%(Bern01)s, -Holes@n0*%(Bern01)s, -Holes@n0*%(vdiff)s)" % tdict
+    CreateEdgeModel(device, region, HoleCurrent, Jp)
+    for i in ("Holes", "Potential", "Electrons"):
+        CreateEdgeModelDerivatives(device, region, HoleCurrent, Jp, i)
+
+def CreateAroraMobilityLF(device, region):
+    '''
+      Creates node mobility models and then averages them on edge
+      Uses model from Muller and Kamins
+      Add T derivative dependence later
+    '''
+    models = (
+        ('Tn', 'T/300'),
+      ('mu_arora_n_node',
+            'MUMN * pow(Tn, MUMEN) + (MU0N * pow(T, MU0EN))/(1 + pow((NTOT/(NREFN*pow(Tn, NREFNE))), ALPHA0N*pow(Tn, ALPHAEN)))'),
+      ('mu_arora_p_node',
+            'MUMP * pow(Tn, MUMEP) + (MU0P * pow(T, MU0EP))/(1 + pow((NTOT/(NREFP*pow(Tn, NREFPE))), ALPHA0P*pow(Tn, ALPHAEP)))')
+    )
+
+    for k, v in models:
+        CreateNodeModel(device, region, k, v)
+    CreateArithmeticMean(device, region, 'mu_arora_n_node', 'mu_arora_n_lf')
+    CreateArithmeticMean(device, region, 'mu_arora_p_node', 'mu_arora_p_lf')
+    CreateElectronCurrent(device, region, mu_n = 'mu_arora_n_lf', Potential="Potential", sign=-1, ElectronCurrent="Jn_arora_lf", V_t="V_t_edge")
+    CreateHoleCurrent(device, region, mu_p = 'mu_arora_p_lf', Potential="Potential", sign=-1, HoleCurrent="Jp_arora_lf", V_t="V_t_edge")
+    return {
+        'mu_n' : 'mu_arora_n_lf',
+      'mu_p' : 'mu_arora_p_lf',
+      'Jn'   : 'Jn_arora_lf',
+      'Jp'   : 'Jp_arora_lf',
+    }
+
+
+def CreateHFMobility(device, region, mu_n, mu_p, Jn, Jp):
+    '''
+      Add T derivatives when debugged
+      use parameters to set model flags
+      Caughey Thomas
+    '''
+
+    tdict = {
+        'Jn' : Jn,
+      'mu_n' : mu_n,
+      'Jp' : Jp,
+      'mu_p' : mu_p
+    }
+    tlist = (
+        ("vsat_n", "VSATN0 * pow(T, VSATNE)" % tdict, ('T')),
+      ("beta_n", "BETAN0 * pow(T, BETANE)" % tdict, ('T')),
+      ("Epar_n",
+            "ifelse((%(Jn)s * EField) > 0, abs(EField), 1e-15)" % tdict, ('Potential')),
+      ("mu_n", "%(mu_n)s * pow(1 + pow((%(mu_n)s*Epar_n/vsat_n), beta_n), -1/beta_n)"
+            % tdict, ('Electrons', 'Holes', 'Potential', 'T')),
+      ("vsat_p", "VSATP0 * pow(T, VSATPE)" % tdict, ('T')),
+      ("beta_p", "BETAP0 * pow(T, BETAPE)" % tdict, ('T')),
+      ("Epar_p",
+            "ifelse((%(Jp)s * EField) > 0, abs(EField), 1e-15)" % tdict, ('Potential')),
+      ("mu_p", "%(mu_p)s * pow(1 + pow(%(mu_p)s*Epar_p/vsat_p, beta_p), -1/beta_p)"
+            % tdict, ('Electrons', 'Holes', 'Potential', 'T')),
+    )
+
+    variable_list = ('Electrons', 'Holes', 'Potential')
+    for (model, equation, variables) in tlist:
+        CreateEdgeModel(device, region, model, equation)
+        for v in variable_list:
+            if v in variables:
+                CreateEdgeModelDerivatives(device, region, model, equation, v)
+
+    # This create derivatives automatically
+    CreateElectronCurrent(device, region, mu_n='mu_n', Potential="Potential", sign=-1, ElectronCurrent="Jn", V_t="V_t_edge")
+    CreateHoleCurrent(    device, region, mu_p='mu_p', Potential="Potential", sign=-1, HoleCurrent="Jp", V_t="V_t_edge")
+    return {
+        'mu_n' : 'mu_n',
+      'mu_p' : 'mu_p',
+      'Jn'   : 'Jn',
+      'Jp'   : 'Jp',
+    }
+
+
+
+
@@ -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          清除所有生成檔案"
@@ -143,13 +139,13 @@ refine2d:
 # sim2d: 執行 2D 模擬 (使用細化後的網格)
 # 指令: python3 opto_simplified_run.py opto_refined.msh
 # 輸入: opto_refined.msh
-# 輸出: opto_simplified_result.msh, opto_simplified_result.tec
+# 輸出: opto_simplified_result.tec
 # -----------------------------------------------------------------------------
 sim2d: $(REFINED_MESH_2D)
 	@echo ">>> [sim2d] 執行 2D 模擬..."
 	@echo "    指令: $(PYTHON) opto_simplified_run.py $(REFINED_MESH_2D) $(RESULT_2D)"
 	$(PYTHON) opto_simplified_run.py $(REFINED_MESH_2D) $(RESULT_2D)
-	@echo ">>> [sim2d] 完成: $(RESULT_2D).msh, $(RESULT_2D).tec"
+	@echo ">>> [sim2d] 完成: $(RESULT_2D).tec"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -162,7 +158,7 @@ sim2d: $(REFINED_MESH_2D)
 #   - $(MESH_2D)         - 初始網格
 #   - $(BG_POS_2D)       - 背景場
 #   - $(REFINED_MESH_2D) - 細化後網格
-#   - $(RESULT_2D).msh/.tec - 模擬結果
+#   - $(RESULT_2D).tec      - 模擬結果
 # -----------------------------------------------------------------------------
 2d:
 	@echo "=========================================="
@@ -193,26 +189,9 @@ sim2d: $(REFINED_MESH_2D)
 	@echo "    初始網格:   $(MESH_2D)"
 	@echo "    細化網格:   $(REFINED_MESH_2D)"
 	@echo "    背景場:     $(BG_POS_2D)"
-	@echo "    結果檔案:   $(RESULT_2D).msh, $(RESULT_2D).tec"
+	@echo "    結果檔案:   $(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 模擬流程
@@ -275,13 +254,13 @@ refine3d:
 # sim3d: 執行 3D 模擬 (使用細化後的網格)
 # 指令: python3 opto_3d_run.py opto_3d_refined.msh
 # 輸入: opto_3d_refined.msh
-# 輸出: opto_3d_result.msh, opto_3d_result.tec
+# 輸出: opto_3d_result.tec
 # -----------------------------------------------------------------------------
 sim3d: $(REFINED_MESH_3D)
 	@echo ">>> [sim3d] 執行 3D 模擬..."
 	@echo "    指令: $(PYTHON) opto_3d_run.py $(REFINED_MESH_3D) $(RESULT_3D)"
 	$(PYTHON) opto_3d_run.py $(REFINED_MESH_3D) $(RESULT_3D)
-	@echo ">>> [sim3d] 完成: $(RESULT_3D).msh, $(RESULT_3D).tec"
+	@echo ">>> [sim3d] 完成: $(RESULT_3D).tec"
 	@echo ""

 # -----------------------------------------------------------------------------
@@ -294,7 +273,7 @@ sim3d: $(REFINED_MESH_3D)
 #   - $(MESH_3D)         - 初始網格
 #   - $(BG_POS_3D)       - 背景場
 #   - $(REFINED_MESH_3D) - 細化後網格
-#   - $(RESULT_3D).msh/.tec - 模擬結果
+#   - $(RESULT_3D).tec      - 模擬結果
 # -----------------------------------------------------------------------------
 3d:
 	@echo "=========================================="
@@ -325,32 +304,17 @@ sim3d: $(REFINED_MESH_3D)
 	@echo "    初始網格:   $(MESH_3D)"
 	@echo "    細化網格:   $(REFINED_MESH_3D)"
 	@echo "    背景場:     $(BG_POS_3D)"
-	@echo "    結果檔案:   $(RESULT_3D).msh, $(RESULT_3D).tec"
+	@echo "    結果檔案:   $(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 - 清除所有生成檔案
 # =============================================================================
 clean:
 	@echo ">>> [clean] 清除所有生成檔案..."
-	@echo "    rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh"
-	rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh
+	@echo "    rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh *.log"
+	rm -f *.msh *.tec *.pos *.vtu *.vtm *.vtk current_mesh*.msh *.log
+	@echo "    rm -rf __pycache__ *.pyc"
+	rm -rf __pycache__ *.pyc
 	@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,38 @@ 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 = 10.0e-4   # 10 um
+MESH_SIZE_MAX_3D = 2.0e-3   # 200 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
+# E-field thresholds (2D)
+E_VERY_HIGH_2D = 1.0e4   # 10 kV/cm  -> 1x
+E_HIGH_2D = 5.0e3        # 5 kV/cm   -> 2x
+E_MEDIUM_2D = 1.0e3      # 1 kV/cm   -> 4x
+E_LOW_2D = 5.0e2         # 500 V/cm  -> 8x
+E_VERYLOW_2D = 1.0e2     # 100 V/cm  -> 16x

-
-# --- 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 (3D)
+E_VERY_HIGH_3D = 2.5e3   # 2.5 kV/cm -> 1x
+E_HIGH_3D = 2.0e3        # 2 kV/cm   -> 2x
+E_MEDIUM_3D = 1.5e3      # 1.5 kV/cm -> 4x
+E_LOW_3D = 1.0e3         # 1 kV/cm   -> 8x
+E_VERYLOW_3D = 5.0e2     # 500 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,24 +126,33 @@ 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)."""
-        if self.dimension == 3:
-            print("--- Creating 3D Element Models ---")
-            for region in self.regions:
-                try:
-                    element_from_edge_model(device=self.name, region=region, edge_model="ElectricField")
-                except:
-                    pass
+        """Create element models for visualization (Emag, etc.)."""
+        from devsim import element_model
+        
+        print(f"--- Creating Element Models ({self.dimension}D) ---")
+        for region in self.regions:
+            try:
+                # Convert edge ElectricField to element model
+                element_from_edge_model(device=self.name, region=region, edge_model="ElectricField")
+                
+                # Calculate Emag (electric field magnitude)
+                if self.dimension == 2:
+                    element_model(device=self.name, region=region, name="Emag",
+                                  equation="(ElectricField_x^2 + ElectricField_y^2)^(0.5)")
+                else:  # 3D
+                    element_model(device=self.name, region=region, name="Emag",
+                                  equation="(ElectricField_x^2 + ElectricField_y^2 + ElectricField_z^2)^(0.5)")
+                
+                print(f"  {region}: Emag created")
+            except Exception as e:
+                print(f"  {region}: Failed - {e}")
                    
    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}.tec", type="tecplot")
+        print(f"Done: {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,56 +1,42 @@
 #!/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,
-            ifelse(Emag > {opto_common.E_MEDIUM}, 4.0,
-              ifelse(Emag > {opto_common.E_LOW}, 8.0,
-                if(Emag > {opto_common.E_VERYLOW}, 16)))))
+        ifelse(Emag > {opto_common.E_VERY_HIGH_2D}, 1.0,
+          ifelse(Emag > {opto_common.E_HIGH_2D}, 2.0,
+            ifelse(Emag > {opto_common.E_MEDIUM_2D}, 4.0,
+              ifelse(Emag > {opto_common.E_LOW_2D}, 8.0,
+                if(Emag > {opto_common.E_VERYLOW_2D}, 16)))))
    ''')
    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,77 +1,57 @@
 #!/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,
-            ifelse(Emag > {opto_common.E_MEDIUM}, 4.0,
-              ifelse(Emag > {opto_common.E_LOW}, 8.0,
-                if(Emag > {opto_common.E_VERYLOW}, 16)))))
+        ifelse(Emag > {opto_common.E_VERY_HIGH_3D}, 1.0,
+          ifelse(Emag > {opto_common.E_HIGH_3D}, 2.0,
+            ifelse(Emag > {opto_common.E_MEDIUM_3D}, 4.0,
+              ifelse(Emag > {opto_common.E_LOW_3D}, 8.0,
+                if(Emag > {opto_common.E_VERYLOW_3D}, 16)))))
    ''')
    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()
@@ -10,8 +10,7 @@ Mesh.Algorithm = 6;
 /* Scale & Mesh */
 sf = 0.1;                         // 1 mm = 0.1 cm
 um = 1.0e-4;                      // 1 µm = 1e-4 cm
-lc_fine = 5.0 * um;
-lc_base = 50.0 * um;
+lc_coarse = 50.0 * um;

 /* Dimensions */
 W_cu = 3.5*sf;    H_cu = 0.2*sf;      cu_ext = 1.0*sf;
@@ -62,12 +61,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 +135,8 @@ 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 === */
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthMax = lc_coarse;

-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;
@@ -17,8 +17,7 @@ Mesh.Algorithm3D = 10;
 /* Scale & Mesh */
 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_coarse = 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 +212,6 @@ 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 === */
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthMax = lc_coarse;
@@ -77,6 +77,7 @@ def run_simulation():
    # Solve and output
    device.solve()
    device.print_results()
+    device.create_element_models()
    device.output(output_prefix)
Author	SHA1	Message	Date
Dean Huang	be5107731d	feat: Stabilizing High-Voltage BJT Simulation	2026-03-27 16:13:23 +08:00
Dean Huang	90c0da0194	feat: Remove .msh file generation for Paraview	2026-03-19 16:30:19 +08:00
Dean Huang	e0f941ad13	feat: Add PROJECT GUIDE.md.	2026-03-16 14:48:02 +08:00
Dean Huang	0d65a7a0d1	refactor: 更新程式碼以及Makefile	2026-03-16 09:51:13 +08:00
Dean Huang	00afb79e86	docs: 更新註解與README。	2026-03-13 16:01:52 +08:00
Dean Huang	5df4b3432b	feat: 新增靜態直流偏壓分析，更新 README 和 Makefile。	2026-03-13 15:13:51 +08:00
Dean Huang	dd44ab08c4	feat: 新增並更新了多個3D繪圖相關圖檔。	2025-12-29 16:33:20 +08:00
Dean Huang	fe66120592	feat: 新增 P 與 E 的 3D 繪圖	2025-12-29 16:25:41 +08:00
Dean Huang	be38fce34a	Update	2025-12-29 15:07:45 +08:00
Dean Huang	3be76542f6	Update	2025-12-29 12:31:47 +08:00
Dean Huang	d75b08dcc2	Update	2025-12-29 12:23:32 +08:00
Dean Huang	b92e0ae19d	Update	2025-12-29 12:09:32 +08:00
Dean Huang	a664fca506	Update	2025-12-29 11:29:22 +08:00
Dean Huang	eb849b3e39	Update	2025-12-29 10:58:29 +08:00
Dean Huang	b767114c1d	Update	2025-12-29 10:28:21 +08:00
Dean Huang	f1299dadeb	feat: 調整 3D 網格最大尺寸並區分 2D/3D 電場閾值。	2025-12-26 22:02:32 +08:00
Dean Huang	8f12dc4bc1	Merge branch 'refacter'	2025-12-26 16:14:16 +08:00
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