!     Contributed by Noem T
#include <petsc/finclude/petscdmplex.h>
program test
  use petscdm
  use petscdmplex
  implicit none
  DM       :: dm
  PetscInt :: numCells = 1
  PetscInt :: cStart
  PetscInt :: numVertices = 3, numCorners = 3
  PetscErrorCode :: ierr
  PetscInt, parameter :: numDim = 2
  PetscReal, parameter :: eps = 5.0*epsilon(1.0)
  PetscCallA(PetscInitialize(ierr))
  triangle: block

!
! Single triangle element
! Corner coordinates (1, 1), (5, 5), (3, 6)
! Using a 3-point quadrature rule
!
    PetscInt :: i
    PetscInt :: zero = 0

! Number of quadrature points
    PetscInt, parameter :: tria_qpts = 3
! Quadrature order
    PetscInt, parameter :: tria_qorder = 2
! Mapped quadrature points
    PetscReal, parameter :: tria_v(tria_qpts*numDim) = [2.0, 2.5, 3.0, 5.0, 4.0, 4.5]
! Jacobian (constant, repeated tria_qpts times)
    PetscReal, parameter :: tria_J(tria_qpts*numDim**2) = [([2.0, 1.0, 2.0, 2.5], i=1, tria_qpts)]

    PetscQuadrature :: q
    PetscReal :: J(tria_qpts*numDim**2), invJ(tria_qpts*numDim**2), v(tria_qpts*numDim), detJ(tria_qpts)
    PetscReal :: vertexCoords(6) = 1.0*[1.0, 1.0, 5.0, 5.0, 3.0, 6.0]
    PetscInt  :: cells(3) = [0, 1, 2]

    PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, numDim, numCells, numVertices, numCorners, PETSC_FALSE, cells, numDim, vertexCoords, dm, ierr))
    PetscCallA(PetscDTSimplexQuadrature(numDim, tria_qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, ierr))
    PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, PETSC_NULL_INTEGER, ierr))
    PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr))
    PetscCheckA(all(abs(v - tria_v) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong mapped quadrature points (triangle)')
    PetscCheckA(all(abs(J - tria_J) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong jacobian (triangle)')
    PetscCallA(PetscQuadratureDestroy(q, ierr))
    PetscCallA(DMDestroy(dm, ierr))
  end block triangle

  quadrilateral: block

!
! Single quadrilateral element
! Corner coordinates (-3, -2), (3, -1), (2, 4), (-2, 3)
! Using a 4-point (2x2) Gauss quadrature rule
!

! Number of quadrature points
    PetscInt, parameter :: quad_qpts = 4

    PetscQuadrature :: q
    PetscReal :: vertexCoords(8) = [-3.0, -2.0, 3.0, -1.0, 2.0, 4.0, -2.0, 3.0]
    PetscReal :: J(quad_qpts*numDim**2), invJ(quad_qpts*numDim**2), v(quad_qpts*numDim), detJ(quad_qpts)
    PetscReal :: a = -1.0, b = 1.0
    PetscInt  :: cells(4) = [0, 1, 2, 3]
    PetscInt  :: nc = 1, npoints = 2
    PetscInt  :: k
    PetscInt :: zero = 0

    numCells = 1
    numCorners = 4
    numVertices = 4

    PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, numDim, numCells, numVertices, numCorners, PETSC_FALSE, cells, numDim, vertexCoords, dm, ierr))
    PetscCallA(PetscDTGaussTensorQuadrature(numDim, nc, npoints, a, b, q, ierr))
    PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, PETSC_NULL_INTEGER, ierr))
    PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr))
    do k = 1, quad_qpts
      PetscCheckA(all(abs(v((k - 1)*numDim + 1:k*numDim) - quad_v(Gauss_rs(k))) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong mapped quadrature points (quadrilateral)')
      PetscCheckA(all(abs(J((k - 1)*numDim**2 + 1:k*numDim**2) - quad_J(Gauss_rs(k))) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong jacobian (quadrilateral)')
    end do
    PetscCallA(PetscQuadratureDestroy(q, ierr))
    PetscCallA(DMDestroy(dm, ierr))
  end block quadrilateral
  PetscCallA(PetscFinalize(ierr))

contains
! Quadrature points in a quadrilateral in [-1,+1]
  function Gauss_rs(n)
    PetscInt, intent(in) :: n
    PetscReal :: Gauss_rs(2)

    PetscReal, parameter :: p = 1.0/sqrt(3.0)

    select case (n)
    case (1)
      Gauss_rs = [-p, -p]
    case (2)
      Gauss_rs = [-p, +p]
    case (3)
      Gauss_rs = [+p, -p]
    case (4)
      Gauss_rs = [+p, +p]
    end select
  end function Gauss_rs
! Mapped quadrature points
  function quad_v(rs)
    PetscReal, intent(in) :: rs(2)
    PetscReal :: quad_v(2)

    PetscReal :: r, s

    r = rs(1)
    s = rs(2)
    quad_v(1) = -0.5*r*(s - 5)       ! sum N_i * x_i
    quad_v(2) = 0.5*(r + 5*s + 2)      ! sum N_i * y_i

  end function quad_v
! Jacobian
  function quad_J(rs)
    PetscReal, intent(in) :: rs(2)
    PetscReal :: quad_J(4)

    PetscReal :: r, s
    PetscReal :: pfive = .5, twopfive = 2.5

    r = rs(1)
    s = rs(2)
    quad_J = [-0.5*(s - 5), -0.5*r, pfive, twopfive]
  end function quad_J
end program test

! /*TEST
!
! test:
!   output_file: output/empty.out
!
! TEST*/