Sphere.cpp

Create a Grid object that describes sphere.

Create a Grid object that describes sphere

Author

"Bidur Bohara" bboha.nosp@m.r1@l.nosp@m.su.ed.nosp@m.u

 
#include <stdio.h>
#include <math.h>
 
#include <ocean/plankton/VInputCreator.hpp>
 
#include <ocean/shrimp/VStringList.hpp>
 
//#include <eagle/rotor3.hpp>
#include <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 <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
{
        struct edge
        {
                unsigned start,midpoint,end;
        };
}
 
namespace Fiber
{
 
template <>
string Fiber::element_to_string(const edge&)
{
        return {};
}
 
} // Fiber
 
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;
 
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);
  }
 
  unsigned searchMidpoint(const unsigned&indexStart, const unsigned&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,
                      std::vector<bool>&vertex_pairs)
  {
 
    unsigned a = face[0];
    unsigned b = face[1];
    unsigned c = face[2];
    unsigned 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;
 
    Ref<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;
 
    std::vector<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++)
      {
    unsigned a = triangle[i][0];
    unsigned b = triangle[i][1];
    unsigned c = triangle[i][2];
 
    unsigned ab = searchMidpoint(b,a,*edgeArray,newCrds);
    unsigned bc = searchMidpoint(c,b,*edgeArray,newCrds);
    unsigned 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);
        faces++;
 
        newTriangle[faces][0] = ca;
        newTriangle[faces][1] = ab;
        newTriangle[faces][2] = bc;
        createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, vertex_pairs);
        faces++;
 
        newTriangle[faces][0] = ca;
        newTriangle[faces][1] = bc;
        newTriangle[faces][2] = c;
        createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, vertex_pairs);
        faces++;
 
        newTriangle[faces][0] = ab;
        newTriangle[faces][1] = b;
        newTriangle[faces][2] = bc;
        createEdgeArray(edgeCount, newTriangle[faces], edge2Vertices, vertex_pairs );
        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(this, "sphereGrid")
    , ColorField(this, "colors")
  {
    BundlePtr B;
 
    SubDivision.setProperty("max", 3);
    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());
  }
 
 
  bool update( VRequest&R, double precision ) override;
};
 
 
 
// 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;
}
 
} // namespace
 
namespace
{
using namespace Panthalassa;
 
static VFilter< point, SpherePoints, Grid >
MyCreator(      Category( "CreateGrid" ) + VIdentifier( "Sphere" ) + Application( "General;HydroVish" ) +
                Description( "Create a sphere as triangular surface." ) +
                Help( "todo" ) +
                Url( "http://ahm.co.at/wiki/en" ),
                ObjectQuality::EXPERIMENTAL );
}