您可以通过从 Gmsh 读取 "physical names" 在 FiPy 中创建单独的网格吗？

Question

我有以下代码代表 Gmsh 的 .geo 文件，我想在 FiPy 中加载它以进行过程模拟。它包含：(i) 氧化物，(ii) 硅和 (iii) 气体作为需要识别为单独网格的物理实体。

SetFactory("OpenCASCADE");
//
ns = 1e-1;
ns2 = 1e-2;
//
x1 = 0;
x2 = 1;
y1 = 0;
y2 = 0.5;
Point(1) = {x1, y1, 0, ns};
Point(2) = {x2, y1, 0, ns};
Point(3) = {x2, y2, 0, ns2};
Point(4) = {x1, y2, 0, ns2};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Physical Surface("Oxide") = {1};
//
y1 = 0.5;
y2 = 1.0;
shiftX = 0.4;
Point(51) = {x1, y1, 0, ns2};
Point(61) = {x1+shiftX, y1, 0, ns2};
Point(71) = {x1+shiftX, y2, 0, ns2};
Point(81) = {x1, y2, 0, ns2};
Line(51) = {51, 61};
Line(61) = {61, 71};
Line(71) = {71, 81};
Line(81) = {81, 51};
Curve Loop(21) = {51, 61, 71, 81};
Plane Surface(21) = {21};
Physical Surface("Silicon") = {21};
//
Point(52) = {x2-shiftX, y1, 0, ns2};
Point(62) = {x2, y1, 0, ns2};
Point(72) = {x2, y2, 0, ns2};
Point(82) = {x2-shiftX, y2, 0, ns2};
Line(52) = {52, 62};
Line(62) = {62, 72};
Line(72) = {72, 82};
Line(82) = {82, 52};
Curve Loop(22) = {52, 62, 72, 82};
Plane Surface(22) = {22};
Physical Surface("Silicon") += {22};
//
y1 = 1.0;
y2 = 1.5;
shiftX = 0.4;
shiftY = 0.5;
Point(9) = {x1, y1, 0, ns2};
Point(91) = {x1+shiftX, y1, 0, ns2};
Point(92) = {x1+shiftX, y1-shiftY, 0, ns2};
Point(101) = {x2-shiftX, y1-shiftY, 0, ns2};
Point(102) = {x2-shiftX, y1, 0, ns2};
Point(10) = {x2, y1, 0, ns2};
Point(11) = {x2, y2, 0, ns};
Point(12) = {x1, y2, 0, ns};
Line(9) = {9, 91};
Line(91) = {91, 92};
Line(92) = {92, 101};
Line(101) = {101, 102};
Line(102) = {102, 10};
Line(10) = {10, 11};
Line(11) = {11, 12};
Line(12) = {12, 9};
Curve Loop(3) = {9, 91, 92, 101, 102, 10, 11, 12};
Plane Surface(3) = {3};
Physical Surface("Gas") = {3};
//

有没有办法在 FiPy 中读取它并创建 3 个具有给定物理名称的网格？

Answer 1

这里是对我最初提供的蛮力解决方案的改进：

import fipy as fp

