1

我用 Gmsh 创建了一个网格（一个有孔的表面，然后挤压它）。现在，我想在使用 MatplotlibViewer （Mayavi 在我的两台计算机上都不起作用）进行模拟之后将模型绘制在单独的切片中。我曾希望可以使用mesh.physicalFaces 定义一个新网络，但如果可以的话，我还没有想通。我的第二次尝试是使用 Gmsh 再次将网格应用到 Extrude 命令。但网格与 3D 版本的网格不对应。有人可以给我一个线索吗？也像另一种表示方法。

我正在开发 Win10、Fipy 3.1.3、Python 3.6

import numpy as np
from fipy import *
#%%
def func_mesh():
    mesh = Gmsh3D('''
    Geometry.OCCAutoFix = 0;
    SetFactory("OpenCASCADE");

    x = 1.;
    bseg = 0.08;
    bs= bseg*x;
    ls = 2.1; 
    cl = 0.01;
    radius = 0.006;

    // Exterior (bounding box) of mesh
    Point(1) = {0, 0, 0, cl};
    Point(2) = {0, bs, 0, cl};
    Point(4) = { bs,  0, 0, cl};
    Point(3) = {bs,  bs, 0, cl};
    Line(1) = {1, 2};
    Line(2) = {2, 3};
    Line(3) = {3, 4};
    Line(4) = {4, 1};
    Line Loop (21) = {1,2,3,4};


    //Circle
    Point(5) = {bseg/2 - radius, bseg/2, 0, cl};
    Point(6) = {bseg/2, bseg/2 + radius, 0, cl};
    Point(7) = { bseg/2 + radius, bseg/2, 0, cl};
    Point(8) = {bseg/2, bseg/2 - radius, 0, cl};
    Point(9) =  {bseg/2, bseg/2, 0, cl};

    Circle(10) = {5,9,6};
    Circle(11) = {6,9,7};
    Circle(12) = {7,9,8};
    Circle(13) = {8,9,5};
    Line Loop(22) = {10,11,12,13};
    Plane Surface(40) = {22}; //cycl

    Plane Surface(15) = {21, 22}; //Surface with a hole


    id[] = Extrude {0, 0, ls} {Surface{15}; Layers{210}; Recombine;};

    Surface Loop(2) = {46, 45, 48, 47, 49, 41, 44, 43, 42, 15};

    Physical Volume("Vol") = {id[]};

    Physical Surface("surf_ges") = {41, 42, 43, 44, 49, 47, 45, 48, 46, 15};

    Physical Surface("HX") = {45, 46, 48, 47};

    Physical Surface("Extr") = {15};

    ''')
    return mesh

    mesh = func_mesh()

    x,y,z = mesh.cellCenters 
    X,Y,Z = mesh.faceCenters

    tS = CellVariable(name="storage", 
                      mesh=mesh, 
                      value=367., 
                      hasOld=True)

submesh = mesh.physicalFaces['Extr']
xsub, ysub = submesh.cellCenters
tSslice = CellVariable(name = 'tSsclice',
                 mesh = submesh,
                 value = tS[z== z[0]])
viewer = MatplotlibViewer(vars = tSslice)

此尝试的错误消息是：AttributeError: 'binOp' object has no attribute 'cellCenters'。如果我只重新定义网格，直到挤出顺序我得到：“ValueError：解包的值太多（预期为 2）”，因为 tSslice 的形状。我很感激任何帮助

4

1 回答 1

2

submesh不是一个Mesh；它是一个掩码，用于标识表面mesh中包含哪些面Extr。

FiPy 无法从另一个网格中提取网格。Mesh2D使用submesh掩码、mesh.vertexCoords和来创建一个应该是可行的mesh.faceVertexIDs，但这并非易事。

从理论上讲，您可以调用Gmsh2D最多 的所有内容Plane Surface(15) = {21, 22}; //Surface with a hole，但我发现生成的元素数量与您的 3D 切片的数量不同z == z[0]。

啊哈，我看到了这个问题。我以为Extrude手术会产生棱柱形细胞，但事实并非如此。细胞呈四面体。由于 的单元格mesh并非都具有相同的四面体几何形状，因此其基点位于 的单元格Extr并不都保证其中心位于 处z == z[0]。更好的方法是使用 FiPy 的CellVariable插值来提取 的tS坐标 的值tSslice：

from fipy import *

geo = '''
Geometry.OCCAutoFix = 0;
SetFactory("OpenCASCADE");

x = 1.;
bseg = 0.08;
bs= bseg*x;
ls = 2.1; 
cl = 0.01;
radius = 0.006;

// Exterior (bounding box) of mesh
Point(1) = {0, 0, 0, cl};
Point(2) = {0, bs, 0, cl};
Point(4) = { bs,  0, 0, cl};
Point(3) = {bs,  bs, 0, cl};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop (21) = {1,2,3,4};


//Circle
Point(5) = {bseg/2 - radius, bseg/2, 0, cl};
Point(6) = {bseg/2, bseg/2 + radius, 0, cl};
Point(7) = { bseg/2 + radius, bseg/2, 0, cl};
Point(8) = {bseg/2, bseg/2 - radius, 0, cl};
Point(9) =  {bseg/2, bseg/2, 0, cl};

Circle(10) = {5,9,6};
Circle(11) = {6,9,7};
Circle(12) = {7,9,8};
Circle(13) = {8,9,5};
Line Loop(22) = {10,11,12,13};
Plane Surface(40) = {22}; //cycl

Plane Surface(15) = {21, 22}; //Surface with a hole


id[] = Extrude {0, 0, ls} {Surface{15}; Layers{210}; Recombine;};

Surface Loop(2) = {46, 45, 48, 47, 49, 41, 44, 43, 42, 15};

Physical Volume("Vol") = {id[]};

Physical Surface("Extr") = {15};
'''

mesh = Gmsh3D(geo + '''
Physical Surface("surf_ges") = {41, 42, 43, 44, 49, 47, 45, 48, 46, 15};

Physical Surface("HX") = {45, 46, 48, 47};
'''
)
submesh = Gmsh2D(geo)

x,y,z = mesh.cellCenters 
X,Y = submesh.cellCenters[..., submesh.physicalCells['Extr']]
Z = numerix.ones(X.shape) * z[0]

tS = CellVariable(name="storage", 
                  mesh=mesh, 
                  value=mesh.x * mesh.y * mesh.z, 
                  hasOld=True)

tSslice = CellVariable(name = 'tSsclice',
                 mesh = submesh)

# interpolate values of tS at positions of tSslice
tSslice[..., submesh.physicalCells['Extr']] = tS(numerix.vstack([X, Y, Z]))

viewer = MatplotlibViewer(vars = tSslice)

在这里，我使用相同的.geo脚本来定义mesh和submesh。我只将surf_ges和HX物理表面添加到 的定义中mesh，否则所有这些面也将被导入submesh，尽管有效z值为 0，因此它们会掩盖您感兴趣的面。

坦率地说，我认为通过这样的 3D 数据实现切片的更好方法是使用自定义MayaviClient（参见Cahn-Hilliard 球体和sphereDaemon.py示例）或使用VTKViewer导出，然后使用类似工具渲染数据切片ParaView、VisIt或Mayavi。

于 2019-06-14T13:58:44.637 回答