//REFERENCE: https://github.com/FreeFem/FreeFem-sources/blob/master/plugin/mpi/PETSc-code.hpp
#include <petsc.h>
static char help[] = "Creates the element stiffness matrix.\n"
"Option prefix = opt_.\n";

int main(int argc, char **args){
  //DEFINITIONS OF PETSC's DMPLEX LINGO:
  //POINT: A topology element (cell, face, edge, or vertex).
  //CHART: It an interval from 0 to the number of "points." (the range of admissible linear indices)
  //STRATUM: A subset of the "chart" which corresponds to all "points" at a given "level."
  //LEVEL: This is either a "depth" or a "height".
  //HEIGHT: Dimensionality of an element measured from 0D to 3D. Heights: cell = 0, face = 1, edge = 2, vertex = 3.
  //DEPTH: Dimensionality of an element measured from 3D to 0D. Depths: cell = 3, face = 2, edge = 1, vertex = 0;
  //CLOSURE: *of an element is the collection of all other elements that define it.I.e., the closure of a surface is the collection of vertices and edges that make it up.
  //STANDARD LABELS: These are default tags that DMPlex has for its topology. ("depth")
  PetscErrorCode ierr;//Error tracking variable.
  DM dm;
  PetscSection s;
  //DMLabel Dirichlet1,
  //        Dirichlet2;
  Mat K; //Global stiffness matrix (Square, assume unsymmetric).
      //KE;//Element stiffnes matrix (Square, assume unsymmetric).
  Vec XYZ, //Coordinate vector, contains spatial locations of mesh's vertices (NOTE: This vector self-destroys!).
      U; //Displacement vector.
  PetscBool interpolate = PETSC_TRUE,//Instructs Gmsh importer whether to generate faces and edges (Needed when using P2 or higher elements).
            useCone = PETSC_TRUE;//Needed to extract closure of cells.
  PetscInt nc,//number of cells. (PETSc lingo for "elements")
           nv,//number of vertices. (PETSc lingo for "nodes")
           nf,//number of faces. (PETSc lingo for "surfaces")
           ne,//number of edges. (PETSc lingo for "lines")
           pStart,//starting LI of global elements.
           pEnd,//ending LI of all elements.
           cStart,//starting LI for cells global arrangement.
           cEnd,//ending LI for cells in global arrangement.
           vStart,//starting LI for vertices in global arrangement.
           vEnd,//ending LI for vertices in global arrangement.
           fStart,//starting LI for faces in global arrangement.
           fEnd,//ending LI for faces in global arrangement.
           eStart,//starting LI for edges in global arrangement.
           eEnd,//ending LI for edges in global arrangement.
           sizeK,//Size of the element stiffness matrix.
           ii,jj,//Dedicated looping variables.
           *closure,//Pointer to the closure elements of a cell.
           size_closure,//Size of the closure of a cell.
           dim;//Dimension of the mesh.
  const PetscScalar *xyz_el;//Pointer to array of global XYZ nodal data.
  PetscReal *x_el,//Pointer to a vector that will store the x-coordinates of an element's vertices.
            *y_el,//Pointer to a vector that will store the y-coordinates of an element's vertices.
            *z_el;//Pointer to a vector that will store the z-coordinates of an element's vertices.
  const char* gmshfile = "TopOptmesh2.msh";//Name of gmsh file to import.

  ierr = PetscInitialize(&argc,&args,NULL,help); if(ierr) return ierr; //And the machine shall work....

  //MESH IMPORT=================================================================
  ierr = DMCreate(PETSC_COMM_WORLD,&dm); CHKERRQ(ierr);//Create distributed memory object.
  //ierr = DMSetType(dm,DMPLEX); CHKERRQ(ierr);//Create distributed memory object.

  //IMPORTANT NOTE: Gmsh only creates "cells" and "vertices," it does not create the "faces" or "edges."
  //Gmsh probably can generate them, must figure out how to.
  ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,gmshfile,interpolate,&dm); CHKERRQ(ierr);//Read Gmsh file and generate the DMPlex.
  ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);//1-D, 2-D, or 3-D
  ierr = DMPlexGetChart(dm, &pStart, &pEnd); CHKERRQ(ierr);//Extracts linear indices of cells, vertices, faces, and edges.
  ierr = DMGetCoordinatesLocal(dm,&XYZ); CHKERRQ(ierr);//Extracts coordinates from mesh.(Contiguous storage: [x0,y0,z0,x1,y1,z1,...])
  //PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
  //ierr = VecGetArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's data.
  //PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
  ierr = VecGetArrayRead(XYZ,&xyz_el);
  //PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
  //ierr = VecRestoreArray(XYZ,&xyz_el); CHKERRQ(ierr);
  //ierr = VecView(XYZ,PETSC_VIEWER_STDOUT_WORLD);


  //IMPORTANT NOTE: PETSc assumes that vertex-cell meshes are 2D even if they were 3D, so its ordering changes.
  //Cells remain at height 0, but vertices move to height 1 from height 3. To prevent this from becoming an issue
  //the "interpolate" variable is set to PETSC_TRUE to tell the mesh importer to generate faces and edges.
  //PETSc, therefore, technically does additional meshing. Gotta figure out how to get this from Gmsh directly.
  //PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
  ierr = DMPlexGetDepthStratum(dm,3, &cStart, &cEnd);//Get LI of cells.
  ierr = DMPlexGetDepthStratum(dm,2, &fStart, &fEnd);//Get LI of faces
  ierr = DMPlexGetDepthStratum(dm,1, &eStart, &eEnd);//Get LI of edges.
  ierr = DMPlexGetDepthStratum(dm,0, &vStart, &vEnd);//Get LI of vertices.
  ierr = DMGetStratumSize(dm,"depth", 3, &nc);//Get number of cells.
  ierr = DMGetStratumSize(dm,"depth", 2, &nf);//Get number of faces.
  ierr = DMGetStratumSize(dm,"depth", 1, &ne);//Get number of edges.
  ierr = DMGetStratumSize(dm,"depth", 0, &nv);//Get number of vertices.
  DMView(dm,PETSC_VIEWER_STDOUT_WORLD);

  PetscPrintf(PETSC_COMM_WORLD,"global start = %10d\t global end = %10d\n",pStart,pEnd);
  PetscPrintf(PETSC_COMM_WORLD,"#cells    = %10d\t i = %10d\t i < %10d\n",nc,cStart,cEnd);
  PetscPrintf(PETSC_COMM_WORLD,"#faces    = %10d\t i = %10d\t i < %10d\n",nf,fStart,fEnd);
  PetscPrintf(PETSC_COMM_WORLD,"#edges    = %10d\t i = %10d\t i < %10d\n",ne,eStart,eEnd);
  PetscPrintf(PETSC_COMM_WORLD,"#vertices = %10d\t i = %10d\t i < %10d\n",nv,vStart,vEnd);
  //DMLabel label;
  //ierr = DMGetLabel(dm, "X-Pin", &label);CHKERRQ(ierr);

  //MESH IMPORT=================================================================

  //NOTE: This section extremely hardcoded right now.
  //Current setup would only support P1 meshes.
  //MEMORY ALLOCATION ==========================================================
  ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s); CHKERRQ(ierr);
  //The chart is akin to a contiguous memory storage allocation. Each chart entry is associated
  //with a "thing," could be a vertex, face, cell, or edge, or anything else.
  ierr = PetscSectionSetChart(s, pStart, pEnd); CHKERRQ(ierr);
  //For each "thing" in the chart, additional room can be made. This is helpful for associating
  //nodes to multiple degrees of freedom. These commands help associate nodes with
  for(ii = cStart; ii < cEnd; ii++){//Cell loop.
    ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with cells.
  for(ii = fStart; ii < fEnd; ii++){//Face loop.
    ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with faces.
  for(ii = vStart; ii < vEnd; ii++){//Vertex loop.
    ierr = PetscSectionSetDof(s, ii, 3);CHKERRQ(ierr);}//Sets x, y, and z displacements as dofs.
  for(ii = eStart; ii < eEnd; ii++){//Edge loop
    ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with edges.
  ierr = PetscSectionSetUp(s); CHKERRQ(ierr);
  ierr = PetscSectionGetStorageSize(s,&sizeK);CHKERRQ(ierr);//Determine the size of the global stiffness matrix.
  ierr = DMSetLocalSection(dm,s); CHKERRQ(ierr);//Associate the PetscSection with the DM object.
  PetscPrintf(PETSC_COMM_WORLD,"sizeK = %10d\n",sizeK);

  //OBJECT SETUP================================================================
  //Global stiffness matrix.
  ierr = MatCreate(PETSC_COMM_WORLD,&K); CHKERRQ(ierr);
  ierr = MatSetType(K,MATAIJ); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.
  ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,sizeK,sizeK); CHKERRQ(ierr);

  //Displacement vector.
  ierr = VecCreate(PETSC_COMM_WORLD,&U); CHKERRQ(ierr);
  ierr = VecSetType(U,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.
  ierr = VecSetSizes(U,PETSC_DECIDE,sizeK); CHKERRQ(ierr);
  //OBJECT SETUP================================================================

for(ii=cStart;ii<cEnd;ii++){//loop through cells.
  ierr = DMPlexGetTransitiveClosure(dm, ii, useCone, &size_closure, &closure); CHKERRQ(ierr);//NOTE: NO need to restore the arrays.
  PetscInt cells=0, edges=0, vertices=0, faces=0;
  PetscPrintf(PETSC_COMM_WORLD,"cell # %5d\t",ii);

  for(jj=0;jj<(2*size_closure);jj+=2){//Scan the closure of the current cell.
    //Use information from the DM's strata to determine composition of cell_ii.
    if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices.
      vertices++;
      PetscPrintf(PETSC_COMM_WORLD,"%5d\t",closure[jj]);
    }//end if
    else if(eStart <= closure[jj] && closure[jj]< eEnd){//Check for edge ID's
      edges++;
    }//end else if
    else if(fStart <= closure[jj] && closure[jj]< fEnd){//Check for face ID's
      faces++;
    }//end else if
    else if(cStart <= closure[jj] && closure[jj]< cEnd){//Check for cell ID's
      cells++;
    }//end else if
    //else{//If index is out of bounds
      //const char *message = "Closure index out of bounds!";
      //ierr =  SETERRQ (PETSC_COMM_SELF, ierr,message);
    //}//end else
    //PetscPrintf(PETSC_COMM_WORLD,"%5d\n",closure[jj]);

  }//end "jj" loop.
  PetscPrintf(PETSC_COMM_WORLD,"%\n");
}//end "ii" loop.






  //CLEANUP=====================================================================
  DMDestroy(&dm);
  MatDestroy(&K);
  VecDestroy(&U);
  PetscSectionDestroy(&s);
  //CLEANUP=====================================================================
  //PetscErrorCode  PetscMallocDump(FILE *fp)
  //ierr = PetscMallocDump(NULL);
  return PetscFinalize();//And the machine shall rest....
  }