1! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and 2! DMPlexComputeCellGeometryFVM() 3! Contributed by Adrian Croucher <a.croucher@auckland.ac.nz> 4 program main 5 use petscsys 6 use petscdmplex 7#include <petsc/finclude/petscsys.h> 8#include <petsc/finclude/petscdmplex.h> 9 implicit none 10 DM :: dm, dmi 11 PetscFV :: fvm 12 PetscInt, parameter :: dim = 3 13 14 PetscInt :: depth = 1 15 PetscErrorCode :: ierr 16 PetscInt, dimension(2) :: numPoints 17 PetscInt, dimension(14) :: coneSize 18 PetscInt, dimension(16) :: cones 19 PetscInt, dimension(16) :: coneOrientations 20 PetscScalar, dimension(36) :: vertexCoords 21 PetscReal :: vol = 0. 22 PetscReal, target, dimension(dim) :: centroid 23 PetscReal, target, dimension(dim) :: normal 24 PetscReal, pointer :: pcentroid(:) 25 PetscReal, pointer :: pnormal(:) 26 27 PetscReal, target, dimension(dim) :: v0 28 PetscReal, target, dimension(dim*dim) :: J 29 PetscReal, target, dimension(dim*dim) :: invJ 30 PetscReal, pointer :: pv0(:) 31 PetscReal, pointer :: pJ(:) 32 PetscReal, pointer :: pinvJ(:) 33 PetscReal :: detJ 34 35 PetscInt :: i 36 37 numPoints = [12, 2] 38 coneSize = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] 39 cones = [2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8] 40 coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0] 41 vertexCoords = [ & 42 & -0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, & 43 & -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, & 44 & 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0 ] 45 pcentroid => centroid 46 pnormal => normal 47 48 pv0 => v0 49 pJ => J 50 pinvJ => invJ 51 52 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 53 if (ierr .ne. 0) then 54 print*,'Unable to initialize PETSc' 55 stop 56 endif 57 58 call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr) 59 call PetscObjectSetName(dm, 'testplex', ierr);CHKERRA(ierr) 60 call DMSetDimension(dm, dim, ierr);CHKERRA(ierr) 61 62 call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones,coneOrientations, vertexCoords, ierr);CHKERRA(ierr) 63 64 call DMPlexInterpolate(dm, dmi, ierr);CHKERRA(ierr) 65 call DMPlexCopyCoordinates(dm, dmi, ierr);CHKERRA(ierr) 66 call DMDestroy(dm, ierr);CHKERRA(ierr) 67 dm = dmi 68 69 call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr) 70 71 do i = 0, 1 72 call DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr);CHKERRA(ierr) 73 write(*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ',pcentroid(1), pcentroid(2), pcentroid(3) 74 call DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ,detJ, ierr);CHKERRA(ierr) 75 end do 76 77 call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr);CHKERRA(ierr) 78 call PetscFVSetUp(fvm, ierr);CHKERRA(ierr) 79 call PetscFVDestroy(fvm, ierr);CHKERRA(ierr) 80 81 call DMDestroy(dm, ierr);CHKERRA(ierr) 82 call PetscFinalize(ierr) 83 end program main 84 85!/*TEST 86! 87! test: 88! suffix: 0 89! 90!TEST*/ 91