HelicalLines.cpp

A VObject which provides a set of lines.

A VObject which provides a set of lines.

Author

werner

See also

Fiber::LineSet

#include <ocean/shrimp/TimeDependent.hpp>
 
#include <bundle/Bundle.hpp>
#include <bone/BundleInput.hpp>
 
#include <grid/CartesianChart.hpp>
#include <eagle/PhysicalSpace.hpp>
#include <bone/GridObject.hpp>
#include <grid/types/LineSet.hpp>
#include <bone/FishSaver.hpp>
#include <field/ArrayRef.hpp>
 
#include <eagle/rotor3.hpp>
 
using namespace Wizt;
using namespace Fiber;
using namespace Eagle;
using namespace Eagle::PhysicalSpace;
 
namespace
{
 
class   HelicalLines : public virtual VObject, public FishGridSavable
{
public:
        VOutput<Grid>           LineGrid;
 
        TypedSlot<int>          tsNumberOfLines;
        TypedSlot<int>          TimeSteps;
        TypedSlot<double>       TMax;
 
 
        /*
          Create helical lines on the given Grid object
         */
        void    makeGrid(Grid&G, double time, index_t NumberOfLines)
        {
        index_t VerticesPerLine = 512; 
 
        index_t TotalVertices = VerticesPerLine*NumberOfLines;
 
        //
        // Define the vertices of the lines
        //
        Ref<LineSet::CoordsMemArray_t> CoordsArray( TotalVertices );
        std::vector<point>&Crds = CoordsArray->myChunk()->get_vector(); 
 
        ArrayRef<tvector3>      Velocity( TotalVertices ); 
        ArrayRef<tvector3>      Normals( TotalVertices );
        ArrayRef<tvector3>      dT_ds( TotalVertices );
        ArrayRef<double>        CurvatureCross( TotalVertices );
        ArrayRef<double>        TorsionCross( TotalVertices );
        ArrayRef<double>        CurvatureGA( TotalVertices );
        ArrayRef<double>        TorsionGA( TotalVertices );
 
        ArrayRef<double>        VertexEnumeration( TotalVertices );
 
        /*
          A data array that contains 128 values per vertex of the lines.
         */
        int     nVariant = 128;
        ArrayRef<unsigned char,2> VariantData( MIndex( TotalVertices, nVariant ) );
 
                {for(index_t j=0; j<NumberOfLines; j++)
                 for(index_t i=0; i<VerticesPerLine; i++)
                 {
                 index_t n( i + j*VerticesPerLine ); 
 
                         VertexEnumeration[n] = double(n)/ TotalVertices;
 
                 double omega = 1.0/40.0*(j+0.5); 
 
                 double si = sin(i*omega), 
                        co = cos(i*omega);
 
                 double A = (.1+time)*18, 
                        B = (.1+time)*16,
                        C = 0.04/(time*time+.2);
 
                 double Scale = 1.0 / 20.0;
 
                 double x = A*si        *Scale,
                        y = B*co        *Scale,
                        z = C*(i+j)     *Scale - 0.5; 
 
                         Crds[ n ] =  {x, y, z};
 
                        // velocity is derivative along curve parameter, 
                        // which is "i" in this case: 
 
                 double Vx = A*co*omega,
                        Vy = -B*si*omega,
                        Vz = C; 
 
                         Velocity[ n ] = {Vx, Vy, Vz}; 
 
 
                        // second derivative, the curve normal
                double Nx = -A*si*omega*omega,
                       Ny = -B*co*omega*omega,
                       Nz = 0; 
 
                        Normals[ n ] = {Nx, Ny, Nz}; 
 
                //||dot c <x> ddot c|| / ||dot c||^3 
                tvector Cross( Vy*Nz - Ny*Vz, -(Vx*Nz-Nx*Vz), Vx*Ny-Nx*Vy ); 
                double s = norm(Velocity[n]); 
                        CurvatureCross[ n ] = norm(Cross)/(s*s*s); 
 
                // derivative of normalized velocity: 
 
                double A2 = A*A,
                       B2 = B*B,
                   omega2 = omega*omega,
                   omega3 = omega2*omega; 
 
                double D = A2*co*co*omega2 + B2*si*si*omega2 + C*C; 
 
                double E = -2*A2*co*omega3*si + 2*B2*si*omega3*co; 
 
                double dTUx = -A*si*omega2 / sqrt(D) - A*co*omega*E / (2*pow(D, 3./2)),
                       dTUy = -B*co*omega2 / sqrt(D) + B*si*omega*E / (2*pow(D, 3./2)),
                       dTUz = -C*E / (2*pow(D, 3./2)); 
 
                tvector dTU(dTUx, dTUy, dTUz ); 
 
                        dT_ds[ n ] = dTU; 
 
                        CurvatureGA[ n ] = norm( dTU );
 
                // second derivative of normalized velocity 
                double omega4 = omega3*omega; 
 
                double F = 2*A2*si*si*omega4-2*A2*co*co*omega4+2*B2*co*co*omega4-2*B2*si*si*omega4; 
 
                double ddTUx = -A*co*omega3/sqrt(D) + A*si*omega2*E / pow(D, 3./2) + 3./4 * A*co*omega*E*E / pow(D, 5./2) - 1./2 * A*co*omega*F / pow(D, 3./2),
                       ddTUy = -B*si*omega3/sqrt(D) + B*co*omega2*E / pow(D, 3./2) - 3./4 * B*si*omega*E*E / pow(D, 5./2) + 1./2 * B*si*omega*F / pow(D, 3./2),
                       ddTUz = 3./4 * C *E*E / pow(D, 5./2) - 1./2 * C * F / pow(D, 3./2); 
 
                tvector ddTU(ddTUx, ddTUy, ddTUz ); 
                        // third derivative 
                tvector U; 
                        U ={ (.1+time)*-18*cos(i*omega)*omega*omega*omega,
                             (.1+time)*16*sin(i*omega)*omega*omega*omega,
                             0}; 
 
                        TorsionCross[ n ] = dot(Cross, U) / norm2(Cross); 
 
 
                //double tau = ( ( U ^ Normals[n] ) * Velocity[n] ).trivec().r(); 
 
//      double  tau = ( ( U ^ Normals[n] ) ^ Velocity[n] ).r() / norm2( Normals[n] ^ Velocity[n] ); 
 
                double  tau = ( ( ddTU ^ dTU ) ^ Velocity[n] ).r() / norm2( dTU ^ Velocity[n] ); 
 
                        TorsionGA[ n ] = tau;
 
                        /*
                          setup a funny variant field
                         */
                        for(int vidx = 0; vidx < nVariant; vidx++)
                        {
                        double  V = 1.0*vidx / nVariant,
                                N = 1.0 * i/ VerticesPerLine; 
                        //      J = 1.0*j/NumberOfLines;
                        //      V *= (1.0-V); 
 
                                V = 0.5+0.5*(sin(V*3.141592*(1+4*N) ));
 
//                      V = N;
 
                                VariantData[ MIndex( n, vidx) ] = static_cast<unsigned char>(255.0*V);
                        }
         }}
 
        // 
        // Define the indices that form a line 
        // (this could actually be done in the first loop as well)
        // This is actually an array of arrays, since each line could be 
        // of different length. Just in this case here each of them consists 
        // of 1000 vertices each.
        // 
        //
        Ref<LineSet::LinesetArray_t>     LinesetArray(NumberOfLines);
        MultiArray<1, LineSet::LineIndices_t>&LineIndices = *LinesetArray; 
                {for(index_t j=0; j<NumberOfLines; j++)
                {
                LineSet::LineIndices_t&E = LineIndices[ j ];
                        E.resize(VerticesPerLine);
                        for(index_t i=0; i<VerticesPerLine; i++)
                        {
                                E[i] = i + j*VerticesPerLine;
                        }
                }}
 
        LineSet theLines(G.self(), CoordsArray, LinesetArray ); 
                (*theLines.CartesianVertices)[ LineSet::PredefinedFieldNames::Velocity ]->setPersistentData( Velocity ); 
                (*theLines.CartesianVertices)[ "Analytic CrossCurvature" ]->setPersistentData( CurvatureCross );
                (*theLines.CartesianVertices)[ "Analytic GACurvature" ]->setPersistentData( CurvatureGA ); 
                (*theLines.CartesianVertices)[ "Analytic Normals" ]->setPersistentData( Normals ); 
                (*theLines.CartesianVertices)[ "Analytic CrossTorsion" ]->setPersistentData( TorsionCross ); 
                (*theLines.CartesianVertices)[ "Analytic GATorsion" ]->setPersistentData( TorsionGA ); 
                (*theLines.CartesianVertices)[ "Analytic dT_ds" ]->setPersistentData( dT_ds ); 
 
                (*theLines.CartesianVertices)[ "Variant Data" ]->setPersistentData( VariantData ); 
                (*theLines.CartesianVertices)[ "LocalVertexNumber" ]->setPersistentData( VertexEnumeration );
 
                theLines.setupStandardFields(G.self() );
        }
 
 
        HelicalLines(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VObject(name, p, VP)
        , FishGridSavable(this, LineGrid)
        , LineGrid( this, "grid")
        , tsNumberOfLines(this, "NumberOfLines", 4, 0)
        , TimeSteps( this, "timesteps", 30,2)
        , TMax(this, "tmax", 2.0, 1)
        {
                TMax.setProperty("max", 100.0); 
                tsNumberOfLines.setProperty("max", 20); 
        }
 
        bool update(VRequest&Context, double precision) override
        {
//      puts("bool update(VRequest&Context, double precision)");
        BundlePtr B; // create new bundle each time on update
        double  t_min = 0.0,
                t_max = 1.0;
 
        int     nTimeSteps = 10;
                TimeSteps << Context >> nTimeSteps; 
                TMax      << Context >> t_max;
 
#define gridname        "helicalgrid" 
 
        int     NumberOfLines = 4;
                tsNumberOfLines << Context >> NumberOfLines;
 
                for(int t = 0; t<nTimeSteps; t++)
                {
                double  time = t*(t_max - t_min) / nTimeSteps + t_min;
                Grid&G = B[ time ][ gridname ]; 
                        makeGrid( G, time, NumberOfLines);
                }
 
        GridSelector GS(gridname, B);
                LineGrid << Context << GS;
 
                return true;
        }
};
 
} // namespace
 
namespace
{
using namespace Panthalassa;
 
static VSource< HelicalLines, Grid >
MyCreator(      Category( "CreateGrid" ) + VIdentifier( "HelicalLines" ) + Application( "General;HydroVish" ) +
                Description( "Create a time series of lines." ) +
                Help( "todo" ) +
                Url( "http://ahm.co.at/wiki/en" ),
                ObjectQuality::EXPERIMENTAL );
}