1! Contributed by Noem T 2#include <petsc/finclude/petscdmplex.h> 3program 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 16 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 24 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)] 33 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] 38 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 48 49 quadrilateral: block 50 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! 56 57! Number of quadrature points 58 PetscInt, parameter :: quad_qpts = 4 59 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 68 69 numCells = 1 70 numCorners = 4 71 numVertices = 4 72 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)) 85 86contains 87! Quadrature points in a quadrilateral in [-1,+1] 88 function Gauss_rs(n) 89 PetscInt, intent(in) :: n 90 PetscReal :: Gauss_rs(2) 91 92 PetscReal, parameter :: p = 1.0/sqrt(3.0) 93 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) 109 110 PetscReal :: r, s 111 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 116 117 end function quad_v 118! Jacobian 119 function quad_J(rs) 120 PetscReal, intent(in) :: rs(2) 121 PetscReal :: quad_J(4) 122 123 PetscReal :: r, s 124 PetscReal :: pfive = .5, twopfive = 2.5 125 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 130end program test 131 132! /*TEST 133! 134! test: 135! output_file: output/empty.out 136! 137! TEST*/ 138