1: !     Contributed by Noem T
  2: #include <petsc/finclude/petscdmplex.h>
  3: module dmplexcomputecellgeometryfemmodule
  4:   use petscdm
  5:   use petscdmplex
  6:   implicit none

  8: contains
  9: ! Quadrature points in a quadrilateral in [-1,+1]
 10:   pure function Gauss_rs(n)
 11:     PetscInt, intent(in) :: n
 12:     PetscReal :: Gauss_rs(2)

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

 16:     select case (n)
 17:     case (1)
 18:       Gauss_rs = [-p, -p]
 19:     case (2)
 20:       Gauss_rs = [-p, +p]
 21:     case (3)
 22:       Gauss_rs = [+p, -p]
 23:     case (4)
 24:       Gauss_rs = [+p, +p]
 25:     end select
 26:   end function Gauss_rs

 28: ! Mapped quadrature points
 29:   pure function quad_v(rs)
 30:     PetscReal, intent(in) :: rs(2)
 31:     PetscReal :: quad_v(2)

 33:     PetscReal :: r, s

 35:     r = rs(1)
 36:     s = rs(2)
 37:     quad_v(1) = -0.5*r*(s - 5)       ! sum N_i * x_i
 38:     quad_v(2) = 0.5*(r + 5*s + 2)    ! sum N_i * y_i
 39:   end function quad_v

 41: ! Jacobian
 42:   pure function quad_J(rs)
 43:     PetscReal, intent(in) :: rs(2)
 44:     PetscReal :: quad_J(4)

 46:     PetscReal :: r, s

 48:     r = rs(1)
 49:     s = rs(2)
 50:     quad_J = [-0.5*(s - 5), -0.5*r, .5_PETSC_REAL_KIND, 2.5_PETSC_REAL_KIND]
 51:   end function quad_J
 52: end module dmplexcomputecellgeometryfemmodule

 54: program test
 55:   use petscdm
 56:   use petscdmplex
 57:   use dmplexcomputecellgeometryfemmodule

 59:   implicit none
 60:   DM       :: dm
 61:   PetscInt, parameter :: numCells = 1
 62:   PetscInt :: cStart
 63:   PetscErrorCode :: ierr
 64:   PetscInt, parameter :: numDim = 2
 65:   PetscReal, parameter :: eps = 5.0*epsilon(1.0)
 66:   PetscQuadrature :: q
 67:   PetscInt :: i

 69:   PetscCallA(PetscInitialize(ierr))

 71:   triangle: block
 72: !
 73: ! Single triangle element
 74: ! Corner coordinates (1, 1), (5, 5), (3, 6)
 75: ! Using a 3-point quadrature rule
 76: !
 77:     PetscInt, parameter :: numVertices = 3, numCorners = 3

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

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

 92:     PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, numDim, numCells, numVertices, numCorners, PETSC_FALSE, cells, numDim, vertexCoords, dm, ierr))
 93:     PetscCallA(PetscDTSimplexQuadrature(numDim, tria_qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, ierr))
 94:     PetscCallA(DMPlexGetHeightStratum(dm, 0_PETSC_INT_KIND, cStart, PETSC_NULL_INTEGER, ierr))
 95:     PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr))
 96:     PetscCheckA(all(abs(v - tria_v) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong mapped quadrature points (triangle)')
 97:     PetscCheckA(all(abs(J - tria_J) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong jacobian (triangle)')
 98:     PetscCallA(PetscQuadratureDestroy(q, ierr))
 99:     PetscCallA(DMDestroy(dm, ierr))
100:   end block triangle

102:   quadrilateral: block
103: !
104: ! Single quadrilateral element
105: ! Corner coordinates (-3, -2), (3, -1), (2, 4), (-2, 3)
106: ! Using a 4-point (2x2) Gauss quadrature rule
107: !
108:     PetscInt, parameter :: numVertices = 4, numCorners = 4

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

113:     PetscReal :: J(quad_qpts*numDim**2), invJ(quad_qpts*numDim**2), v(quad_qpts*numDim), detJ(quad_qpts)
114:     PetscReal, parameter :: vertexCoords(8) = [-3.0, -2.0, 3.0, -1.0, 2.0, 4.0, -2.0, 3.0]
115:     PetscInt, parameter  :: cells(4) = [0, 1, 2, 3]
116:     PetscReal, parameter :: a = -1.0, b = 1.0
117:     PetscInt, parameter :: nc = 1, npoints = 2

119:     PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, numDim, numCells, numVertices, numCorners, PETSC_FALSE, cells, numDim, vertexCoords, dm, ierr))
120:     PetscCallA(PetscDTGaussTensorQuadrature(numDim, nc, npoints, a, b, q, ierr))
121:     PetscCallA(DMPlexGetHeightStratum(dm, 0_PETSC_INT_KIND, cStart, PETSC_NULL_INTEGER, ierr))
122:     PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr))
123:     do i = 1, quad_qpts
124:       PetscCheckA(all(abs(v((i - 1)*numDim + 1:i*numDim) - quad_v(Gauss_rs(i))) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong mapped quadrature points (quadrilateral)')
125:       PetscCheckA(all(abs(J((i - 1)*numDim**2 + 1:i*numDim**2) - quad_J(Gauss_rs(i))) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong jacobian (quadrilateral)')
126:     end do
127:     PetscCallA(PetscQuadratureDestroy(q, ierr))
128:     PetscCallA(DMDestroy(dm, ierr))
129:   end block quadrilateral
130:   PetscCallA(PetscFinalize(ierr))

132: end program test

134: ! /*TEST
135: !
136: ! test:
137: !   output_file: output/empty.out
138: !
139: ! TEST*/