3D : Total Curvature

#!/usr/bin/env python3

import os
import vedo

import numpy as np
import pygranite as pg  # trajectory library
import pgt  # pygranite toolkit

from tabulate import tabulate


def heightmap(csvf):
    return np.genfromtxt(csvf, delimiter=",")


if __name__ == "__main__":
    root = os.path.abspath(os.path.dirname(__file__))

    # simulation properties
    settings = pg.IntegratorSettings()
    settings.Space = pg.Space3D
    settings.MinimumAliveParticles = 16
    settings.WindfieldMode = pg.WindfieldMode.Dynamic
    settings.Integrator = pg.Integrator.ClassicRungeKutta
    settings.DeltaT = 1
    settings.DataTimeDistance = 60
    settings.MaximumSimulationTime = 1000
    settings.GridScale = [25, 25, 25]
    settings.SaveInterval = 0
    settings.CurvatureMode = pg.CurvatureMode.TotalCurvature

    # load heightmap
    hm = heightmap(root + "pyr_1000.csv")
    settings.Topography = hm

    # generate start positions
    spawn_min = np.array([2000, 2500, 150])
    spawn_max = np.array([2300, 2600, 170])

    inp = np.random.rand(1024, 3)
    for i in range(inp.shape[0]):
        inp[i] = ((spawn_max - spawn_min) * inp[i] + spawn_min)

    start_set = pg.TrajectorySet(inp)

    with pgt.Timer("Angle Curvature"):
        # initialize windfield loader
        directory = root + "fields_scaled/"
        loader = pgt.netCDFLoader3D([
            directory + f"{i}.nc" for i in range(len(os.listdir(directory)))
        ])
        loader.PrintProgress = False
        settings.CurvatureMode = pg.CurvatureMode.FastTotalCurvature

        # computation
        integrator = pg.TrajectoryIntegrator(
            settings, loader, start_set)
        result_set: pg.TrajectorySet = integrator.compute()

        # visualization using vedo
        angle_curvatures = result_set.totalCurvatures()

    with pgt.Timer("Cross Curvature"):
        # initialize windfield loader
        directory = root + "fields_scaled/"
        loader = pgt.netCDFLoader3D([
            directory + f"{i}.nc" for i in range(len(os.listdir(directory)))
        ])
        loader.PrintProgress = False
        settings.CurvatureMode = pg.CurvatureMode.TotalCurvature

        # computation
        integrator = pg.TrajectoryIntegrator(
            settings, loader, start_set)
        result_set: pg.TrajectorySet = integrator.compute()

        # visualization using vedo
        cross_curvatures = result_set.totalCurvatures()

    print("\n### Computed Values:\n")
    print(tabulate(np.transpose(np.array([angle_curvatures, cross_curvatures, abs(
        cross_curvatures - angle_curvatures), angle_curvatures / cross_curvatures])), headers=["ANGLE", "CROSS", "ABS-DIFF", "RATIO"], tablefmt='orgtbl'))