example3d_0.py - master - chukharev/geometry3d

example3d_0.py
526 строк24 KB
chukharev
Introduce end-of-line normalization as in git-doc/gitattributes.html
20 июл 2026, 08:57
20 июл 2026, 08:57a56a0a1
#!/usr/bin/env python
# -*- coding: utf8 -*-

"""A usage example for 3D modeling library."""

from __future__ import annotations

import mayavi
import numpy as np
from mayavi.modules.surface import Surface
from mayavi.scripts import mayavi2
from mayavi.sources.poly_data_reader import PolyDataReader
from stl import mesh

from FormPARAMS import (B_width, C_width, F_max_x, F_max_y, Fname,
                        Half, Height, L1_cx, L1_cy, L1_r, L_length,
                        N1_x0, N1_x1, N1_y0, N1_y1, N2_x0, N2_x1, N2_y0, N2_y1,
                        P2_x, P2_y, P3_x, P3_y, Parity, S_length, Xx, Yy, show_params, x_step, y_step)
from geometry3d import BumpPit as BuP
from geometry3d import Collect as Col
from geometry3d import HalfPipes as HaP
from geometry3d import LipUpDown as Lip
from geometry3d import LogError3D as LoE
from geometry3d import Point3D as Pnt
from geometry3d import SideNose as SiN
from geometry3d import Vector3D as V3d

__version__ = "$Id$"
__all__ = ["main"]


def in_limits(xx_: float, yy_: float, siz: float = 0.0) -> bool:
    """Check if the coordinates are in limits"""
    return (-F_max_x + siz <= xx_ <= F_max_x - siz) and (-F_max_y + siz <= yy_ <= F_max_y - siz) and (
            (abs(xx_) - L1_cx) ** 2 + (abs(yy_) - L1_cy) ** 2 >= (L1_r + siz) ** 2)


