可以通过读取Gmsh中的“物理名称”在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

谢谢。我能够使用physicalCell标识符创建单独的网格。关于找到气体/氧化物界面的顶点，我能够使用以下代码中的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)]

Answer 2

以下是我最初提供的暴力解决方案的改进：

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仍然无法满足您的要求，因为气体网格的绝大多数面（外部或其他）都不在氧化物网格附近。您只需要重合的面。要了解这有多复杂，您需要大约98%的this code

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