Sphere.cpp

Create a Grid object that describes sphere

Author:

"Bidur Bohara" <bbohar1@lsu.edu>

#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

---