Matematické Fórum

Zdravím,

mám problém s výpočtem působiště vztlakové síly. Naše zadání zní takto:
Vytvořit animaci Bézierovy plochy (vč. řídicích bodů), která několikrát změní svůj tvar, nakonec vytvoří trup lodi.

Konkrétní zadání:

- Vypočítat polohu působiště vztlakové síly působící na trup.

- Odevzdat animaci, zdrojový kód, výpočetní zprávu.

Animaci mám, ale to působiště je problém. Nejprve jsme počítali souřadnice jako jednotlivé váhy, ale toto řešení je podle učitele špatně. Máme to dělat přes momenty a že se z-ová souřadnice nedá vypočítat. Kód je napsán v Pythonu. Přijde mi, že to těžiště má být ve težišti ponořené části lodě. Mohl by mi s tím někdo poradit?
Kód: 
from math import comb
from typing import TypeAlias
import matplotlib.pyplot as plt
import numpy as np
import numpy.typing as npt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter

# Typ alias pro kontrolní body Bézierových ploch
ControlPoints: TypeAlias = npt.NDArray[np.float64]

# Konstanty
rho_water = 1000  # Hustota vody v kg/m^3
g = 9.81  # Gravitační zrychlení v m/s^2

# Funkce pro výpočet Bernsteinových polynomů
def bernstein_polynom(i, n, t):
    return comb(n, i) * (t ** i) * ((1 - t) ** (n - i))

# Funkce pro výpočet Bézierovy plochy
def bezier_surface(control_points, u_steps=50, v_steps=50):
    n = control_points.shape[0] - 1
    m = control_points.shape[1] - 1
    u_values = np.linspace(0, 1, u_steps)
    v_values = np.linspace(0, 1, v_steps)
    surface = []

    for u in u_values:
        for v in v_values:
            point = np.zeros(3)
            for i in range(n + 1):
                for j in range(m + 1):
                    b_ij = bernstein_polynom(i, n, u) * bernstein_polynom(j, m, v)
                    point += b_ij * control_points[i, j]
            surface.append(point)

    surface = np.array(surface)
    return surface[:, 0].reshape(u_steps, v_steps), surface[:, 1].reshape(u_steps, v_steps), surface[:, 2].reshape(u_steps, v_steps)

# Výpočet těžiště lodi
def calculate_centroid(control_points: ControlPoints):
    total_points = control_points.shape[0] * control_points.shape[1]
    centroid = np.sum(control_points, axis=(0, 1)) / total_points
    print(f"Těžiště lodi: {centroid}")
    return centroid

# Výpočet působiště vztlakové síly s ohledem na těžiště lodi
def calculate_buoyant_force_center(control_points: ControlPoints, water_level: float):
    centroid = calculate_centroid(control_points)
    u_values = np.linspace(0, 1, 50)
    v_values = np.linspace(0, 1, 50)

    total_buoyant_force = 0
    moment_x, moment_y, moment_z = 0, 0, 0

    for u in u_values:
        for v in v_values:
            point = np.zeros(3)
            for i in range(control_points.shape[0]):
                for j in range(control_points.shape[1]):
                    b_ij = bernstein_polynom(i, control_points.shape[0] - 1, u) * bernstein_polynom(j, control_points.shape[1] - 1, v)
                    point += b_ij * control_points[i][j]
           
            # Výška zvlněné hladiny vody na daném bodě (x, y)
            water_level = 2
            z_wave = water_level + 0.05 * np.sin(point[0] * 2 + point[1] * 2)
           
            # Kontrola, zda bod leží pod zvlněnou hladinou vody
            if point[2] < z_wave:
                buoyant_force = rho_water * g * abs(point[2] - z_wave)  # Vztlaková síla závisí na ponořené hloubce
                total_buoyant_force += buoyant_force
               
                # Přispění k momentům vzhledem k těžišti
                dx, dy, dz = point - centroid
                moment_x += buoyant_force * dx
                moment_y += buoyant_force * dy
                moment_z += buoyant_force * dz

    if total_buoyant_force == 0:
        return None

    # Výpočet působiště vztlakové síly vzhledem k těžišti
    center_of_buoyancy = centroid + np.array([
        moment_x / total_buoyant_force,
        moment_y / total_buoyant_force,
        moment_z / total_buoyant_force
    ])
    print(f"Působiště vztlakové síly: {center_of_buoyancy}")
    return center_of_buoyancy


# Kontrolní body pro různé tvary
control_points_start = np.array([
    [[-1, -2, 0], [-1, -1, 0], [-1, 0, 0], [-1, 1, 0], [-1, 2, 0]],
    [[0, -2, 0], [0, -1, 0], [0, 0, 0], [0, 1, 0], [0, 2, 0]],
    [[1, -2, 0], [1, -1, 0], [1, 0, 0], [1, 1, 0], [1, 2, 0]],
    [[2, -2, 0], [2, -1, 0], [2, 0, 0], [2, 1, 0], [2, 2, 0]],
    [[3, -2, 0], [3, -1, 0], [3, 0, 0], [3, 1, 0], [3, 2, 0]],
])

