program fe_test
#include <petsc/finclude/petsc.h>

  use petsc
  implicit none


  PetscInt, parameter :: dim = 3
  PetscInt :: depth = 1
  
  PetscInt, dimension(2) :: numPoints
  PetscInt, dimension(14) :: coneSize
  PetscInt, dimension(16) :: cones
  PetscInt, dimension(16) :: coneOrientations
  PetscScalar, dimension(36) :: vertexCoords
  
  DM :: dm
  PetscFE :: fem
  PetscInt :: start_cell, end_cell, num_cells
  PetscQuadrature :: quad
  PetscInt :: q_dim, q_num, q_nc
  PetscReal, target, dimension(8) :: q_weights
  PetscReal, target, dimension(24) :: q_points
  PetscReal, pointer :: pq_points(:)
  PetscReal, pointer :: pq_weights(:)

  PetscInt :: i,j
  
  PetscErrorCode :: ierr

  pq_points => q_points
  pq_weights => q_weights
  
  numPoints = [12, 2]
  coneSize  = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  cones = [2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8]
  coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0]
  vertexCoords = [                                                  &
       &  -0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,     &
       &  -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,     &
       &   0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0 ]

  call PetscInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)
  
  call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
  CHKERRQ(ierr)
  call PetscObjectSetName(dm, 'testplex', ierr)
  CHKERRQ(ierr)
  call DMSetDimension(dm, dim, ierr)
  CHKERRQ(ierr)

  call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, &
       cones,coneOrientations, vertexCoords, ierr)
  CHKERRQ(ierr)
  call DMPlexGetHeightStratum(dm, 0, start_cell, end_cell, ierr)
  CHKERRQ(ierr)

  num_cells = end_cell-start_cell

  ! -----------------------------------------------------------------------
  ! FE routine stuff
  ! -----------------------------------------------------------------------
  call PetscFECreateDefault(dm,dim,num_cells,PETSC_FALSE,"fe",1,fem,ierr)
  CHKERRQ(ierr)
  call PetscFEGetQuadrature(fem,quad,ierr); CHKERRQ(ierr)
  !call PetscQuadratureView(quad,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr)
  call PetscQuadratureGetData(quad,q_nc,q_dim,q_num,pq_points,pq_weights,ierr);
  CHKERRQ(ierr)
  write(*,*) 'weights:',pq_weights
  write(*,*) 'points:',pq_points

  !CHKERRQ(ierr)
  call PetscFEDestroy(fem,ierr); CHKERRQ(ierr)
  ! -----------------------------------------------------------------------
  
  call DMView(dm,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr)
  call DMDestroy(dm, ierr);CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program fe_test