mesh = fp.Gmsh2D('''
SetFactory("OpenCASCADE");
//
ns = 1e-1;
ns2 = 1e-2;
//
x1 = 0;
x2 = 1;
y1 = 0;
y2 = 0.5;
Point(1) = {x1, y1, 0, ns};
Point(2) = {x2, y1, 0, ns};
Point(3) = {x2, y2, 0, ns2};
Point(4) = {x1, y2, 0, ns2};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Physical Surface("Oxide") = {1};
//
y1 = 0.5;
y2 = 1.0;
shiftX = 0.4;
Point(51) = {x1, y1, 0, ns2};
Point(61) = {x1+shiftX, y1, 0, ns2};
Point(71) = {x1+shiftX, y2, 0, ns2};
Point(81) = {x1, y2, 0, ns2};
Line(51) = {51, 61};
Line(61) = {61, 71};
Line(71) = {71, 81};
Line(81) = {81, 51};
Curve Loop(21) = {51, 61, 71, 81};
Plane Surface(21) = {21};
Physical Surface("Silicon") = {21};
//
Point(52) = {x2-shiftX, y1, 0, ns2};
Point(62) = {x2, y1, 0, ns2};
Point(72) = {x2, y2, 0, ns2};
Point(82) = {x2-shiftX, y2, 0, ns2};
Line(52) = {52, 62};
Line(62) = {62, 72};
Line(72) = {72, 82};
Line(82) = {82, 52};
Curve Loop(22) = {52, 62, 72, 82};
Plane Surface(22) = {22};
Physical Surface("Silicon") += {22};
//
y1 = 1.0;
y2 = 1.5;
shiftX = 0.4;
shiftY = 0.5;
Point(9) = {x1, y1, 0, ns2};
Point(91) = {x1+shiftX, y1, 0, ns2};
Point(92) = {x1+shiftX, y1-shiftY, 0, ns2};
Point(101) = {x2-shiftX, y1-shiftY, 0, ns2};
Point(102) = {x2-shiftX, y1, 0, ns2};
Point(10) = {x2, y1, 0, ns2};
Point(11) = {x2, y2, 0, ns};
Point(12) = {x1, y2, 0, ns};
Line(9) = {9, 91};
Line(91) = {91, 92};
Line(92) = {92, 101};
Line(101) = {101, 102};
Line(102) = {102, 10};
Line(10) = {10, 11};
Line(11) = {11, 12};
Line(12) = {12, 9};
Curve Loop(3) = {9, 91, 92, 101, 102, 10, 11, 12};
Plane Surface(3) = {3};
Physical Surface("Gas") = {3};
//''')

from fipy.meshes.mesh2D import Mesh2D

def extract_mesh(mesh, mask):
    cellFaceIDs = mesh.cellFaceIDs[..., mask]
    faceIDs = numerix.unique(cellFaceIDs.flatten())
    facemap = numerix.zeros(mesh.faceVertexIDs.shape[1], dtype=int)
    facemap[faceIDs] = faceIDs.argsort()
    
    faceVertexIDs = mesh.faceVertexIDs[..., faceIDs]
    vertIDs = numerix.unique(faceVertexIDs.flatten())
    vertmap = numerix.zeros(mesh.vertexCoords.shape[1], dtype=int)
    vertmap[vertIDs] = vertIDs.argsort()

    return Mesh2D(mesh.vertexCoords[..., vertIDs],
                  vertmap[faceVertexIDs],
                  facemap[cellFaceIDs])
    
mesh_gas = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Gas"])
mesh_silicon = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Silicon"])
mesh_oxide = extract_mesh(mesh=mesh, mask=mesh.physicalCells["Oxide"])

派生的网格具有紧凑的顶点和面列表以及合理数量的外表面。

nearest 仍然不会做您想要的，因为气体网格的绝大多数面（外部或其他）都不在氧化物网格附近。您只需要重合的面孔。要了解这有多复杂，您需要 this code.

的大约 98%

我认为更好的答案是定义您感兴趣的边界，使用 Gmsh GEO 脚本，Physical Line 定义。然后可以将这些提取为 mesh.physicalFaces（需要对 extract_mesh() 做一些进一步的工作，但我认为它很容易处理）。

Answer 2

谢谢jeguyer。我能够使用 physicalCell 标识符创建单独的网格。关于在 gas/oxide 界面找到顶点，我能够使用以下代码中的 surfCoords 追踪界面轮廓：

gasFaces = gasMesh.exteriorFaces
gasVertices = numerix.unique(gasMesh.faceVertexIDs[..., gasFaces].flatten())
gasVertexCoords = gasMesh.vertexCoords[..., gasVertices]

oxideFaces = oxideMesh.exteriorFaces
oxideVertices = numerix.unique(oxideMesh.faceVertexIDs[..., oxideFaces].flatten())
oxideVertexCoords = oxideMesh.vertexCoords[..., oxideVertices]

surfVertexID = numerix.nearest(gasVertexCoords, oxideVertexCoords)
surfCoords = gasVertexCoords[..., numerix.unique(surfVertexID)]

您可以通过从 Gmsh 读取 "physical names" 在 FiPy 中创建单独的网格吗？

Can you create separate meshes in FiPy by reading "physical names" from Gmsh?

python

fipy

scikit-fmm