6. DMPlex 和 Mesh DMPlex and Mesh#
DMPlex 的基本用法参见 DMPlex.
Overview
How Firedrake meshes are represented by PETSc DMPlex: creating a DMPlex, PETSc.Section, the relation between DMPlex and Firedrake Mesh/FunctionSpace, and locating geometric entities and their DoF numbering.
准备并行环境
import ipyparallel as ipp
cluster = ipp.Cluster(engines="mpi", n=2)
client = cluster.start_and_connect_sync()
Starting 2 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>
%%px --block
from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
import numpy as np
在串行环境中导入必要的包
from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
import numpy as np
6.1. 创建 DMPlex#
6.1.1. 底层创建方法#
plex = PETSc.DMPlex().create()
plex.setDimension(2)
plex.setChart(0, 11)
# plex.setConeSize(point, number of points that cover the point)
plex.setConeSize(0, 3)
plex.setConeSize(1, 3)
plex.setConeSize(6, 2)
plex.setConeSize(7, 2)
plex.setConeSize(8, 2)
plex.setConeSize(9, 2)
plex.setConeSize(10, 2)
plex = plex.setUp() # plex.setUp() return self
# plex.setCone(point, [points that cover the point])
plex.setCone(0, [6, 7, 8])
plex.setCone(1, [7, 9, 10])
plex.setCone(6, [2, 3])
plex.setCone(7, [3, 4])
plex.setCone(8, [4, 2])
plex.setCone(9, [4, 5])
plex.setCone(10, [5, 3])
plex.symmetrize()
plex.stratify()
plex.view()
DM Object: 1 MPI process
type: plex
DM_0x436b70c0_0 in 2 dimensions:
Number of 0-cells per rank: 4
Number of 1-cells per rank: 5
Number of 2-cells per rank: 2
Labels:
depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))
celltype: 3 strata with value/size (0 (4), 1 (5), 3 (2))
6.1.2. 使用高级接口#
cells = [
[0, 1, 2],
[1, 3, 2]
]
coords = [
[-1, 0],
[0, -1],
[0, 1],
[1, 0]
]
plex = PETSc.DMPlex().createFromCellList(
dim=2, cells=cells, coords=coords, interpolate=True, comm=None)
plex.view()
DM Object: DM_0x436b70c0_1 1 MPI process
type: plex
DM_0x436b70c0_1 in 2 dimensions:
Number of 0-cells per rank: 4
Number of 1-cells per rank: 5
Number of 2-cells per rank: 2
Labels:
celltype: 3 strata with value/size (0 (4), 1 (5), 3 (2))
depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))
6.2. PETSc.Section#
section1 = PETSc.Section().create()
section1.setChart(*plex.getChart())
section2 = PETSc.Section().create()
section2.setChart(*plex.getChart())
old_to_new = plex.getOrdering(PETSc.Mat.OrderingType.RCM).indices
reordering = np.empty_like(old_to_new) # reordering[new] -> old
reordering[old_to_new] = np.arange(old_to_new.size, dtype=old_to_new.dtype)
perm = PETSc.IS().createGeneral(reordering)
section2.setPermutation(perm)
ps, pe = plex.getDepthStratum(0)
for p in range(ps, pe):
section1.addDof(p, 1)
section2.addDof(p, 1)
section1.setUp()
section2.setUp()
section1.view()
section2.view()
PetscSection Object: 1 MPI process
type not yet set
Process 0:
( 0) dof 0 offset 0
( 1) dof 0 offset 0
( 2) dof 1 offset 0
( 3) dof 1 offset 1
( 4) dof 1 offset 2
( 5) dof 1 offset 3
( 6) dof 0 offset 4
( 7) dof 0 offset 4
( 8) dof 0 offset 4
( 9) dof 0 offset 4
( 10) dof 0 offset 4
PetscSection Object: 1 MPI process
type not yet set
Process 0:
( 0) dof 0 offset 0
( 1) dof 0 offset 0
( 2) dof 1 offset 3
( 3) dof 1 offset 0
( 4) dof 1 offset 2
( 5) dof 1 offset 1
( 6) dof 0 offset 4
( 7) dof 0 offset 4
( 8) dof 0 offset 4
( 9) dof 0 offset 4
( 10) dof 0 offset 4
6.3. Firedrake 中的 Mesh#
mesh = RectangleMesh(4, 4, 1, 1)
plex = mesh.topology_dm
rank, size = mesh.comm.rank, mesh.comm.size
%%px --block
from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
mesh = RectangleMesh(4, 4, 1, 1)
plex = mesh.topology_dm
rank, size = mesh.comm.rank, mesh.comm.size
%%px --block
mesh._dm_renumbering.indices
Out[1:3]:
array([ 0, 39, 57, 47, 21, 23, 20, 1, 61, 49, 24, 2, 40, 48, 4, 41, 58,
26, 3, 50, 5, 62, 51, 27, 6, 42, 8, 43, 59, 29, 7, 52, 9, 63,
53, 30, 10, 44, 12, 45, 60, 32, 11, 54, 13, 64, 55, 33, 14, 46, 15,
56, 22, 65, 25, 66, 28, 67, 31, 68, 34, 16, 69, 73, 35, 17, 70, 74,
36, 18, 71, 75, 37, 19, 72, 76, 38], dtype=int32)
Out[0:3]:
array([15, 46, 64, 56, 27, 29, 14, 54, 11, 44, 63, 25, 10, 52, 7, 42, 62,
23, 6, 50, 3, 40, 61, 21, 2, 48, 28, 60, 26, 13, 45, 55, 12, 53,
59, 24, 9, 43, 8, 51, 58, 22, 5, 41, 4, 49, 57, 20, 1, 39, 0,
47, 38, 19, 68, 76, 72, 37, 36, 18, 67, 75, 71, 35, 34, 17, 66, 74,
70, 33, 32, 73, 30, 16, 65, 69, 31], dtype=int32)
和重新编号相关的函数
%%px --block
# https://petsc.org/release/manualpages/DMPlex/DMPlexGetOrdering
old_to_new = plex.getOrdering(PETSc.Mat.OrderingType.RCM).indices
reordering = np.empty_like(old_to_new)
reordering[old_to_new] = np.arange(old_to_new.size, dtype=old_to_new.dtype)
mesh._default_reordering, reordering
Out[0:4]:
(array([15, 14, 13, 19, 11, 12, 10, 9, 18, 7, 8, 6, 5, 17, 3, 4, 2,
1, 0, 16, 27, 28, 29, 26, 38, 36, 37, 25, 24, 34, 35, 23, 22, 32,
33, 21, 20, 30, 31, 46, 64, 56, 60, 54, 45, 55, 68, 76, 72, 44, 63,
53, 59, 52, 43, 67, 75, 71, 42, 62, 51, 58, 50, 41, 66, 74, 70, 40,
61, 49, 57, 48, 39, 73, 47, 65, 69], dtype=int32),
array([15, 14, 13, 19, 11, 12, 10, 9, 18, 7, 8, 6, 5, 17, 3, 4, 2,
1, 0, 16, 27, 28, 29, 26, 38, 36, 37, 25, 24, 34, 35, 23, 22, 32,
33, 21, 20, 30, 31, 46, 64, 56, 60, 54, 45, 55, 68, 76, 72, 44, 63,
53, 59, 52, 43, 67, 75, 71, 42, 62, 51, 58, 50, 41, 66, 74, 70, 40,
61, 49, 57, 48, 39, 73, 47, 65, 69], dtype=int32))
Out[1:4]:
(array([ 0, 1, 2, 16, 4, 3, 5, 6, 17, 8, 7, 9, 10, 18, 12, 11, 13,
14, 15, 19, 21, 23, 20, 24, 22, 35, 25, 26, 27, 36, 28, 29, 30, 37,
31, 32, 33, 34, 38, 39, 57, 47, 61, 49, 40, 48, 69, 65, 73, 41, 58,
50, 62, 51, 42, 70, 66, 74, 43, 59, 52, 63, 53, 44, 71, 67, 75, 45,
60, 54, 64, 55, 46, 68, 56, 72, 76], dtype=int32),
array([ 0, 1, 2, 16, 4, 3, 5, 6, 17, 8, 7, 9, 10, 18, 12, 11, 13,
14, 15, 19, 21, 23, 20, 24, 22, 35, 25, 26, 27, 36, 28, 29, 30, 37,
31, 32, 33, 34, 38, 39, 57, 47, 61, 49, 40, 48, 69, 65, 73, 41, 58,
50, 62, 51, 42, 70, 66, 74, 43, 59, 52, 63, 53, 44, 71, 67, 75, 45,
60, 54, 64, 55, 46, 68, 56, 72, 76], dtype=int32))
%%px --block
mesh.topology_dm.view()
entity_classes = mesh._entity_classes
ec_msg = ', '.join([f'{i}-cells: {ec}' for i, ec in enumerate(entity_classes)])
PETSc.Sys.syncPrint(f'[{rank}/{size}] mesh._entity_classes (core, owned, ghost): {ec_msg}')
PETSc.Sys.syncFlush()
[stdout:0] DM Object: DM_0x25fc9210_0 2 MPI processes
type: plex
DM_0x25fc9210_0 in 2 dimensions:
Number of 0-cells per rank: 19 19
Number of 1-cells per rank: 38 38
Number of 2-cells per rank: 20 20
Labels:
depth: 3 strata with value/size (0 (19), 1 (38), 2 (20))
celltype: 3 strata with value/size (0 (19), 1 (38), 3 (20))
marker: 1 strata with value/size (1 (19))
Face Sets: 3 strata with value/size (2 (9), 3 (5), 4 (7))
exterior_facets: 1 strata with value/size (1 (19))
interior_facets: 1 strata with value/size (1 (47))
pyop2_core: 1 strata with value/size (1 (26))
pyop2_owned: 1 strata with value/size (1 (26))
pyop2_ghost: 1 strata with value/size (1 (25))
[0/2] mesh._entity_classes (core, owned, ghost): 0-cells: [ 5 10 19], 1-cells: [13 26 38], 2-cells: [ 8 16 20]
[1/2] mesh._entity_classes (core, owned, ghost): 0-cells: [10 15 19], 1-cells: [26 30 38], 2-cells: [16 16 20]
计算 node_set 相关的函数
计算 node_set 相关的概念
V1 = FunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'CG', 2)
PETSc.Sys.Print(f'V1 entity_dofs: {V1.finat_element.entity_dofs()}')
PETSc.Sys.Print(f'V2 entity_dofs: {V2.finat_element.entity_dofs()}')
nodes_per_entity_V1 = tuple(mesh.make_dofs_per_plex_entity(V1.finat_element.entity_dofs()))
nodes_per_entity_V2 = tuple(mesh.make_dofs_per_plex_entity(V2.finat_element.entity_dofs()))
PETSc.Sys.Print(f'V1 nodes_per_entity: {nodes_per_entity_V1}')
PETSc.Sys.Print(f'V2 nodes_per_entity: {nodes_per_entity_V2}')
V1 entity_dofs: {0: {0: [0], 1: [1], 2: [2]}, 1: {2: [], 1: [], 0: []}, 2: {0: []}}
V2 entity_dofs: {0: {0: [0], 1: [3], 2: [5]}, 1: {2: [1], 1: [2], 0: [4]}, 2: {0: []}}
V1 nodes_per_entity: (1, 0, 0)
V2 nodes_per_entity: (1, 1, 0)
%%px --block
V1 = FunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'CG', 2)
PETSc.Sys.Print(f'V1 entity_dofs: {V1.finat_element.entity_dofs()}')
PETSc.Sys.Print(f'V2 entity_dofs: {V2.finat_element.entity_dofs()}')
nodes_per_entity_V1 = tuple(mesh.make_dofs_per_plex_entity(V1.finat_element.entity_dofs()))
nodes_per_entity_V2 = tuple(mesh.make_dofs_per_plex_entity(V2.finat_element.entity_dofs()))
PETSc.Sys.Print(f'V1 nodes_per_entity: {nodes_per_entity_V1}')
PETSc.Sys.Print(f'V2 nodes_per_entity: {nodes_per_entity_V2}')
[stdout:0] V1 entity_dofs: {0: {0: [0], 1: [1], 2: [2]}, 1: {2: [], 1: [], 0: []}, 2: {0: []}}
V2 entity_dofs: {0: {0: [0], 1: [3], 2: [5]}, 1: {2: [1], 1: [2], 0: [4]}, 2: {0: []}}
V1 nodes_per_entity: (1, 0, 0)
V2 nodes_per_entity: (1, 1, 0)
6.4. FunctionSpace#
下图使用 Sphinx / Jupyter Book 图形插件 plantuml 制作.
Fig. 6.1 class FunctionSpace#
6.5. 示例: 寻找特定的几何实体以及节点编号#
6.5.1. 二维示例#
使用 global_numbering 寻找某条边界上端点, 以及在这条边界上与它相邻的点 (2D mesh)
from firedrake import *
from py.intro_utils import triplot
import matplotlib.pyplot as plt
import numpy as np
rectangle = Mesh("gmsh/rectangle.msh")
fig, axes = plt.subplots(figsize=[4, 3])
triplot(rectangle, axes=axes)
axes.set_aspect('equal')
rectangle.topology_dm.view()
DM Object: firedrake_default_topology 1 MPI process
type: plex
firedrake_default_topology in 2 dimensions:
Number of 0-cells per rank: 98
Number of 1-cells per rank: 259
Number of 2-cells per rank: 162
Labels:
celltype: 3 strata with value/size (0 (98), 1 (259), 3 (162))
depth: 3 strata with value/size (0 (98), 1 (259), 2 (162))
Cell Sets: 1 strata with value/size (1 (162))
Face Sets: 4 strata with value/size (1 (17), 2 (17), 3 (17), 4 (17))
exterior_facets: 1 strata with value/size (1 (64))
interior_facets: 1 strata with value/size (1 (325))
pyop2_core: 1 strata with value/size (1 (519))
pyop2_owned: 0 strata with value/size ()
pyop2_ghost: 0 strata with value/size ()
%%px --block
from firedrake import *
from py.intro_utils import triplot
import matplotlib.pyplot as plt
import numpy as np
rectangle = Mesh("gmsh/rectangle.msh")
fig, axes = plt.subplots(figsize=[4, 3])
triplot(rectangle, axes=axes)
axes.set_aspect('equal')
axes.set_xlim([-0.1, 1.1])
axes.set_ylim([-0.1, 1.1])
rectangle.topology_dm.view()
[stdout:0] DM Object: firedrake_default_topology 2 MPI processes
type: plex
firedrake_default_topology in 2 dimensions:
Number of 0-cells per rank: 61 60
Number of 1-cells per rank: 150 149
Number of 2-cells per rank: 90 90
Labels:
depth: 3 strata with value/size (0 (61), 1 (150), 2 (90))
celltype: 3 strata with value/size (0 (61), 1 (150), 3 (90))
Cell Sets: 1 strata with value/size (1 (90))
Face Sets: 3 strata with value/size (2 (17), 3 (11), 4 (9))
exterior_facets: 1 strata with value/size (1 (35))
interior_facets: 1 strata with value/size (1 (194))
pyop2_core: 1 strata with value/size (1 (190))
pyop2_owned: 1 strata with value/size (1 (60))
pyop2_ghost: 1 strata with value/size (1 (51))
[output:1]
[output:0]
def get_interface_element_with_contact_point(mesh, V, interface_tag, adj_line_tag):
dm = mesh.topology_dm
edge_label = dm.getLabel("Face Sets") # 2D Face is Edge
core_label = dm.getLabel("pyop2_core")
owned_label = dm.getLabel("pyop2_owned")
edge_label_values = edge_label.getValueIS().indices
interface_indices = []
if interface_tag in edge_label_values:
interface_indices = edge_label.getStratumIS(interface_tag).indices
adj_line = []
if adj_line_tag in edge_label_values:
adj_line = edge_label.getStratumIS(adj_line_tag).indices
points = np.intersect1d(interface_indices, adj_line)
plex_element = []
if len(points) > 0:
point = points[0]
support = dm.getSupport(point)
for edge in support:
if edge_label.getValue(edge) == interface_tag and \
(core_label.getValue(edge) == 1 or owned_label.getValue(edge) == 1):
cone = dm.getCone(edge)
adj_point = cone[1] if cone[0] == point else cone[0]
plex_element = [point, adj_point]
break
local_section = V.global_numbering # global_numbering is a local section
element = [local_section.getOffset(_) for _ in plex_element]
return element
# interface_tag, adj_line_tag = 2, 4 # 2: left boundary, 4: upper boundary
V = FunctionSpace(rectangle, 'CG', 1)
element = get_interface_element_with_contact_point(rectangle, V, interface_tag=1, adj_line_tag=3)
coords_data = rectangle.coordinates.dat.data_ro_with_halos # This must be outside the if condition (mpi collective)
rank, size = rectangle.comm.rank, rectangle.comm.size
if len(element) > 0:
coords = [coords_data[_] for _ in element]
PETSc.Sys.syncPrint(f"[{rank}/{size}] node {element[0]}: {coords[0]}), node {element[1]}: {coords[1]}")
PETSc.Sys.syncFlush()
[0/1] node 41: [0. 0.]), node 47: [0.125 0. ]
%%px --block
def get_interface_element_with_contact_point(mesh, V, interface_tag, adj_line_tag):
dm = mesh.topology_dm
edge_label = dm.getLabel("Face Sets") # 2D Face is Edge
core_label = dm.getLabel("pyop2_core")
owned_label = dm.getLabel("pyop2_owned")
edge_label_values = edge_label.getValueIS().indices
interface_indices = []
if interface_tag in edge_label_values:
interface_indices = edge_label.getStratumIS(interface_tag).indices
adj_line = []
if adj_line_tag in edge_label_values:
adj_line = edge_label.getStratumIS(adj_line_tag).indices
points = np.intersect1d(interface_indices, adj_line)
plex_element = []
if len(points) > 0:
point = points[0]
support = dm.getSupport(point)
for edge in support:
if edge_label.getValue(edge) == interface_tag and \
(core_label.getValue(edge) == 1 or owned_label.getValue(edge) == 1):
cone = dm.getCone(edge)
adj_point = cone[1] if cone[0] == point else cone[0]
plex_element = [point, adj_point]
break
local_section = V.global_numbering # global_numbering is a local section
element = [local_section.getOffset(_) for _ in plex_element]
return element
# interface_tag, adj_line_tag = 2, 4 # 2: left boundary, 4: upper boundary
V = FunctionSpace(rectangle, 'CG', 1)
element = get_interface_element_with_contact_point(rectangle, V, interface_tag=2, adj_line_tag=4)
coords_data = rectangle.coordinates.dat.data_ro_with_halos # This must be outside the if condition (mpi collective)
rank, size = rectangle.comm.rank, rectangle.comm.size
if len(element) > 0:
coords = [coords_data[_] for _ in element]
PETSc.Sys.syncPrint(f"[{rank}/{size}] node {element[0]}: {coords[0]}), node {element[1]}: {coords[1]}")
PETSc.Sys.syncFlush()
[stdout:0] [0/2] node 32: [1. 1.]), node 29: [0.875 1. ]
6.5.2. 三维示例#
使用 global_numbering 寻找界面上与接触线相邻的三角形 (3D mesh)
本示例网格文件 cylinder.msh 由几何文件 cylinder.geo 生成,
如图: 
import matplotlib.pyplot as plt
import numpy as np
from py.intro_utils import triplot
cylinder = Mesh("gmsh/cylinder.msh")
fig, axes = plt.subplots(figsize=[4, 3], subplot_kw={'projection': '3d'})
triplot(cylinder, axes=axes)
axes.set_aspect('equal')
cylinder.topology_dm.view()
DM Object: firedrake_default_topology 1 MPI process
type: plex
firedrake_default_topology in 3 dimensions:
Number of 0-cells per rank: 157
Number of 1-cells per rank: 774
Number of 2-cells per rank: 1106
Number of 3-cells per rank: 488
Labels:
celltype: 4 strata with value/size (0 (157), 1 (774), 3 (1106), 6 (488))
depth: 4 strata with value/size (0 (157), 1 (774), 2 (1106), 3 (488))
Cell Sets: 2 strata with value/size (1 (244), 2 (244))
Face Sets: 5 strata with value/size (3 (209), 4 (209), 5 (209), 6 (230), 7 (230))
Edge Sets: 1 strata with value/size (1 (16))
exterior_facets: 1 strata with value/size (1 (782))
interior_facets: 1 strata with value/size (1 (1745))
pyop2_core: 1 strata with value/size (1 (2525))
pyop2_owned: 0 strata with value/size ()
pyop2_ghost: 0 strata with value/size ()
%%px --block
import matplotlib.pyplot as plt
import numpy as np
from py.intro_utils import triplot
cylinder = Mesh("gmsh/cylinder.msh")
fig, axes = plt.subplots(figsize=[4, 3], subplot_kw={'projection': '3d'})
triplot(cylinder, axes=axes)
axes.set_aspect('equal')
cylinder.topology_dm.view()
[output:1]
[stdout:0] DM Object: firedrake_default_topology 2 MPI processes
type: plex
firedrake_default_topology in 3 dimensions:
Number of 0-cells per rank: 104 109
Number of 1-cells per rank: 468 476
Number of 2-cells per rank: 637 640
Number of 3-cells per rank: 272 272
Labels:
depth: 4 strata with value/size (0 (104), 1 (468), 2 (637), 3 (272))
celltype: 4 strata with value/size (0 (104), 1 (468), 3 (637), 6 (272))
Cell Sets: 2 strata with value/size (1 (136), 2 (136))
Face Sets: 5 strata with value/size (3 (125), 4 (129), 5 (129), 6 (129), 7 (129))
Edge Sets: 1 strata with value/size (1 (9))
exterior_facets: 1 strata with value/size (1 (459))
interior_facets: 1 strata with value/size (1 (1047))
pyop2_core: 1 strata with value/size (1 (742))
pyop2_owned: 1 strata with value/size (1 (470))
pyop2_ghost: 1 strata with value/size (1 (269))
[output:0]
def get_interface_element_include_contact_line(mesh, V, interface_tag, contact_line_tag):
dm = mesh.topology_dm
edge_label = dm.getLabel("Edge Sets")
face_label = dm.getLabel("Face Sets")
core_label = dm.getLabel("pyop2_core")
owned_label = dm.getLabel("pyop2_owned")
edge_label_values = edge_label.getValueIS().indices
contact_line = []
if contact_line_tag in edge_label_values:
contact_line = edge_label.getStratumIS(contact_line_tag).indices
faces = []
for seg in contact_line:
for face in dm.getSupport(seg):
if face_label.getValue(face) == interface_tag and \
(core_label.getValue(face) == 1 or owned_label.getValue(face) == 1):
faces.append(int(face))
break
plex_cell_node_map = np.zeros((len(faces), 3), dtype=np.int32)
for i, (seg, face) in enumerate(zip(contact_line, faces)):
seg_nodes = dm.getCone(seg)
plex_cell_node_map[i, :2] = seg_nodes
plex_cell_node_map[i, 2:] = np.setdiff1d(
np.unique(np.array([dm.getCone(_) for _ in dm.getCone(face)]).flatten()),
seg_nodes)
local_section = V.global_numbering # global_numbering is a local section
cell_node_map = np.zeros_like(plex_cell_node_map)
for i, cell in enumerate(plex_cell_node_map):
assert np.all(np.array([local_section.getDof(_) for _ in cell]) > 0)
cell_node_map[i, :] = [local_section.getOffset(_) for _ in cell]
return cell_node_map
CONTACT_LINE = 1
INTERFACE = 3
V = FunctionSpace(cylinder, 'CG', 1)
cell_node_map = get_interface_element_include_contact_line(cylinder, V, INTERFACE, CONTACT_LINE)
%%px --block
def get_interface_element_include_contact_line(mesh, V, interface_tag, contact_line_tag):
dm = mesh.topology_dm
edge_label = dm.getLabel("Edge Sets")
face_label = dm.getLabel("Face Sets")
core_label = dm.getLabel("pyop2_core")
owned_label = dm.getLabel("pyop2_owned")
edge_label_values = edge_label.getValueIS().indices
contact_line = []
if contact_line_tag in edge_label_values:
contact_line = edge_label.getStratumIS(contact_line_tag).indices
faces = []
for seg in contact_line:
for face in dm.getSupport(seg):
if face_label.getValue(face) == interface_tag and \
(core_label.getValue(face) == 1 or owned_label.getValue(face) == 1):
faces.append(int(face))
break
plex_cell_node_map = np.zeros((len(faces), 3), dtype=np.int32)
for i, (seg, face) in enumerate(zip(contact_line, faces)):
seg_nodes = dm.getCone(seg)
plex_cell_node_map[i, :2] = seg_nodes # set the seg nodes first
plex_cell_node_map[i, 2:] = np.setdiff1d(
np.unique(np.array([dm.getCone(_) for _ in dm.getCone(face)]).flatten()),
seg_nodes)
local_section = V.global_numbering # global_numbering is a local section
cell_node_map = np.zeros_like(plex_cell_node_map)
for i, cell in enumerate(plex_cell_node_map):
assert np.all(np.array([local_section.getDof(_) for _ in cell]) > 0)
cell_node_map[i, :] = [local_section.getOffset(_) for _ in cell]
return cell_node_map
CONTACT_LINE = 1
INTERFACE = 3
V = FunctionSpace(cylinder, 'CG', 1)
cell_node_map = get_interface_element_include_contact_line(cylinder, V, INTERFACE, CONTACT_LINE)
# plot the triangle to check if they are on the interface
coords = cylinder.coordinates.dat.data_ro_with_halos
if len(cell_node_map) > 0:
assert np.allclose(coords[:, 2][cell_node_map], 0.25)
fig, axes = plt.subplots(figsize=[4, 3])
c = axes.triplot(coords[:, 0], coords[:, 1], triangles=cell_node_map)
lines = [[(coords[_, 0], coords[_, 1]) for _ in __[:2] ] for __ in cell_node_map]
from matplotlib.collections import LineCollection
line_collection = LineCollection(lines, colors='k', linestyles=':')
axes.add_collection(line_collection)
axes.set_xlim([-0.52, 0.52])
axes.set_ylim([-0.52, 0.52])
axes.set_aspect("equal")
axes.grid("on")
%%px --block
# plot the triangle to check if they are on the interface
coords = cylinder.coordinates.dat.data_ro_with_halos
if len(cell_node_map) > 0:
assert np.allclose(coords[:, 2][cell_node_map], 0.25)
fig, axes = plt.subplots(figsize=[4, 3])
c = axes.triplot(coords[:, 0], coords[:, 1], triangles=cell_node_map)
lines = [[(coords[_, 0], coords[_, 1]) for _ in __[:2] ] for __ in cell_node_map]
from matplotlib.collections import LineCollection
line_collection = LineCollection(lines, colors='k', linestyles=':')
axes.add_collection(line_collection)
axes.set_xlim([-0.52, 0.52])
axes.set_ylim([-0.52, 0.52])
axes.set_aspect("equal")
axes.grid("on")
[output:1]
[output:0]
