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