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