Create a Grid object that describes sphere
#include <stdio.h> #include <math.h> #include <ocean/plankton/VCreator.hpp> #include <ocean/plankton/VInputCreator.hpp> #include <ocean/shrimp/VStringList.hpp> #include <ocean/eagle/rotor3.hpp> #include <ocean/eagle/Matrix.hpp> #include <ocean/shrimp/EaglePhysicalSpacePoint.hpp> #include <ocean/shrimp/EaglePhysicalSpaceRotor.hpp> #include <ocean/shrimp/VEnum.hpp> #include <ocean/shrimp/VObjectStatus.hpp> #include <bone/FishField.hpp> #include <bone/FishSaver.hpp> #include <fish/fiber/bundle/Bundle.hpp> #include <fish/lakeview/bone/FishGrid.hpp> #include <fish/fiber/grid/Grid.hpp> #include <fish/fiber/grid/CartesianChart.hpp> #include <fish/fiber/field/Cell.hpp> #include <fish/fiber/bundle/GridSelector.hpp> namespace { using namespace Wizt; using namespace Eagle::PhysicalSpace; using Eagle::Matrix; using namespace Fiber; #define GRIDNAME "spheregrid" class SpherePoints : public virtual VObject , public StatusIndicator , public FishGridSavable { int n_vertices, n_faces, n_edges, edge_walk; double radius, frac, shiftDist; struct edge { int start,midpoint,end; }; public: TypedSlot<int> SubDivision; TypedSlot<double> Radius; TypedSlot<double> Separation; TypedSlot<point> Center; VOutput<Fiber::Grid> SphereGrid; VOutput<Fiber::Field> ColorField; typedef MemArray<1, point> CoordsArray_t; typedef MemArray<1, TriangleCell> CellArray_t; typedef MemArray<1, EdgeCell> EdgeList_t; typedef MemArray<1, edge> EdgeArray_t; typedef pair<RefPtr<CoordsArray_t>,RefPtr<CellArray_t> > pairArray_t; // Contains the information about the associativity // of edge to vertices and also to the triangles. RefPtr<EdgeList_t> edgeVerticesData; RefPtr<EdgeList_t> edgeFacesData; //RefPtr<Field> Coords; //RefPtr<Field> TriangleCells; point NormalizeRadius(const point&p1) { point result; double length = sqrt(p1[0]*p1[0]+ p1[1]*p1[1] + p1[2]*p1[2]); length = radius/length; result[0] = p1[0]*length; result[1] = p1[1]*length; result[2] = p1[2]*length; return result; } pairArray_t initTetra(const double&rad) { printf("\nInitializing the Tetrahedron...\n"); n_vertices = n_faces = 4; n_edges = 6; frac = 1/10; RefPtr<CoordsArray_t> coordsData = new CoordsArray_t(n_vertices); RefPtr<CellArray_t> cellsData = new CellArray_t(n_faces); MultiArray<1, point>coordinates = *coordsData; MultiArray<1, TriangleCell>triangles = *cellsData; coordinates[0] = NormalizeRadius(point(rad, rad, rad)); coordinates[1] = NormalizeRadius(point(-rad, -rad, rad)); coordinates[2] = NormalizeRadius(point(-rad, rad, -rad)); coordinates[3] = NormalizeRadius(point(rad, -rad, -rad)); triangles[0] = 0,2,1; triangles[1] = 0,1,3; triangles[2] = 2,3,1; triangles[3] = 3,2,0; //TriangleCells->setPersistentData(cells); //Coords->setPersistentData( CoordsData); return pairArray_t(coordsData, cellsData); } pairArray_t initOcta(const double&rad) { printf("\nInitializing the Tetrahedron...\n"); n_vertices = 6; n_faces = 8; n_edges = 12; frac = 1/10; RefPtr<CoordsArray_t> coordsData = new CoordsArray_t(n_vertices); RefPtr<CellArray_t> cellsData = new CellArray_t(n_faces); edgeVerticesData = new EdgeList_t(n_edges); edgeFacesData = new EdgeList_t(n_edges); MultiArray<1, point>coordinates = *coordsData; MultiArray<1, TriangleCell>triangles = *cellsData; MultiArray<1, EdgeCell>edge2Vertices = *edgeVerticesData; MultiArray<1, EdgeCell>edge2Faces = *edgeFacesData; coordinates[0] = NormalizeRadius(point(0.0, 0.0, -rad)); coordinates[1] = NormalizeRadius(point(rad, 0.0, 0.0)); coordinates[2] = NormalizeRadius(point(0.0, -rad, 0.0)); coordinates[3] = NormalizeRadius(point(-rad, 0.0, 0.0)); coordinates[4] =NormalizeRadius(point( 0.0, rad, 0.0)); coordinates[5] =NormalizeRadius(point( 0.0, 0.0, rad)); triangles[0] = 0,1,2; triangles[1] = 0,2,3; triangles[2] = 0,3,4; triangles[3] = 0,4,1; triangles[4] = 5,2,1; triangles[5] = 5,3,2; triangles[6] = 5,4,3; triangles[7] = 5,1,4; // Contains the index of Vertices as in array 'coordinates' edge2Vertices[0] = 0,1; edge2Vertices[1] = 0,2; edge2Vertices[2] = 0,3; edge2Vertices[3] = 0,4; edge2Vertices[4] = 1,4; edge2Vertices[5] = 1,2; edge2Vertices[6] = 2,3; edge2Vertices[7] = 3,4; edge2Vertices[8] = 1,5; edge2Vertices[9] = 2,5; edge2Vertices[10] = 3,5; edge2Vertices[11] = 4,5; // Contains the index of Trinagle as in arry 'triangles' edge2Faces[0] = 0,3; edge2Faces[1] = 0,1; edge2Faces[2] = 1,2; edge2Faces[3] = 2,3; edge2Faces[4] = 3,7; edge2Faces[5] = 0,4; edge2Faces[6] = 1,5; edge2Faces[7] = 2,6; edge2Faces[8] = 4,7; edge2Faces[9] = 4,5; edge2Faces[10] = 5,6; edge2Faces[11] = 6,7; //TriangleCells->setPersistentData(cells); //Coords->setPersistentData( CoordsData); return pairArray_t(coordsData, cellsData); } pairArray_t initIcosa(const double&rad) { printf("\nInitializing the Icosahedron...\n"); n_vertices = 12; n_faces = 20; n_edges = 30; frac = 1/10; double t, tau, one; t=rad;//(1/sqrt(5))/2; tau = (t/sqrt(1+t*t)); one = (1/sqrt(1+t*t)); RefPtr<CoordsArray_t> coordsData = new CoordsArray_t(n_vertices); RefPtr<CellArray_t> cellsData = new CellArray_t(n_faces); MultiArray<1, point>coordinates = *coordsData; MultiArray<1, TriangleCell>triangles = *cellsData; coordinates[0] =NormalizeRadius(point( tau, one, 0.0)); coordinates[1] =NormalizeRadius(point( -tau, one, 0.0)); coordinates[2] =NormalizeRadius(point( -tau, -one, 0.0)); coordinates[3] =NormalizeRadius(point( tau, -one, 0.0)); coordinates[4] =NormalizeRadius(point( one, 0.0, tau)); coordinates[5] =NormalizeRadius(point( one, 0.0, -tau)); coordinates[6] =NormalizeRadius(point( -one, 0.0,-tau)); coordinates[7] =NormalizeRadius(point( -one, 0.0, tau)); coordinates[8] =NormalizeRadius(point( 0.0, tau, one)); coordinates[9] =NormalizeRadius(point( 0.0, -tau, one)); coordinates[10] =NormalizeRadius(point( 0.0, -tau, -one)); coordinates[11] =NormalizeRadius(point( 0.0, tau, -one)); triangles[0] = 4,8,7; triangles[1] = 4,7,9; triangles[2] = 5,6,11; triangles[3] = 5,10,6; triangles[4] = 0,4,3; triangles[5] = 0,3,5; triangles[6] = 2,7,1; triangles[7] = 2,1,6; triangles[8] = 8,0,11; triangles[9] = 8,11,1; triangles[10] = 9,10,3; triangles[11] = 9,2,10; triangles[12] = 8,4,0; triangles[13] = 11,0,5; triangles[14] = 4,9,3; triangles[15] = 5,3,10; triangles[16] = 7,8,1; triangles[17] = 6,1,11; triangles[18] = 7,2,9; triangles[19] = 6,10,2; //TriangleCells->setPersistentData(cells); //Coords->setPersistentData( CoordsData); return pairArray_t(coordsData, cellsData); } pairArray_t initSphere(const double&rad) { printf("\nInitializing the sphere...\n"); n_vertices = 6; n_faces = 8; n_edges = 12; frac = 1/10; RefPtr<CoordsArray_t> coordsData = new CoordsArray_t(n_vertices); RefPtr<CellArray_t> cellsData = new CellArray_t(n_faces); MultiArray<1, point>coordinates = *coordsData; MultiArray<1, TriangleCell>triangles = *cellsData; coordinates[0] = rad, 0.0, 0.0; //a coordinates[1] = 0.0, 0.0, -rad; //b coordinates[2] = -rad, 0.0, 0.0; //c coordinates[3] = 0.0, 0.0, rad; //d coordinates[4] = 0.0, rad, 0.0; //e coordinates[5] = 0.0, - rad, 0.0; //f triangles[0] = 3,0,4; triangles[1] = 0,1,4; triangles[2] = 1,2,4; triangles[3] = 2,3,4; triangles[4] = 0,3,5; triangles[5] = 1,0,5; triangles[6] = 2,1,5; triangles[7] = 3,2,5; return pairArray_t(coordsData, cellsData); } int searchMidpoint(const int&indexStart, const int&indexEnd, EdgeArray_t&edgeArray, MultiArray<1,point>&coords) { for(int i = 0; i < edge_walk; i++) { if((edgeArray[i].start == indexStart && edgeArray[i].end == indexEnd) || (edgeArray[i].start == indexEnd && edgeArray[i].end == indexStart)) { int result = edgeArray[i].midpoint; edgeArray[i].start = edgeArray[edge_walk -1].start; edgeArray[i].end = edgeArray[edge_walk-1].end; edgeArray[i].midpoint = edgeArray[edge_walk-1].midpoint; edge_walk--; return result; } } //vertex not in the list so we add it edgeArray[edge_walk].start = indexStart; edgeArray[edge_walk].end = indexEnd; edgeArray[edge_walk].midpoint = n_vertices; // create a new vertex coords[n_vertices][0] = (coords[indexStart][0]+coords[indexEnd][0])/2.0; coords[n_vertices][1] = (coords[indexStart][1]+coords[indexEnd][1])/2.0; coords[n_vertices][2] = (coords[indexStart][2]+coords[indexEnd][2])/2.0; double length = sqrt(coords[n_vertices][0]*coords[n_vertices][0] + coords[n_vertices][1]*coords[n_vertices][1] + coords[n_vertices][2]*coords[n_vertices][2]); length = radius/length; //length = 1/length; coords[n_vertices][0] *= length; coords[n_vertices][1] *= length; coords[n_vertices][2] *= length; n_vertices++; edge_walk++; return edgeArray[edge_walk -1].midpoint; } int createEdgeArray(int&edgeCount, const TriangleCell&face, MultiArray<1, EdgeCell>&edges, bool *vertex_pairs) { int a = face[0]; int b = face[1]; int c = face[2]; int start, end; //Edge 1 if(a<b) { start = a; end = b; } else { start = b; end = a; } if(vertex_pairs[start*n_edges + end] != 1) { edges[edgeCount++] = start, end; vertex_pairs[start*n_edges + end] = 1; } //Edge 2 if(b<c) { start = b; end = c; } else { start = c; end = b; } if(vertex_pairs[start*n_edges + end] != 1) { edges[edgeCount++] = start, end; vertex_pairs[start*n_edges + end] = 1; } //Edge 3 if(a<c) { start = a; end = c; } else { start = c; end = a; } if(vertex_pairs[start*n_edges+ end] != 1) { edges[edgeCount++] = start, end; vertex_pairs[start*n_edges + end] = 1; } //printf("\nDelete this...EDGECOUNT %d\n", edgeCount); return edgeCount; } pairArray_t subDivideTriangles(const pairArray_t&input) //const MultiArray<1,point>&crds, const MultiArray<1,TriangleCell>&triangle) { int n_vertices_new; int faces, n_faces_new = 4*n_faces; int edgeCount = 0; edge_walk = faces = 0; n_vertices_new = n_vertices + n_edges; n_edges = 4*n_edges; //2*n_vertices + 3*n_faces; //4*n_edges; const MultiArray<1,point>&crds = *input.first; const MultiArray<1,TriangleCell>&triangle = *input.second; EdgeArray_t edgeArray(n_edges); //EdgeArray_t end = new EdgeArray_t(n_edges); //EdgeArray_t midPoint = new EdgeArray_t(n_edges); // Infomation about coordinates and triangle faces. RefPtr<CoordsArray_t> coordsData = new CoordsArray_t(n_vertices_new*2); RefPtr<CellArray_t> cellsData = new CellArray_t(n_faces_new*2); MultiArray<1, point> newCrds = *coordsData; MultiArray<1, TriangleCell> newTriangle = *cellsData; // Information about edges-vertices and egde-triangles. edgeVerticesData = new EdgeList_t(n_edges); edgeFacesData = new EdgeList_t(n_edges); MultiArray<1, EdgeCell>edge2Vertices = *edgeVerticesData; MultiArray<1, EdgeCell>edge2Faces = *edgeFacesData; bool vertex_pairs[n_edges][n_edges]; /*for(int s=0; s < n_edges; s++) for(int e=0; e < n_edges; e++) { vertex_pairs[s][e] = 0; } */ // Copy existing vertices into new arrray. //printf("\n****Vertices %d\n", n_vertices); for(int j = 0; j<n_vertices; j++) { newCrds[j] = crds[j]; } for(int i = 0; i < n_faces; i++) { int a = triangle[i][0]; int b = triangle[i][1]; int c = triangle[i][2]; int ab = searchMidpoint(b,a,edgeArray,newCrds); int bc = searchMidpoint(c,b,edgeArray,newCrds); int ca = searchMidpoint(a,c,edgeArray,newCrds); newTriangle[faces][0] = a; newTriangle[faces][1] = ab; newTriangle[faces][2] = ca; createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, &vertex_pairs[0][0]); faces++; newTriangle[faces][0] = ca; newTriangle[faces][1] = ab; newTriangle[faces][2] = bc; createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, &vertex_pairs[0][0]); faces++; newTriangle[faces][0] = ca; newTriangle[faces][1] = bc; newTriangle[faces][2] = c; createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, &vertex_pairs[0][0]); faces++; newTriangle[faces][0] = ab; newTriangle[faces][1] = b; newTriangle[faces][2] = bc; createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, &vertex_pairs[0][0]); faces++; } n_faces = faces; RefPtr<CoordsArray_t> newCoordsData = new CoordsArray_t(n_vertices*2); MultiArray<1, point> newCoords = *newCoordsData; for(int j = 0; j<n_vertices; j++) { newCoords[j] = newCrds[j]; newCoords[j+n_vertices] = newCrds[j]; } //Copying connectivity information for(int i = 0; i < n_faces; i++) { newTriangle[n_faces+i][0] = newTriangle[i][0] + n_vertices; newTriangle[n_faces+i][1] = newTriangle[i][1] + n_vertices; newTriangle[n_faces+i][2] = newTriangle[i][2] + n_vertices; } printf("\nn_vertice = %d .... n_vertice_new = %d\n", n_vertices, n_vertices_new); return pairArray_t(newCoordsData, cellsData); } void CreateSphere(Grid&G, unsigned subDiv, double sphereRadius, const point&Center) { Skeleton&Vertices = G.makeVertices(3); Skeleton&Cells = G[ SkeletonID(2,1) ]; Skeleton&Edges = G[ SkeletonID(1,1) ]; RefPtr<Chart> myCartesian = G.makeChart(typeid(Fiber::CartesianChart3D)); Representation&CartesianVertices = Vertices[myCartesian]; Representation&CellsVertices = Cells[Vertices]; Representation&EdgesVertices = Edges[Vertices]; Representation&EdgesFaces = Edges[Cells]; RefPtr<Field>Coords = CartesianVertices[FIBER_POSITIONS]; RefPtr<Field>TriangleCells = CellsVertices[FIBER_POSITIONS]; RefPtr<Field>SphereColor = CartesianVertices["COLORS"]; RefPtr<Field>EdgesPerVertices = EdgesVertices["EDGEVERTEX"]; RefPtr<Field>EdgesPerFaces = EdgesFaces["EDGEFACE"]; //pairArray_t surface = initTetra(sphereRadius); //( CoordsData, cells); //pairArray_t surface = initIcosa(sphereRadius); pairArray_t surface = initOcta(sphereRadius); for(index_t i = 0; i < subDiv; i++) { printf("\nWe are in %d level of sub division\n", subDiv); surface = subDivideTriangles(surface); } MultiArray<1,point>&crds = *surface.first; Ref<MemArray<1, double> > Colors = new MemArray<1, double>(crds.nElements());//n_vertices*2); MultiArray<1, double> vertexColor = *Colors; //MultiArray<1,point>&crds = *surface.first; int size = crds.nElements()/2; tvector shift = Center - point(0,0,0); for(int j = 0; j<size; j++) { crds[j] += shift; vertexColor[j] = 0.0; crds[j+size] += shift; crds[j+size][0] -= shiftDist;//2*1.3;//0.115; vertexColor[j+size] = 1.0; //printf("Point %d = %f %f %f\n",j, crds[j][0],crds[j][1],crds[j][2]); } Coords->setPersistentData(surface.first); TriangleCells->setPersistentData(surface.second); SphereColor->setPersistentData(Colors); EdgesPerVertices->setPersistentData(edgeVerticesData); EdgesPerFaces->setPersistentData(edgeFacesData); } SpherePoints(const string&name, int p, const RefPtr<VCreationPreferences>&VP) : VObject(name, p, VP) , StatusIndicator(this, "Testing the Error.") , FishGridSavable(this, SphereGrid) , SubDivision(this, "subdivision", 3) , Radius(this, "radius", 1/sqrt(3)) //(1+sqrt(5))/2) , Separation(this, "shiftdistance", 0.0) , Center(this, "center", point(0,0,0) ) , SphereGrid( self(), "sphereGrid") , ColorField(self(), "colors") { BundlePtr B; SubDivision.setProperty("max", 5); Radius.setProperty("max", 1.8); Separation.setProperty("max", 2.0); CreateSphere(B[0.0][GRIDNAME], 0, 1/sqrt(3.0), point(0,0,0) ); GridSelector GS(GRIDNAME, B); SphereGrid() = GS; //printf("\nSphere Points Object.....%s\n", SphereGrid().Gridname().c_str()); } override bool update( VRequest&R, double precision ); }; // Update function of the Sphere class. bool SpherePoints::update(VRequest&R, double precision) { GridSelector GS; int subDivision; double rad, sphereShift; SphereGrid << R >> GS; SubDivision << R >> subDivision; Radius << R >> rad; Separation << R >> sphereShift; point C; Center << R >> C; //Grid&G = GS.theSourceBundle[0.0][ GRIDNAME ]; RefPtr<Grid> G = GS.theSourceBundle->newGrid(); radius = rad; shiftDist = sphereShift; setStatusInfo(R, "Creating sphere now!!!"); if(rad <= 0.0) return true; CreateSphere(*G, subDivision, rad, C); double T = 0.0; GS.theSourceBundle[T].insert( GRIDNAME, G); GS.theGridname = GRIDNAME; Info<Grid> gridInfo(0.0, (*GS.theSourceBundle)(T), G); FieldSelector FS(GS, gridInfo); FS.selectField("COLORS"); SphereGrid << R << GS; ColorField << R << FS; puts("bool SpherePoints::update(VRequest&R, double precision)"); return true; } static Ref<VCreator<SpherePoints> > mySphere("CreateGrid/Sphere", ObjectQuality::BETA); static Ref<VInputCreator<Grid, VCreator<SpherePoints, AcceptList<point> > > > mySphereIO("Compute/Sphere", ObjectQuality::BETA); } // namespace