КТМИАД 1
Построить треугольную сетку по варианту (рис. 8)
Подготовка
class Mesh:
def __init__(self, vert, nobsx, nobsy, sratiox, sratioy):
self.points = self.__quadrangulate(nobsx, nobsy)
self.points_ids = self.__triangulate(nobsx, nobsy)
self.xs, self.ys = self.points[:, 0], self.points[:, 1]
self.triangulation = tri.Triangulation(self.points[:, 0], self.points[:, 1], triangles=self.points_ids)
self.triangles = self.triangulation.triangles
def __quadrangulate(self, nobsx: np.ndarray, nobsy: np.ndarray) -> np.ndarray:
xpoint_num = np.sum(nobsx) + 1 - nobsx.shape[0]
ypoint_num = np.sum(nobsy) + 1 - nobsy.shape[0]
points_shape = (ypoint_num * xpoint_num, 2)
points = np.zeros(points_shape)
for y_region_num, y_region in enumerate(nobsy):
for x_region_num, x_region in enumerate(nobsx):
verts = vert[x_region_num : 2 + x_region_num, y_region_num : y_region_num + 2]
vert_directions = (verts[:, 1] - verts[:, 0]) / (y_region - 1)
left_vert_pts = verts[0, 0] + np.arange(y_region)[:, None] * vert_directions[0]
right_vert_pts = verts[1, 0] + np.arange(y_region)[:, None] * vert_directions[1]
horiz_directions = (right_vert_pts - left_vert_pts) / (x_region - 1)
raw_points = left_vert_pts[:, :, None] + horiz_directions[:, :, None] * np.arange(x_region)
local_points = np.zeros((y_region, x_region, 2))
local_points[..., 0] = raw_points[:, 0, :]
local_points[..., 1] = raw_points[:, 1, :]
x_start_id = np.sum(nobsx[:x_region_num]) - x_region_num
x_end_id = x_start_id + x_region
y_start_id = np.sum(nobsy[:y_region_num]) - y_region_num
y_end_id = y_start_id + y_region
x_ids = np.arange(x_start_id, x_end_id)
y_ids = np.arange(y_start_id, y_end_id) * xpoint_num
local_points_ids = x_ids + y_ids[:, None]
points[local_points_ids, :] = local_points
return points
def __triangulate(self, nobsx: np.ndarray, nobsy: np.ndarray) -> np.ndarray:
xpoint_num = np.sum(nobsx) + 1 - nobsx.shape[0]
ypoint_num = np.sum(nobsy) + 1 - nobsy.shape[0]
i = np.arange(ypoint_num - 1)
j = np.arange(xpoint_num - 1)
points_left_bottom_ids = (j + i[:, None] * xpoint_num).flatten()
points_left_top_ids = (j + (i[:, None] + 1) * xpoint_num).flatten()
points_ids = np.zeros((2 * (xpoint_num - 1) * (ypoint_num - 1), 3))
even_elems = np.arange(points_ids.shape[0]) % 2 == 0
points_ids[even_elems, 0] = points_left_bottom_ids
points_ids[even_elems, 1] = points_left_top_ids
points_ids[even_elems, 2] = points_left_top_ids + 1
odd_elems = np.arange(points_ids.shape[0]) % 2 == 1
points_ids[odd_elems, 0] = points_left_bottom_ids
points_ids[odd_elems, 1] = points_left_bottom_ids + 1
points_ids[odd_elems, 2] = points_left_top_ids + 1
return points_ids
def __repr__(self) -> str:
repr = f"<Mesh{PAD}triangles: {lo(self.triangles.shape[0])}{PAD}points: {lo(self.points)}{PAD}xs: {lo(self.xs)}{PAD}ys: {lo(self.ys)}>"
return repr
def __str__(self) -> str:
return repr(self)
def show_calc_area(mesh: Mesh, figsize: Tuple[int, int], show_ticks: bool = True):
fig, ax = plt.subplots(figsize=figsize)
ax.set_aspect('equal')
if show_ticks:
plt.xticks(np.unique(mesh.xs))
plt.yticks(np.unique(mesh.ys))
ax.triplot(mesh.triangulation, 'bo-', lw=1)
fig.tight_layout()
plt.show()
Построение сетки
Задание вертикальных линий (координатами XY)
# from left to right
vert = np.array([
[[0, 0], [0, 0.2], [0, 1]],
[[1e-5, 0], [1e-5, 0.2], [1e-5, 1]],
[[3, 0], [3, 1e-4], [3, 1]],
[[7, 0], [7, 1e-4], [7, 1]],
[[10-(1e-4), 0], [10-(1e-4), 0.2], [10-(1e-4), 1]],
[[10, 0], [10, 0.2], [10, 1]]
], dtype=np.float64)
nobsx, sratiox = np.array([4, 3, 2, 3, 4]), np.array([1, 1, 1, 1, 1])
nobsy, sratioy = np.array([5, 3]), np.array([1, 1])
mesh = Mesh(vert, nobsx, nobsy, sratiox, sratioy)
mesh
Визуализация сетки
show_calc_area(mesh, figsize=(15, 8))