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 PetscCallA(PetscInitialize(ierr)) 53 PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)) 54 PetscCallA(PetscObjectSetName(dm, 'testplex', ierr)) 55 PetscCallA(DMSetDimension(dm, dim, ierr)) 56 57 PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones,coneOrientations, vertexCoords, ierr)) 58 59 PetscCallA(DMPlexInterpolate(dm, dmi, ierr)) 60 PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr)) 61 PetscCallA(DMDestroy(dm, ierr)) 62 dm = dmi 63 64 PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)) 65 66 do i = 0, 1 67 PetscCallA(DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr)) 68 write(*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ',pcentroid(1), pcentroid(2), pcentroid(3) 69 PetscCallA(DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ,detJ, ierr)) 70 end do 71 72 PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr)) 73 PetscCallA(PetscFVSetUp(fvm, ierr)) 74 PetscCallA(PetscFVDestroy(fvm, ierr)) 75 76 PetscCallA(DMDestroy(dm, ierr)) 77 PetscCallA(PetscFinalize(ierr)) 78 end program main 79 80!/*TEST 81! 82! test: 83! suffix: 0 84! 85!TEST*/ 86