def fill_dispenser(x0: float) -> Col.Collect:
    """Create the main part of the surface"""
    dx_, dy_ = 2 * V3d.Vector(y_step, 0, 0), 2 * V3d.Vector(0, x_step, 0)
    assert isinstance(dx_, V3d.Vector) and isinstance(dy_, V3d.Vector)
    siz: float = B_width + L_length  # sum of Bump/Pit and Lip sizes
    pit: BuP.Pit | BuP.Bump
    bump: BuP.Pit | BuP.Bump
    if Parity == "Even":
        pit, bump = BuP.Pit(Pnt.Point(0, 0, 0)), BuP.Bump(Pnt.Point(0, 0, 0))
    else:
        pit, bump = BuP.Bump(Pnt.Point(0, 0, 0)), BuP.Pit(Pnt.Point(0, 0, 0))

    xx_: float
    yy_: float
    xx_, yy_ = x0 + dx_.x / 2.0, 0
    row0 = Col.Collect(pit.shifted_copy(V3d.Vector(xx_, yy_, 0)))
    row1 = Col.Collect(pit.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
    yy_ += dy_.y
    while in_limits(xx_, yy_, siz):
        row0.append(bump.shifted_copy(V3d.Vector(xx_, yy_, 0)))
        row0.prepend(bump.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
        row1.append(bump.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
        row1.prepend(bump.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
        yy_ += dy_.y
        if in_limits(xx_, yy_, siz):
            row0.append(pit.shifted_copy(V3d.Vector(xx_, yy_, 0)))
            row0.prepend(pit.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
            row1.append(pit.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
            row1.prepend(pit.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
            yy_ += dy_.y
    parts = Col.Collect(row0, row1)

    xx_, yy_ = xx_ + dx_.x, dy_.y / 2
    while in_limits(xx_, yy_, siz):
        row0 = Col.Collect(bump.shifted_copy(V3d.Vector(xx_, yy_, 0)),
                           bump.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
        row1 = Col.Collect(bump.shifted_copy(V3d.Vector(-xx_, yy_, 0)),
                           bump.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
        yy_ += dy_.y
        while in_limits(xx_, yy_, siz):
            row0.append(pit.shifted_copy(V3d.Vector(xx_, yy_, 0)))
            row0.prepend(pit.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
            row1.append(pit.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
            row1.prepend(pit.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
            yy_ += dy_.y
            if in_limits(xx_, yy_, siz):
                row0.append(bump.shifted_copy(V3d.Vector(xx_, yy_, 0)))
                row0.prepend(bump.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
                row1.append(bump.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
                row1.prepend(bump.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
                yy_ += dy_.y
        parts.append(row0)
        parts.prepend(row1)
        yy_ = 0
        xx_ += dx_.x

        if in_limits(xx_, yy_, siz):
            row0 = Col.Collect(pit.shifted_copy(V3d.Vector(xx_, yy_, 0)))
            row1 = Col.Collect(pit.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
            yy_ += dy_.y
            while in_limits(xx_, yy_, siz):
                row0.append(bump.shifted_copy(V3d.Vector(xx_, yy_, 0)))
                row0.prepend(bump.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
                row1.append(bump.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
                row1.prepend(bump.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
                yy_ += dy_.y
                if in_limits(xx_, yy_, siz):
                    row0.append(pit.shifted_copy(V3d.Vector(xx_, yy_, 0)))
                    row0.prepend(pit.shifted_copy(V3d.Vector(xx_, -yy_, 0)))
                    row1.append(pit.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
                    row1.prepend(pit.shifted_copy(V3d.Vector(-xx_, -yy_, 0)))
                    yy_ += dy_.y
            parts.append(row0)
            parts.prepend(row1)
            yy_ = dy_.y / 2
            xx_ += dx_.x
    xx_ += dy_.y - dx_.x
    yy_ = 0
    if in_limits(xx_, yy_, siz):
        parts.append(bump.shifted_copy(V3d.Vector(xx_, yy_, 0)))
        parts.append(bump.shifted_copy(V3d.Vector(-xx_, yy_, 0)))
    return parts


def make_part() -> Col.Collect:
    """Create a surface of a part"""
    # parts = fill_dispenser(0.0)
    # parts = fill_dispenser(x_step/2.0)
    parts = Col.Collect(HaP.HalfCylinderDown(Pnt.O_point))
    for i in range(1, 8, 2):
        parts.append(HaP.HalfCylinderUp(Pnt.Point(0, i * C_width, 0)))
        parts.append(HaP.HalfCylinderUp(Pnt.Point(0, -i * C_width, 0)))
        parts.append(HaP.HalfCylinderDown(Pnt.Point(0, (i + 1) * C_width, 0)))
        parts.append(HaP.HalfCylinderDown(Pnt.Point(0, -(i + 1) * C_width, 0)))
    parts.append(fill_dispenser(parts.containing_box_xyz()[1].x))
    parts.append(SiN.Side(Pnt.Point(P3_x - L_length, (P3_y - S_length), 0), Pnt.Point(P2_x + L_length, P2_y, 0)))
    parts.prepend(SiN.Side(Pnt.Point(P3_x - L_length, -P3_y, 0), Pnt.Point(P2_x + L_length, S_length - P2_y, 0)))
    parts.append(Lip.LipUp(Pnt.Point(L1_cx, L1_cy, 0), L1_r + L_length, L1_r))
    parts.append(Lip.LipUp(Pnt.Point(-L1_cx, -L1_cy, 0), L1_r + L_length, L1_r))
    parts.append(Lip.LipDown(Pnt.Point(-L1_cx, L1_cy, 0), L1_r + L_length, L1_r))
    parts.append(Lip.LipDown(Pnt.Point(L1_cx, -L1_cy, 0), L1_r + L_length, L1_r))
    parts.append(SiN.Nose(Pnt.Point(N1_x0, N1_y0, 0), Pnt.Point(N1_x1, N1_y1, 0)))
    parts.append(SiN.Nose(Pnt.Point(N2_x0, N2_y0, 0), Pnt.Point(N2_x1, N2_y1, 0)))
    return parts


def tessellate_ribbon(a_points: np.ndarray, b_points: np.ndarray, norm: str = "down") -> np.ndarray:
    """Generate a (part of) mesh from two arrays of points. The points should be arranged so that any line connecting
    a_points[i] to any of a_points[i+1], b_points[i] or b_points[i+1] as well as connecting b_points[i] to b_points[i+1]
    completely belongs to the part described by the mesh. Additionally, the triangles should not overlap.
    """
    assert norm in ("up", "down")
    if norm == 'up':
        a_points, b_points = b_points, a_points
    mm, nn = a_points.shape, b_points.shape
    assert mm[1] == nn[1] == 3, f"{mm=}, {nn=}."
    m0: int = mm[0]
    n0: int = nn[0]
    if m0 != n0:
        LoE.log.info(f"{m0=}, {n0=}, {norm=}\n{a_points=}\n{b_points=}")
    assert m0 == n0 or m0 == n0 + 1 or m0 + 1 == n0, f"{m0=}, {n0=}."
    facet_num: int = 2 * (m0 - 1)  # at maximum this number of facets
    data: np.ndarray = np.zeros(facet_num, dtype=mesh.Mesh.dtype)
    facet_num = 0  # count actual number of facets
    for ii in range(min(m0, n0) - 1):
        variant = 0
        if ii + 1 < n0 and not np.isnan((a_points[ii, 2], b_points[ii, 2], b_points[ii + 1, 2])).any():
            data['vectors'][facet_num] = np.array((a_points[ii], b_points[ii], b_points[ii + 1]))
            facet_num += 1
            variant = 1
        if variant != 1 and ii + 1 < m0 and not np.isnan((a_points[ii, 2], b_points[ii, 2], a_points[ii + 1, 2])).any():
            data['vectors'][facet_num] = np.array((a_points[ii], b_points[ii], a_points[ii + 1]))
            facet_num += 1
            variant = 2
        if ii + 1 >= min(m0, n0):
            continue
        if variant != 2 and not np.isnan((a_points[ii, 2], b_points[ii + 1, 2], a_points[ii + 1, 2])).any():
            data['vectors'][facet_num] = np.array((a_points[ii], b_points[ii + 1], a_points[ii + 1]))
            facet_num += 1
        if variant != 1 and not np.isnan((a_points[ii + 1, 2], b_points[ii, 2], b_points[ii + 1, 2])).any():
            data['vectors'][facet_num] = np.array((a_points[ii + 1], b_points[ii], b_points[ii + 1]))
            facet_num += 1
    LoE.log.info(f"{data.shape=}, {facet_num=} in {__name__}.")
    return data[:facet_num]


def tessellate_array(x_arr: np.ndarray, y_arr: np.ndarray, z_arr: np.ndarray, norm: str = "down") -> np.ndarray:
    """Generate a mesh from a rectangle 2D array. The coordinates of a point at indexes (i,j) are
    (x_arr[i], y_arr[j], z_arr[i, j]).
    """
    assert norm in ("up", "down")
    data: list[np.ndarray] = []
    a_line: np.ndarray = np.array([(x_arr[iii], y_arr[0], z_arr[iii, 0]) for iii in range(x_arr.size - 1)])
    for jjj in range(1, y_arr.size - 1):
        b_line: np.ndarray
        b_line = np.array([(x_arr[iii], y_arr[jjj], z_arr[iii, jjj]) for iii in range(x_arr.size - 1)])
        data.append(tessellate_ribbon(a_line, b_line, norm=norm))
        a_line = b_line
    return np.concatenate(data)


def is_on_edge(z_arr: np.ndarray, ind: np.ndarray) -> bool:  # ind contains 2 ints
    """Return True is the point belongs to body and a neighboring point does not."""
    if not np.isfinite(z_arr[ind[0]][ind[1]]):
        return False  # outside the part body, not on its edge
    max_ix: int
    max_iy: int
    max_ix, max_iy = z_arr.shape
    if ind[0] == max_ix - 1 or ind[0] == 0 or ind[1] == max_iy - 1 or ind[1] == 0:
        return True  # edge of the mesh is edge of the part body
    inner: int
    outer: int
    inner, outer = 0, 0  # count number of neighboring points of these types
    for i_step, j_step in ((0, 1), (1, 0), (1, 1), (0, -1), (-1, 0), (-1, -1), (1, -1), (-1, 1)):
        if np.isfinite(z_arr[ind[0] + i_step][ind[1] + j_step]):
            inner += 1
        else:
            outer += 1
        if outer > 1 and inner > 2:
            return True  # return early, if possible
    return False


def clockwise45(vec: tuple[int, int]) -> tuple[int, int]:
    """Return turned 45 degrees clockwise unit vector on square grid"""
    if vec[0] == 1 and vec[1] == 0:
        return 1, -1
    if vec[0] == 1 and vec[1] == 1:
        return 1, 0
    if vec[0] == 0 and vec[1] == 1:
        return 1, 1
    if vec[0] == -1 and vec[1] == 1:
        return 0, 1
    if vec[0] == -1 and vec[1] == 0:
        return -1, 1
    if vec[0] == -1 and vec[1] == -1:
        return -1, 0
    if vec[0] == 0 and vec[1] == -1:
        return -1, -1
    if vec[0] == 1 and vec[1] == -1:
        return 0, -1
    raise ValueError(f"Should not happen, {vec=}.")


def edge_points(z_arr: np.ndarray, start: np.ndarray | None = None) -> list[np.ndarray]:
    """Return a list of index tuples of points neighboring object body limits
    We go around the body counterclockwise. Always first try to turn away from the body, i.e. 45 degrees clockwise
    from the direction on the last step. Then try the direction of the last step, and then all other directions."""
    if start is None:
        start = np.array((0, 0))
    max_ix: int
    max_iy: int
    max_ix, max_iy = z_arr.shape
    while not is_on_edge(z_arr, start):
        LoE.log.info(f"z[{start[0]:d},{start[1]:d}]: {z_arr[start[0]][start[1]]:g}")
        start += 1, 1
        assert start[0] < max_ix and start[1] < max_iy
    edge: list[np.ndarray] = [start]
    front: list[np.ndarray] = [start]
    beam: tuple[int, int] = (1, -1)
    for step in ((1, -1), (1, 1), (-1, -1), (0, -1), (1, 0)):
        if is_on_edge(z_arr, start + step):
            edge.append(start + step)
            front.append(start + step)
            beam = step
            break
    LoE.log.info(f"{front=}")
    assert len(front) > 1
    while len(front) > 1:
        point: np.ndarray = front.pop()

        # We go around the body counterclockwise. Always first try to turn away from the body, i.e. 45 degrees clockwise
        # from the direction on the last step. Then try the direction of the last step, and then all other directions.
        for step in clockwise45(beam), beam, (0, 1), (1, 0), (1, 1), (0, -1), (-1, 0), (-1, -1), (1, -1), (-1, 1):
            pnt: np.ndarray = point + step
            tpl = tuple(pnt)
            if (any((ex == pnt[0] and ey == pnt[1] for ex, ey in edge)) or
                    not (0 <= tpl[0] < max_ix) or not (0 <= tpl[1] < max_iy)):
                continue
            if is_on_edge(z_arr, pnt):
                edge.append(pnt)
                front.append(pnt)
                beam = step
                break
    edge.append(edge[0])
    LoE.log.info(f"edge_points: {len(edge):d}: {edge!r}")
    return edge


def segmentate(points: list[np.ndarray]) -> list[list[np.ndarray]]:
    """Divide given list into list of lists of linear (horizontal, vertical, or diagonal) parts"""
    ind: int = 0
    res: list[list[np.ndarray]] = []
    while ind < len(points) - 1:
        ind += 1
        vec: np.ndarray = np.array(points[ind]) - points[ind - 1]
        segment: list[np.ndarray] = [points[ind - 1], points[ind]]
        while (ind < len(points) - 1) and (vec + segment[-1] == np.array(points[ind + 1])).all():
            segment.append(points[ind + 1])
            ind += 1
        res.append(segment)
    vec0: np.ndarray = np.array(res[0][1]) - res[0][0]
    vec_l: np.ndarray = np.array(res[-1][1]) - res[-1][0]
    if res[0][0][0] == res[-1][-1][0] and res[0][0][1] == res[-1][-1][1] and vec_l.all() == vec0.all():
        segment = res.pop()
        res[0] = segment[:-1] + res[0]
    return res


def tessellate_2_n(two_points: np.ndarray, n_points: np.ndarray, norm: str = "down", show: bool = False) -> np.ndarray:
    """Generate a (part of) mesh from two arrays of points. The points should be arranged so that any line connecting
    a_points[i] to any of a_points[i+1] or b_points[i] or b_points[i+1] as well as connecting b_points[i] to
    b_points[i+1] completely belongs to the part described by the mesh. Additionally, the triangles should not overlap.
    """
    assert norm in ("up", "down")
    mm: tuple[int, int] = two_points.shape
    nn: tuple[int, int] = n_points.shape
    assert (mm[0] == 2) and (nn[0] >= 2) and (mm[1] == nn[1] == 3)
    data: np.ndarray = np.zeros(nn[0], dtype=mesh.Mesh.dtype)
    if norm == 'up':
        two_points, n_points = two_points[::-1], n_points[::-1]
    mid: int = nn[0] // 2
    facet_num: int = 0  # count actual number of facets
    for ind in range(mid):
        data['vectors'][facet_num] = np.array((two_points[0], n_points[ind], n_points[ind + 1]))
        facet_num += 1
    for ind in range(mid, nn[0] - 1):
        data['vectors'][facet_num] = np.array((two_points[1], n_points[ind], n_points[ind + 1]))
        facet_num += 1
    data['vectors'][facet_num] = np.array((two_points[0], n_points[mid], two_points[1]))
    facet_num += 1
    if show:
        LoE.log.info(f"t2n: {data['vectors'].shape}, {facet_num}, {data}.")
    return data[:]


def monotone_parts(ij_array: list[np.ndarray]) -> tuple[tuple[list[np.ndarray], list[np.ndarray]],
                                                        tuple[list[np.ndarray], list[np.ndarray]],
                                                        tuple[list[np.ndarray], list[np.ndarray]],
                                                        tuple[list[np.ndarray], list[np.ndarray]]]:
    """Collect monotone parts of the bottom edge line. Variable names come from where in average the part is aimed:
    nne - towards North-North-East, i.e. strictly monotone towards +Y (North) and non-negative towards +X (East);
    wws - towards West-West-South, etc.
    The borderline starts at the left end, at middle height, and goes counterclockwise around the part.
    The corners are concave, so it first goes ESE, then SSE, NNE, ENE, WNW, NNW, SSW, WSW."""
    nne: list[np.ndarray]; ene: list[np.ndarray]; wnw: list[np.ndarray]; nnw: list[np.ndarray]
    ssw: list[np.ndarray]; wsw: list[np.ndarray]; ese: list[np.ndarray]; sse: list[np.ndarray]
    nne, ene, wnw, nnw = [], [], [], []
    ssw, wsw, ese, sse = [], [], [], []

    ind: int
    ind, max_ind = 0, len(ij_array)
    di, dj = ij_array[1][0] - ij_array[0][0], ij_array[1][1] - ij_array[0][1]
    assert di > 0 >= dj
    while ind < max_ind and di > 0 >= dj:
        ese.append(ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"ese done:, {ind}, {len(ese)}, {ese}")

    sse.insert(0, ij_array[ind - 1])
    while ind < max_ind and di >= 0 > dj:
        sse.insert(0, ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"sse done:, {ind}, {len(sse)}, {sse}")

    nne.append(ij_array[ind])
    ind += 1
    di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    while ind < max_ind and di >= 0 and dj > 0:
        nne.append(ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"nne done:, {ind}, {len(nne)}, {nne}")

    ene.insert(0, ij_array[ind - 1])
    while ind < max_ind and di > 0 and dj >= 0:
        ene.insert(0, ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"ene done:, {ind}, {len(ene)}, {ene}")

    wnw.append(ij_array[ind])
    ind += 1
    di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    while ind < max_ind and di < 0 <= dj:
        wnw.append(ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"wnw done:, {ind}, {len(wnw)}, {wnw}")

    nnw.insert(0, ij_array[ind - 1])
    while ind < max_ind and di <= 0 < dj:
        nnw.insert(0, ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"nnw done:, {ind}, {len(nnw)}, {nnw}")

    ssw.append(ij_array[ind])
    ind += 1
    di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    while ind < max_ind and di <= 0 and dj < 0:
        ssw.append(ij_array[ind])
        ind += 1
        di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"ssw done:, {ind}, {len(ssw)}, {ssw}")

    wsw.insert(0, ij_array[ind - 1])
    while ind < max_ind and di < 0 and dj <= 0:
        wsw.insert(0, ij_array[ind])
        ind += 1
        if ind < max_ind:
            di, dj = ij_array[ind][0] - ij_array[ind - 1][0], ij_array[ind][1] - ij_array[ind - 1][1]
    LoE.log.info(f"wsw done:, {ind}, {len(wsw)}, {wsw}")

    # Move a few points between neighboring monotone chains
    while len(ese) > len(wsw) and len(wsw) < len(ssw):
        assert all(ssw[-1] == wsw[-1]), f"{ssw[-1]}, {wsw[-1]}."
        del ssw[-1]
        wsw.append(ssw[-1])
    while len(wnw) > len(ene) and len(ene) < len(nne):
        assert all(nne[-1] == ene[-1]), f"{nne[-1]}, {ene[-1]}."
        del nne[-1]
        ene.append(nne[-1])

    return (wsw, ese), (sse, nne), (ene, wnw), (nnw, ssw)


def numpy2mesh(xx: np.ndarray, yy: np.ndarray, zz: np.ndarray) -> mesh.Mesh:
    """Generate a mesh structure from m-element xx, n-element yy and m*n-element z-array."""

    data0: np.ndarray = tessellate_array(xx, yy, zz)

    # Select indexes of edge points
    ih_points: list[np.ndarray] = edge_points(zz, start=np.array((0, yy.size // 2)))
    segments: list[list[np.ndarray]] = segmentate(ih_points)
    LoE.log.info(f"segments: {len(segments)}, {segments}.")
    data1: list[np.ndarray] = []
    bottom: list[np.ndarray] = []
    aa: int = 0
    for segment in segments:
        high_points = np.array([(xx[x_ind], yy[y_ind], zz[x_ind][y_ind]) for x_ind, y_ind in segment])
        low_points = np.array([(xx[x_ind], yy[y_ind], 0.0) for x_ind, y_ind in (segment[0], segment[-1])])
        if aa < 3:
            aa += 1
            LoE.log.info(f"segment: {len(segment)=}: {segment}.")
            LoE.log.info(f"high: {len(high_points)=}: {high_points}.")
            LoE.log.info(f"low: {len(low_points)=}: {low_points}.")
            data1.extend(tessellate_2_n(low_points, high_points, show=True))
        else:
            data1.extend(tessellate_2_n(low_points, high_points))
            # break
        bottom.append(segment[-1])
    LoE.log.info(f"bottom: {len(bottom)}, {bottom}.")
    bottom_segments = segmentate(bottom)
    LoE.log.info(f"bottom_segments: {len(bottom_segments)}, {bottom_segments}.")
    len_bottom_segments = len(bottom_segments)
    assert len(bottom) - 1 == len_bottom_segments, f"{len(bottom)=} {len_bottom_segments=:d}"

    mono_lines = monotone_parts(bottom)
    data2: list[np.ndarray] = []
    for ii, (a_line, b_line) in enumerate(mono_lines):
        directions = (("wsw", "ese"), ("sse", "nne"), ("ene", "wnw"), ("nnw", "ssw"))
        LoE.log.info(f'{directions[ii][0]} a_line: {len(a_line)}, {a_line}')
        LoE.log.info(f'{directions[ii][1]} b_line: {len(b_line)}, {b_line}')
        a_points = np.array([(xx[x_ind], yy[y_ind], 0.0) for x_ind, y_ind in a_line])
        b_points = np.array([(xx[x_ind], yy[y_ind], 0.0) for x_ind, y_ind in b_line])
        data2.extend(tessellate_ribbon(a_points, b_points))
    b_points = np.array([(xx[x_ind], yy[y_ind], 0.0) for x_ind, y_ind in (mono_lines[0][1][-1], mono_lines[2][0][-1])])
    a_points = np.array([(xx[x_ind], yy[y_ind], 0.0) for x_ind, y_ind in (mono_lines[0][0][-1], mono_lines[2][1][-1])])
    data2.extend(tessellate_ribbon(a_points, b_points))
    data: np.ndarray = data0
    np.concatenate((data, data1, data2))
    # np.append(data, data2)
    return mesh.Mesh(data)
    # return mesh.Mesh(np.concatenate([data0, data1]))


@mayavi2.standalone
def view(fname: str | None = None) -> None:
    """Show an STL file on the screen"""
    if fname is None:
        fname = Fname
    mayavi.new_scene()
    source: PolyDataReader = PolyDataReader()
    source.initialize(fname)
    mayavi.add_source(source)
    s: Surface = Surface()
    mayavi.add_module(s)


def main() -> None:
    """Create few touching surfaces and show them"""

    LoE.log.info(f"Start example3d_0.py.")
    LoE.log.info(f"{show_params()}\n\n\n")
    part: Col.Collect = make_part()
    zzz: np.ndarray = np.zeros((Xx.size, Yy.size))
    iii: int
    jjj: int
    xxx: float
    yyy: float
    if Half == "Low":
        for jjj, yyy in enumerate(Yy):
            for iii, xxx in enumerate(Xx):
                zzz[iii][jjj] = Height + part(xxx, yyy) if in_limits(xxx, yyy) else np.nan
    else:  # "Top"
        for jjj, yyy in enumerate(Yy[::-1]):
            for iii, xxx in enumerate(Xx):
                zzz[iii][jjj] = Height - part(xxx, yyy) if in_limits(xxx, yyy) else np.nan
    LoE.log.info(f"{zzz.shape=}")

    mesh1: mesh.Mesh = numpy2mesh(Xx, Yy, zzz)
    LoE.log.info(f"{mesh1.dtype=}, {mesh1.data.shape=}.")
    mesh1.save(Fname)
    LoE.log.info(f"File {Fname} saved, starting view().")

    # mesh1.save(Fname, mode=stl.Mode.ASCII)

    # stl.numpy2stl(zzz, Fname, scale=1, force_python=True, solid=True,
    #               max_width=Length, max_depth=Width, max_height=Height)
    # mask_val=-10.0, max_width=Length, max_depth=Width, max_height=Height, solid=True)
    #    view()
    LoE.log.info(f"The end.")


if __name__ == "__main__":
    main()