1! Contributed by Noem T 2program test 3#include <petsc/finclude/petscdmplex.h> 4#include <petsc/finclude/petscdm.h> 5use PETScDM 6use PETScDMplex 7implicit none 8DM :: dm 9PetscInt :: numCells = 1 10PetscInt :: cStart 11PetscInt :: numVertices = 3, numCorners = 3 12PetscErrorCode :: ierr 13PetscInt, parameter :: numDim = 2 14PetscReal, parameter :: eps = 5.0 * epsilon(1.0) 15PetscCallA(PetscInitialize(ierr)) 16triangle: block 17 18! 19! Single triangle element 20! Corner coordinates (1, 1), (5, 5), (3, 6) 21! Using a 3-point quadrature rule 22! 23PetscInt :: i 24PetscInt :: zero = 0 25 26! Number of quadrature points 27PetscInt, parameter :: tria_qpts = 3 28! Quadrature order 29PetscInt, parameter :: tria_qorder = 2 30! Mapped quadrature points 31PetscReal, parameter :: tria_v(tria_qpts*numDim) = [2.0, 2.5, 3.0, 5.0, 4.0, 4.5] 32! Jacobian (constant, repeated tria_qpts times) 33PetscReal, parameter :: tria_J(tria_qpts*numDim**2) = [([2.0, 1.0, 2.0, 2.5], i = 1, tria_qpts)] 34 35PetscQuadrature :: q 36PetscReal :: J(tria_qpts*numDim**2), invJ(tria_qpts*numDim**2), v(tria_qpts*numDim), detJ(tria_qpts) 37PetscReal :: vertexCoords(6) = 1.0 * [1.0, 1.0, 5.0, 5.0, 3.0, 6.0] 38PetscInt :: cells(3) = [0, 1, 2] 39 40PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD,numDim,numCells,numVertices,numCorners,PETSC_FALSE,cells,numDim,vertexCoords,dm,ierr)) 41PetscCallA(PetscDTSimplexQuadrature(numDim, tria_qorder, PETSCDTSIMPLEXQUAD_DEFAULT, q, ierr)) 42PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, PETSC_NULL_INTEGER, ierr)) 43PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr)) 44PetscCheckA(all(abs(v - tria_v) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong mapped quadrature points (triangle)') 45PetscCheckA(all(abs(J - tria_J) < eps), PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Wrong jacobian (triangle)') 46PetscCallA(PetscQuadratureDestroy(q, ierr)) 47PetscCallA(DMDestroy(dm,ierr)) 48end block triangle 49 50quadrilateral: block 51 52! 53! Single quadrilateral element 54! Corner coordinates (-3, -2), (3, -1), (2, 4), (-2, 3) 55! Using a 4-point (2x2) Gauss quadrature rule 56! 57 58! Number of quadrature points 59PetscInt, parameter :: quad_qpts = 4 60 61PetscQuadrature :: q 62PetscReal :: vertexCoords(8) = [-3.0, -2.0, 3.0, -1.0, 2.0, 4.0, -2.0, 3.0] 63PetscReal :: J(quad_qpts*numDim**2), invJ(quad_qpts*numDim**2), v(quad_qpts*numDim), detJ(quad_qpts) 64PetscReal :: a = -1.0, b = 1.0 65PetscInt :: cells(4) = [0, 1, 2, 3] 66PetscInt :: nc = 1, npoints = 2 67PetscInt :: k 68PetscInt :: zero = 0 69 70numCells = 1 71numCorners = 4 72numVertices = 4 73 74PetscCallA(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD,numDim,numCells,numVertices, numCorners,PETSC_FALSE,cells,numDim,vertexCoords,dm,ierr)) 75PetscCallA(PetscDTGaussTensorQuadrature(numDim, nc, npoints, a, b, q, ierr)) 76PetscCallA(DMPlexGetHeightStratum(dm, zero, cStart, PETSC_NULL_INTEGER, ierr)) 77PetscCallA(DMPlexComputeCellGeometryFEM(dm, cStart, q, v, J, invJ, detJ, ierr)) 78do k = 1, quad_qpts 79 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)') 80 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)') 81end do 82PetscCallA(PetscQuadratureDestroy(q, ierr)) 83PetscCallA(DMDestroy(dm,ierr)) 84end block quadrilateral 85PetscCallA(PetscFinalize(ierr)) 86 87contains 88! Quadrature points in a quadrilateral in [-1,+1] 89function Gauss_rs(n) 90PetscInt, intent(in) :: n 91PetscReal :: Gauss_rs(2) 92 93 PetscReal, parameter :: p = 1.0/sqrt(3.0) 94 95 select case(n) 96 case(1) 97 Gauss_rs = [-p, -p] 98 case(2) 99 Gauss_rs = [-p, +p] 100 case(3) 101 Gauss_rs = [+p, -p] 102 case(4) 103 Gauss_rs = [+p, +p] 104end select 105end function Gauss_rs 106! Mapped quadrature points 107function quad_v(rs) 108PetscReal, intent(in) :: rs(2) 109PetscReal :: quad_v(2) 110 111 PetscReal :: r, s 112 113 r = rs(1) 114 s = rs(2) 115 quad_v(1) = -0.5*r*(s-5) ! sum N_i * x_i 116 quad_v(2) = 0.5*(r+5*s+2) ! sum N_i * y_i 117 118end function quad_v 119! Jacobian 120function quad_J(rs) 121PetscReal, intent(in) :: rs(2) 122PetscReal :: quad_J(4) 123 124 PetscReal :: r, s 125 PetscReal :: pfive = .5, twopfive = 2.5 126 127 r = rs(1) 128 s = rs(2) 129 quad_J = [-0.5*(s-5), -0.5*r, pfive, twopfive] 130end function quad_J 131end program test 132 133! /*TEST 134! 135! test: 136! output_file: output/empty.out 137! 138! TEST*/ 139