Coverage for dragonfly_sim/utils/meshwarping_utils.py: 21%

1import numpy as np 

2 

3def inner_bdry_function_from_corner_list(ffd_block_corner_list, overset_margin=1e-4, rel_tol=1e-6): 

4 """ 

5 Build an inner-boundary indicator function out of an FFD block corner list. 

6 

7 ffd_block_corner_list carries the same information an inner_bdry_function 

8 encodes: the inner boundary facets are the ones contained in the FFD block, 

9 so the block itself can serve as the containment test. The returned closure 

10 has the signature Mesh expects: it takes x of shape (n_dimensions, n_points) 

11 and returns a boolean mask of the points inside the block. 

12 

13 The block is reconstructed exactly as in 

14 shape_parameterization.ffd_corners_from_corner_list -- two surfaces offset 

15 along a single connecting direction, linearly interpolated in between -- so a swept 

16 or tapered block is tested against its actual tapered extent rather than its 

17 bounding box. 

18 

19 overset_margin matches the one the control point construction applies, so the 

20 test region is the FFD block itself rather than the bare corner list; rel_tol 

21 widens it further by that fraction of the block extent in each direction, so 

22 that geometry sitting exactly on the block boundary still counts as inside. 

23 

24 NOTE the test is purely geometric: any OTHER boundary passing through the FFD 

25 block (e.g. the symmetry plane at the root of a wing) is flagged as well, 

26 whereas a hand-written inner_bdry_function can exclude it by construction. 

27 """ 

28 num_surfaces = len(ffd_block_corner_list) 

29 assert num_surfaces == 2, "Expected exactly 2 surface corner pairs, got {}".format(num_surfaces) 

30 ndim = len(ffd_block_corner_list[0][0]) 

31 

32 # Per-surface min/max in each coordinate, as in the control point construction 

33 surface_mins = np.zeros((num_surfaces, ndim)) 

34 surface_maxs = np.zeros((num_surfaces, ndim)) 

35 for s, (corner_a, corner_b) in enumerate(ffd_block_corner_list): 

36 corner_a, corner_b = np.array(corner_a), np.array(corner_b) 

37 surface_mins[s] = np.minimum(corner_a, corner_b) 

38 surface_maxs[s] = np.maximum(corner_a, corner_b) 

39 

40 surface_midpoints = 0.5 * (surface_mins + surface_maxs) 

41 coord_diffs = np.abs(surface_midpoints[1] - surface_midpoints[0]) 

42 connecting_dir = int(np.argmax(coord_diffs)) 

43 span_dirs = [d for d in range(ndim) if d != connecting_dir] 

44 

45 # Same overset margin the control points are laid out with, so that the 

46 # containment test covers the whole FFD block 

47 for s in range(num_surfaces): 

48 for d in range(ndim): 

49 delta = surface_maxs[s, d] - surface_mins[s, d] 

50 if delta > 0: 

51 surface_mins[s, d] -= overset_margin * delta 

52 surface_maxs[s, d] += overset_margin * delta 

53 

54 # Tolerances, scaled per direction by the block's bounding box extent. A 

55 # direction the block is collapsed in falls back to the largest extent, so 

56 # its tolerance stays a real (finite) length. 

57 block_mins = np.minimum(surface_mins[0], surface_mins[1]) 

58 block_maxs = np.maximum(surface_maxs[0], surface_maxs[1]) 

59 extents = block_maxs - block_mins 

60 tol = rel_tol * np.where(extents > 0., extents, extents.max()) 

61 

62 connect_lo = surface_midpoints[0][connecting_dir] 

63 connect_hi = surface_midpoints[1][connecting_dir] 

64 connect_span = connect_hi - connect_lo 

65 t_tol = tol[connecting_dir] / abs(connect_span) 

66 

67 print("inner_bdry_function_from_corner_list: FFD block used as inner boundary test") 

68 print(" connecting direction: {}".format(connecting_dir)) 

69 print(" block min corner: {}".format(list(block_mins))) 

70 print(" block max corner: {}".format(list(block_maxs))) 

71 

72 def inner_bdry_function(x): 

73 """ 

74 Input: x is an array of size (n_dimensions, n_points) 

75 """ 

76 coords = np.asarray(x)[:ndim, :] 

77 

78 # Parametric position along the connecting direction; points outside 

79 # [0, 1] are beyond the two surfaces bounding the block 

80 t = (coords[connecting_dir, :] - connect_lo) / connect_span 

81 inside = (t >= -t_tol) & (t <= 1. + t_tol) 

82 

83 # Within the block, the span-direction bounds are the linear 

84 # interpolation of the two surfaces' bounds 

85 t_clipped = np.clip(t, 0., 1.) 

86 for d in span_dirs: 

87 lower = (1. - t_clipped) * surface_mins[0, d] + t_clipped * surface_mins[1, d] 

88 upper = (1. - t_clipped) * surface_maxs[0, d] + t_clipped * surface_maxs[1, d] 

89 inside &= (coords[d, :] >= lower - tol[d]) & (coords[d, :] <= upper + tol[d]) 

90 

91 return inside 

92 

93 return inner_bdry_function 

94 

95def find_local_surface_nodes(msh, facet_tags, surface_tag: int) -> np.ndarray: 

96 """Return *local* geometry indices of owned surface nodes (ghosts excluded). 

97 

98 Parameters 

99 ---------- 

100 msh : dolfinx.mesh.Mesh 

101 facet_tags : dolfinx.mesh.MeshTags | None 

102 surface_tag : int 

103 

104 Returns 

105 ------- 

106 np.ndarray, int64, shape (M_local,) 

107 """ 

108 

109 tdim = msh.topology.dim 

110 fdim = tdim - 1 

111 msh.topology.create_connectivity(fdim, 0) 

112 f_to_v = msh.topology.connectivity(fdim, 0) 

113 

114 node_set: set[int] = set() 

115 if facet_tags is not None: 

116 for f in facet_tags.find(surface_tag): 

117 node_set.update(f_to_v.links(int(f)).tolist()) 

118 else: 

119 msh.topology.create_entities(fdim) 

120 msh.topology.create_connectivity(fdim, tdim) 

121 f_to_c = msh.topology.connectivity(fdim, tdim) 

122 for f in range(msh.topology.index_map(fdim).size_local): 

123 if len(f_to_c.links(f)) == 1: 

124 node_set.update(f_to_v.links(f).tolist()) 

125 

126 n_local_geom = msh.geometry.index_map().size_local # excludes ghosts 

127 return np.array( 

128 sorted(idx for idx in node_set if idx < n_local_geom), 

129 dtype=np.int64, 

130 )