control_points_triangle = np.array([
    [[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0]],
    [[0, 1, 0], [1, 1, 0], [2, 1, 0], [3, 1, 0], [4, 1, 0]],
    [[0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [4, 2, 0]],
    [[0, 3, 0], [1, 3, 0], [2, 3, 0], [3, 3, 0], [4, 3, 0]],
    [[0, 4, 0], [1, 4, 0], [2, 4, 0], [3, 4, 0], [4, 4, 0]],
])

control_points_cylinder = np.array([
    [[1, 0, 0], [0.309, 0.951, 0], [-0.809, 0.588, 0], [-0.809, -0.588, 0], [0.309, -0.951, 0]],
    [[1, 0, 1.25], [0.309, 0.951, 1.25], [-0.809, 0.588, 1.25], [-0.809, -0.588, 1.25], [0.309, -0.951, 1.25]],
    [[1, 0, 2.5], [0.309, 0.951, 2.5], [-0.809, 0.588, 2.5], [-0.809, -0.588, 2.5], [0.309, -0.951, 2.5]],
    [[1, 0, 3.75], [0.309, 0.951, 3.75], [-0.809, 0.588, 3.75], [-0.809, -0.588, 3.75], [0.309, -0.951, 3.75]],
    [[1, 0, 5.0], [0.309, 0.951, 5.0], [-0.809, 0.588, 5.0], [-0.809, -0.588, 5.0], [0.309, -0.951, 5.0]],
])

control_points_hull = np.array([
    [[1, -3, 3], [1, -2, 1], [1, -1, 0], [1, -2, 1], [1, -3, 3]],
    [[0, -1, 3], [0, -1, 1], [1, -1, 0], [2, -1, 1], [2, -1, 3]],
    [[0, 0, 3], [0, 0, 1], [1, 0, 0], [2, 0, 1], [2, 0, 3]],
    [[0, 1, 3], [0, 1, 1], [1, 1, 0], [2, 1, 1], [2, 1, 3]],
    [[1, 3, 3], [1, 2, 1], [1, 1, 0], [1, 2, 1], [1, 3, 3]],
])

control_points_flow = [control_points_start, control_points_triangle, control_points_cylinder, control_points_hull]

# Nastavení proměnných pro animaci
OBJECT_FRAME_WINDOW = 15
CURR_OBJECT_FRAME_WINDOW = OBJECT_FRAME_WINDOW
NUM_OF_SHAPES = len(control_points_flow) - 1
CURR_SHAPE_INDEX = 1

# Funkce pro aktualizaci animace
def update(frame: int):
    global CURR_SHAPE_INDEX, CURR_OBJECT_FRAME_WINDOW
    ax.clear()

    # Nastavení os
    ax.set_xlim(-1, 3)
    ax.set_ylim(-3, 3)
    ax.set_zlim(-3, 3)
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")

    t = (OBJECT_FRAME_WINDOW - CURR_OBJECT_FRAME_WINDOW) / OBJECT_FRAME_WINDOW
    control_points = (1 - t) * np.array(control_points_flow[CURR_SHAPE_INDEX - 1]) + t * np.array(control_points_flow[CURR_SHAPE_INDEX])

    u_vals = np.linspace(0, 1, 30)
    v_vals = np.linspace(0, 1, 30)

    # Výpočet Bézierovy plochy
    X, Y, Z = bezier_surface(control_points, len(u_vals), len(v_vals))

    # Vykreslení kontrolních bodů
    for i in range(len(control_points)):
        for j in range(len(control_points[0])):
            pos = control_points[i][j]
            ax.scatter(pos[0], pos[1], pos[2], color='orange', s=20)

    # Vykreslení plochy
    ax.plot_surface(X, Y, Z, color="green", alpha=0.4, edgecolor="k")

    # Vykreslení hladiny vody
    water_level = 2
    x, y = np.linspace(-1, 3, 30), np.linspace(-3, 3, 30)
    x, y = np.meshgrid(x, y)
    z = water_level + 0.05 * np.sin(x * 2 + y * 2)
    ax.plot_surface(x, y, z, color='blue', alpha=0.3)

    # Na konci animace vykresli šipku a zobraz souřadnice
    if CURR_SHAPE_INDEX == NUM_OF_SHAPES and CURR_OBJECT_FRAME_WINDOW == 0:
        center_of_buoyancy = calculate_buoyant_force_center(control_points, water_level)
        if center_of_buoyancy is not None:
            # Vykreslení šipky
            ax.quiver(
                center_of_buoyancy[0], center_of_buoyancy[1], center_of_buoyancy[2],
                0, 0, 1, color='red', length=1.5, normalize=True, linewidth=2.5
            )

            # Přidání textu s aktuálními souřadnicemi
            coordinates_text = f"Působiště vztlakové síly: ({center_of_buoyancy[0]:.2f}, {center_of_buoyancy[1]:.2f}, {center_of_buoyancy[2]:.2f})"
            fig.text(0.5, 0.95, coordinates_text, ha="center", va="top", fontsize=14, color="red")

    # Přepnutí na další tvar
    if CURR_OBJECT_FRAME_WINDOW == 0:
        if CURR_SHAPE_INDEX == NUM_OF_SHAPES:
            print("Animation complete.")
            anim.event_source.stop()
        else:
            CURR_OBJECT_FRAME_WINDOW = OBJECT_FRAME_WINDOW
            CURR_SHAPE_INDEX += 1
    else:
        CURR_OBJECT_FRAME_WINDOW -= 1

# Nastavení grafu a animace
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# Inicializace animace
anim = FuncAnimation(fig, update, frames=100, repeat=False)

# Uložení animace jako GIF
#output_gif_path = "bezier_animation.gif"
#writer = PillowWriter(fps=20)
#anim.save(output_gif_path, writer=writer)

#print(f"Animace byla uložena jako {output_gif_path}")

# Zobrazení animace
plt.show()