8. DMPlex 和 Mesh#
准备并行环境
import ipyparallel as ipp
cluster = ipp.Cluster(profile="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
在串行环境中导入必要的包
from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance
8.1. 创建 DMPlex#
8.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_0x84000002_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 (3 (2), 0 (4), 1 (5))
8.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_0x84000002_1 1 MPI process
type: plex
DM_0x84000002_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), 3 (2), 1 (5))
depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))
8.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) dim 0 offset 0
( 1) dim 0 offset 0
( 2) dim 1 offset 0
( 3) dim 1 offset 1
( 4) dim 1 offset 2
( 5) dim 1 offset 3
( 6) dim 0 offset 4
( 7) dim 0 offset 4
( 8) dim 0 offset 4
( 9) dim 0 offset 4
( 10) dim 0 offset 4
PetscSection Object: 1 MPI process
type not yet set
Process 0:
( 0) dim 0 offset 0
( 1) dim 0 offset 0
( 2) dim 1 offset 3
( 3) dim 1 offset 0
( 4) dim 1 offset 2
( 5) dim 1 offset 1
( 6) dim 0 offset 4
( 7) dim 0 offset 4
( 8) dim 0 offset 4
( 9) dim 0 offset 4
( 10) dim 0 offset 4
8.3. Firedrake 中的 Mesh#
mesh = RectangleMesh(4, 4, 1, 1)
mesh.init()
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)
mesh.init()
plex = mesh.topology_dm
rank, size = mesh.comm.rank, mesh.comm.size
%%px --block
mesh._dm_renumbering.indices
Out[0:3]:
array([15, 63, 64, 62, 27, 29, 14, 58, 11, 57, 56, 25, 10, 52, 7, 51, 50,
23, 6, 46, 3, 45, 44, 21, 2, 43, 28, 60, 26, 13, 61, 59, 12, 55,
54, 24, 9, 53, 8, 49, 48, 22, 5, 47, 4, 42, 41, 20, 1, 40, 0,
39, 38, 19, 75, 76, 74, 36, 37, 18, 72, 73, 71, 34, 35, 17, 69, 70,
68, 32, 33, 66, 30, 16, 67, 65, 31], dtype=int32)
Out[1:3]:
array([ 0, 39, 40, 41, 20, 21, 23, 1, 42, 43, 24, 2, 44, 45, 4, 48, 49,
26, 3, 47, 5, 50, 51, 27, 6, 52, 8, 55, 56, 29, 7, 54, 9, 57,
58, 30, 10, 59, 12, 62, 63, 32, 11, 61, 13, 64, 65, 33, 14, 66, 15,
68, 22, 46, 25, 53, 28, 60, 31, 67, 34, 16, 69, 70, 35, 17, 71, 72,
36, 18, 73, 74, 37, 19, 75, 76, 38], 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[1:4]:
(array([ 0, 1, 2, 16, 4, 3, 5, 6, 17, 8, 7, 9, 10, 18, 12, 11, 13,
14, 15, 19, 20, 21, 23, 24, 22, 35, 25, 26, 27, 36, 28, 29, 30, 37,
31, 32, 33, 34, 38, 39, 40, 41, 42, 43, 44, 45, 69, 70, 46, 48, 49,
47, 50, 51, 52, 71, 72, 53, 55, 56, 54, 57, 58, 59, 73, 74, 60, 62,
63, 61, 64, 65, 66, 67, 68, 75, 76], dtype=int32),
array([ 0, 1, 2, 16, 4, 3, 5, 6, 17, 8, 7, 9, 10, 18, 12, 11, 13,
14, 15, 19, 20, 21, 23, 24, 22, 35, 25, 26, 27, 36, 28, 29, 30, 37,
31, 32, 33, 34, 38, 39, 40, 41, 42, 43, 44, 45, 69, 70, 46, 48, 49,
47, 50, 51, 52, 71, 72, 53, 55, 56, 54, 57, 58, 59, 73, 74, 60, 62,
63, 61, 64, 65, 66, 67, 68, 75, 76], dtype=int32))
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, 29, 28, 26, 38, 36, 37, 25, 24, 34, 35, 23, 22, 32,
33, 21, 20, 30, 31, 63, 64, 62, 58, 60, 61, 59, 75, 76, 74, 57, 56,
55, 52, 54, 53, 72, 73, 71, 51, 50, 49, 46, 48, 47, 69, 70, 68, 45,
44, 42, 43, 41, 40, 39, 66, 67, 65], dtype=int32),
array([15, 14, 13, 19, 11, 12, 10, 9, 18, 7, 8, 6, 5, 17, 3, 4, 2,
1, 0, 16, 27, 29, 28, 26, 38, 36, 37, 25, 24, 34, 35, 23, 22, 32,
33, 21, 20, 30, 31, 63, 64, 62, 58, 60, 61, 59, 75, 76, 74, 57, 56,
55, 52, 54, 53, 72, 73, 71, 51, 50, 49, 46, 48, 47, 69, 70, 68, 45,
44, 42, 43, 41, 40, 39, 66, 67, 65], 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: firedrake_default_topology 2 MPI processes
type: plex
firedrake_default_topology 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))
Face Sets: 3 strata with value/size (1 (5), 2 (7), 4 (9))
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)
8.4. FunctionSpace#
下图使用 Sphinx / Jupyter Book 图形插件 plantuml 制作.
Fig. 8.1 class FunctionSpace#
8.5. 示例: 寻找特定的几何实体以及节点编号#
8.5.1. 二维示例#
使用 global_numbering 寻找某条边界上端点, 以及在这条边界上与它相邻的点 (2D mesh)
from firedrake import *
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), 3 (162), 1 (259))
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), 4 (17), 2 (17), 3 (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 *
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 61
Number of 1-cells per rank: 150 150
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 (7), 4 (11))
exterior_facets: 1 strata with value/size (1 (33))
interior_facets: 1 strata with value/size (1 (195))
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 24: [1. 1.]), node 21: [0.875 1. ]
8.5.2. 三维示例#
使用 global_numbering 寻找界面上与接触线相邻的三角形 (3D mesh)
本示例网格文件 cylinder.msh 由几何文件 cylinder.geo 生成,
如图: 
import matplotlib.pyplot as plt
import numpy as np
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), 6 (488), 3 (1106), 1 (774))
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 (7 (230), 3 (209), 5 (209), 6 (230), 4 (209))
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
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()
[stdout:0] DM Object: firedrake_default_topology 2 MPI processes
type: plex
firedrake_default_topology in 3 dimensions:
Number of 0-cells per rank: 106 104
Number of 1-cells per rank: 471 468
Number of 2-cells per rank: 639 638
Number of 3-cells per rank: 273 273
Labels:
depth: 4 strata with value/size (0 (106), 1 (471), 2 (639), 3 (273))
celltype: 4 strata with value/size (0 (106), 1 (471), 3 (639), 6 (273))
Cell Sets: 2 strata with value/size (1 (140), 2 (133))
Face Sets: 5 strata with value/size (3 (125), 4 (133), 5 (123), 6 (137), 7 (127))
Edge Sets: 1 strata with value/size (1 (8))
exterior_facets: 1 strata with value/size (1 (467))
interior_facets: 1 strata with value/size (1 (1053))
pyop2_core: 1 strata with value/size (1 (744))
pyop2_owned: 1 strata with value/size (1 (470))
pyop2_ghost: 1 strata with value/size (1 (275))
[output:1]
[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:0]
[output:1]
