Извлечь PhysicalFace как 2D-сетку в фипы - PullRequest
Новая
       48

Извлечь PhysicalFace как 2D-сетку в фипы

1 голос
/ 14 июня 2019

Я создал сетку с помощью 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'не имеет атрибута "cellCenters".Если я переопределю сетку только до порядка выдавливания, я получу: «ValueError: слишком много значений для распаковки (ожидается 2)» из-за формы tSslice.Я благодарен за любую помощь

1 Ответ

1 голос
/ 14 июня 2019

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, поэтому они затемните лица, которые вас интересуют.

Честно говоря, я думаю, что лучший способ добиться среза с помощью трехмерных данных, как это, это использовать либо MayaviClient (см. сфера Кан-Хиллиарда и sphereDaemon.py для примера) или экспортируйте с помощью VTKViewer , а затем визуализируйте срезы данных с помощью инструмента, такого как ParaView , VisIt или MayaVi .

...