A VObject which provides a set of lines.
#include <ocean/plankton/VCreator.hpp> #include <ocean/shrimp/TimeDependent.hpp> #include <bundle/Bundle.hpp> #include <bone/BundleInput.hpp> #include <grid/CartesianChart.hpp> #include <ocean/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> TimeSteps; TypedSlot<double> TMax; /* Create helical lines on the given Grid object */ void makeGrid(Grid&G, double time) { index_t NumberOfLines = 4; 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.1/(time*time+.2); double x = A*si, y = B*co, z = C*(i+j); 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) ] = (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( self(), "grid") , TimeSteps( this, "timesteps", 30,2) , TMax(this, "tmax", 2.0, 1) { TMax.setProperty("max", 100.0); } override bool update(VRequest&Context, double precision) { puts("override 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" 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); } GridSelector GS(gridname, B); LineGrid << Context << GS; return true; } }; static Ref<VCreator<HelicalLines> > myCreator("CreateGrid/HelicalLines", ObjectQuality::DEMO); } // namespace
1.6.2