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

 17: !
 18: ! Single triangle element
 19: ! Corner coordinates (1, 1), (5, 5), (3, 6)
 20: ! Using a 3-point quadrature rule
 21: !
 22:     PetscInt :: i
 23:     PetscInt :: zero = 0

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

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

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

 49:   quadrilateral: block

 51: !
 52: ! Single quadrilateral element
 53: ! Corner coordinates (-3, -2), (3, -1), (2, 4), (-2, 3)
 54: ! Using a 4-point (2x2) Gauss quadrature rule
 55: !

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

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

 69:     numCells = 1
 70:     numCorners = 4
 71:     numVertices = 4

 73:     PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, numDim, numCells, numVertices, numCorners, PETSC_FALSE, cells, numDim, vertexCoords, dm, ierr))
 74:     PetscCallA(PetscDTGaussTensorQuadrature(numDim, nc, npoints, a, b, q, ierr))
 75:     PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, PETSC_NULL_INTEGER, ierr))
 76:     PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr))
 77:     do k = 1, quad_qpts
 78:       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)')
 79:       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)')
 80:     end do
 81:     PetscCallA(PetscQuadratureDestroy(q, ierr))
 82:     PetscCallA(DMDestroy(dm, ierr))
 83:   end block quadrilateral
 84:   PetscCallA(PetscFinalize(ierr))

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

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

 94:     select case (n)
 95:     case (1)
 96:       Gauss_rs = [-p, -p]
 97:     case (2)
 98:       Gauss_rs = [-p, +p]
 99:     case (3)
100:       Gauss_rs = [+p, -p]
101:     case (4)
102:       Gauss_rs = [+p, +p]
103:     end select
104:   end function Gauss_rs
105: ! Mapped quadrature points
106:   function quad_v(rs)
107:     PetscReal, intent(in) :: rs(2)
108:     PetscReal :: quad_v(2)

110:     PetscReal :: r, s

112:     r = rs(1)
113:     s = rs(2)
114:     quad_v(1) = -0.5*r*(s - 5)       ! sum N_i * x_i
115:     quad_v(2) = 0.5*(r + 5*s + 2)      ! sum N_i * y_i

117:   end function quad_v
118: ! Jacobian
119:   function quad_J(rs)
120:     PetscReal, intent(in) :: rs(2)
121:     PetscReal :: quad_J(4)

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

126:     r = rs(1)
127:     s = rs(2)
128:     quad_J = [-0.5*(s - 5), -0.5*r, pfive, twopfive]
129:   end function quad_J
130: end program test

132: ! /*TEST
133: !
134: ! test:
135: !   output_file: output/empty.out
136: !
137: ! TEST*/